Introduction

Mobility at Low Field Strength

The movement of ions through a collision gas at elevated pressures induced by an electrical field is of great interest in mass spectrometry and related areas: Not only is ion mobility important in modern atmospheric pressure ionization sources [1,2,3], but in ion mobility spectrometry (IMS), it is also used to separate ions according to their size prior to mass analysis [4,5,6]. The mobility of ions has also been subject to theoretical investigations [7,8,9,10] since it offers insights into fundamental physics of collisions and ion-neutral interactions, e.g., kinetic theory in the context of the Boltzmann transport equation [11].

The acceleration of the ion ensemble through the electrical field is countered by collisions with the background gas, leading to a constant drift velocity vD proportional to the applied field strength, E:

$$ {v}_{\mathrm{D}}=K(E)\cdot E $$
(1)

The proportionality constant K is called the ion mobility and, for low field strengthsFootnote 1 (as applied in IMS), can be regarded as constant. In this zero-field limit, the mobility can accurately be described by the Mason-Schamp equation [7]:

$$ K(0)=\frac{3}{16}{\left(\frac{2\pi }{\mu {k}_{\mathrm{B}}T}\right)}^{1/2}\frac{ze}{N\Omega (T)} $$
(2)

Here, μ is the reduced mass of the ion-neutral pair, kB is the Boltzmann constant, T the gas temperature, ze the ions charge, N is the neutral particle density, and Ω is the collision cross section (CCS) of the ion in the particular collision gas. This equation is used, for example, to determine the CCS from mobility measurements, opening the opportunity to identify compounds when compared to theoretical determination of the CCS (see, for example, [13,14,15,16]). The latter can be performed by numerically solving the following equation [17, 18]:

$$ {\displaystyle \begin{array}{l}{\Omega}^{\left(1,1\right)}(T)=\frac{1}{8{\pi}^2}\underset{0}{\overset{2\pi }{\int }}\mathrm{d}\theta \underset{0}{\overset{\pi }{\int }}\sin \varphi \kern.5em \mathrm{d}\varphi \underset{0}{\overset{2\pi }{\int }}\mathrm{d}\gamma \frac{\pi }{8}{\left(\frac{\mu }{k_{\mathrm{B}}T}\right)}^3\\ {}{\int}_0^{\infty}\exp \kern.5em \left(-\frac{\mu {v}_r^2}{2{k}_{\mathrm{B}}T}\right){v}_r^5\mathrm{d}{v}_r{\int}_0^{\infty }2b\left[1-\cos \kern.5em \left(\chi \left(\theta, \varphi, \gamma, {v}_r,b\right)\right)\right]\mathrm{d}b\end{array}} $$
(3)

Here, the angles θ, φ, and γ describe the orientation of the impact between ion and collision gas; b is the impact parameter; vr the relative velocity; and χ the resulting scattering angle, which depends on all these parameters and the interaction potential between the collision partners. We note that this equation provides only a first order approximation to the collision cross section (hence the (1,1) superscript), but it is commonly used and yields sufficiently accurate results [17].

Mobility at Higher Field Strengths

If the field strength is further increased, the gas-phase ion mobility starts to change [7]. From a theoretical standpoint, the change of the mobility becomes significant when the velocity distribution of the ion ensemble starts to differ from the velocity distribution of the collision gas, i.e., when the acceleration through the electrical field significantly increases the collision frequency over the thermal collision rate [12]. One popular approach to describe the two different velocity distributions is the so-called two-temperature theory [19,20,21] in which the velocity distribution of the ions (in the direction of the electric field) is described by a different temperature than the distribution of the collision gas. While still assumed to be a 1-D Maxwellian distribution, the higher temperature, Tion, broadens the distribution, allowing for higher ion velocities. The collisions between ions and neutrals are then determined by the relative velocity distribution, which depends on the reduced mass of the ion-neutral pair and an effective temperature, Teff, which can be expressed as follows [19,20,21]:

$$ {T}_{\mathrm{eff}}={T}_{\mathrm{bath}}+\frac{m_{\mathrm{bath}}{v_{\mathrm{D}}}^2}{3{k}_{\mathrm{B}}} $$
(4)

Here, mbath is the mass of the collision gas. Although this formula is used for any type of potential with good accuracy [7 p. 154], it is only exact when a repulsive r−4 potential between ion and neutral is assumed [22].

Applying this effective temperature instead of the background gas temperature to the Mason-Schamp equation (Eq. (2)) greatly improves the accuracy of the equation at higher field strengths [22]. It can also be applied to the theoretical description of the CCS by calculating Ω(1,1) at Teff instead of T. With that, the Boltzmann weighting of the relative velocities vr (exponential in fourth integral of Eq. (3)) then accounts for the correct shape of the distribution by applying not only the reduced mass but also the effective temperature. If the ions are polyatomic, the high impact collisions will also deposit energy into the internal degrees of freedom of the molecule. While there is some theoretical work on the influence of inelastic collision on the mobility [7 chapter 6-6], their contributions are often ignored for simplicity reasons.

Differential Mobility Spectrometry

The electric field dependency of the mobility is exploited in differential mobility spectrometry (DMS) [23,24,25,26,27], also known as field-asymmetric waveform ion mobility spectrometry (FAIMS). Briefly, the ions are introduced to a small gap between two parallel plate electrodes at atmospheric pressure. On the plates, an asymmetric waveform is applied to produce an oscillating electric field orthogonal to the gas flow stream. The so-called separation field (peak-to-peak voltage in the order of 1–4 kV) creates a high-field phase in one and a low-field phase in the opposite direction such that the integral over one waveform cycle is zero [28]. Ideally, this would be a rectangular waveform, where the low field portion has 1/nth amplitude and lasts n times as long. However, in practice, other, more practical waveforms are applied, e.g., a two-harmonics or double sine waveform (see below). Because of the field dependency of the ion mobility, the overall displacement towards the electrodes is not canceled and ions are dragged towards one of the plates, eventually being neutralized and precluding their detection. This drift is compensated by a constant voltage (called compensation voltage; 0–100 V in magnitude) applied on the electrodes. By scanning this compensation voltage, different ions are able to pass through the DMS cell. Thus, while in traditional IMS, ions are mainly separated by their CCS; in DMS, the difference in high- and low-field mobility is the characteristic which separates the ions.

Adding solvent vapor in low concentrations to the gas phase often increases the differential mobility effect and thus the separation capability [25]. Three main trends are observed in experiments [25, 29]: the ions mobility appear larger in high- than in low-field (type A), larger in the low-field than in the high-field (type C), and an intermediate case, where the high-field mobility is larger at first, but eventually, the low-field mobility exceeds that of the high-field condition (type B). These types can be easily identified by varying the separation field amplitude and recording the compensation voltage required for optimal transmission, yielding a so-called dispersion plot. It should be noted that these types are not intrinsically different. For example, any type A case would become a type B case for high enough field strengths since the low-field mobility eventually overcomes the high-field one (see below). Also, type C can be viewed as limiting case of type B, where the initial excess of the high-field mobility over the low-field mobility vanishes.

