Design of Nonideal Eutectic Mixtures Based on Correlations with Molecular Properties.

In this work, a statistical analysis was performed to reveal how the molecular properties are correlated with the nonideal behavior observed in eutectic mixtures. From this, a statistical model, combined with theory and experimental results, was developed to predict the nonideal behavior of a specific set of eutectic mixtures, consisting of quaternary ammonium bromides with dicarboxylic acids and polyols. The combination of this analysis and this model can be considered as a first step toward the a priori design of eutectic mixtures. The analysis performed is based on principal components. The descriptors used for this are molecular properties of the constituents of these mixtures. The molecular properties are a combination of experimental, theoretical, and computed properties. The analysis reveals that there are strong correlations between the nonideality of the mixtures and a measure of the acidity of the hydrogen bond donating protons, the displacement of the bromide anion, and the bulkiness of the quaternary ammonium salt. Our analysis highlights the design rules of deep eutectic systems (DES), enabling control over the extent of the liquid window. Our model enables prediction of the eutectic temperature for a range of related mixtures.

Introduction

Since the first mentioning of eutectic mixtures as deep eutectic solvents[1] (DESs), they have been compared to ionic liquids (ILs). Ionic liquids are eulogized for being designer solvents. This means that the physical properties of ILs can be tailored based on the choice of the anion and cation. Similar to ILs, DESs have the potential of being designer solvents, tailored toward specific applications based on the choice of the mixture’s constituents. However, the understanding of the relation between the molecular properties of DES constituents and the properties of the mixture is still far from complete. Without understanding which properties of the constituents connect to the physical and chemical properties of the mixtures, designer mixtures will not materialize. Research on these relations has until now mostly focused on the dependence of physicochemical properties, such as density and viscosity, on the hydrogen bond donor bonding capacity,[2,3] the water content,[4,5] or temperature.[6,7] Recently, Abranches et al. reported results on similar salt mixtures and observed unexpected nonidealities.[8,9] It is however difficult to draw conclusions from these types of mixtures, as the ideal reference state is unknown due to a lack of fusion property data of the pure compounds. Therefore, this reference state can only be estimated.[10,11] Recently, the importance of differences in acidities was pointed out for the preparation of nonionic, nonideal eutectic mixtures.[12]

García and co-workers aimed at relating the electronic energy, a computed molecular property, of several mixtures to their observed eutectic melting point.[13] Their study is of particular interest when screening new mixtures for applicability in the desired temperature range. However, the method employed still requires the electronic energy to be computed for every mixture under consideration. For ILs, progress was made with a priori screening for certain properties of not yet synthesized ILs via standardized computed quantities, for example, predictions of the Kamlet–Taft parameters.[14]

In this paper, a step toward a priori design of the liquid window of nonideal eutectic mixtures is endeavored. For this, common physical properties and estimations of the molecular interactions obtained via density functional theory (DFT) of the individual components were gathered. It is hypothesized that these properties contribute to the attractions between the different components and to be responsible for the nonideal behavior of the eutectic mixtures. These properties were examined using a principal component analysis to identify their relations to interaction parameters obtained from phase diagrams of 18 systems. Using these relations, the nonideality and thus the eutectic temperatures for the 22 other systems were predicted. With this approach, we developed a first step toward a theoretical framework for the quantification, the description, and the prediction of the phase behavior of DESs.

Thermodynamic Quantification of the Nonideality

The nonideal contributions to the melting point depression can be quantified by obtaining interaction parameters via a description of the molar Gibbs energy g of mixing and by differentiating between contributions from (molar) entropy s and (molar) enthalpy h:Here g can be related directly to experimentally obtained solid–liquid phase diagrams via the chemical potential μwhere n is the total number of moles, n is the number of moles of component i, x is the mole fraction, and γ is the activity coefficient. Here μ – μ* = Δμ is related to the melting point depression of the mixture and to pure component properties through the Schroeder–van Laar equationwhere ΔH is the enthalpy of fusion, T* the melting point of the pure component i, and T the melting point of the mixture.

The entropy model used here is the volume- and surface-based Staverman–Guggenheim entropy with s given byHere ϕ denotes the volume fraction, defined byand θ represents the surface fractionwith Vm, and Am, being, respectively, the molar volume and surface for component i. Here Q is a direct function of ϕ and θ:To characterize the enthalpic contributions, regular solution (RS) theory[15,16] is used. RS theory quantifies the enthalpic contributions using one interaction parameter, χ, and was shown to accurately describe nonidealities typically found in DESs:[17]This expression was first proposed by Margules,[18] without invoking lattice approaches introduced later, which allow the interaction parameter χ to take both positive and negative values, depending on the relative interactions between the components. It should also be stressed that, while for apolar systems the χ parameter may be predicted from (Hansen) solubility parameters, this is not straightforward for polar systems where hydrogen bonding plays a role.[19]

