Revealing intrinsic electric double layer viscoelasticity in ionic liquid solutions via quartz crystal microbalance

Abstract

Ionic liquids (IL) are room-temperature molten salts that function as complex fluids with tunable interfacial and bulk properties, making them attractive for applications ranging from electrochemical energy storage to lubrication. In such systems, the viscoelasticity of the electric double layer (EDL) at charged interfaces can strongly influence performance, yet its characterization remains challenging due to the nanometric EDL thickness. Herein, we use quartz crystal microbalance (QCM) to measure changes in the resonant frequency and energy dissipation of a gold-coated quartz crystal upon deposition of IL solutions. Since the gold surface of the QCM is negatively charged at an open-circuit potential, we can estimate the loss modulus of the EDL near the charged surface through a wave propagation model under non-confining conditions. Using 1-butyl-3-methylimidazolium (Bmim)-based ILs with three distinct anions-bis(trifluoromethanesulfonyl)imide (TFSI), trifluoromethanesulfonate (TfO), and tetrafluoroborate (BF₄), we find that the EDL loss modulus increases sharply with increasing IL concentrations in the low concentration regime, eventually reaching values up to three orders of magnitude higher than that of the bulk solution and saturating at high concentrations. Notably, this concentration-dependent scaling is consistent across the three anion types tested, in contrast to reports for nanoconfined ILs where ion identity markedly affects this behavior. Our results demonstrate that bulk viscoelastic properties can be used to infer the EDL loss modulus under non-confining conditions, providing a practical framework for engineering soft, ion-rich interfaces in electrochemical and tribological systems.

---

1 Introduction

Ionic liquids (ILs) are molten salts that consist of weakly coordinated bulky cations and anions.¹ Owing to their low lattice energy,² ILs can remain in a liquid state at room temperature and exhibit unique physicochemical properties, such as non-flammability due to negligible vapor pressure, high ionic conductivity, and high thermal and electrochemical stabilities. These attributes have motivated applications across diverse fields, such as energy storage and conversion,³ sensing,⁴ water treatment,⁵ lubrication,⁶ and drug delivery.⁷ The performance of IL-based materials can be tailored by appropriate selection of cation and anion structures^8–10 and by adjusting IL concentration, which influences properties, such as viscosity and ionic conductivity, in mixtures with non-ionic organic solvents.^11–15 Understanding how ionic structures and concentrations govern IL material properties is therefore central to optimizing IL-based system performance.

Ionic liquids often have direct contact with solid surfaces during material synthesis and processing-related applications (e.g., batteries and lubrication).^16–18 When the solid surfaces are electrified, electrostatic attraction drives the assembly of ions at the interface to screen the surface charge, forming an electric double layer (EDL) that extends over several nanometers, depending on the ion concentration.¹⁹ In many of these systems, the structure and dynamics of the EDL can affect interfacial charge and mass transport,^20–25 energy conversion and storage,^26–29 fluid lubrication,^30–32 and stability of colloids under electrical or mechanical perturbations,^33–37 thereby strongly influencing overall system performance. The EDL structure in IL solutions has been extensively investigated using surface force apparatus (SFA)^38–43 and atomic force microscopy (AFM).^44–48 At low IL concentrations, the EDL structure is well described using the Gouy–Chapman–Stern (GCS) model.¹⁹ In the GCS framework, a compact layer of oppositely charged ions is first formed directly on the charge surface to partially reduce the surface potential. Beyond this region, a diffuse layer is formed, and the potential decays exponentially toward its bulk level. On the other hand, at high IL concentrations, the EDL structure deviates significantly from the GCS model due to strong ion correlations. Specifically, EDLs in IL solutions exhibit overscreening, where the first ion layer overcompensates the surface charge, followed by underscreening, in which the charge-density profile decays oscillatory toward the bulk with a screening length that increases with concentration.^49–51 This anomalous behavior at high IL concentrations has been confirmed through fluorescence imaging⁵² and scattering⁵³ experiments, and is consistent with theoretical predictions.^54–61 However, directly probing EDL viscoelasticity remains difficult because of its nanometer-scale dimensions involved and the complex coupling between electrostatic and hydrodynamic effects. A deeper understanding of EDL viscoelasticity may offer valuable insights into the selection and design of ions for optimizing material processing via controlled fluid transport, as well as for enhancing the performance of electric double layer capacitors, lubrication systems, and electrochemical reactions.

To address this challenge, the viscoelastic response of confined EDLs formed in pure ILs has been investigated using SFA- and AFM-based techniques, showing that their shear viscosity depends strongly on ion identity, confinement geometry, and surface chemistry.^62–67 For example, Kurihara and co-workers^62,63 found that the shear viscosity of pure ILs confined between charged surfaces was 1–3 orders of magnitude greater than the bulk value and decayed exponentially with distance from the surface. Such differences have been attributed to confinement-induced ionic ordering,⁶⁸ corroborated by X-ray diffraction measurements.⁶⁹ While SFA and AFM studies have shown dramatic confinement-induced changes in IL viscosity, such measurements inherently obscure the intrinsic viscoelasticity of EDLs. In addition, it remains unclear how these results can be translated to open, non-confining geometries, more relevant to practical electrochemical and tribological applications.^16–18

In this work, we use quartz crystal microbalance (QCM), a non-confining, surface-sensitive technique, to isolate and quantify the loss modulus of EDLs at charged gold surfaces. QCM probes EDL dynamics by coupling oscillatory shear motion at a solid–liquid interface with frequency-shift analysis under non-confined conditions.⁷⁰ With suitable acoustic wave propagation models, QCM can provide storage and loss moduli of interfacial layers. While this approach has been applied to aqueous electrolytes,^71–73 polymer solutions,⁷⁴ and cells,^75,76 it has rarely been used for ILs, particularly for quantifying the EDL complex modulus without nanoconfinements.

We combine QCM measurements with a shear-wave propagation model⁷⁷ to quantify the EDL loss modulus in IL solutions under non-confining conditions. 1-Butyl-3-methylimidazolium(Bmim)-based ILs with three representative anions, bis(trifluoromethanesulfonyl)imide (TFSI), trifluoromethanesulfonate (TfO), and tetrafluoroborate (BF₄), are selected for their differences in size and symmetry (see Fig. 1). By systematically varying IL concentration, we identify scaling behaviors and saturation trends in the EDL loss modulus, and assess their dependence on anion type. This work advances understanding of ion-rich interface rheology and establishes a framework for linking bulk viscoelastic properties to EDL mechanics in IL-based electrochemical and tribological systems.