Besides the theoretical improvements with regard to the mobility at higher field strengths mentioned above, no thorough model yet exists to explain differential mobility under various experimental conditions. The effect of solvent vapor in particular is still a matter of debate, while its effect through the change of the interaction potential and subsequent change in CCS has been discussed [30], other literature reports suggest a dynamic clustering/declustering mechanism [30,31,32]. In the high field, where the ion internal energy increases due to collisional heating, firmly bound solvent molecules will boil off, decreasing the overall average CCS and thus increasing the ensemble mobility (as observed for type A ions). Even in IMS, it has be shown that (dynamic) clustering has to be considered when modeling the reactant ion peak (RIP), e.g., the proton bound water cluster system [H+(H2O)n]+ [33].

In this contribution, we present a thorough model for the calculation of an ion’s mobility under different conditions and over a wide range of field strengths, incorporating two-temperature theory, the temperature dependence of the CCS, and the effect of ion-solvent clustering. To test our model, we simulate dispersion plots of tetramethylammonium (Me4N+) in a pure nitrogen atmosphere and seeded with methanol (MeOH), acetonitrile (ACN), and acetone (ACE). Experimental DMS data for those conditions are available in the literature [34] and show a range of different behavior, i.e., from type C, very weak type B (no modifier added), strong type B (MeOH added), and type A (ACN, ACE) behavior. We reproduce these data by applying the same formalism for each system studied and having the results solely determined by the nature of the system (i.e., not through fitting to experimental data or introduction of system specific parameters). Our model strongly suggests the importance of clustering/declustering reactions for differential mobility and yields the first thorough method for the estimation of trends in of dispersion plots from first principles.

Computational Methods

Ab Initio Calculations

To explore the large configurational space accessible to larger cluster sizes (different binding sites/configurations of the solvents to the ion), we performed basin hopping (BH) simulations to find candidate structures for further optimization. Details of the basin hopping algorithm can be found elsewhere [35, 36]. Briefly, to generate random input structures, different moieties (ion and n solvents) were each rotated by a random angle of −10° ≤ ϕ ≤ 10° around their body-fixed x, y, or z axis and randomly translated along the overall x, y, or z axis by \( -0.5{\mathrm{A}}^{\circ}\le l\le 0.5{\mathrm{A}}^{\circ } \). If dihedral angles were suspected to influence the potential energy (e.g., methyl group rotation in MeOH compared to a fixed ion-HO framework), they were also randomly altered with an angle −5° ≤ φ ≤ 5°. These random structures were optimized using the AMBER force field [37] accompanied by partial charges calculated with the Merz-Singh-Kollman scheme [38, 39] for a “guess structure” optimized at the B3LYP-GD3/6-31++G(d,p) [40,41,42,43] level of theory. This procedure yielded up to 20,000 structures per cluster. These structures were subsequently combined into energy bins of a certain size, and the lowest energy structure in each bin was used for geometry comparison. Describing the geometry of a molecule by a vector of mass-weighted distances to the center of mass, two geometries can be compared by cosine similarity, i.e., measuring the angle between these two vectors. If the angle was smaller than a predefined threshold, the structures were categorized as being identical. The BH algorithm would typically find 10 to 50 unique candidate structures.

Each of the unique candidates was then re-optimized at the B3LYP-GD3/6-31++G(d,p) level of theory, as suggested for the MobCal-MPI code [44] (see below). Many candidate structures identified by molecular mechanics yielded the same structure following density functional theory (DFT) treatment. DFT optimizations typically produced 1–4 structures for which ESP charges [38, 39] and frequency calculations were also performed to obtain input for partition function calculations (see below). Additional single-point energy calculations were performed at the B2PLYP-GD3BJ/def2-TZVPP [45,46,47] level of theory for a better description of the electronic energy. Each of those structures is defined by the number of solvent molecules n = 0, 1, 2, … and an index to name different conformers α = a, b, c, …, e.g., [Me4N+(ACN)2]b+.

For some test cases, molecular dynamic (MD) simulations were performed to test the influence of vibrational broadening (due to elevated ion temperatures) on the CCS. Because simple force fields neglect the anharmonic nature of loosely bound clusters, we conducted ab initio MD simulations with the atom-centered density matrix propagation (ADMP) [48,49,50] model using, again, the B3LYP-GD3/6-31++G(d,p) level of theory. This version of a Car-Parinello [51] density matrix propagation is much faster than to actually solve the SCF equations at each time step, allowing us to perform 5000 × 0.2 fs trajectories on a full quantum mechanical (QM) potential energy surface (PES). To model the temperature effect, a kinetic energy of (3N − 6)kBTion, with N being the number of atoms, kB the Boltzmann constant and Tion the ions temperature, was applied, randomly distributed over the nuclei. This yields a mean kinetic energy of \( \frac{1}{2}{k}_{\mathrm{B}}{T}_{\mathrm{ion}} \) per vibrational degree of freedom, as should be the case for a thermalized multidimensional harmonic oscillator.

All ab initio and molecular mechanics calculations were conducted with the Gaussian16 program package [52].

CCS Calculations

Collision cross section (CCS) calculations (Eq. (3)) were performed with the recently developed MobCal-MPI code [44] using the trajectory method [17, 18], which allows for fast calculations of the CCS through parallelization of the different trajectories. We slightly modified the code to be able to adjust the temperature which is used for the Boltzmann weighting of the relative velocities (see Eq. (3)) and set it to the effective temperature studied (Eq. (4)). All calculations were performed in molecular nitrogen and using 10 cycles, 48 velocity parameter integrations, and 512 impact parameter integrations as suggested by the authors [44]. We used either the equilibrium geometry or 100 random samples from the MD simulations (conducted at the corresponding temperature), all accompanied by the determined ESP charges and MMFF94 van der Waals parameters [53, 54], to calculate the CCS over a range of effective temperatures. These data points were subsequently fitted to a simple function of the form a Teffb + c. We have no complete physical meaning to this fitting function. It is stated in the literature [7, 30] that the CCS is proportional to Teff−1/2 for low temperatures/field strengths and an ion/neutral potential proportional to r−4. However, the corresponding fitting function a Teff−0.5 + c showed systematic deviation from the data, probably due to elevated (effective) temperatures. Thus, we chose to also adjust the exponent for a more accurate interpolation. In terms of physical meaning, the offset c can be seen as CCS in the infinite temperature limit, when all long-range interactions become negligible and only the hard-sphere core is contributing to the CCS [55].

Mobility Calculations