Previously, it was demonstrated that a Redlich–Kister-like expansion can be used to more accurately describe the solid–liquid phase boundaries without changing the value, or physical meaning, of χ due to the orthogonal nature of Legendre polynomials.[20] Therefore, χ can be determined accurately, based only on the fusion properties of both components and the eutectic temperature measured for the mixture using only RS theory.

Combining RS theory and Staverman–Guggenheim[21−24] results in the following expression for Δμ:

Model Systems

The full set of systems consists of eight hydrogen bond donors (HBDs) and five hydrogen bond acceptors (HBAs), resulting in 40 possible combinations. New DESs were developed rather than using previously published data to ensure that components have well-defined physical properties and do not react with themselves or each other; e.g., combinations of carboxylic acids and alcohols, which can react to esters, were avoided.[25] Since already the selection of compounds that are suitable for determining the phase behavior is not straightforward, we therefore restricted ourselves to these 40 combinations. This selection resulted in the hydrogen bond donors and acceptors shown in Figure . These components differ systematically in specific features. For example, for the HBD in the dicarboxylic acid series, ranging from succinic acid to suberic acid, the alkyl chain length increases stepwise by a single carbon atom. This causes the molecule to become slightly less hydrophilic and to adopt more conformations to alternate between cis and trans. This alternation for the odd and even dicarbolyxic acids is known to influence the melting point and enthalpy of fusion.[26,27] The polyols used range from erythritol to sorbitol, also increasing in chain length, but here the number of hydroxyl groups increases simultaneously. The alkylammonium bromides (HBAs) range from tetraethyl- to tetrahexylammonium bromide, which we will refer to as Et4NBr, Pr4NBr, Bu4NBr, Pe4NBr, and He4NBr.

Figure 1

Molecular structures of the individual components of the binary mixtures studied here. The H-bond donating components are displayed individually. The H-bond accepting components are depicted via a general structure, where R represents the alkyl groups ethyl, propyl, butyl, pentyl, and hexyl.

Table shows the fusion properties, molar volume, and molecular surface for the considered components measured by us and compared to literature values. The values obtained for the components used for this work are in most cases in good agreement with the literature values. The fusion properties for Et4NBr and Pr4NBr could not be obtained experimentally, as these salts degrade before melting. For these components, an alternate procedure has been adopted. For these systems, the liquidus curves for the mixtures were measured at molar ratios accessible by our techniques. The data points, of all mixtures containing the salt of interest, were then simultaneously fitted to eq , from which the unknown fusion properties were derived as global fit parameters for all mixtures and χ as an individual fit parameter for each mixture.

Table 1

Experimentally Obtained Melting Point T, Enthalpy of Fusion ΔH, Molar Volume Vm, and Molecular Surface Am of the Individual Components of the Binary Mixtures Studied as well as the Literature Values (ref) of Several Compounds

	T (K)		ΔH (kJ·mol–1)
component	meas.	lit. (ref.)	meas.	lit. (ref.)	Vmb (cm3·mol–1)	Amb (cm2·mol–1)
succinic acid	460.0	457.1 ([28])	37.1	32.6 ([28])	41.9	10.6
glutaric acid	371.6	371.0 ([28])	21.5	20.9 ([28])	52.2	13.1
adipic acid	426.2	426.4 ([28])	37.4	34.9 ([28])	62.2	15.5
pimelic acid	378.5	377.5 ([28])	26.1	27.6 ([28])	66.5	16.5
suberic acid	415.3	417.3 ([28])	30.2	28.8 ([28])	75.4	18.6
erythritol	394.7	365.9 ([29])	39.3	36.5 ([29])	46.4	11.7
xylitol	370.4	391.9 ([29])	36.6	40.7 ([29])	55.6	13.9
sorbitol	369.6	366.4 ([29])	22.5	27.7 ([29])	66.9	16.6
Et4NBr	633a		13.9a		83.3	20.5
Pr4NBr	563a		7.8a		121.6	29.7
Bu4NBr	390.8	395.2 ([30])	10.0	15.5 ([30])	150.0	36.5
Pe4NBr	375.9	376.2 ([30])	40.1	41.4 ([30])	190.7	46.2
He4NBr	370.6	377.2 ([30])	11.6	15.9 ([30])	220.0	53.2

Estimated from solid–liquid equilibrium data.[31]