Fig. 1 Chemical structure of Bmim–TFSI, Bmim–TfO, and Bmim–BF₄.

2 Experimental methods

2.1 Materials

1-Butyl-3-methylimidazolium-based ionic liquids having three different anions, i.e., Bmim–TFSI (IL-0029-UP), Bmim–BF₄ (IL-0013-HP), and Bmim–TfO (IL-0012-HP), were purchased from Ionic Liquid Technologies (Germany). The water content in these ionic liquids, provided by the manufacturer, was less than 100 ppm for Bmim–TFSI, 250 ppm for Bmim–BF₄, and 500 ppm for Bmim–TfO. Ethylene glycol (EG) was purchased from Sigma-Aldrich (Japan) and used as a reference fluid, as it contains no ions and therefore does not form an EDL. Super dehydrated N,N-dimethylformamide (DMF) was purchased from FUJIFILM Wako Chemicals (Japan) and used as a solvent.

2.2 Sample preparation

The tested solutions were prepared by mixing the components directly in a glass vial under ambient conditions (room temperature: approx. 25 °C; humidity: approx. 40%). Specifically, the desired amount of solute was first added to a glass vial, and its mass m_s was measured using an electronic microbalance (AP225W, SHIMADZU; accuracy: ±0.1 mg). The solute was then diluted with DMF, and the total mass of the resulting mixture, m_tot, was measured. The prepared solution was hermetically sealed into the glass vial with a Teflon tape and stored at room temperature before use.

The concentration of solutes, c_s, is commonly calculated as the ratio of the molar amount of solutes to the volume of the solvent. However, in our tested solutions, the volume of the resulting mixture may deviate from the volume of the added solvent, especially at high solute concentrations. On the other hand, the total mass of the components is conserved before and after mixing. Therefore, we calculated c_s as

where ρ_B is the density of the prepared solution, and M_s is the molar mass of the solute, calculated as the sum of the atomic molar masses comprising the solute. This calculation yields c_s = 3.42 M for pure Bmim–TFSI, c_s = 4.49 M for pure Bmim–TfO, c_s = 5.31 M for pure Bmim–BF₄, and c_s = 17.1 M for pure EG, since in these cases c_s = ρ_B/M_s. The number of significant digits for c_s was chosen based on the measurement accuracy of the solution density, which leads to a larger uncertainty than the mass measurements.

2.3 Density measurements

The density ρ_B of the bulk solution, (i.e., the mixture of DMF and an ionic liquid or EG), was measured by using a density meter (DMA 35 Basic, Anton Paar; accuracy: ±0.001 g cm⁻³) at room temperature. Fig. 2 plots the measured values of ρ_B as a function of c_s, displaying a monotonic increase of ρ_B with increasing c_s for both IL and EG solutions. The measured density of pure DMF was ρ_DMF = (944 ± 0.4) kg m⁻³, in good agreement with the literature value.⁷⁸ We fitted the measured ρ_B using a quadratic function of ρ_B = ρ_DMF + Ac_s + Bc_s², and obtained the best fit curve for each solution (dashed lines in Fig. 2). The values of the corresponding fitting coefficients are given in Table 1 and used to estimate the dependence of the density ρ_EDL of the EDL on c_s, given by eqn (2) in Section 2.6.2.

Fig. 2 The density ρ_B of the bulk solution at 25 °C is plotted as a function of the solute concentration c_s for Bmim–TFSI (black squares), Bmim–TfO (red circles), Bmim–BF₄ (blue triangles), and EG (green stars) solutions in DMF. The dashed line represents the best fit curve to the measured ρ_B for each solution, obtained using a quadratic function of ρ_B = ρ_DMF + Ac_s + Bc_s². The corresponding fitting coefficients providing the best fit are provided in Table 1.

Table 1 Fitting coefficients providing the best fit to the empirical polynomial functions for the density (A, B) and shear viscosity (D, E, F) of Bmim–TFSI, Bmim–TfO, Bmim–BF₄, and EG solutions in DMF, and the dielectric constant (H, I) of the Bmim–TFSI, Bmim–TfO, and Bmim–BF₄ solutions in DMF

	Bmim–TFSI	Bmim–TfO	Bmim–BF4	EG
A (10^−2 kg mol^−1)	14.5 ± 0.1	8.37 ± 0.04	5.74 ± 0.33	1.15 ± 0.18
B (10^−3 kg mol^−2 L)	0	−1.39 ± 0.19	−2.32 ± 0.15	−0.546 ± 0.089
D (10^−4 Pa s mol^−1 L)	0	4.52 ± 0.13	5.24 ± 0.14	−1.01 ± 0.23
E (10^−4 Pa s mol^−2 L^2)	8.17 ± 0.69	0	0	1.64 ± 0.26
F (10^−4 Pa s mol^−3 L^3)	1.28 ± 0.32	2.17 ± 0.03	1.42 ± 0.03	−0.27 ± 0.07
H (mol^−1 L)	−8.95 ± 0.13	−6.72 ± 0.17	−6.05 ± 0.06	N/A
I (10^−1 mol^−2 L^2)	5.82 ± 0.51	2.77 ± 0.66	2.62 ± 0.22	N/A

---

2.4 Shear viscosity measurements

The shear viscosity η_B of the bulk solution was measured using a strain-controlled rheometer (ARES-G2, TA Instruments) with shear rates ([small gamma, Greek, dot above]) ranging from 0.01 to 1000 s⁻¹. A stainless steel cone and plate geometry with 50 mm diameter and 1° cone angle was used as an upper geometry element, while a stainless steel flat plate with 60 mm diameter was attached to an advanced Peltier system (TA Instruments) as a lower geometry element to regulate the temperature with an accuracy of ±0.1 °C. The temperature was set at 25 °C. The value of η_B was determined by averaging the measured shear viscosities over shear rates in the Newtonian regime where η_B is independent of [small gamma, Greek, dot above]. The shear viscosity curves of the tested solutions are provided in Fig. S3 of the SI.