The individual drift velocity vD of each structure (n, α) considered is calculated with Eq. (1) using the Mason-Schamp expression for the mobility, K (Eq. (2)); the applied field strength, E; and two-temperature theory, i.e., setting T → Teff (Eq. (4)). Since the effective temperature is not known a priori, we start by using the background gas temperature, Tbath, use the obtained drift velocity to calculate the effective temperature with Eq. (4) and, again, calculate the drift velocity with Eqs. (1) and (2). This is done iteratively until the effective temperature obtained converges to a predefined threshold of 10−4 K. In this self-consistent calculation, Ω(T) is also updated at every iteration step to the effective temperature applied using the analytical fitting function determined beforehand. Having an analytical function for Ω(T) dramatically speeds up the calculations since otherwise the CCS would have to be computed at every step of this self-consistent calculation explicitly using MobCal-MPI. As a consequence of the different CCS values (viz. their fitting functions), each cluster has a different effective temperature and a different drift velocity.

Boltzmann Weighting

For simplicity, we assume complete equilibration of the system at any given field strength. A detailed discussion of this assumption is given below. To obtain the population distribution of a thermalized system with different local minima, e.g., different conformers, the quantum-harmonic superposition approximation (QHSA) is typically used [56,57,58]. In this formalism, the density of states in each minimum is assumed to be independent from each other and, thus, the relative population Pi of minimum i can be calculated from the vibrational partition function Qvib,i and the zero-point corrected electronic energy ε0,i, both readily available from ab initio calculations:

$$ {P}_i=\frac{Z_i}{\sum_j{Z}_j}\kern2em \mathrm{with}\kern1.50em {Z}_i={Q}_{\mathrm{vib},i}\exp \kern.5em \left(-\frac{\varepsilon_{0,i}-{\varepsilon}_{0,\mathrm{ref}}}{k_{\mathrm{B}}T}\right) $$
(5)

where the sum in the denominator runs over all minima considered.

We have to modify this description in two ways. First, we are considering dissociation equilibria (viz. solvent molecules evaporating from the ion) and thus the total partition function has to be used (including rotational and translational contributions). Also, the concentration of the solvent [M] has to be considered, because it determines the collision frequency between ions and solvent molecules which, in turn, determines the equilibrium constant. In particular, if the solvent concentration is 0, the population of any cluster has to be 0 as well, independent from any energy differences. This leads to

$$ {Z}_i={\left(\frac{\left[M\right]}{N}\right)}^{n_i}\kern.1em {Q}_{\mathrm{tot},i}\cdot {Q_{\mathrm{tot},M}}^{\overset{\sim }{n}-{n}_i}\cdot \exp \kern.1em \left(-\frac{\varepsilon_{0,i}-{\varepsilon}_{0,\mathrm{ref}}}{k_{\mathrm{B}}T}\right) $$
(6)

where N is the total particle density and ni the number of solvent molecules firmly attached to the ion in structure i. The partition function for the system is thus a product of the (total) partition function of the ionic cluster i and the (total) partition function of each unbound solvent, with \( \overset{\sim }{n} \) being the largest cluster size considered. For example, if \( \overset{\sim }{n}=1 \), the partition function of the bound cluster is just Qtot,i whereas the partition function of the dissociated system is Qtot,i−1 · Qtot,M. This ensures that the partition functions are being calculated in the same configurational space, treating the bound and dissociated states as two “local minima” on the same PES. The contributions of the additional translational degrees of freedom in QM under dissociation are describing an increase in entropy when forming two particles out of one. A similar, although not as thorough, formalism was already used in the literature for the determination of field dependent mobilities under consideration of clustering [30].

The second modification is necessary due to the different ion temperatures of each structure considered, assuming that the internal degrees of freedom are also heated and thus also have a temperature of Tion. Initially, we want to deduce the appropriate formula for a simple system containing two local minima, A and B, as shown in Figure 1, and then move on to the generalization of multiple minima. For simplicity, we drop the index on the energy (ε0 → ε). The total partition function of the system can then be written as an integral over the total density of states ρ(ε) and the Boltzmann weighting function f(ε):

$$ Z=\underset{0}{\overset{\infty }{\int }}\rho \left(\varepsilon \right)f\left(\varepsilon \right)\kern.1em \mathrm{d}\varepsilon $$
(7)
Figure 1
figure 1

Quantum harmonic superposition approximation (QHSA). The population density in minimum a is according to exp(−ε/kTA) (red) and minimum B is exp(−(ε − εB)/kTB) (green) scaled by the constant value of the exponential of minimum a at the energy of minimum B, i.e., exp(−εB/kTA). See text for details

The normal superposition approximation now states:

$$ Z=\underset{\varepsilon_{\mathrm{A}}=0}{\overset{\infty }{\int }}{\rho}_{\mathrm{A}}\left(\varepsilon \right)f\left(\varepsilon \right)\kern.1em \mathrm{d}\varepsilon +\underset{\varepsilon_{\mathrm{B}}}{\overset{\infty }{\int }}{\rho}_{\mathrm{B}}\left(\varepsilon \right)f\left(\varepsilon \right)\kern.1em \mathrm{d}\varepsilon ={Z}_{\mathrm{A}}+{Z}_{\mathrm{B}} $$
(8)

While the first term, ZA, is just the partition function QA of minimum A, the second term can be written in a different way:

$$ \underset{\varepsilon_{\mathrm{B}}}{\overset{\infty }{\int }}{\rho}_{\mathrm{B}}\left(\varepsilon \right)\exp \kern.3em \left(-\frac{\varepsilon }{k_{\mathrm{B}}T}\right)\kern.3em \mathrm{d}\varepsilon =\underset{\varepsilon_{\mathrm{B}}}{\overset{\infty }{\int }}{\rho}_{\mathrm{B}}\left(\varepsilon \right)\exp \kern.3em \left(-\frac{\varepsilon -{\varepsilon}_{\mathrm{B}}}{k_{\mathrm{B}}T}\right)\exp \kern.3em \left(-\frac{\varepsilon_{\mathrm{B}}}{k_{\mathrm{B}}T}\right)\kern.3em \mathrm{d}\varepsilon =\exp \kern.3em \left(-\frac{\varepsilon_{\mathrm{B}}}{k_{\mathrm{B}}T}\right)\cdot \underset{\varepsilon_B}{\overset{\infty }{\int }}{\rho}_{\mathrm{B}}\left(\varepsilon \right)\exp \kern.3em \left(-\frac{\varepsilon -{\varepsilon}_{\mathrm{B}}}{k_{\mathrm{B}}T}\right)\mathrm{d}\varepsilon $$
(9)

If we now set \( \varepsilon -{\varepsilon}_{\mathrm{B}}=\overset{\sim }{\varepsilon } \) and recognize that ρB(ε) = 0 for ε < εB, the integral becomes the partition function of B:

$$ \exp \left(-\frac{\varepsilon_{\mathrm{B}}}{k_{\mathrm{B}}T}\right)\cdot \underset{0}{\overset{\infty }{\int }}{\rho}_{\mathrm{B}}\left(\overset{\sim }{\varepsilon}\right)\exp \kern.1em \left(-\frac{\overset{\sim }{\varepsilon }}{k_{\mathrm{B}}T}\right)\mathrm{d}\overset{\sim }{\varepsilon }={Q}_{\mathrm{B}}\cdot \exp \kern.1em \left(-\frac{\varepsilon_{\mathrm{B}}}{k_{\mathrm{B}}T}\right) $$
(10)