Generated by MolModPro software.

Table lists the obtained eutectic temperatures and interaction parameters χ for the mixtures considered. All of the phase diagrams obtained for this work are presented in Figures S1–S3. When comparing the interaction parameters for the different mixtures, a few trends can be observed. For example, binary mixtures of adipic acid with three different quaternary ammonium salts were analyzed, in which χ becomes less negative with increasing length of the side chains of these salts, indicating stronger attraction with smaller quaternary ammonium salt. This might suggest that the bromide anion is more accessible when the side chains are shorter. However, this trend does not extend to succinic acid, of which also binary mixtures were studied with three salts. Here Pe4NBr shows a larger negative deviation than Bu4NBr.

Table 2

Measured Eutectic Temperatures Te and Extracted Interaction Parameters χ of the Eutectic Mixtures Studied with Hydrogen Bond Donors and Acceptors Indicated

HBD	HBA	Te (K)	χ
succinic acid	Et4NBr	352.6	–7.92
adipic acid	Et4NBr	352.5	–5.53
glutaric acid	Et4NBr	303.3	–5.44
erythritol	Et4NBr	346.1	–4.79
pimelic acid	Et4NBr	323.5	–4.29
adipic acid	Pr4NBr	344.9	–3.67
pimelic acid	Pr4NBr	321.0	–2.26
sorbitol	Pr4NBr	336.3	–0.64
suberic acid	Bu4NBr	315.9	–3.15
xylitol	Bu4NBr	313.2	–2.24
succinic acid	Bu4NBr	355.8	–1.30
sorbitol	Bu4NBr	312.5	–1.03
pimelic acid	Pe4NBr	309.6	–5.95
succinic acid	Pe4NBr	349.3	–4.01
sorbitol	Pe4NBr	319.7	–3.43
xylitol	Pe4NBr	332.4	–2.77
erythritol	Pe4NBr	346.4	–1.90
adipic acid	He4NBr	332.0	–2.16

Another observation from the interaction parameters from these mixtures comes from the mixtures with Et4NBr as HBA. Here the interaction parameter becomes more negative with decreasing molecular size, while it follows the odd–even effect described earlier for these dicarboxylic acids. This effect explains why χ for adipic acid with Et4NBr is slightly more negative then one would expect and very close to the χ obtained for glutaric acid with Et4NBr. When considering these trends, an oddity is observed in the mixture of pimelic acid with Pe4NBr. As pimelic acid is larger than succinic acid, it would be expected to have a χ closer to zero than succinic acid with Pe4NBr. Looking further into the trend of pimelic acid with the different quaternary ammonium salts, again an oddity is observed. Also, here the previously observed trend of stronger interactions with smaller cations is broken.

Properties as Descriptors

Trends in nonideality in eutectic mixtures observed in earlier publications[17,20] suggest a relation with hydrogen bonding, where nonideality is supposed to increase with stronger hydrogen bonds. However, no straightforward and general rule of thumb is yet established to support such a relationship, as discussed in the previous section. Here, it is hypothesized that the nonideality of mixtures can be qualitatively predicted based on their corresponding H-bond characteristics. Achieving these predictions will result in a better understanding of the interactions responsible for the dramatic melting depressions from which deep eutectics gained their names.

In order to include a wide range of properties in this exploration, common and new properties of the individual components were used. The properties chosen are expected to indicate the strength of the H-bond interactions, such as pKa, molecular dipole moments, and number of H-bonding moieties. Furthermore, H2O solubility as well as the Hansen solubility parameters of the components of the binary mixtures were considered. Quantitative information on hydrogen bonds can also be computed,[32] for example, by using density functional theory.

These properties will be used as descriptors in a principal component analysis and for a regression analysis. The goal is to identify which properties, or combination of properties, can explain the trends in nonideality, as observed in a selected set of eutectic mixtures, which later also will be used to predict eutectic properties of the mixtures.

Properties

It can be rationalized that all properties used are candidates for explaining the variety in nonideality, quantified as interaction parameters, observed in the eutectic mixtures studied here. Logically, the fusion properties, the enthalpy of fusion ΔH and melting point T*, are already considered as part of the entropy of fusion and may not contribute to the enthalpy of mixing. However, a synergistic effect between the fusion properties and another property should not be excluded from the correlations.