Fig. 3 shows the dependence of bulk shear viscosity η_B on c_s for Bmim–TFSI, Bmim–TfO, Bmim–BF₄, and EG solutions in DMF at 25 °C. The measured shear viscosity of pure DMF was η_DMF = 0.917 ± 0.047 mPa s, consistent with the literature.⁷⁹ The value of η_B increased monotonically with increasing c_s, and the noticeable curvature was captured by fitting the measured η_B(c_s) with the lowest-order polynomial that avoided overfitting: η_B = η_DMF + Dc_s + Ec_s² + Fc_s³. This polynomial form is a common empirical representation of mixture viscosities at fixed temperature and aligns with the low-order series expansions that approximate electrolyte and liquid-mixture models within limited concentration ranges.^80,81 The best fit curves are shown as the dashed lines in Fig. 3, and the obtained corresponding fitting coefficients are given in Table 1. These fitting coefficients were then used to estimate the dependence of the loss modulus of the EDL on c_s, given by eqn (4) in Section 2.6.2. Note that the loss modulus of the bulk solution can be expressed as where f₀ is the fundamental frequency of the quartz employed in our QCM measurements.⁸²

Fig. 3 The shear viscosity η_B of the bulk solution at 25 °C is plotted as a function of the solute concentration c_s for Bmim–TFSI (black squares), Bmim–TfO (red circles), Bmim–BF₄ (blue triangles), and EG (green stars) solutions in DMF. The dashed line represents the best fit curve to the measured η_B for each solution, obtained using a cubic function of η_B = η_DMF + Dc_s + Ec_s² + Fc_s³. The corresponding fitting coefficients providing the best fit are provided in Table 1.

2.5 Dielectric constant measurements

The dielectric constant ε_B of the bulk solution was measured using a RF LCR meter (4287A, Agilent Technologies) equipped with a homemade stainless steel electrode cell. Specifically, the complex permittivity ε* = ε′ − iε″ was measured with varying frequencies from 1 MHz to 3 GHz. The instrument was calibrated using an ultrapure reference water sample (Wako Pure Chemicals). The vacant electric capacitance C₀ was measured without samples and obtained as C₀ = 0.23 pF. After loading a test solution in the cell, the experiment was performed at room temperature (∼25 °C), and the electric capacitance C and conductance K of the solution as function of frequency were determined using a parallel connected equivalent circuit model. The real and imaginary parts of the complex permittivity were then calculated as ε′ = CC₀⁻¹ and ε″ = K(2πC₀f)⁻¹, respectively. The value of ε_B was determined by averaging ε′ values over frequencies where the measured ε′ was independent of frequency.

Fig. 4 shows the dependence of ε_B on c_s for Bmim–TFSI, Bmim–TfO, and Bmim–BF₄ solutions in DMF at 25 °C. In the plot, the values of ε_B for pure Bmim–TfO (filled circle) and pure Bmim–BF₄ (filled triangle) were taken from the literature.⁸³ The dielectric constant of the tested IL solution decreased monotonically with increasing c_s over the measured c_s range. We fitted the values of ε_B by using a quadratic function of ε_B = ε_DMF + Hc_s + Ic_s², and the obtained best fits are shown as the dashed curves in Fig. 4. Here, ε_DMF = 38.3 is the dielectric constant of pure DMF measured using the RF LCR meter. The values of the corresponding fitting coefficients are provided in Table 1 and used to estimate the Debye screening length λ_Debye, given by eqn (7) in Section 2.6.2.

Fig. 4 The dielectric constant ε_B of the bulk solution at 25 °C is plotted as a function of the solute concentration c_s for Bmim–TFSI (black squares), Bmim–TfO (red circles), and Bmim–BF₄ (blue triangles) solutions in DMF. The dashed line represents the best fit curve to the measured ε_B for each solution, obtained using a quadratic function of ε_B = ε_DMF + Hc_s + Ic_s². The corresponding fitting coefficients providing the best fit are provided in Table 1.

2.6 QCM measurements and data analysis

2.6.1 QCM-D experimental procedure. The complex modulus of the EDL was evaluated using a QCM-D (openQCM Q-1, Novaetech, S.r.l. Italy) equipped with a custom made fused silica open cell. An oscillating potential of 0.7 V amplitude was applied to an AT-CUT gold-coated quartz (Industria Elettronica Varese, S.r.l.) in the QCM device, enabling the quartz to resonate at a fundamental frequency f₀ of ≈10 MHz in air. Although higher overtones are commonly employed in QCM-D studies to minimize the contribution from the bulk solution and to emphasize the near-surface response,^70,84,85 we performed measurements at the fundamental frequency because the quartz crystal exhibits its maximum resonance amplitude in air at this mode. Under these conditions, the measurable solute concentration range is broader than that obtained at overtone frequencies. As shown in Fig. 5, the resonance amplitude decreased with increasing c_s, reflecting the rise in fluid viscoelasticity. A higher starting amplitude at the fundamental frequency therefore extends the accessible concentration range before the signal becomes too attenuated.

Fig. 5 The amplitude–frequency curve of the mixture of Bmim–TFSI and DMF at various representative Bmim–TFSI concentrations, measured using QCM-D at room temperature (∼25 °C). The dotted arrow guides the eye to see how the QCM spectrum changes with the increasing concentration of the Bmim–TFSI.

After the resonance of the quartz stabilized over time, a sample solution with a volume of 200 µL was loaded on the quartz resonator to measure the change of the resonant frequency and energy dissipation over time. Since the gold-coated quartz carries a net negative charge on its surface,⁸⁶ we expect that an EDL is formed when the tested solution contains ions. Finally, we stopped recording when the frequency and dissipation signals became stable. After each QCM measurement, the quartz resonator was cleaned by sonication in acetone, isopropyl alcohol, and Milli-Q water for 5 min, and prepared for the follow-up QCM experiments at different c_s. Measurements were performed in triplicate at each c_s to obtain an estimation of the standard deviation. All experiments were conducted at room temperature (∼25 °C).