In total, this yields the following:

$$ Z={Z}_{\mathrm{A}}+{Z}_{\mathrm{B}}={Q}_{\mathrm{A}}\cdot \underset{=1,\mathrm{since}\kern.1em {\varepsilon}_{\mathrm{A}}=0}{\underbrace{\exp \kern.1em \left(-\frac{\varepsilon_{\mathrm{A}}}{k_{\mathrm{B}}T}\right)}}+{Q}_{\mathrm{B}}\cdot \exp \kern.1em \left(-\frac{\varepsilon_{\mathrm{B}}}{k_{\mathrm{B}}T}\right) $$
(11)

which is the normal superposition approximation described in Eq. (5).

If we consider Figure 1, it can be seen that the constant \( \exp \left(-\frac{\varepsilon_{\mathrm{B}}}{k_{\mathrm{B}}T}\right) \) is a scaling factor for the population density function for minimum B with respect to minimum A because the population density of minimum A declines to this value at that particular energy. For the case of different temperatures in minima A and B, their Boltzmann functions will become \( {f}_{\mathrm{A}}\left(\varepsilon \right)=\exp \left(-\frac{\varepsilon }{k_{\mathrm{B}}{T}_{\mathrm{A}}}\right) \) and \( {f}_{\mathrm{B}}\left(\varepsilon \right)=\exp \left(-\frac{\varepsilon -{\varepsilon}_{\mathrm{B}}}{k_{\mathrm{B}}{T}_{\mathrm{B}}}\right) \). However, one must still scale the population distribution in minimum B with the exponential of minimum A since that describes the overall population density at the given energy. Hence, we set

$$ Z=\underset{0}{\overset{\infty }{\int }}{\rho}_{\mathrm{A}}\left(\varepsilon \right){f}_{\mathrm{A}}\left(\varepsilon \right)\kern.5em \mathrm{d}\varepsilon +\exp \kern.3em \left(-\frac{\varepsilon_{\mathrm{B}}}{k_{\mathrm{B}}{T}_{\mathrm{A}}}\right)\cdot \underset{\varepsilon_{\mathrm{B}}}{\overset{\infty }{\int }}{\rho}_{\mathrm{B}}\left(\varepsilon \right){f}_{\mathrm{B}}\left(\varepsilon \right)\kern.1em \mathrm{d}\varepsilon ={Q}_{\mathrm{A}}\left({T}_{\mathrm{A}}\right)+{Q}_{\mathrm{B}}\left({T}_{\mathrm{B}}\right)\cdot \exp \kern.1em \left(-\frac{\varepsilon_{\mathrm{B}}}{k_{\mathrm{B}}{T}_{\mathrm{A}}}\right) $$
(12)

Thus, the partition function of each minimum is calculated at their respective temperature; however, the scaling factor of the partition function in minimum B is determined by the exponential decay of the population density in the next lower energy minimum. For any system with multiple, energy sorted minima (i.e., ε0,i+1 > ε0,i for all i), the scaling factors of all lower energy minima have to be taken into account. With the lowest energy set to 0 (ε0,1 = 0), it follows that:

$$ {Z}_1={Q}_1\left({T}_1\right) $$
(13a)
$$ {Z}_i={Q}_i\left({T}_i\right)\cdot \prod \limits_{j=2}^i\exp \kern.1em \left(-\frac{\varepsilon_{0,j}-{\varepsilon}_{0,j-1}}{k_{\mathrm{B}}{T}_{j-1}}\right)={Q}_i\left({T}_i\right)\cdot \exp \kern.1em \left(-\sum \limits_{j=2}^i\frac{\varepsilon_{0,j}-{\varepsilon}_{0,j-1}}{k_{\mathrm{B}}{T}_{j-1}}\right)\kern2em \mathrm{for}\kern1em i\ge 2 $$
(13b)

It is worth mentioning that this weighting reduces back to the standard superposition approximation if all Ti are equal.

Assuming that the dissociated solvent molecules are not heated by the electric field, their partition function is calculated at the background gas temperature. Hence, we set Qi · QM → Qi(Ti) · QM(Tbath). Including the concentration dependency, this leads to

$$ {Z}_1={\left(\frac{\left[M\right]}{N}\right)}^{n_1}{Q}_1\left({T}_1\right){Q}_M{(T)}^{\overset{\sim }{n}-{n}_1} $$
(14a)
$$ {Z}_i={\left(\frac{\left[M\right]}{N}\right)}^{n_i}{Q}_i\left({T}_i\right){Q}_M{(T)}^{\overset{\sim }{n}-{n}_i}\cdot \exp \kern.1em \left(-\sum \limits_{j=2}^i\frac{\varepsilon_{0,j}-{\varepsilon}_{0,j-1}}{k_{\mathrm{B}}{T}_{j-1}}\right)\kern0.75em \mathrm{for}\kern0.5em i\ge 2 $$
(14b)

We calculate all partition functions using standard formula (see SI S1) since they have a relatively simple functional dependence of the temperature, and all other input data (e.g., mass, rotational constants, normal frequencies; read in from the ab initio output) have to be determined only once, independent from the temperature.

For the determination of the ensemble mobility, we just weight the individual mobilities by the now accessible relative populations:

$$ {\left\langle K\right\rangle}_{\mathrm{ens}}=\sum \limits_i{K}_i{P}_i $$
(15)

We want to point out that this approach assumes fast thermal equilibrium compared to any change in the field strength. This means that we assume an ionic cluster, as soon as it adopts a new structure through dissociation or change of conformation, equilibrates to the new effective temperature determined by its CCS instantly. For fast varying fields (compared to collision frequencies), this may not hold true and kinetic simulations with RRKM theory determined rate constants would be needed to accurately determine the relative population of each cluster structure. This would also circumvent the rather strange assumption of thermal equilibrium between systems having different temperatures, which clearly is in contrast to the zeroth law of thermodynamics. One would then calculate the rate constants for association A+ + M → [A + M]+ with the collision frequency applying Teff of A+, since this temperature describes the relative velocity distribution between ion and neutral. The rate constant for dissociation \( {\left[A+M\right]}_a^{+}\to {A}^{+}+M \) or rearrangement \( {\left[A+M\right]}_a^{+}\to {\left[A+M\right]}_b^{+} \) would be calculated applying Tion of \( {\left[A+M\right]}_a^{+} \), since this temperature determines the internal energy and thus state density at the dissociation or transition state energy. This approach will be subject to future work.

Dispersion Plot Calculation