Hydrogen bonds are directly related to the Brønsted acidity of components.[35,36] Hence, the acid dissociation constant Ka, expressed as pKa, is added as a descriptor to quantify the hydrogen bond capabilities. The proton affinity[37] is closely related to the pKa and is calculated here using DFT, as described in the Methods section. Next to the proton affinity, the solubility of the components in H2O as well as their molecular dipole moments were used as descriptors. The dipole moments again show the odd/even behavior for the dicarboxylic acids. These properties are listed in the Supporting Information, Table S3. The enthalpic interactions expected to be responsible for the observed nonideality will not only depend on the hydrogen bond capable moiety but also depend on the rest of the molecule. It is observed in the trends of Table that smaller dicarboxylic acids usually result in stronger attractions.

Similarly, for the quaternary ammonium bromides, longer side chains could make them less accessible, as the interactions are assumed to be taking place through the bromide anion. It is also expected that the amount of hydrogen bond capable groups will influence the nonideal behavior in the mixtures. A selection of these parameters is presented in Table . The complete set of parameters can be found in the Supporting Information.

Table 3

Physical Properties of the Individual Components Related to the Molecular Structures: the Molecular Mass M, Density of the Solid Crystalline Components ρ, and Number (#) of −CH2 Groupsa

component	M (g·mol–1)	ρb (g·cm–3)	# of −CH2	proton affinity (eV)	molecular dipole moment (D)
succinic acid	118.09	1.56	2	1.159	2.14[33]
glutaric acid	132.12	1.402	3	1.448	2.37[33]
adipic acid	152.1	1.355	4	1.513	2.30[33]
pimelic acid	160.17	1.335	5	1.619	2.36[33]
suberic acid	174.2	1.28	6	1.71	2.34[33]
erythritol	122.12	1.459	0	1.41	3.20[34]
xylitol	152.15	1.516	0	1.745	3.60[34]
sorbitol	182.17	1.509	0	1.578	3.90[34]
Et4NBr	210.16	1.397	4	1.218
PR4NBr	266.26	1.213	8	1.213
Bu4NBr	322.37	1.191	12	1.127
Pe4NBr	378.47	1.100	16	1.204
He4NBr	434.59	1.094	20	1.206

Properties of the individual components related to the acidity of the hydrogen bonding moieties: the proton affinity (PA) and the molecular dipole moment.

Experimentally obtained using a Micromeritics AccuPyc II 1340 Gas pycnometer.

The Hansen solubility parameters are based on the “like dissolves like” concept and are partial solubility parameters which differentiate between dispersive forces, polarity, and hydrogen bonding ability. The Hansen solubility parameters,[38] listed in the Supporting Information (Table S4), were estimated using a group contribution method.[39]

Other interesting features may result from interactions of the component of interest with a probe molecule. Their difference in behavior toward the same molecule may yield insight into their behavior in mixtures with each other. To obtain these interactions, density functional theory (DFT) computations were used. Details on the calculations are described in the Methods section. The obtained properties from these calculations are the ground state energies and the optimized distance between the oxygen atom and proton in the −OH groups used to form the hydrogen bond. These properties are calculated for the isolated molecules and in the presence of another probe molecule. Water and formic acid were taken as probe molecules. These probe molecules were chosen for their small size and dual nature of being both donors as well as acceptors of hydrogen bonds. For the acceptors, the optimized distance between the oxygen and the proton in the probe molecule was calculated additionally. Figure presents an overview of the situations considered for the computed properties. The labels are the notations used to refer to specific distances in the molecules. The quantity h refers to a distance between an oxygen atom and a proton. The variable b is used to denote the distance between the N+ cation and Br– anion. The sub- and superscripts refer to the molecule considered and the molecule inducing an interaction, respectively. For example, hFA refers to the distance between the oxygen atom and the proton in component i while it interacts with formic acid. The annotations of the differences in ground state energies of the molecules are defined in Figure and Table . Here the sub- and superscripts refer to the final ground state and the initial ground state(s), respectively.

Figure 2

Visualization of the characteristic distances obtained from DFT calculations; h denotes the distance between an oxygen atom and a proton, indicated by black arrows; b represents the distance between the N+ cation and Br– anion, indicated by the dotted yellow arrows. The subscripts refer to the central molecule; the superscripts refer to the molecule inducing intermolecular interactions.

Table 4

Properties Resulting from the Differences in Electronic Energy with the Probe Molecules Calculated via ΔE = Electronic Energy of the Final State – Electronic Energy of the Initial State 1 – Electronic Energy of the Initial State 2 for Component i

final ground state energy	initial ground state energy 1	initial ground state energy 2	abbreviation
2i	i	i	ΔE2ii
i + H2O	i	H2O
i + FA	i	FA	ΔEFAi
ΔE2ii			ΔE2iH2O
ΔE2ii	ΔEFAi		ΔE2iFA

The resulting ground state energies ΔE are listed in the Supporting Information (Table S5) for all studied components. The optimized distances of each component i, h and b, in the different environments are listed in the Supporting Information (Tables S6 and S7).

Principal Component Analysis

In order to thoroughly evaluate the variety of and correlations between the properties discussed, a principal component analysis (PCA) is performed. The goal of this PCA is to extract information from the data and to express this information as a set of new orthogonal variables called principal components (PCs). These PCs are obtained as linear combinations of the original variables and reduce the dimensionality of the data set by keeping only the most distinguishing information.[40,41] The first principal component is the linear combination of the variables that shows the largest variance. The second principal component is the combination that gives the largest variance while orthogonal to PC1.

Since the eutectic mixtures are combinations of the hydrogen bond donors and the quaternary ammonium salts, every mixture will be characterized by 48 descriptors, 25 of the hydrogen bond donor and 23 of the hydrogen bond acceptor; see the Supporting Information (Tables S9 and S10). These 48 variables will be reduced to principal components (PCs). The eigenvalues of the principal components, of which the total sum equals unity, indicate the contribution of the respective PC to the variation in the data. Table lists the eigenvalues of the first six principal components. These eigenvalues show that PC1 and PC2 contribute equally while the eigenvalues decrease from PC3 and further on. The sum of the eigenvalues reveals that the first four principal components can account for about 80% of the variation in the data and the first six for about 90%. However, in order to minimize the number of free parameters, only the first four PCs will be considered further.

Table 5

Eigenvalues of the First Six Principal Components, PC1, PC2, PC3, PC4, PC5, and PC6

principal component	eigenvalues	cumulative eigenvalues	donor or acceptor
PC1	0.315	0.315	acceptor
PC2	0.311	0.626	donor
PC3	0.093	0.719	donor
PC4	0.074	0.793	acceptor
PC5	0.054	0.847	acceptor
PC6	0.048	0.895	donor

The so-called scores of the mixtures are the new coordinates describing the mixtures generated by the PCs. By plotting these scores for all mixtures, one can visualize whether the PCs are able to distinguish the mixtures, as the goal of the PCA is to generate unique coordinates for each mixture. The first four principal components are plotted in Figure . Here every marker represents a mixture with a specific acceptor and a specific donor. It is evident that the mixtures are well separated in this representation. The coordinates for each specific mixture can be found in the Supporting Information, Table S8.

Figure 3

Score plot of the first four principal components of the eutectic mixtures.

The principal components describe either the hydrogen bond donating component or the hydrogen bond accepting component, since the properties we use of HBD and HBA are independent. The properties were not divided beforehand. For this reason, it is evident that the principal components are to be used in pairs, both accounting for the accepting components as their donating counterparts. This is also intuitive from a physical–chemical point of view.

Next to generating unique coordinates for all mixtures, the principal components contain information about the significance of the variables used. The eigenvectors of the principal components, usually referred to as the loadings, directly describe the composition of the PC by the variables and can be visualized in a loadings plot. The loadings plot of PC1 and PC4, both describing the hydrogen bond acceptors, is presented in Figure (left panel). Variables with similar loadings in the PCs are highly correlated or contain similar information. The length of these vectors represents their contribution to the principal component. Analyzing the loadings plot of PC1 and PC4, two main clusters can be observed. These clusters contain variables which are closely related to each other and contain more or less the same information. Similar clusters are observed in PC2 and PC3; however, here the right cluster is less dense, thus containing less variables.

Figure 4

Loadings plots of the descriptors contributions to the principal components. Colors indicate the distance to the origin, where lighter is closer to the origin.

The distances of the loadings to the origin, which indicate their contribution to the principal components, are listed in Table . The distances are determined using the two most important PCs per category.

Table 6

Loadings of the Principal Componentsa

descriptor donors	length	descriptor acceptors	length
dipole moment	0.993	ΔE2iFA	1.000
δH	0.992	ΔE2ii	0.999
δP	0.989		0.998
# of −COH	0.986	bi	0.993
# of −CH2	0.982	Vm	0.988
hii	0.982	biH2O	0.987
δD	0.981	ρ	0.983
# of =CO	0.977	δD	0.982
pKa	0.971	δH	0.982
ΔE2ii	0.926	# of −CH2	0.982
Vm	0.917	δP	0.982
ΔE2iFA	0.915	M	0.982
M	0.905	biFA	0.979
ρ	0.889	bii	0.967
hi	0.872	T*	0.938
hiH2O	0.871	hFAi	0.926
H2O solubility	0.862	H2O solubility	0.878
	0.836	ΔS	0.781
proton affinity	0.824	ΔH	0.760
T*	0.815	proton affinity	0.704
ΔH	0.728	hiH2O	0.602
	0.714		0.536
hiFA	0.698	ΔEiFA	0.430
ΔS	0.649
ΔEiFA	0.543

The values of the loadings are the distance to the origin and are determined using the two most important PCs per category.

Results and Discussion

What can be noted directly from the loadings of the descriptors of the donors in Table in the left column is that the dipole moment, the Hansen solubility parameters δH and δP, and the number of −COH groups in the molecule are highly related and contribute most to the variation observed in the donors. It is not surprising that these properties are highly related, as the dipole moment strongly depends on the different moieties related to hydrogen bonding. Also, the Hansen solubility parameters δH and δP describe the hydrogen bonding ability and the polarity of a molecule. This cluster of properties is closely followed by the remaining structure descriptors, being # of CO, # of CH2, the bond length between the oxygen and the proton of the acid group when interacting with itself h, and δD. Also, here the findings are according to expectations. It is well-known that the availability of a proton for hydrogen bonding depends not only on the direct neighbors but also on the entire molecular structure for such small molecules. Descending with the length of the vectors means that these properties do not contribute as much to the variety in the information provided for these components. The lowest contributing properties for the HBD are the fusion properties, as here the lengths suddenly drop significantly, and several computed properties related to the probe molecules.

For the loadings of the acceptors, in Table in the right column, a similar grouping of descriptor clusters can be observed. Here the electronic energy differences with the probe molecules and two identical molecules all contribute highly. It is somewhat surprising that the response to the probe molecule compared to the self-interaction and an identical molecule all contain the same information, making the use any of the probe molecules here good candidates for testing interaction abilities of quaternary ammonium salts. Another interesting observation is the apparent difference in information provided by the bromide anion distance between the anion and the cation, b, and the length of the side chains of the cation. These quantities were expected to be highly related, as longer side chains are assumed to make the charged cation less accessible. Again, here the fusion properties score very low as contributors. Together with the proton affinity and the reaction of a single molecule to a probe molecule, they contribute least.

Based on these results, a selection was made of the properties which contribute most to the variation of the data, and are not highly correlated. These properties are listed in Table and are used to develop a statistical model to predict the nonideality of the remaining mixtures which were not measured in this work.

Table 7

Descriptors of Both Donors and Acceptors Used for the Regression and Their Aliases

descriptor	donor/acceptor	symbol
ΔE2iFA (eV)	acceptor	A1
bi (Å)	acceptor	A2
ρ (g·cm–3)	acceptor	A3
# of −CH2	acceptor	A4
δH	acceptor	A5
dipole moment (D)	donor	D1
# of −COH	donor	D2
ΔE2iFA (eV)	donor	D3
ρ (g·cm–3)	donor	D4
proton affinity (eV)	donor	D5

One may question why the solid phase density could be a fundamental descriptor for both hydrogen bond donors and acceptors. The density gives information on how well molecules are packed into lattices. The packing of crystal lattices has a great influence on the melting properties. This should not affect the nonideality of course, but it does influence the melting temperatures. Indeed, we do see more deviations from nonideality when glutaric acid and pimelic acid (the less densely packed ones) are involved, so one could postulate about how this relates to the interactions in the liquid phase due the stereochemistry of the molecule. Perhaps the density is a good indicator of this.

Regression Analysis

Regression was performed to correlate the differences in nonideal behavior observed in the eutectic mixtures and quantified by the interaction parameter χ, as listed in Table . The principal component analysis resulted in four principal components, which together account for 80% of the variation in the data. These principal components were used for multiple linear regression, and subsequently to identify the most important properties to use as descriptors. Table lists the descriptors used for the regression and their aliases.

The regression was carried out using JMP software. This software was used for the principal component analysis, for the visualization of results, and for the regression analysis, leading to the various modeling strategies described below.

Models

The regression was performed for various situations: linear regression with the principal components, linear regression with descriptors directly chosen from the physical properties, and a nonlinear version of the latter regression analysis. These descriptors were chosen based on the loadings resulting from the principal components analysis. The results are listed in Table . The principal component regression resulted in a model with a correlation coefficient of R2 = 0.69.

Table 8

Descriptors Used for the Regressiona

	descriptors
regression	linear terms	quadratic terms	cross terms	R2
linear	PC1			0.69
	PC2
	PC4
nonlinear	PC2	PC4		0.71
	PC4
linear	A2			0.73
	A5
	D1
nonlinear	A2	D1	D5·A4	0.93
	A4	A2	D5·A2
	D1
	D5

The last column lists the correlation coefficients.

The second method of modeling employs the separate dominant donor and acceptor individual descriptors. They are selected based on the PCs. The subsequent models are independent from those based on PCA, and resulting R2 values can either be poorer or better, depending on the number of allowed terms. Interestingly, the linear model with singular descriptors resulted in a model with a similar accuracy of R2 = 0.73. This means that using a combination of all of the descriptors, in the form of four principal components, does not yield a better description. The reason is that the principle component analysis extracts the greatest variance in the data, but there is no guarantee that this variance will correlate strongly to another quantity (in this case the interaction parameter χ). Therefore, even when using only three singular descriptors, two for the acceptor and one for the donor, a slightly better description is obtained. The descriptors used here are the distance between cation and anion b, the Hansen solubility parameter δH which is related to the hydrogen bonding of the acceptor, and the dipole moment of the donor.

Both linear models can already be considered as statistically relevant; however, the nonlinear model with singular descriptors is significantly better with R2 = 0.93. The fact that the nonlinear model provides a more accurate description suggests the importance of synergistic effects. The prediction of the interaction parameter χ and the subsequently calculated eutectic temperature Te and eutectic composition xe of the remaining 22 combinations of the hydrogen bond donors and acceptors, using the model equation (eq S1) presented in the Supporting Information, can be found in Table S2, based on this nonlinear model.

The accuracy of the model is visualized in Figure , where the predicted values of the interaction parameter χ and the subsequently resulting eutectic temperature Te and eutectic composition xe are plotted versus their experimentally obtained counterparts. The data show close agreement between predictions and observations.

Figure 5

Comparison of the nonideality parameter χ and the eutectic temperature Te as predicted by the nonlinear model and the measured values.

Concluding Remarks

In this work, the nonideal behavior of 18 eutectic mixtures based on 5 hydrogen bond acceptors and 8 hydrogen bond donors was quantified in terms of the interaction parameter χ from regular solution theory with experimentally determined phase diagrams. For the 22 other possible combinations of acceptor and donor, the nonideal behavior was subsequently predicted using a statistical model developed using the measured 18 mixtures. The theoretical approach used also allowed estimation of the fusion properties of two quaternary ammonium salts which decompose before melting.

In order to model the trends observed in these 18 interaction parameters, a principal component analysis (PCA) was performed on 48 descriptors quantifying various properties, including the hydrogen bonding capabilities, of the components in the mixtures. This resulted in four principal components which were used to identify the relevant properties to perform regression with the interaction parameters.

From the correlations of the variables with the interaction parameter χ, it follows that the bonding strength of the acidic proton of the hydrogen bond donating groups is important, where a larger dipole moment seems to indicate larger attraction. This bonding strength is measured by the proton affinity and the dipole moment. For the proton affinity, the effect is not as straightforward. In the case of the hydrogen bond acceptor, it shows that the ionic bond length between the bromide anion and the ammonium cation is an important indicator for interactions. Here it seems that a smaller ionic bond leads to stronger attractions; however, this trend is not observed in all cases. Next to this, it follows that the length of the aliphatic side chains of the quaternary ammonium salts is of high importance; also, here the trend is intuitive but again does not apply everywhere. One way to explain these observations is that, due to the shorter side chains, the distance between N+ and Br− decreases, and the Br− anion will experience a larger effective attraction with the −OH of the carboxylic acid or the polyol as compared to itself, which is in agreement with the literature.[8]

From the properties which are present in the model, some generalizations can be made, since these could be regarded as contributors to the enthalpy of mixing. However, as the model comprises nonlinearities, the exact nature of the effect of the properties cannot be derived from this. What is remarkable is that for each donor and acceptor two descriptors are present in the model. The descriptors used for the acceptor here are the distance between the anion and cation of the hydrogen bond acceptor and the length of the aliphatic side chain, thus providing information about the degree of hydrophilicity as well as the bulkiness. For the hydrogen bond donor, both the dipole moment as well as the proton affinity, both being well-known indicators for hydrogen bond donating abilities, are of importance to consider when looking for new nonideal eutectic mixtures. The fact that both HBAs and HBDs are represented in the model does indicate that the nonideality observed in these eutectic mixtures can not easily be ascribed to single properties of the components, for example, to the hydrogen bond donating or accepting ability of the components. It is more likely that the obtained liquid mixtures are the result of a synergistic combination of the properties identified here.

Methods

Experimental Methods

Determination of the Eutectic Temperature and Liquidus Curve

The eutectic temperatures of the different mixtures were measured using melting point capillaries. The DES compositions used were prepared inside a glovebox with a dry nitrogen atmosphere, yielding DES mixtures with moisture levels below 10 ppm. The temperature of the first liquid visible at a heating rate of 5 K·min–1 was taken for the solidus line (eutectic temperature). The liquidus curve was (if obtained) measured using DSC on the acid-rich side and using the cloud-point method on the salt-rich side.

DSC Method

A TA Q2000 Dynamic Scanning Calorimeter (DSC) was used to determine ΔH. This equipment has an accuracy of ±0.1 K for T and <1% for ΔH and a precision for T measured of ±0.4 K. Samples were measured in Tzero aluminum hermetic pans. All samples were prepared and pressed inside the glovebox using the mixtures prepared as described in the previous sections. The melting point was obtained by measuring samples of various compositions at three different heating rates (dT/dt), 10, 5, and 1 K·min–1. Prior to every heating cycle, the sample was heated isothermally for at least 10 min at 288 K in order to avoid any metastable intermediate phases.

Cloud-Point Method

A sample was prepared by mixing a known amount of both individual components in a glass vial using a magnetic stirring bar at a certain temperature. The sample was heated in an aluminum heating block, while stirring at 250 rpm, on an IKA RCT basic hot plate with an ETS D-5 temperature controller (accuracy ±0.5 K, precision 0.1 K). The temperature sensor used to control the temperature of the aluminum heating block was immersed in a stirred glass vial with glycerol placed in the aluminum heating block.

At the first chosen temperature, it was evaluated whether enough liquid had formed in order to sample only liquid. Prior to sampling, stirring was turned off and left off for at least 1 h to allow the solids remaining in the vial to settle. Sampling was performed with a micropipette, and care was taken to only sample liquid. After this, the temperature was raised in steps of approximately 4 K until 398 K, or until the melting point of the pure component has been reached, and the sampling procedure was repeated for each step.

The composition of the samples were analyzed using 1H NMR spectroscopy. For this, the samples were dissolved in approximately 1 mL of D2O (Cambridge Isotope Laboratories Inc., deuterium oxide 99.9%, DLM-4-100, CAS 7789-20-0). Measurements were carried out in 5 mm thin-walled Wilmad NMR tubes using a Bruker BZH 400/52 spectrometer equipped with an autosampler carousel. The spectra were autoshimmed and autophased and recorded in 16 scans with a relaxation time of 5 s by the TopSpin software. The obtained spectra were edited by applying a Withaker Smoother baseline correction and autophase correction using MestReNova 12.0.0-20080.

In order to determine the composition of the sample, the peaks corresponding to specific groups from both components were integrated and normalized to the amount of protons represented by the first peak appearing on the spectrum on the left.

Computational Methods (DFT)

All calculations were performed using the Vienna ab initio simulation package (VASP).[42−45] The Perdew–Burke–Ernzerhof exchange-correlation functional was used.[46] For the plane wave basis set, a cutoff energy of 400 eV was chosen. All molecules were placed at the center of a sufficiently large unit cell in such a way that the vacuum distance between the compounds in the neighboring super cells was at least 15 Å. It was carefully checked that the charge density approached zero at the border of the unit cell. For all calculations, only the gamma point was employed. Structure optimization was done using the conjugate-gradient algorithm as implemented in VASP. It was verified that after structure optimization the forces on the atoms were less than 0.02 eV/Å in any Cartesian coordinate. Subsequently, ground state energies as well as all of the distances between the available protons and their oxygen for the donors, and cation to anion distances for acceptors, were extracted using a Python script as available within the Atomic Simulation Environment (ASE) package. After the geometry optimization, the probe molecule was added. The probe molecule was either another molecule of the same component or a different molecule, namely, water or formic acid. The positioning of the probe molecule was done identically for all components: The accepting oxygen atom was placed 1.5 Å in front of the donor hydrogen. For HBAs, the donating hydrogen was positioned 2.5 Å from the accepting halide ion. After positioning, another geometry optimization was performed at the same level of theory and using the same methodology as described above. Proton affinities were calculated as the difference in the electronic energy between the components with and without a proton. For the donors, this was achieved by removing a proton from the optimized molecule while keeping the number of valence electrons the same. For the acceptors, a slightly different approach was taken. Since the quaternary ammonium salts do not have a free proton to remove and since the assumption was made that the hydrogen bonding takes place via the bromide anion, a proton was placed at 1.2 Å from the bromide anion in line with the nitrogen cation. To avoid reorganization of the alkyl side chains of the cation, the atoms corresponding to the alkyl side chain were frozen, whereas all other atoms remained unconstrained.