The resonant frequency f and energy dissipation D were determined as the peak frequency of the resonant peak and the full width at half maximum divided by the peak frequency respectively, shown in the amplitude vs. frequency plot. Fig. 5 displays representative amplitude–frequency curves at various c_s observed for the mixture of DMF and Bmim–TFSI at different concentrations. The exposure of the quartz crystal to DMF at c_s = 0 M produced 2 major effects on the resonance peak of the piezoelectric material. First, the sharp and intense peak observed in air was damped and broadened by the exposure to the solvent. Second, the peak position shifted toward lower frequencies. Both these spectral changes are due to variations in the mass of the fluid, caused by the variation of ρ_B (Fig. 2), as well as due to the viscous damping of the quartz oscillation, associated mainly with the variation of η_B (Fig. 3). For the mixture of Bmim–TFSI and DMF, increasing c_s led to further downshifts in frequency and increased peak broadening, reflecting the combined effects of higher density and viscosity.

2.6.2 QCM data analysis. We evaluated the viscoelastic properties of the EDL by using a wave propagation model proposed previously to predict the QCM response for ordinary electrolyte solutions.⁷⁷ In this model, we assumed a viscoelastic layer of the EDL with density ρ_EDL and complex modulus decaying exponentially over distance toward the complex modulus of the bulk solution:where ρ_EDL,max, , and are the limiting values of density and storage and loss moduli of the EDL, i.e., the density and complex modulus near the surface. Here, the bulk solution is assumed to be a Newtonian fluid, i.e., and . Since ions are condensed near the charged surface due to the electrostatic attraction even at relatively low c_s, we also assume ρ_EDL,max to be the density at the saturation concentration, i.e., ρ_B of the pure IL. In eqn (2)–(4), the characteristic length of the exponential decay for ρ_EDL and is represented by the charge screening length λ_EDL for IL solutions.

Here, we recall the recent experimental report⁴² based on the surface force apparatus that the Debye–Hückel theory was not applicable to predict λ_EDL for IL solutions at high c_s, and the value of λ_EDL exhibited a non-monotonic dependence on c_s. In order to capture the observed non-monotonic trend of λ_EDL against c_s, we employ a model proposed by Lee et al.⁸⁷ as,

λ_EDL = λ_Debye for λ_Debye ≥ a (5)

where the Debye screening length λ_Debye is given by¹⁹

Here, a, ε₀, k_B, N_A, T, and e are the ion diameter, the dielectric constant of the vacuum, the Boltzmann constant, the Avogadro constant, the absolute temperature, and the elementary charge. Following Lee's method, we estimate the ion diameter of the pure ionic liquid as . Fig. 6 shows the predicted non-monotonic dependence of λ_EDL at 25 °C for each IL solution, estimated using eqn (5) and (6). The value of λ_EDL for our IL solutions decreased significantly with increasing c_s for c_s < 0.25 M, followed by a gradual increase in λ_EDL at higher c_s.

Fig. 6 The non-monotonic dependence of the charge screening length λ_EDL at 25 °C on c_s for Bmim–TFSI (black solid), Bmim–TfO (red dashed), and Bmim–BF₄ (blue dotted) solutions in DMF.

When a viscoelastic EDL layer is formed adjacent to the bulk solution, the change in magnitude of the resonant frequency f_exp(c_s) and dissipation D_exp(c_s) of the quartz with respect to those in air can be expressed using a wave propagation model proposed by Voinova et al.⁸⁸ Voinova's model accounts for the oscillation of a quartz covered by viscoelastic layers in a semi-infinite Newtonian fluid, and predicts f_exp(c_s) and D_exp(c_s) as:

where the contributions of the bulk solution (denoted by the subscript B) and the EDL are given bywhere ρ_q = 2.648 g cm⁻³ and G_q = 2.947 × 10¹¹ g cm⁻¹ s⁻² are the density and the storage shear modulus of the AT-cut quartz crystal, respectively. Here, it should be noted that eqn (8) may not hold when the resonant frequency shift is estimated as the change in magnitude of frequencies in fluids with respect to those in air.⁸⁹ The mechanism behind this breakdown (eqn (8)) is still unknown but has been attributed to compressional waves generated by the flexural deformation of quartz crystals.⁸⁹ Since the compressional wave effects in fluids are more significant than those in air, they can never be eliminated and thus contribute to f_exp(c_s). In this case, assuming that the compressional-wave effects on the resonant frequency are independent of c_s, eqn (8) can be modified as f_exp(c_s) = f_B(c_s) + f_EDL(c_s) + f_com, where f_com is the concentration-independent offset arising from compressional-wave coupling.

We also note that the viscoelastic relationship described above, based on QCM measurements, is restricted to solutions at relatively low c_s with relatively low viscosities and densities. We assume that the contribution of the bulk solution to the measured frequency and dissipation shifts can be eliminated by using eqn (10) and (11). However, as shown in Fig. S1 of the SI, we observed f_exp(c_s) ≠ f_B(c_s) for EG solutions at c_s > 8 M in the absence of EDL (i.e., f_EDL(c_s) = 0). These results suggest that f_B(c_s) is no longer a simple function of η_B and ρ_B, as expressed by eqn (10). Although the origin of the discrepancy between the experimental results and the model prediction is not clear, we speculate that a strong damping of the resonating quartz has caused the oscillation of the quartz insensitive to the variation of η_B and ρ_B with c_s in the high c_s regime. In fact, the amplitude of the resonant peak at c_s = 2.61 M in Fig. 5 is almost close to the level of the baseline around 2 dB. To mitigate this problem, we limited our studies with c_s ≤ 2.6 M. More detailed discussion can be found in Section S1 of the SI.

3 Results and discussion

3.1 Frequency and dissipation shifts

Fig. 7(a) shows the resonant frequency f changing over time when the QCM sensitive surface is exposed to Bmim–TFSI in DMF at concentrations ranging from 0 to 2.12 M. The value of f was normalized by the corresponding fundamental frequency in air (i.e., f₀ = 10 MHz). Moreover, the obtained curves are shifted horizontally for the sake of clarity in presenting the observed results. In all the experiments, the initial phase involves the stabilization of the quartz oscillation frequency in air. For instance, following the black line in Fig. 7(a), the normalized resonant frequency remains constant, indicating that the quartz oscillation has stabilized in air. After the injection of the DMF at 500 s, the normalized resonant frequency dropped down significantly due to the increase in the inertial resistance for the oscillating quartz, followed by the stabilization of the quartz resonance frequency. When the QCM was exposed to a mixture of Bmim–TFSI and DMF at c_s = 0.655 M (red line in Fig. 7(a)), a similar frequency response was observed, but the magnitude of the normalized resonant frequency reduction was larger than that at c_s = 0 M. This effect became more evident as the Bmim–TFSI concentration increased (i.e., blue, green and yellow lines in Fig. 7(a)).