With the described formalism, the ensemble mobility of a cluster-solvent system (Eq. (15)) can be calculated for any given field strength E, which is introduced in the Mason-Schamp equation for the self-consistent determination of the drift velocity and effective temperature. By varying the field strength E, the field dependency of the (ensemble) mobility can be calculated, typically described by the α function, defined as

$$ \alpha (E)=\frac{{\left\langle K(E)\right\rangle}_{\mathrm{ens}}}{{\left\langle K(0)\right\rangle}_{\mathrm{ens}}}-1 $$
(16)

Just to repeat, at any given field strength, the CCS, the mobility, the partition function, and the Boltzmann weight for all structures are re-calculated at the effective temperature determined at this field strength. In case of the zero-field mobility 〈K(0)〉ens, this reduces to the background gas temperature.

The α function and its derivative with respect to the field strength, α′ = dα/dE, can now be used to calculate the compensation voltage (CV) for a known separation field. Given the peak-to-peak separation voltage SVpp and using a double sine waveform with maximum amplitude D = 2/3 · SVpp [26]:

$$ E(t)=\frac{D}{d}\cdot \left(\frac{2}{3}\sin \kern.1em \left(\omega t\right)+\frac{1}{3}\sin \kern.1em \left(2\omega t-\frac{\pi }{2}\right)\right) $$
(17)

in which d is the gap size in the DMS cell, the CV can be calculated according to [23, 30]:

$$ CV=-\frac{{\left\langle \alpha \left(E(t)\right)\cdot E(t)\right\rangle}_{\mathrm{wf}}}{1+{\left\langle \alpha \left(E(t)\right)\right\rangle}_{\mathrm{wf}}+{\left\langle {\alpha}^{\prime}\left(E(t)\right)\cdot E(t)\right\rangle}_{\mathrm{wf}}}\cdot d $$
(18)

In this equation, the averages 〈⋯〉wf are taken over one cycle of the waveform. Doing this for multiple values of SVpp, a full dispersion plot can be calculated from first principles. We want to point out that through our assumption of fast equilibration and thus “equilibrium” α functions, the actual waveform frequency, ω, does not appear in our calculations. Our results thus correspond to the limiting case where ω is small compared to all chemical processes. Although the calculated α curves could be compared to experimentally determined α functions (available from experimental dispersion plots through a least square fit procedure also relying on Eq. (18) [23, 25]), we wanted to provide a method to produce dispersion plots since they are more readily available from experiment and thus easier for comparison.

Results and Discussion

We applied the above formalism to tetramethylammonium (Me4N+) in gaseous molecular nitrogen seeded with either methanol (MeOH), acetonitrile (ACN), or acetone (ACE) vapor. These systems are a good test of theory, because experimental DMS data is readily available [34]; the measured dispersion plots show very different behavior (type A and type B, depending on the solvent-modified environment), and the system size is relatively small, allowing for faster computation.

Cluster Structures

Figure 2 shows all structures n, α for the [Me4N+(ACN)n]+ ion-solvent cluster systems found by the BH search outlined above. Structures for the ACE and MeOH systems, as well as Cartesian coordinates and energies (electronic and thermochemical), are found in the Supporting Information (S2 and Figs. S2.1 and S2.2). For n = 1, only one unique conformer was found, i.e., where the solvent binds to the triangular face of the tetramethyl ammonium ion. Multiple structures were found for n = 2 where the two solvents either bind to two different faces of the tetrahedron (α = a structures) or interact with each other. For example, MeOH shows a structure where the two solvents are forming a hydrogen bond chain ([Me4N+(MeOH)2]c+). Since we found [Me4N+(ACN)2]+ to be highly populated for low field strengths (see below, especially Figure 5), we also added the lowest-energy conformer n = 3 cluster to the formalism. Hence, for all solvents, we ensured that the highest cluster numbers considered are only slightly populated at low field strengths.

Figure 2
figure 2

Geometries of all [Me4N+(ACN)n]+ clusters optimized at the B3LYP-GD3/6-31++G(d,p) level of theory. Zero-point corrected binding energies, Δε0, are given in electron volt relative to the fully dissociated state. Solvents prefer interacting via the triangular faces of the tetramethyl ammonium ion

Taking a look at the calculated zero-point energy corrected binding energies Δε0, it becomes clear that MeOH is a rather weak binding solvent (Δε0 = 0.45 eV), whereas ACN and ACE interactions with Me4N+ are much stronger (Δε0 = 0.63 eV and Δε0 = 0.64 eV, respectively). This drives the cluster stability and the field dependency of the cluster size distribution, ultimately leading to different α functions and dispersion plots (see below).

The cluster geometries reveal an inherent problem with the described formalism; since loosely bound, the cluster structures are not very rigid and describing them as isolated states whose partition functions are calculated in the rigid-rotor/harmonic-oscillator approximation is a source of error. For example, internal rotations may have small barriers and should thus be treated as a free rotation rather than a vibration. While this could be done in principle—although it would take some work—other problems are more difficult to deal with. The inherent anharmonicity of the dissociation coordinates, the number of low frequency modes, and the (presumably) low barrier transition between two cluster conformers of the same size all introduce error within the harmonic-oscillator approximation to the PES and the superposition approximation. We chose our approach for the sake of having closed formula for the partition function in which the (ion) temperature can be applied in a straight forward way and we expect that this should be qualitatively correct for relatively rigid systems. However, our approach is not expected to perform well for fluxional, highly anharmonic species.

We should also note that there is always the possibility that some local minima were missed with the BH workflow and hence were not included in the formalism. We tried to minimize this source of error by optimizing additional manually generated cluster structures not provided by BH (e.g., solvents binding to the edges of the tetrahedron), but none of those was found to be stable minima.

CCS Fits: MD vs. Rigid Samples

The collision cross section (CCS) of a molecule is temperature dependent in two ways: (1) the effective temperature (following two-temperature theory), which describes the relative velocity distribution between ions and collision gas, enters Eq. (3) directly, and (2) the internal energy of a molecular ion, corresponding to the ion temperature, will increase the vibrational motion and thus the occupied configurational space. Point (2) is not taken into account by Eq. (3). Thus, performing MD simulations at a certain ion temperature and taking random sample geometries to calculate the CCS for each geometry with Eq. (3) at the respective effective temperature should incorporate both effects. We note, though, that inelastic collisions are not modeled this way because the geometries employed for MobCal-MPI calculations are treated as rigid. New approaches to this problem have been published recently [59], and it would be interesting to compare the results of the method presented here when modeling the inelasticity of collisions.

To conduct our treatment, we took 100 random geometries from each MD simulation, hence performing 100 CCS calculations for each temperature. We assume Tion ≈ Teff as should be true for an atomic collision gas and molecular ions [7 chapter 6-6 p. 358]. Although our calculations are performed in molecular nitrogen, we assume this relation holds true, rationalized by the rather stiff nature of the N2 bond.

Figure 3 shows the results for Me4N+, [Me4N+ACE]+, and [Me4N+(ACE)2]a+. We only studied the temperature range up to 700 K because of the time-consuming MD simulations. Generally, it can be noted that \( \overline{\Omega_{\mathrm{MD}}}>\overline{\Omega_{\mathrm{rigid}}} \), i.e., the rigid CCS are always within the MD distributions but towards the lower end, hinting at the anharmonic nature of the PES. For a completely harmonic oscillator, these averages should be equal. Taking a look at the results for Me4N+ (Figure 3a), it is safe to say that the contribution due to vibrational extension is negligible compared to the error of the CCS calculation, i.e., the distributions are very narrow and easily within the error of the rigid calculations. To some degree, this still holds true for the [Me4N+ACE]+ cluster (Figure 3b), although a significant increase in the width of the distributions is visible. In both cases, though, the difference between the fitting functions is negligible. In the case of the [Me4N+(ACE)2]a+ cluster (Figure 3c), the width of the CCS distributions from MD sampling becomes relatively large. We conducted an in-depth analysis of the MD simulation performed at 450 K, where the distribution covers a range of more than \( 10\kern.1em {{\mathrm{A}}^{\circ}}^2 \). During the trajectory, one of the ACE molecules moved from one face of the tetramethyl ammonium to another. When ACE passes over the edge, the ion-solvent cluster increases in size, explaining the broad distribution (skewed towards higher CCS) observed. This, again, hints that the harmonic oscillator treatment for the partition function could be improved by explicitly correcting for anharmonicity. While the mean values still agree within error, the error of the rigid structure CCS does not encompass the entire MD distribution.

Figure 3
figure 3

Calculated CCSs and respective fitting functions for the MD simulations (blue) and the rigid ion/cluster structures (red) over a range of effective temperatures. (a) Me4N+. (b) [Me4N+ACE]+. (c) [Me4N+(ACE)2]+a. The CCS distributions for each MD simulation are also shown (right) as histogram and smooth Gaussian kernel distribution

We chose to perform rigid CCS calculation only and ignore the vibrational broadening for two reasons. First, a single fitting function, even when considering all MD data, cannot account for the width of the CCS distribution, and we are only using the average value of the CCS in the Mason-Schamp equation. Insufficient sampling of the accessible configurational space due to missed ergodic mixing is also a major problem of the MD approach and is likely the reason for the nonphysical set of parameters obtained for [Me4N+(ACE)2]a+ (c should be positive, b is very different compared to the other functions; see Figure 3c). Second, there are two counteracting effects regarding the temperature. Although the actual size of the molecule increases with temperature due to vibrational broadening, the overall CCS decreases because the long-range interactions with the background gas become less important for higher relative velocities [55]. The fact that the two fitting functions in the [Me4N+(ACE)2]a+ case seem to meet again towards higher temperatures could be interpreted in this manner. This might not be the case for all systems: Especially, peptides or proteins can show significant unfolding and thus overall increase of CCS with (effective) temperature [55]. Relying only on a single geometry optimization (for CCS calculations at all temperatures) in contrast to a set of rather time-consuming MD simulations (for each different temperature) dramatically speeds up the whole formalism. The use of MD simulations for the CCS fit and an opportunity to calculate the partition function from the trajectories are discussed below.

We now focus on the fitting functions themselves, which are given for all molecules in Table 1. Functions in Table 1 differ from those shown in Figure 3 owing to the extended temperature range. The ranges were extended because mobility calculation for some species yielded higher temperatures than 700 K; extended fits ensure that the CCS at a specific temperature was always interpolated and never extrapolated. Because the exponent can vary, the fitting parameters are rather insensitive, i.e., two different sets of parameters can yield a very similar fitting function in the temperature range studied. However, those two functions can differ significantly outside of that temperature range, rendering extrapolation of the CCS rather difficult. It would be beneficial to fix the exponent to a value derived from theoretical considerations, since it might well be a constant. While this is beyond the scope of this work, we do note that the average exponent is around −0.72 and no trend regarding cluster size is visible. Also, it is always lower than the value of −0.5, which was previously described in the literature [30].

Table 1 Parameters for CCS Fitting Functions (Ω(1,1) = a Teffb + c ). For Each Ion, a Different Temperature Range Is Considered, Depending on the Observed Maximum Effective Temperature. The Smaller Ions Can Move Faster and Thus Can Have a Higher Effective Temperature

Regarding the other two parameters, it can be seen that there is a positive correlation of the offset c with size of the cluster geometry. This is expected since the offset can be interpreted as CCS in the infinite temperature limit, i.e., when long-range interactions are minimized and the CCS is determined only by hard-sphere collisions [55]. The proportionality constant a seems to decrease with size of the ion. Another benefit of knowing the exponent, b, from theory would be that a and c could be determined from just two CCS calculations, covering the needed temperature range.

We should note that the CCS calculations were performed in nitrogen gas only. All interaction with the solvent vapor is modeled only by the stable clusters. However, having about 1.5 mol% solvent in the gas phase should have an influence on the CCS not only because they are larger than N2 but also because they have a permanent dipole moment and thus different interaction potentials (V ∝ r−2 for ion-permanent dipole compared to V ∝ r−4 for ion-induced dipole). Blanc’s law [60] could be used to account for this additional effect by calculating the CCS in N2 and in the appropriate solvent vapor, then averaging the two according to their respective mole percent. However, we have no means to calculate the CCS in a pure, e.g., MeOH environment with comparable accuracy as for N2. Moreover, Blanc’s law would be a rather rough approximation for the mixture at higher field strengths. It is therefore rather difficult to estimate the error associated with modifier collisions, but we expect it to be relatively small owing to the relatively low partial pressure. Note also that the interaction potential parameters implemented in the MobCal-MPI code were trained to reproduce experimental CCS at 298.15 K. Although these parameters could be different at different temperatures, we assumed them to be constant at all temperatures. Since long-range interactions become less important with increasing effective temperature, the error introduced by using fixed parameters should vanish for higher temperatures. [61]

Mobility Calculations

Having determined the geometries, CCS fits, partition function parameters, and energies, we can now conduct the mobility calculations. Figure 4 shows exemplary results for the [Me4N+(ACE)n]+ system for two different field strengths, i.e., 25 Td and 76 Td. A background gas temperature of 373 K and a solvent concentration of 1.5 mol% was applied.

Figure 4
figure 4

(a) Drift velocities. (b) Effective temperatures. (c) Collision cross sections. (d) Boltzmann weights. (e) Partition functions. (f) Cluster size distributions for 25 Td (500 V/mm at 1 atm and 373 K, blue data) and 76 Td (1500 V/mm at 1 atm and 373 K, red data) for the [Me4N+(ACE)n]+ system. As expected, the increase in field strength, and thus effective temperature, shifts the cluster size distribution to smaller cluster sizes, increasing the overall mobility of the ensemble (type A). For a detailed explanation of the panels, see text