Fig. 7 The dependence of (a) the resonant frequency f and (b) the energy dissipation D on time for the mixture of Bmim–TFSI and DMF at various c_s. The value of f is normalized by the fundamental frequency f₀ employed, while that of D is normalized by the dissipation value D₀ in air after the quartz resonator is stabilized. For comparison purpose, the data for pure DMF at c_s = 0 M are also included in the plot (i.e., black solid line). All experiments are performed at room temperature (∼25 °C).

Fig. 7(b) shows the corresponding time dependence of the normalized energy dissipation, D/D₀, obtained for the mixture of Bmim–TFSI and DMF at the same IL concentrations as those used for the resonant frequency response in Fig. 7(a). Here, D₀ represents the energy dissipation D in air at a time where the quartz resonator stabilized. In contrast to the behavior of the resonant frequency, the normalized energy dissipation increased by the injection of the tested solution due to the damping of the propagating wave excited in the interfacial solid–liquid region. However, the magnitude of the increase in the normalized energy dissipation increased with increasing c_s, similar to the observed trend for the normalized resonant frequency. We observed a similar dependence of f/f₀ and D/D₀ on time for Bmim–TfO and Bmim–BF₄ as well as EG solutions in DMF (see Fig. S3 of the SI). In analysis described next, we used the values of f and D at the stabilization time as the resonant frequency f_exp(c_s) and the energy dissipation D_exp(c_s) for solutions at a given c_s.

In order to visualize the change in magnitude of f and D observed in Fig. 7, we chose the resonant frequency and dissipation at c_s = 0 M, corresponding to pure DMF, as the reference signals, and plotted the values of Δf_exp = f_exp(c_s) − f_exp(0) and ΔD_exp = D_exp(c_s) − D_exp(0) as a function of c_s in Fig. 8(a) and (b), respectively. By doing so, the effects of the compressional-wave on f_exp(c_s) can be naturally canceled out in our data analysis since Δf_exp = {f_exp(c_s) + f_comp} − {f_exp(0) + f_comp}. Here, although the frequency shift Δf_exp was negative because |f_exp(c_s)| < |f_exp(0)|, its absolute value |Δf_exp| is plotted in Fig. 8(a) for simplicity to highlight the change in the resonant frequency. We observed that both |Δf_exp| and ΔD_exp increased monotonically with increasing c_s. Moreover, the magnitude of |Δf_exp| and ΔD_exp was large for Bmim–TFSI, intermediate for Bmim–TfO, and small for Bmim–BF₄ among IL solutions tested. The mixture of EG and DMF exhibited |Δf_exp| and ΔD_exp values much smaller than those for the IL solution at high c_s.

Fig. 8 The value of (a) the resonant frequency |Δf_exp| and (b) the dissipation ΔD_exp shifts at room temperature (∼25 °C) is plotted as a function of the concentration c_s of Bmim–TFSI (black squares), Bmim–TfO (red circles), Bmim–BF₄ (blue triangles), and EG (green stars). For the Δf_exp data, its absolute value |Δf_exp| is shown for simplicity to highlight the change in the resonant frequency. The dashed line represents the predicted curve for the contribution of the bulk solution, given by eqn (10) and (11).

3.2 Estimation of the complex modulus of the EDL

To extract the contribution of the EDL to the total frequency and dissipation shifts based on the wave propagation model discussed in Section 2.6.2, we first calculated the frequency |Δf_B| and dissipation ΔD_B shifts from the contribution of the bulk solution by using eqn (10) and (11). To ensure that these equations are applicable, we verified that the tested fluids behave as Newtonian fluids at 10 MHz within the measured c_s range (see Section S2 for more details). The obtained values of |Δf_B| and ΔD_B (dashed lines) were compared with the measured |Δf_exp| and ΔD_exp values (symbols), shown in Fig. 8. Fig. 8(a) shows that the measured |Δf_exp| for IL solutions was substantially larger than |Δf_B| over the measured c_s range, while the value of |Δf_exp| for EG solutions agreed well with the model prediction of |Δf_B|. Since ethylene glycol does not contain ions but ionic liquids does, it is reasonable to attribute the observed excess frequency shift to the presence of an EDL in the tested IL solution. Accordingly, the contribution of the EDL to the total frequency shift, Δf_EDL (= f_EDL(c_s) − f_EDL(0)), was estimated by subtracting the bulk contribution from the experimentally measured total frequency shift, i.e., Δf_EDL = Δf_exp − Δf_B. The resulting values of |Δf_EDL| are plotted as a function of c_s in Fig. 9, revealing a monotonic increase of |Δf_EDL| with increasing c_s in the tested IL solutions. In contrast, |Δf_EDL| for the EG solutions remained nearly zero over the entire concentration range, confirming the absence of an EDL. On the other hand, Fig. 8(b) illustrates that the dependence of the measured ΔD_exp on c_s described well by that of ΔD_B, regardless of the type of solutes (i.e., ionic vs. non-ionic). These results indicate that the dissipation shift is dominated by the change of the density and shear viscosity of the bulk solution. According to the expression of ΔD_EDL, given by eqn (13), the observed agreement between ΔD_exp and ΔD_B suggests that the loss modulus of the EDL is much larger than the storage modulus of the EDL.

Fig. 9 The contribution of the EDL to the total frequency shift, |Δf_EDL|, in Bmim–TFSI (black squares), Bmim–TfO (red circles), and Bmim–BF₄ (blue triangles) increases monotonically with the increasing ionic liquid concentration, c_s. In contrast, |Δf_EDL| in EG solutions (green stars) remained nearly zero with respect to c_s.