In Figure 4a–c, which shows the self-consistent mobility results, it can be seen that a larger CCS correlates with a lower effective temperature and also a lower drift velocity, as one would expect. During the low-field condition, the effective temperature is only slightly increased compared to the background gas, whereas substantial heating occurs for the small Me4N+ ion under high-field conditions. This behavior arises due to the increase in drift velocity owing to the higher electric field and concomitant decrease in CCS (see Figures 3 and 4c). The effective temperatures of the larger clusters also increase with field strength, but not as much as it does for the small bare ion.

The data from Figure 4b, c gives the opportunity to calculate collision frequencies according to two-temperature theory [7, p., 275]. For the smallest ion, i.e., Me4N+, we estimate the collision frequency with N2 to be 1.5 × 1010 s−1 at 25 Td. While the larger clusters have a higher CCS, their relative velocities are smaller due to the lower effective temperature and the higher reduced mass. The resulting collision frequency is in the same order of magnitude, e.g., for the [Me4N+(ACE)2]a+ cluster at 25 Td about 2.1 × 1010 s−1. Typical DMS devices operate with a waveform frequency in the megahertz regime; thus, the collision frequency is four orders of magnitude higher. This validates our assumption of instant thermalization since there is a sufficient number of collision per cycle of the waveform.

Figure 4d shows the Boltzmann weighing factors derived from our implementation of the superposition approximation, again assuming Tion ≈ Teff [7 chapter 6-6 p. 359]. Since the [Me4N+(ACE)2]a+ cluster is energetically lowest, its weight factor is always 1, independent of temperature (cf. Eqs. (13a) and (14a)). In contrast, the bare ion has the lowest weighting because it is the highest in energy (no contributions from bound solvents). Multiplying these weights with the partition functions plotted in Figure 4e result in the population distribution shown in Figure 4f. Although the bare ion has the lowest weight, its partition function is orders of magnitude larger due to the additional translational degrees of freedom (or in other words, the entropy gain due to complete dissociation). Changing the field strength results in a change in the Boltzmann weights and the associated partition functions. The careful interaction of those two values ultimately determines the relative population of each cluster structure. Since the variation in weights and partition functions are not the same for all structures, the resulting population distribution changes with field strength. As expected, the mean cluster size decreases with higher field strength, ultimately leading to a higher ensemble mobility (type A behavior).

It should be noted that, while the self-consistent calculation is an elegant way of calculating drift velocity and effective temperature from just the CCS, the use of the Mason-Schamp equation is still a source of some error. Even when applying two-temperature theory, there are known deviations in the order of 10% at higher field strengths [22]. Although there are higher order corrections possible to two-temperature theory, as well as a different approach called momentum-transfer theory [12, 22], to date there is no closed formula available which correctly predicts the mobility of an ion for high field strengths.

We also want to stress that the error in the partition function calculations increases with temperature. This is due to the following: (1) neglecting anharmonic contributions, which become more important at higher temperatures (as apparent in the CCS distributions; see Figure 3) and (2) neglecting ro-vibrational coupling, which increases with temperature and is especially significant for low frequency dissociative modes.

Hence, both factors in Eq. (15), i.e., the individual mobility and the relative population are source to an increasing error with increasing field strength.

Figure 5a shows the α function and its derivative, as determined with our method, for the [Me4N+(ACN)n]+ system in N2 with Tbath = 373 K and a solvent vapor concentration of 1.5 mol%. Figure 5b shows the corresponding relative population of each cluster structure. The small changes of the distribution up to 20 Td lead to an almost constant mobility (or α function), which can be interpreted as the well-known field independency of the mobility in the low-field limit. Further increasing the field strength leads to a gradual decrease of the mean cluster size and thus mean CCS, which in turn increases the mobility of the ensemble (type A). At very high field strengths, when there is no larger change in cluster size distribution (Pn = 0 → 1) and the individual CCS has almost converged to its asymptote c, the mobility actually decreases due to the fact that \( K\propto 1/\sqrt{T_{\mathrm{eff}}} \) (cf. Eq. (2)). Thus, three major factors contribute to the overall mobility of an ensemble: (1) the general \( K\propto 1/\sqrt{T_{\mathrm{eff}}} \) effect (sometimes called “hard-sphere effect” [25]), (2) the change of CCS when the background gas is polarizable, and (3) the change in the cluster distribution when there are polar solvents able to form stable clusters with the ion.

Figure 5
figure 5

(a) α curve and its derivative. (b) Relative ion-solvent cluster population as a function of reduced field strength (E/N) for the [Me4N+(ACN)n]+ system. The strong increase in α corresponds to a type A behavior and is reasoned by the overall decrease of CCS due to a smaller average cluster size. Only for high field strengths, where only the bare ion is populated, the hard-sphere effect becomes significant and the α curve returns to smaller values. Similar figures for MeOH and ACE are found in the SI (Figs. S3.1 and S3.3)

Dispersion Plots

Figure 6 shows the calculated and experimental dispersion plots for Me4N+ in pure N2 and seeded with 1.5 mol% of MeOH, ACN, and ACE. Although the authors of the experimental paper report a temperature of 423 K inside the DMS cell [34], using a thermocouple for a coarse measurement of the gas temperature inside the DMS cell, we found the temperature to be only around 373 K. Consequently, we applied Tbath = 373 K in the calculations. Variations of the temperature used for the calculations showed some noticeable change in the dispersion plots within ±20° but no qualitative difference within the aim of this work. A gap size of d = 1 mm is used is Eqs. (17) and (18), and we discretized the waveform into 1000 time steps for the averages in Eq. (18), which was sufficient with regard to convergence.

Figure 6
figure 6

Calculated (blue) and experimental (orange) dispersion plots for Me4N+ in (a) N2 at 1 atm and 373 K and (b) N2 seeded with MeOH, (c) ACN, and (d) ACE as modifiers (1.5 mol% at 373 K). The experimental data are taken from [34]

First, we discuss the dispersion plot of the pure nitrogen environment. Because the CCS of the bare ion decreases with effective temperature, the mobility increases slightly at lower field strengths, thus requiring negative CVs to correct trajectories. Further increasing the separation field strength leads to hard-sphere scattering, which requires increasingly positive CV values for optimal ion transmission. Thus, a very weak type B curve is generated. The calculated dispersion plot nicely reproduces the experimental data and the two curves match qualitatively.

Upon adding a weakly interacting solvent vapor to the system, i.e., MeOH, where only n = 1 is populated even at low field strengths (see Fig. S3.1), the type B curve becomes more pronounced (i.e., the minimum CV becomes more negative and shifts to higher SV values). Our model is able to reproduce this behavior qualitatively, as well as the overall shape of the curve. In particular, the SV, where the CV has a minimum, called SV@CVmin in the literature and reported to correlate with binding energy and other physicochemical parameters such as pKa, pKb, and Hammett parameters [62], is nicely reproduced. Introducing a more strongly interacting solvent vapor (ACN, ACE) yields type A behavior for the range of SV studied. This indicates that the ion mobility under the high-field conditions is larger than under low-field conditions. Again, this behavior is qualitatively reproduced by our model. Note that the formalism is the same in all cases; no system-specific parameters are introduced and every result arises from the (modeled) nature of the system itself.