Based on eqn (12), the observed increase in |Δf_EDL| suggests the variation of the density and loss modulus of the EDL with respect to the ionic liquid concentration. Therefore, we estimated the limiting loss modulus of the EDL near the charged gold surface by using eqn (12) while assuming and ρ_EDL,max = ρ_B,pureIL, where ρ_B,pureIL is the density of the pure ionic liquid. We note that under these assumptions, only is the unknown parameter in eqn (12), and thus its value can be estimated using eqn (12) with the |Δf_EDL| data. Fig. 10 shows the dependence of on c_s for Bmim–TFSI, Bmim–TfO, and Bmim–BF₄ solutions at 25 °C. In the plot, the value of was normalized by the loss modulus of the bulk solution to highlight the difference between the EDL and the bulk solution. The experimental error in is relatively large for c_s < 0.5 M. This behavior is likely due to the fact that the EDL is not yet fully developed, and its influence on the measured resonant frequency is significantly smaller than that of the bulk solution. Consequently, the relative uncertainty in |Δf_EDL| becomes comparable to the experimental noise level. For instance, for the Bmim–TFSI solution at c_s = 0.1 M, the experimental uncertainty for |Δf_EDL| in the Bmim–TFSI solution at c_s = 0.1 M was approximately 90% of its estimated value, whereas it decreased to about 1–3% at higher c_s. The value of increased rapidly with increasing c_s for c_s < 0.75 M and then leveled off at higher c_s. The limiting loss modulus of the EDL was approximately 3 orders of magnitude larger than the loss modulus of the bulk solution for c_s > 1 M, similar to the saturation value reported by using the surface force apparatus.⁶²

Fig. 10 The ratio of the limiting loss modulus of the EDL near the surface to the loss modulus of the bulk solution is plotted as a function of the IL concentration c_s for Bmim–TFSI (black squares), Bmim–TfO (red circles), and Bmim–BF₄ (blue triangles) solutions.

Strikingly, this scaling behavior was independent of anion type over the measured concentration range. This finding stands in sharp contrast to previous SFA and AFM studies under nanoconfinement, where the shear viscosity of ILs was reported to vary strongly with ion identity.^62–64 We attribute the difference to the QCM measurement geometry: the EDL forms on a single charged surface without confinement between opposing interfaces. By minimizing confinement-induced ionic structuring, QCM reveals the intrinsic viscoelastic character of the EDL, unmasking a universal concentration-dependent scaling that is insensitive to the anion structure.

4 Conclusions

Electric double layers (EDLs) are central to the performance of ionic liquid-based materials across electrochemical and tribological applications. Yet, their intrinsic viscoelastic properties have remained difficult to probe because of the nanometric thickness of the EDL and the strong coupling between electrostatics and hydrodynamics.

In this work, we employed quartz crystal microbalance (QCM) to quantify the loss modulus of EDLs formed on charged gold surfaces in Bmim–TFSI, Bmim–TfO, and Bmim–BF₄ solutions in DMF at ∼25 °C. By combining frequency–dissipation measurements with independently determined bulk properties, we isolated the EDL contribution and estimated its limiting loss modulus through our previously established wave propagation model.

Across all systems studied, the EDL exhibited a strongly viscous response, with the loss modulus increasing steeply at low IL concentrations and plateauing above ∼0.75 M. At high concentrations, the EDL loss modulus exceeded that of the bulk by roughly three orders of magnitude (∼800–1000×). Most notably, the scaling and saturation behavior was independent of anion type—a striking departure from confined SFA and AFM studies, where ion identity governs EDL viscosity. This distinction highlights QCM as a non-confining, surface-sensitive probe that suppresses confinement-induced ionic layering and instead reveals the intrinsic rheological character of the EDL.

These results highlight QCM as a low-cost, user-friendly, and quantitative tool for probing EDL rheology—complementary to more complex techniques such as the surface force apparatus—and suggest that within the measured concentration range, EDL loss moduli may be inferred from bulk rheological data. Given the promise of ILs as electrolytes and lubricants, this approach offers a practical route for guiding the design of IL-based systems with improved performance.

Author contributions

Conceptualization: A. M. and R. F.; data collection: R. Y.; data analysis: R. Y., A. M., and R. F.; funding acquisition: A. M., R. F., O. U., T. I., and A. Q. S.; writing – original draft: A. M. and R. F.; writing – review & editing: A. M., R. F., O. U., T. I., and A. Q. S. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability

The data supporting this article are available within the article and the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00877h.

Acknowledgements

The authors thank Marco Mauro from Novaetech, S.r.l. for the energy dissipation analysis using their custom QCM software. The authors acknowledge the support from the Okinawa Institute of Science and Technology Graduate University with subsidy funding from the Cabinet Office, Government of Japan. A. M. acknowledges funding from the Japanese Society for the Promotion of Science (Grants-in-Aid for Early-Career Scientists, Grant No. 21K14686). R. F. acknowledges funding from the Japanese Society for the Promotion of Science (Grants-in-Aid for Early-Career Scientists, Grant No. 20K20237).

Notes and references