At SVpp > 2000 V significant deviations occur between calculated and experimental dispersion curves, especially for ACN. There are several sources that could give rise to these errors. As mentioned previously, the description of the partition functions for high ion temperatures and the use of the Mason-Schamp equation for very high field strengths introduce errors. In addition to those sources of error, our assumption of instant equilibration could break down. While the N2 collision frequency of the order 1010 s−1 should be sufficient to quickly thermalize the ions, the clustering reaction depends on the collision frequency with the solvent molecules rather than with N2. This number is difficult to calculate since we do not have access to \( {\Omega}_{\mathrm{solvent}}^{\left(1,1\right)} \). A rough estimate could be given by just scaling the N2 collision frequency by the mole percent of the solvents which would give a frequency in the order of 108 s−1. This approximation yields only about 100 collisions per cycle of the waveform with solvent molecules. Figure 5 shows that at low-field conditions the n = 2 cluster dominates the population, while at high-field conditions, only the bare ion is populated (SVpp = 4000 V corresponds to a maximum reduced field strength of 135 Td under the given conditions). Thus, having a waveform frequency of ω = 3 MHz, the cluster formation reaction might not be fast enough (due to insufficient collisions) to populate the n = 2 state at its equilibrium value. Additionally, the rate of evaporation given by the unimolecular decay of the heated ions in the high-field portion could also be a limiting factor. A change from n = 2 to n = 0 could be difficult to achieve on these time scales. Indeed, the more negative CV values, compared to the experimental values, suggest a larger dynamic change in clustering in the calculations than in the experiment. Kinetic simulations with theoretically determined rate constants could circumvent this problem and lead to a more accurate population distribution. The deviations in the ACE case are significantly smaller and indeed the change in cluster size between high- and low-field is not as pronounced (see Fig. S3.3), supporting the suggested error due to non-equilibrium conditions. For MeOH, an underestimation instead of an overestimation is visible and the deviations are smaller overall (see Figure 6b). Since MeOH shows only weak clustering compared to ACN and ACE (see Fig. S3.1), we do not expect kinetic effects to play a major role. Thus, the deviations are solely due to other errors already mentioned, e.g., neglecting anharmonicity of the PES and the use of two-temperature theory for very high field strengths.

It is important to stress again the critical role of the CCS fitting function. Applying a \( \Omega \propto 1/\sqrt{{\mathrm{T}}_{\mathrm{eff}}} \) fitting function as stated in the literature [7, 30] would actually cancel the two-temperature effect in the Mason-Schamp equation, as already mentioned in the literature [30]. This would lead to a constant mobility for the pure N2 environment leading to no differential mobility (CV = 0 V for all SVpp). This is clearly not observed in the experimental data (see Figure 6a). The determined exponent of b <  − 0.5 is thus critical for reproducing the correct form of the dispersion plot.

To further test our model, we exploited varying the background gas temperature and the solvent concentration (see Figure 7). Experimental data [23, 27, 30, 63] show that with increasing concentration of the solvent vapor, more negative CV values are needed to transmit the ions through the DMS cell. This results in a pronounced type B behavior that evolves into type A as concentration increases. This can be explained by a higher degree of clustering and thus larger change in CCS upon evaporation since an increase in gas-phase concentration will shift the equilibrium towards more clustered ions. This is modeled by the concentration dependency of the relative population (cf. Eqs. (6) and (14)) and our model nicely reproduces this trend (see Figure 7a).

Figure 7
figure 7

(a) Variations of modifier concentration and (b) background gas temperature for Me4N+ in N2 (1 atm) seeded with MeOH (expt conditions are at 1.5 mol% MeOH and 373 K). Increasing the concentration as well as decreasing the temperature leads to pronounced type B behavior due to higher degree of clustering

Increasing the background gas temperature reduces clustering due to an increased internal energy and subsequent shift of the equilibrium towards higher degree of dissociation. Consequently, less pronounced type B behavior is observed experimentally as Tbath increases [27, 32, 64]. Again, our model reproduces this trend correctly (see Figure 7b). When Tbath is raised to the point at which ion-solvent clustering is prevented even at low-field conditions, the dispersion plot converges to that of the pure N2 atmosphere (cf. Figure 7b dashed line) as expected.

Conclusion

The ability of our model to reproduce dispersion plots for different separation conditions, i.e., different (or no) solvents added, different concentrations and different background gas temperatures, strongly suggests that we capture the most important effects leading to differential mobility, namely clustering, the hard-sphere effect, and the change in individual CCS with effective temperature.

Using two-temperature theory, we model the overall dependence of the ion mobility with field strength. This is important to account for the overall decrease of mobility with field strength, known as the hard-sphere effect. With no solvent added and having a very weakly polarizable gas-like helium, this effect is dominant, leading to type C behavior. While two-temperature theory covers most of this effect, it is only an approximation and known to show deviations for very high field strengths [22]. Higher order terms or a different approach called momentum-transfer theory [12, 22] could improve this description.

Having a polarizable collision gas like N2, it is important to model the temperature dependence of the CCS. The reported dependency of Ω(1,1) ∝ Teff−0.5 [7, 30] not only shows systematic deviations from calculated CCS values, but, in combination with two-temperature theory, it also fails to reproduce type B or type C systems. We model the functional dependency by fitting a more general function to CCS values of the equilibrium structure at different temperatures. Doing this, we neglect vibrational broadening of the molecule and inelasticity of collisions.

If solvents are added to the gas phase, clustering plays an important role and only when this is taken into account, type A dispersion plots can be predicted. We model the cluster size distribution by assuming fast thermal equilibrium and using a special version of the superposition approximation, applying harmonic approximations to the partition functions. While this seems to model the data qualitatively, two main errors are introduced: (1) The equilibrium condition might not hold for fast oscillating fields. Kinetic simulations with theoretically determined rate constants would give a better description of the population distribution. (2) The simple description of the partition functions neglects anharmonicity of the PES, as well as ro-vibrational coupling and easy conversion between cluster structures. This could be partially improved by using, e.g., the quasi-harmonic approximation [65,66,67], where the vibrational partition function is determined from MD simulations on the anharmonic surface. Conducting MD simulations for a set of temperatures (ensuring sufficient sampling of the configurational space) could be used to obtain a more accurate description of the temperature dependence of the vibrational partition function as well as more accurate CCS fits (accounting for vibrational broadening).

Despite the simplifications made, our model is able to predict trends in DMS dispersion plots from first principles for different experimental conditions, i.e., different (or no) solvents added, different solvent concentration, and different background gas temperature. To our knowledge, this is the first time a thorough model is presented that can account for all these effects. The importance of different contributions to overall mobility, especially clustering with solvent molecules, is highlighted.