1. Z. Lei, C. Dai, J. Hallett and M. Shiflett, Chem. Rev., 2024, 124, 7533–7535.
2. I. Krossing, J. M. Slattery, C. Daguenet, P. J. Dyson, A. Oleinikova and H. Weingärtner, J. Am. Chem. Soc., 2006, 128, 13427–13434.
3. Y. Mizutani, T. Watanabe, C. G. Lopez and T. Ono, Soft Matter, 2024, 20, 1611–1619.
4. C. Zhou, Y. Yu, W. Xia, S. Liu, X. Song, Z. Wu and H. Chen, Soft Matter, 2023, 19, 9460–9469.
5. L. Moralez, P. Nakasu and J. Hallett, Soft Matter, 2025, 21, 2480–2492.
6. Y. Shi, S. Yang, X. Zhang and W. Liu, Soft Matter, 2024, 20, 365–374.
7. E. D. Sultanova, A. A. Fedoseeva, A. M. Fatykhova, D. A. Mironova, S. A. Ziganshina, M. A. Ziganshin, V. G. Evtugyn, V. A. Burilov, S. E. Solovieva and I. S. Antipin, Soft Matter, 2024, 20, 7072–7082.
8. A. Matsumoto and A. Q. Shen, Soft Matter, 2022, 18, 4197–4204.
9. D. Rauber, F. Philippi, J. Becker, J. Zapp, B. Morgenstern, B. Kuttich, T. Kraus, R. Hempelmann, P. Hunt, T. Welton and C. W. M. Kay, Phys. Chem. Chem. Phys., 2023, 25, 6436–6453.
10. S. Zhang, T. Xiao, F. Fadaei Tirani, R. Scopelliti, M. K. Nazeeruddin, D. Zhu, P. J. Dyson and Z. Fei, Inorg. Chem., 2022, 61, 5010–5016.
11. D. Kaushik, P. Hitaishi, A. Kumar, D. Sen, S. M. Kamil and S. K. Ghosh, Soft Matter, 2023, 19, 5674–5683.
12. H. Montes-Campos, S. Kondrat, E. Rilo, O. Cabeza and L. M. Varela, J. Phys. Chem. C, 2020, 124, 11754–11759.
13. S. Dongare, M. Zeeshan, A. S. Aydogdu, R. Dikki, S. F. Kurtoğlu-Öztulum, O. K. Coskun, M. Muñoz, A. Banerjee, M. Gautam, R. D. Ross, J. S. Stanley, R. S. Brower, B. Muchharla, R. L. Sacci, J. M. Velázquez, B. Kumar, J. Y. Yang, C. Hahn, S. Keskin, C. G. Morales-Guio, A. Uzun, J. M. Spurgeon and B. Gurkan, Chem. Soc. Rev., 2024, 53, 8563–8631.
14. D. C. Morris, S. W. Prescott and J. B. Harper, Phys. Chem. Chem. Phys., 2021, 23, 9878–9888.
15. A. Matsumoto, R. Yoshizawa, O. Urakawa, T. Inoue and A. Q. Shen, Macromolecules, 2021, 54, 5648–5661.
16. F. Greco, D. Hemmeter, S. Shin, H.-P. Steinrück and F. Maier, J. Ionic Liq., 2022, 2, 100019.
17. F. Wang, J. Ma, D. Duan, X. Liang, K. Zhou, Y. Sun, T. Wang, G. Yang, G. Pei, H. Lin, Y. Shi, Q. Zhu, G. Li and H. Hu, Small, 2024, 20, 2307679.
18. O. J. Curnow, R. Senthooran and A. E. Somers, Phys. Chem. Chem. Phys., 2023, 25, 2401–2410.
19. P. W. Atkins, J. De Paula and J. Keeler, Atkins’physical chemistry, Oxford University Press, 2018.
20. S. Czerny-Holownia, H. R. Boyer, A. J. King, V. Y. Yang, J. Guo, M. J. Liu, J. C. Bui, W. A. Tarpeh and E. W. Lees, EES Catal., 2025, 3, 1272–1284.
21. B. Liu, W. Guo, S. R. Anderson, S. G. Johnstone, S. Wu, M. C. Herrington and M. A. Gebbie, Soft Matter, 2024, 20, 351–364.
22. Z. He, R. Xiao, S. Qu and D.-S. Wu, Soft Matter, 2025, 21, 6070–6077.
23. P. Campione, G. M. Messina, M. Giannetti, C. Mazzuca and A. Palleschi, Soft Matter, 2025, 21, 4442–4454.
24. P. Martin and E. Zussman, Soft Matter, 2024, 20, 5572–5582.
25. P. Bista, A. D. Ratschow, A. Z. Stetten, H.-J. Butt and S. A. Weber, Soft Matter, 2024, 20, 5045–5052.
26. L. Hao, M. I. Khan, Y. Lei, S. Zhou and B. Fan, Lab Chip, 2025, 25, 4161–4220.
27. N. C. X. Stuhlmüller, R. van Roij and M. Dijkstra, Soft Matter, 2025, 21, 6707–6716.
28. A. Choudhury and P. Dayal, Soft Matter, 2024, 20, 6754–6766.
29. A. Barnaveli and R. van Roij, Soft Matter, 2024, 20, 704–716.
30. S. Dhakal, S. Chandrasekar, A. D. Adedeji and S. Morozova, Soft Matter, 2025, 21, 7842–7850.
31. A. M. Bhaskaran, S. Agrawal, K. Sarkar and P. Dhar, Soft Matter, 2024, 20, 5969–5982.
32. S. Laha, D. Dhar, M. Adak, A. Bandopadhyay, S. Das, J. Chatterjee and S. Chakraborty, Soft Matter, 2024, 20, 2610–2623.
33. M. Lattuada, A. Zaccone, H. Wu and M. Morbidelli, Soft Matter, 2016, 12, 5313–5324.
34. K. K. Rangharajan and S. Prakash, Analyst, 2025, 150, 2837–2844.
35. J. K. G. Dhont and W. J. Briels, Soft Matter, 2025, 21, 4498–4516.
36. L. Jansen, T. ter Rele and M. Dijkstra, Soft Matter, 2025, 21, 8284–8299.
37. N. Kolezakis, S. Aime, R. Pastore, V. Guida, G. D’Avino and P. Edera, Soft Matter, 2025, 21, 8130–8139.
38. S. Perkin, T. Albrecht and J. Klein, Phys. Chem. Chem. Phys., 2010, 12, 1243–1247.
39. M. A. Gebbie, M. Valtiner, X. Banquy, E. T. Fox, W. A. Henderson and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 9674–9679.
40. A. M. Smith, K. R. J. Lovelock, N. N. Gosvami, P. Licence, A. Dolan, T. Welton and S. Perkin, J. Phys. Chem. Lett., 2013, 4, 378–382.
41. M. A. Gebbie, H. A. Dobbs, M. Valtiner and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 7432–7437.
42. A. M. Smith, A. A. Lee and S. Perkin, J. Phys. Chem. Lett., 2016, 7, 2157–2163.
43. A. M. Smith, A. A. Lee and S. Perkin, Phys. Rev. Lett., 2017, 118, 096002.
44. R. Atkin, S. Z. El Abedin, R. Hayes, L. H. S. Gasparotto, N. Borisenko and F. Endres, J. Phys. Chem. C, 2009, 113, 13266–13272.
45. R. Hayes, N. Borisenko, M. K. Tam, P. C. Howlett, F. Endres and R. Atkin, J. Phys. Chem. C, 2011, 115, 6855–6863.
46. H. Li, F. Endres and R. Atkin, Phys. Chem. Chem. Phys., 2013, 15, 14624–14633.
47. H.-W. Cheng, J.-N. Dienemann, P. Stock, C. Merola, Y.-J. Chen and M. Valtiner, Sci. Rep., 2016, 6, 1–9.
48. T. Cui, A. Lahiri, T. Carstens, N. Borisenko, G. Pulletikurthi, C. Kuhl and F. Endres, J. Phys. Chem. C, 2016, 120, 9341–9349.
49. M. Z. Bazant, B. D. Storey and A. A. Kornyshev, Phys. Rev. Lett., 2011, 106, 046102.
50. N. Gavish, D. Elad and A. Yochelis, J. Phys. Chem. Lett., 2018, 9, 36–42.
51. A. A. Lee, C. S. Perez-Martinez, A. M. Smith and S. Perkin, Faraday Discuss., 2017, 199, 239–259.
52. P. Gaddam and W. Ducker, Langmuir, 2019, 35, 5719–5727.
53. F. Rondepierre, P.-F. Brevet and J. Duboisset, J. Phys. Chem. Lett., 2025, 16, 2690–2694.
54. S. W. Coles, C. Park, R. Nikam, M. Kanduč, J. Dzubiella and B. Rotenberg, J. Phys. Chem. B, 2020, 124, 1778–1786.
55. N. Anousheh, F. J. Solis and V. Jadhao, AIP Adv., 2020, 10, 125312.
56. R. M. Adar, S. A. Safran, H. Diamant and D. Andelman, Phys. Rev. E, 2019, 100, 042615.
57. B. Rotenberg, O. Bernard and J.-P. Hansen, J. Phys.: Condens. Matter, 2018, 30, 054005.
58. A. Härtel, M. Bültmann and F. Coupette, Phys. Rev. Lett., 2023, 130, 108202.
59. D. Ribar, C. E. Woodward and J. Forsman, Soft Matter, 2025, 21, 5562–5572.
60. S. Wang, H. Tao, J. Yang, J. Cheng, H. Liu and C. Lian, J. Phys. Chem. Lett., 2024, 15, 7147–7153.
61. J. Yang, S. Kondrat, C. Lian, H. Liu, A. Schlaich and C. Holm, Phys. Rev. Lett., 2023, 131, 118201.
62. K. Ueno, M. Kasuya, M. Watanabe, M. Mizukami and K. Kurihara, Phys. Chem. Chem. Phys., 2010, 12, 4066–4071.
63. T. Kamijo, H. Arafune, T. Morinaga, S. Honma, T. Sato, M. Hino, M. Mizukami and K. Kurihara, Langmuir, 2015, 31, 13265–13270.
64. M. Han and R. M. Espinosa-Marzal, J. Phys. Chem. C, 2018, 122, 21344–21355.
65. B. Wu, J. P. Breen, X. Xing and M. D. Fayer, J. Am. Chem. Soc., 2020, 142, 9482–9492.
66. G. Krämer and R. Bennewitz, J. Phys. Chem. C, 2019, 123, 28284–28290.
67. S. Park and J. G. McDaniel, J. Chem. Phys., 2020, 152, 074709.
68. F. Federici Canova, H. Matsubara, M. Mizukami, K. Kurihara and A. L. Shluger, Phys. Chem. Chem. Phys., 2014, 16, 8247–8256.
69. K. Tomita, M. Mizukami, S. Nakano, N. Ohta, N. Yagi and K. Kurihara, Phys. Chem. Chem. Phys., 2018, 20, 13714–13721.
70. D. Johannsmann, A. Langhoff and C. Leppin, Sensors, 2021, 21, 3490.
71. C. Leppin, A. Peschel, F. S. Meyer, A. Langhoff and D. Johannsmann, Analyst, 2021, 146, 2160–2171.
72. R. Etchenique and T. Buhse, Analyst, 2002, 127, 1347–1352.
73. R. Funari, R. Ripa, B. Söderström, U. Skoglund and A. Q. Shen, ACS Sens., 2019, 4, 3023–3033.
74. K. Yamaoka, Y. Fujii and N. Torikai, Langmuir, 2024, 40, 17405–17409.
75. R. Ripa, A. Q. Shen and R. Funari, ACS Omega, 2020, 5, 2295–2302.
76. A. Monemian Esfahani, W. Zhao, J. Y. Chen, C. Huang, N. Xi, J. Xi and R. Yang, Anal. Chem., 2018, 90, 10340–10349.
77. R. Funari, A. Matsumoto, J. R. de Bruyn and A. Q. Shen, Anal. Chem., 2020, 92, 8244–8253.
78. W. M. Haynes, CRC handbook of chemistry and physics, CRC Press, 2014.
79. F. Corradini, G. Franchini, A. Marchetti, M. Tagliazucchi and L. Tassi, Bull. Chem. Soc. Jpn., 1995, 68, 1867–1872.
80. Zhang and Han, J. Chem. Eng. Data, 1996, 41, 516–520.
81. B. R. Breslau and I. F. Miller, J. Phys. Chem., 1970, 74, 1056–1061.
82. C. W. Macosko, Rheology Principles, Measurements, and Ppplications, Wiley-VCH, 1994.
83. M.-M. Huang, Y. Jiang, P. Sasisanker, G. W. Driver and H. Weingärtner, J. Chem. Eng. Data, 2011, 56, 1494–1499.
84. D. Johannsmann, I. Reviakine, E. Rojas and M. Gallego, Anal. Chem., 2008, 80, 8891–8899.
85. A. D. Easley, T. Ma, C. I. Eneh, J. Yun, R. M. Thakur and J. L. Lutkenhaus, J. Polym. Sci., 2022, 60, 1090–1107.
86. S. Donatan, M. Sarikaya, C. Tamerler and M. Urgen, J. R. Soc., Interface, 2012, 9, 2688–2695.
87. A. A. Lee, C. S. Perez-Martinez, A. M. Smith and S. Perkin, Phys. Rev. Lett., 2017, 119, 026002.
88. M. V. Voinova, M. Rodahl, M. Jonson and B. Kasemo, Phys. Scr., 1999, 59, 391.
89. D. Johannsmann, The Quartz Crystal Microbalance in Soft Matter Research: Fundamentals and Modeling, Springer International Publishing, Cham, 2015, pp. 191–204.

---