Mario Špadina

Klemen Bohinc

Multiscale modeling of solvent extraction and the choice of reference state: Mesoscopic modeling as a bridge between nanoscale and chemical engineering

Keywords: solvent extraction, multicomponent systems, reference state, complexation, aggregation, phase structuring

diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The sustainability of the chemical industry demands efficient separations and recycling methods for various metals [1,[START_REF] Binnemans | Recycling of rare earths: A critical review[END_REF][START_REF] Yoon | Solvent extraction, separation and recovery of dysprosium (Dy) and neodymium (Nd) from aqueous solutions: Waste recycling strategies for permanent magnet processing[END_REF]. Novel designs and improvements of already established separation processes are important for the nuclear industry as well, where the goal is to close the fuel cycle in new generations of nuclear systems [START_REF] Modolo | A review of the demonstration of innovative solvent extraction processes for the recovery of trivalent minor actinides from purex raffinate[END_REF][START_REF] Panak | Complexation and Extraction of Trivalent Actinides and Lanthanides by Triazinylpyridine N-Donor Ligands[END_REF][START_REF] Gelis | Actinide Lanthanide Separation Process -ALSEP[END_REF].

One of the most efficient methods used is liquid-liquid extraction (or solvent extraction) [START_REF] Yoon | Solvent extraction, separation and recovery of dysprosium (Dy) and neodymium (Nd) from aqueous solutions: Waste recycling strategies for permanent magnet processing[END_REF][START_REF] Jha | Review on hydrometallurgical recovery of rare earth metals[END_REF]. Solvent extraction works by the principle that the extractant molecule (a weak surfactant-like amphiphilic molecule) facilitates the transfer of ions from the aqueous to the organic phase. This means that even in the simplest laboratory experiment of metals extraction, the system is biphasic and multicomponent. Emerging complexity of solvent extraction systems, makes the optimization of the extraction process in chemical engineering both resource and time consuming [START_REF] Rydberg | Solvent Extraction Principles and Practice, Revised and Expanded[END_REF][START_REF] Nichols | Toward Mechanistic Understanding of Nuclear Reprocessing Chemistries by Quantifying Lanthanide Solvent Extraction Kinetics via Microfluidics with Constant Interfacial Area and Rapid Mixing[END_REF]. The situation is especially problematic since the solubilization and stabilization of ions in the organic phase is achieved in the form of the variety of self-assembled distinct structures and is mediated by supramolecular and colloidal interactions [START_REF] Cote | The supramolecular speciation: a key for improved understanding and modelling of chemical reactivity in complex systems[END_REF][START_REF] Zemb | Recycling metals by controlled transfer of ionic species between complex fluids: en route to "ienaics[END_REF][START_REF] Motokawa | A Telescoping View of Solute Architectures in a Complex Fluid System[END_REF][START_REF] Clark | Amphiphile-Based Complex Fluids: The Self-Assembly Ensemble as Protagonist[END_REF][START_REF] Pleines | A minimal predictive model for better formulations of solvent phases with low viscosity[END_REF]. The diversity in aggregation in the solvent phase (historically called pseudophases) with respect to the change of only two independent system variables is sketched as a simplified phase diagram in Figure 1 and as conceptual view of the nanoscale in Figure 3. A simple change of solutes concentrations in the aqueous or the extractant concentration in the organic phase can set the equilibrium state into the organized complex fluid where the variety of aggregates can coexist such as spherical, ellipsoidal, cylindrical, lamellar, worm-like, or other distinct assemblies [START_REF] Duvail | Mesoscopic modelling of frustration in microemulsions[END_REF][START_REF] Zemb | Recycling metals by controlled transfer of ionic species between complex fluids: en route to "ienaics[END_REF][START_REF] Liang | Nanoscale view of assisted ion transport across the liquid-liquid interface[END_REF] (see Figure 3). Furthermore, the organic phase can undergo an additional phase splitting (see Figure 1, left photograph of dense 'red sphere' inside the flask) [START_REF] Antonio | Structural insights into the multinuclear speciation of tetravalent cerium in the tri-n-butyl phosphate-n-dodecane solvent extraction system[END_REF]. The new third phase is now a dense network highly saturated by cations, whereas the light organic phase at is dilute and Photographs of the third phase, popularly called 'the red sphere' (left), as well as the normal biphasic solution system (right), were kindly provided by Dr. Mark Antonio from Argonne National Laboratory.

gas-like. This phenomenon reflects the complexity of solvent extraction systems since it can occur simply by increasing only the concentration of 'transferable' electrolyte in the aqueous phase. At this point it must be noted that not all electrolytes are readily transferred to the organic phase, or to be more exact, not every type of extractant can efficiently facilitate the transfer of ion to the organic phase. It is the wish of the authors to emphasize that the proper design of the experiment should always consider solute/extractant couple when thinking about forces and effects that govern the extraction, rather than focusing on ion-specific effects alone, or extractant interface partitioning and its affinity towards complexation [START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF].

A good example of how solute should be matched with a proper extractant are highly hydrated alkali cations for which extraction by non-ionic extractants is negligible [START_REF] Meridiano | Aggregation in organic solutions of malonamides: Consequences for water extraction[END_REF], but is far from negligible when cationexchanger extractants are utilized at low acidity conditions [START_REF] Marcus | Extractant class: Phosphoric and pyrophosphoric acids[END_REF]. The practical use, or the application of such salts in the case of extraction with non-ionic extractant is as a background salt (i.e. the salting-out agent) to set the activity of the aqueous phase, whereas for acidic extractant system it can be a target electrolyte, or even be the saponification agent which can lead to changes of the organic phase volume by increased water uptake [START_REF] Qi | Chapter 2 -extractants used in solvent extraction-separation of rare earths: Extraction mechanism, properties, and features[END_REF].

Even though separation chemistry is a mature field of science, with direct applications in various industries, the theoretical knowledge to predict actual extraction efficiency is still scarce. Yet numerous approaches exist, but they are scattered towards elucidating only particular phenomenon which occurs around a single region in high dimensional phase diagram [START_REF] Bauer | Liquid/liquid metal extraction: Phase diagram topology resulting from molecular interactions between extractant, ion, oil and water[END_REF].

To elucidate phenomena which occur at different length scales, the only choice when battling the mentioned complexity of solvent extraction systems is to employ multi-scale modeling. Scaling the sizes of studied systems can quantify average effects which can be then transferred on different theoretical approaches, with the final goal being a successful modeling of the extraction efficiency on the engineering level. The principle of multi-scale modeling is demonstrated in Figure 2.

In fact, the choice of the level of theory one wishes to utilize should follow the extent of the phenomenon, or certain effect that is of major influence on particular scale (both length or time scale of dynamic processes).

In principle, the approach can start at atomic scales where the coordination chemistry of cation and ligands in terms of the geometry is studied (see Quantum mechanics (QM) square in Figure 2). QM can be followed by the molecular dynamics (MD) simulations (see Classical simulations square in Figure 2) to test the colloidal stability of weak self-associated aggregate in explicit solvent. A 'climb' to larger scales where the top is the chemical engineering, leads over mesoscopic modeling where stable aggregates can further be approximated as pseudophases to reach analytical expressions for concentrations of species in the organic phase (see Mesoscopic modeling square in Figure 2). Speciation expressions obtained from mesoscopic modeling in the next step can be incorporated into the chemical engineering multiple equilibria approach to predict the actual extraction process efficiency. It must be noted that although in principle the multi-scale modeling procedure works, in practice this is not the case. Various theoretical studies are scattered and directed towards elucidating distinct experimentally observed phenomena, and the link between them often does not exist.

Therefore, it is only natural to present a paper which deals on multi-scale modeling in solvent extraction by building the model from basic cation/extractant complex and then including higher-order, more delicate soft-matter effects known in colloidal chemistry. Within the paper, we will constantly address current paradigms and emphasize challenges which are present to this date. At the same time we will highlight a much-needed updates in current modeling approaches. An important aspect that needs to be understood is that the choice of the reference state of the solubilized ion in the organic phase dictates the ability of the theoretical approach to predict the extraction efficiency. We will comment on this in detail before proceeding to practical examples of modeling. Note that the special interest will be devoted to the modeling of mostly non-selective collective effects behind solvent extraction, which arise from similarities in chemistry of lanthanides and actinides.

Consequently, this paper will not provide exhaustive literature overview on theoretical work on lanthanides and actinides extraction, but it will rather focus on common aspects of supramolecular and colloidal chemistry of such systems, and on ingredients which are needed to bridge nanoscale description of the organic phase and the chemical engineering.

For the practical systems we will focus on extractants which show great potential for the recycling of spent fuel in IV-generation nuclear systems which assume closed cycles [START_REF] Modolo | A review of the demonstration of innovative solvent extraction processes for the recovery of trivalent minor actinides from purex raffinate[END_REF]. Classes of non-ionic extractants will include mostly malonamides (such as N,N' -dimethyl-N,N' -dioctyl-2-(2-hexyloxy-ethyl)malonamide (DM-DOHEMA) in DIAMEX process, where lantahnides and actinides are coextracted in the first step, prior to the selective extraction like in the example of SANEX process), glycolamides (such as N,N,N',N' -tetraoctyldiglycolamide (TODGA), an alternative to DIAMEX), alkylamides etc [START_REF] Roberts | Treatment and recycling of spent nuclear fuel Actinide partitioning -Application to waste management[END_REF][START_REF] Panak | Complexation and Extraction of Trivalent Actinides and Lanthanides by Triazinylpyridine N-Donor Ligands[END_REF][START_REF] Modolo | A review of the demonstration of innovative solvent extraction processes for the recovery of trivalent minor actinides from purex raffinate[END_REF]. Class of acidic extractants whose progress in multi-scale modeling is reviewed here assumes mostly lipophilic analogues of phosphoric acid such as bis(2-ethyl-hexyl) phosphoric acid (HDEHP) [START_REF] Lumetta | Review: Solvent systems combining neutral and acidic extractants for separating trivalent lanthanides from the transuranic elements[END_REF][START_REF] Modolo | A review of the demonstration of innovative solvent extraction processes for the recovery of trivalent minor actinides from purex raffinate[END_REF]. Also, the colloidal aspect of the famous industrial PUREX process which utilizes tri-n-butyl phosphate (TBP) will be addressed. Especially interesting separation strategies for lanthanides and actinides are mixed non-ionic and acidic extractants systems, which at certain conditions can even exhibit synergistic behavior [START_REF] Lumetta | Review: Solvent systems combining neutral and acidic extractants for separating trivalent lanthanides from the transuranic elements[END_REF][START_REF] Gelis | Actinide Lanthanide Separation Process -ALSEP[END_REF].

Other classes and types of extractants which will not be discussed in this article can be found, among many others, in the following references [START_REF] Wilson | Solvent extraction: the coordination chemistry behind extractive metallurgy[END_REF][START_REF] Jha | Review on hydrometallurgical recovery of rare earth metals[END_REF][START_REF] Rydberg | Solvent Extraction Principles and Practice, Revised and Expanded[END_REF][START_REF] Clark | Coordination of Actinides and the Chemistry Behind Solvent Extraction[END_REF][START_REF] Panak | Complexation and Extraction of Trivalent Actinides and Lanthanides by Triazinylpyridine N-Donor Ligands[END_REF]. Solvent extraction of transition elements and the effect of the organic solvent (diluent) will not be discussed here as well to simplify the paper [START_REF] Qi | Chapter 2 -extractants used in solvent extraction-separation of rare earths: Extraction mechanism, properties, and features[END_REF]. Nonetheless, the concepts that will be discussed throughout the text can be generalized to many hydrometallurgical extraction processes.

The choice of the reference state of solubilized ion in the solvent phase matters

The Gibbs energy of transfer, ∆ tr G -• is the central quantity in separation science since it estimates the magnitude of collective effects responsible for the transfer of an ion from the aqueous to the organic phase. Still, there is ambiguity associated with it. ∆ tr G -• naturally depends on the choice of the reference state of an ion in the organic phase. This means that the calculated extraction efficiency, by means of theoretical chemistry, also depends on the choice of reference state. The authors deem that this small issue generates a barrier among different theoretical approaches in solvent extraction. Let us clarify. [START_REF] Guillaumont | Modeling selectivity in liquid/liquid extraction[END_REF]. Quantum mechanics square: the snapshot of the of density functional theory (DFT) calculations yielding the optimized structure of Nd 3+ coordinated with 6 dimethylphosphoric acid ligands (a truncated version of HDEHP with short alkyl chains) in a vacuum is reproduced from the reference [START_REF] Muller | Understanding the synergistic effect on lanthanides(III) solvent extraction by systems combining a malonamide and a dialkyl phosphoric acid[END_REF] with permission from Elsevier. Classical simulations square: the snapshot of MD simulations of stable Ln(NO 3 ) 3filled aggregate in n-heptane solvent is reprinted from the reference [START_REF] Chen | Stability of reverse micelles in rare-earth separation: a chemical model based on a molecular approach[END_REF] with permission from the Royal Society of Chemistry. The neutral polar core (cation + anions) which contains multiple H 2 O molecules is enclosed by 4 DMDOHEMA non-ionic extractant molecules. Mesoscopic modeling square: idea of the pseudophase-like models where monomeric extractants are in the equilibrium with already formed extractant film. Chemical engineering square: the concept of the chemical engineering multiple equilibrium approach. Every possible equilibria are established and constants are fitted to the experiments [START_REF] Dartiguelongue | Modelling of uranium(VI) extraction by D2EHPA/TOPO from phosphoric acid within a wide range of concentrations[END_REF].

Distances

We can express the transfer of lanthanide cation from the aqueous to the organic phase by chemical reaction M 3+ aq

M 3+ org (1)
where M 3+ aq and M 3+ org represent the hydrated cation in the aqueous and solvated cation in the organic phase, respectively. The corresponding standard Gibbs energy of transfer ∆ tr G -• can be written as the difference between standard states

∆ tr G -• = µ org,-• M 3+ -µ aq,-• M 3+ = -k B T ln a M 3+ org a M 3+ aq ( 2
)
where µ org,-• M 3+ and µ aq,-• M 3+ are standard chemical potentials of solvated cation in the organic and hydrated cation in the aqueous phase, respectively. a M 3+ org and a M 3+ aq are activities of the cation in the organic and aqueous phases, k B is the Boltzmann constant and T is the absolute temperature.

Note that ∆ tr G -• is constant if Eq. 2 is defined on the molarity concentration scale (1:1 ratio of volumes of the two phases is assumed, see analogy with Henry's law). If other concentrations scales are used, due to the differences in densities of the two phases, a volume correction factor should be added to ∆ tr G -• [START_REF] Moeser | The role of the concentration scale in the definition of transfer free energies[END_REF]. A correction factor is entropic in nature and originates from the change of the ion's accessible volume within the transfer to the solvent phase, characterized by different density. A correction factor should be added when, for example, the molality scale is used where same masses of the two phases have different volumes [START_REF] Marcus | The law of mass-action versus non-ideal behaviour in distribution equilibria[END_REF][START_REF] Moeser | The role of the concentration scale in the definition of transfer free energies[END_REF]. The factor can gain even more importance if several extraction/desextraction cycles with different volume ratios are conducted within the experiment [START_REF] Huang | Selective recovery of rare earth elements from ion-adsorption rare earth element ores by stepwise extraction with HEH(EHP) and HDEHP[END_REF].

Based on Eq. 2, to simplify measurement of ∆ tr G -• the following approximations are often made:

1. because of electrostatic interactions, the solvated cation in the organic phase is electroneutral, 2. in dilute aggregation regime (which most efficient extraction formulations are to avoid high viscosity [START_REF] Pleines | A minimal predictive model for better formulations of solvent phases with low viscosity[END_REF]) non-electrostatic inter-aggregate interaction potentials also vanish (e.g. repulsive overlap of hydrophobic films of extractant chains between two interacting aggregates) [START_REF] Ellis | Complexationinduced supramolecular assembly drives metal-ion extraction[END_REF].

With these assumptions, activity coefficients of aggregates and extractants in the bulk solvent are same, thus can be canceled and only the loss of configurational entropy of the Mass action law (MAL) needs to be taken into account. By neglecting activity coefficients, activities can be replaced by equilibrium concentrations. Consequently, Eq. 2 can be rewritten as

∆ tr G -• = -k B T ln c M 3+ org /c M 3+ aq
, where c M 3+ org and c M 3+ org are equilibrium concentrations of ions in both phases. This simplifies the procedure for the experimental estimation of ∆ tr G -• since concentrations are readily measured. The assumption (expressed by Eq. 2) considers that the standard scale for metal activities is the one of a dilute metal ion dispersed in a solvent defined as all the species that are not dilute. Thus for the aqueous phase, the solvent does not only correspond to the water molecules: it also takes into account the other ions (as the ones of the acid). For the organic phase, the solvent is the organic diluent together with the extractant molecules and also the coextracted acid. This choice of solvent defines the standard state and consequently the standard term µ org,-• M 3+ . In the limit of dilute metal ions, the activity ratio is actually constant because of Henry's law and the standard Gibbs energy of transfer ∆ tr G -• is properly defined. It must be noted that experimental and theoretical approaches work in the opposite way: former use equilibrium concentrations to estimate the difference in standard chemical potentials, whereas the latter calculates ∆ tr G -• to estimate concentrations at equilibrium.

While there is the benefit of measuring ∆ tr G -• with respect to the framework of Eq. 2, from the theoretical point of view it is almost impossible to calculate ∆ tr G -• since no information of the solvated ion in the organic phase is given.

Let us imagine that bare trivalent cation of radius a M 3+ , stripped from its hydration sphere, is transferred from the bulk aqueous (ε water ) to the solvent phase (characterized by low dielectric constant ε oil ). We can identify ∆ tr G -• on qualitative level with the Born energy given by the expression ∆G B = (1/ε oil -1/ε water ) q 2 M 3+ /8πε 0 a M 3+ , where q M 3+ is the cation charge number. Since ε water > ε oil , the first term in the expression for ∆G B is always positive, which indicates that the additional work needs to be made on the system in order to maintain the charged particle in low dielectric constant medium. Clearly, the transfer should be unfavorable. Yet, the spontaneous transfer does occur, which means there are other forces in play i.e. the influence of the extractant molecules. While acidic extractants neutralize cation's charge by direct complexation, in the case of non-ionic extractants the charge neutralization is assisted by coextraction of counterions (anions) into the solvent phase. In any case, the effect of the extractant needs to be properly included.

By stressing out the importance of extractant molecules and complex structure, let us rewrite Eq. 2 for the general case of extraction by non-ionic extractant, but this time to include all constituents in MAL:

N M 3+ • M 3+ + N A -• A -+ N L • L + N H2O • H 2 O A N M 3+ ;N A -;N L ;NH 2 O (3)
where M 3+ , A -, L, H 2 O, N M 3+ , N A -, N L , and N H2O are respectively symbols and stoichiometry numbers for the target cation, the anion, the extractant aggregation number, and water molecules.

A N M 3+ ;N A -;N L ;NH 2 O
is the general symbol for any type of the aggregate (or A agg ). Species present in the organic solvent are denoted by overlined symbols. Corresponding ∆ tr G -• can be written in same manner as before

∆ tr G -• = µ org,-• Aagg -N M 3+ µ aq,-• M 3+ -N A -µ aq,-• A --N L µ aq,-• L -N H2O µ aq,-• H2O = -k B T ln   a Aagg a N M 3+ M 3+ a N A - A -a N L L a NH 2 O H2O   (4)
Now, Eqs. 3 and 4 can be work with in terms of modeling. Thermodynamic cycles similar to the Born-Haber cycle can be utilized to roughly estimate whether the transfer of the ion will occur spontaneously or not [START_REF] Cao | Accurate quantum chemical modelling of the separation of eu3+ from am3+/cm3+ by liquid-liquid extraction with cyanex272[END_REF][START_REF] Sun | Quantum chemistry study on the extraction of trivalent lanthanide series by cyanex301: Insights from formation of inner-and outer-sphere complexes[END_REF][START_REF] Clark | Coordination of Actinides and the Chemistry Behind Solvent Extraction[END_REF]. Consequently, the dependency in terms of acid and extractant concentration is taken into account by the activities in Eq. 4 and the corresponding ∆ tr G -• is constant. On the other hand, in Eq. 2, this dependency is taken into account by the standard term because the solvent depends on the acid and extractant concentrations, so that the corresponding ∆ tr G -• depends on the acid and extractant concentrations. This distinction is crucial for any attempt to quantitatively predict the actual extraction efficiency. Note that in the case of ion-exchange extractants, additional activity factors (concentrations in the dilute case) appear in Eq. 4. For HDEHP or similar acidic extractants, pH dependence can be included by modeling H + activity in the aqueous phase, and the same concentration scale should be used for all constituents of MAL. Since ion-exchange occurs at the liquid-liquid interface, it is safe to assume that all H + are released in the aqueous phase. Still, in most occasions, extraction systems cannot be represented successfully with a simple equilibrium reaction. In reality, two major complications arise: Cation/extractant complexes are not isolated but interact via long-range interactions [START_REF] Servis | Square supramolecular assemblies of uranyl complexes in organic solvents[END_REF]. Furthermore, sometimes cation/extractant complexes are part of larger aggregates [START_REF] Bu | Observation of a rare earth ion-extractant complex arrested at the oil-water interface during solvent extraction[END_REF][START_REF] Ellis | In the bottlebrush garden: The structural aspects of coordination polymer phases formed in lanthanide extraction with alkyl phosphoric acids[END_REF][START_REF] Motokawa | A Telescoping View of Solute Architectures in a Complex Fluid System[END_REF]. Both examples interfere with the usual definition of the activity of complex in the solvent phase. In non-dilute state, short-range soft interactions and solvent exclusion effects become more important.

Stoichiometry deduced by experiments varies implying that there is extractants (ligands) exchange between first-and second-coordination sphere, depending on the physical conditions of the extracting system [START_REF] Weßling | Stoichiometry of an(iii)-dmdohema complexes formed during solvent extraction[END_REF][START_REF] Brigham | Trefoil-shaped outer-sphere ion clusters mediate lanthanide(iii) ion transport with diglycolamide ligands[END_REF][START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. The aggregation number is the true stoichiometry (it gives correct speciation in the solvent phase) and it differs from the complexation number (cation's first-coordination sphere) [START_REF] Zemb | Recycling metals by controlled transfer of ionic species between complex fluids: en route to "ienaics[END_REF][START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF]. This was demonstrated in Figure 3. where the complex can be well-defined with few water molecules and ligands in its first coordination sphere (follow the orange arrow), but can also have any other distinct composition. Furthermore, strong experimental evidence suggests that aggregation numbers are larger than the complexation numbers: the number of ligands directly bound to central cation [START_REF] Zemb | Recycling metals by controlled transfer of ionic species between complex fluids: en route to "ienaics[END_REF].

In fact, these discrepancies between the aggregation and the complexation numbers were emphasized in the work by Y. Marcus more than 50 years ago, long before modern scattering techniques were employed in separation science [START_REF] Marcus | Solvent extraction of inorganic species[END_REF]. A valuable discussion concerning the strict definition of the distribution coefficient and its validity in terms of slope method and Job's plots, is given in the aforementioned reference. Even today, theories of activity coefficients in the solvent phase, but also in the case of mixed concentrated electrolytes in the aqueous phase, are not yet mature for scaling to the engineering level [START_REF] Marcus | Solvent extraction of inorganic species[END_REF][START_REF] Marcus | The law of mass-action versus non-ideal behaviour in distribution equilibria[END_REF].

Therefore, neglecting the second-coordination sphere is already a large discrepancy between the predictions of coordination chemistry and actual experiment. Considering the fact that actual industrial processes utilize so-called phase modifiers such as 1-octanol, which in most cases constitute the second-coordination sphere further adds to the complexity [START_REF] Sukhbaatar | Probing the existence of uranyl trisulfate structures in the amex solvent extraction process[END_REF]. Then if we allow the concept to expand to emerging polydispersity by allowing the various cation/ligand complexes as a building units to assemble into larger aggregates where solute fluctuations are energetically less penalized (see cylindrical or lamellar structures in Figure 3.), then simple conclusion can be drawn: choosing the isolated complex as a reference state in solvent phase cannot be used to accurately predict extraction efficiency.

By looking at the conceptual representation of the solvent phase in lanthanide extraction as a complex organized fluid in Figure 3., it becomes unclear what is the actual state (or the environment) of the solubilized ion. Strong experimental evidences do support this conceptual representation, along with some of theoretical approaches that deal on stabilization of ions in apolar medium [START_REF] Zemb | Recycling metals by controlled transfer of ionic species between complex fluids: en route to "ienaics[END_REF][START_REF] Kanduč | Aqueous nanoclusters govern ion partitioning in dense polymer membranes[END_REF]. An important and predominant effect is that all complexes are neutral because of low dielectric constant in a hydrocarbon solvent, thus the Born energy of the transfer is lowered. Another effect can be understood as the formation and swelling of polar nanodomains in solvent phase by increased water uptake. This global effect of ion stabilization in apolar medium is achieved by increasing dielectric constant within the dispersed hydrophilic nanodomains that resemble bulk 'water pools' [START_REF] Kanduč | Aqueous nanoclusters govern ion partitioning in dense polymer membranes[END_REF].

If one chooses the reference state of solubilized cation as a complex then ∆ tr G -• will include mostly the difference in first-coordination spheres between a complex and hydrated cation, but every aspect of the colloidal domain will be neglected. Consequently, quantitative prediction of the extraction efficiency is questionable and often differs substantially from experimental values [START_REF] Vo | Quantum chemistry benchmarking of binding and selectivity for lanthanide extractants[END_REF][START_REF] Jackson | Synergism and aggregation in multi-extractant solvent extraction systems[END_REF].

On the other hand, if the self-assembled aggregate is considered as whole (i.e. the standard state is the aggregate at infinite dilution), additional effects beyond coordination chemistry can be accounted for, such as diluent effects, creation of second coordination sphere by extractants, electrostatic effects, activity of organic phase when completing the MAL, etc. These effects will be discussed in following sections. It must be noted that calculation of the free energy of aggregate at infinite dilution is computationally much more demanding, and at larger scales practically even impossible due to the existence of large structures and their huge co-existing variety. A compromise, or the savior, in this case, is a coarse-grained multi-scale approach with mesoscopic modeling serving as a bridge.

Ion transfer understood as a free energy difference between the complex in the solvent and hydrated ion in the aqueous phases

At the atomic length scales, the assembly made of cation and ligands can be studied as supramolecular complex (see Quantum mechanics simulations square in Figure 2). Lanthanides are highly hydrated complexes in the aqueous phase, with coordination numbers being mostly 9 in dilute solutions [START_REF] Marie | Complexation of lanthanides(iii), americium(iii), and uranium(vi) with bitopic n,o ligands: an experimental and theoretical study[END_REF]. In the organic phase, the cation is usually coordinated by extractants in the first coordination sphere with maximum occupancy (coordination number) and the complex geometry being determined by their sterical aspect. At the same time, the role of water molecules is unclear in the sense of coordination chemistry. The only general rule is that water molecules occupy both first and second coordination spheres [START_REF] Baldwin | Outer-Sphere Water Clusters Tune the Lanthanide Selectivity of Diglycolamides[END_REF]. The extent of water partition in the first coordination sphere is, among others, also the consequence of the steric aspect of the particular extractant molecule. Similar conclusions hold for the case of coextracted anions (counterions) by non-ionic extractants to make the complex electroneutral overall [START_REF] Brigham | Trefoil-shaped outer-sphere ion clusters mediate lanthanide(iii) ion transport with diglycolamide ligands[END_REF][START_REF] Acher | Inner to outer-sphere coordination of plutonium(iv) with n,n-dialkyl amide: influence of nitric acid[END_REF].

Therefore, to calculate the efficiency of the extraction, the first attempt is naturally to calculate ∆ tr G -• between the complexed atom in the organic and the hydrated cation in the aqueous phases (a simplified attempt). The sketch of this procedure is presented in Figure 3 with the orange arrow pointing to the change of the first-coordination sphere composition, where cation's coordinated water molecules in the aqueous phase are replaced by extractant electro-negative atoms in the organic phase.

Transfer energies and selectivities by calculating thermodynamic cycles via QM simulations

The most accurate tool to study the coordination chemistry of solvated complexes at the atomic scales are QM calculations. QM can be used to calculate the transfer energies of certain steps and to complete thermodynamic cycles which give ∆ tr G -• for the transfer of cation. The actual extraction efficiency can then be estimated via Eq. 4. Even if the extraction system can be described by simple chemical equilibrium reaction (as described by Eq. 4, where there are no ionizations of components), the approach requires to undertake multiple steps in calculations. The scheme includes mostly evaporation of solutes to reach gas phase, recombination of species in the gas phase and then solvating these new species from the gas phase. As a result, when all the energies of individual steps within the parent thermodynamic cycles are summed, the binding energies are obtained [START_REF] Vo | Quantum chemistry benchmarking of binding and selectivity for lanthanide extractants[END_REF][START_REF] Sun | Quantum chemistry study on the extraction of trivalent lanthanide series by cyanex301: Insights from formation of inner-and outer-sphere complexes[END_REF].

At this point, it is worth to highlight that QM modeling is intimately bonded with chemical engineering approach [START_REF] Moeyaert | Modeling and speciation study of uranium(vi) and technetium(vii) coextraction with dehiba[END_REF], even though the two approaches are on the exact opposite extremes of multi-scale modeling (see QM in the lower-left corner and Chemical engineering square in the upper right corner in Figure 2). Chemical engineering approaches start by establishing every possible equilibrium reaction, which partitioning in the global MAL is described by associated constants. Not all chemical reactions can be described by means of QM, but some can be, like for example the oligomerization or the dissociation of acidic extractants. Furthermore, the identification of most stable complexes in the organic phase by QM calculations is then transferred to the chemical engineering to shorten the set of chemical reactions included within multicomponent MAL i.e. to include only the most relevant complexes and to reduce the number of constants [START_REF] Moeyaert | Modeling and speciation study of uranium(vi) and technetium(vii) coextraction with dehiba[END_REF].

Formations of homo-and heteromers in the gas phase or in the continuum solvent can be assessed with high accuracy [START_REF] Pecheur | Synergism in a HDEHP/TOPO Liquid-Liquid Extraction System: An Intrinsic Ligands Property?[END_REF][START_REF] Gullekson | Extraction of Water and Speciation of Trivalent Lanthanides and Americium in Organophosphorus Extractants[END_REF][START_REF] Yu | Prediction of solubility properties from transfer energies for acidic phosphorus-containing rare-earth extractants using implicit solvation model[END_REF]. Deprotonation of acidic extractant is a very important aspect as well, for both mesoscopic and later chemical engineering approaches. DFT calculations of common industrial organophosporic acid extractants such as HDEHP and similar molecules yielded an ionization energy trend [START_REF] Jing | Extraction behaviors of heavy rare earths with organophosphoric extractants: The contribution of extractant dimer dissociation, acid ionization, and complexation. a quantum chemistry study[END_REF]. The values of calculated ionization energies in organic solvents are around 16 k B T when a sum of gas phase ionization and corresponding solvation energy (continuum solvent description) is considered. Ionization free energy was identified as a extraction-quenching contribution, along with extractant dimer formation. On the other hand, a strong complex formation between deprotonated acidic extractants and trivalent cations compensates all quenching terms. The conclusion was made that the dominant influence that drives the transfer of ion is the complexation energy, particularly enthalpy driven, but ionization is an important effect as well, and much more pronounced than the effect of extractants oligomerization in the organic phase [START_REF] Jing | Extraction behaviors of heavy rare earths with organophosphoric extractants: The contribution of extractant dimer dissociation, acid ionization, and complexation. a quantum chemistry study[END_REF].

Still, the true power of using QM simulations lies in the modeling of selectivity across lanthanide series, or selectivities of lanthanides versus actinides [START_REF] Cao | Accurate quantum chemical modelling of the separation of eu3+ from am3+/cm3+ by liquid-liquid extraction with cyanex272[END_REF]. The principle exploits cancellation of similar terms in the two thermodynamic cycles, which differ only by the type of the cation of same the valence number [START_REF] Marie | Complexation of lanthanides(iii), americium(iii), and uranium(vi) with bitopic n,o ligands: an experimental and theoretical study[END_REF][START_REF] Clark | Coordination of Actinides and the Chemistry Behind Solvent Extraction[END_REF]. The principle is consistent with statistical thermodynamics, as we will demonstrate by the following argument.

Each MAL is characterized by the corresponding system partition function. In the canonical ensemble, upon factorizing the system partition functions of two MALs, it can be concluded that the selectivity factor is nothing more than the ratio of electronic (or vibronic) partition functions. Practically, selectivity calculated in this way originates from the different cation/ligand interactions in the first coordination sphere. Such calculation supposes that the long-range structure of complexes is the same for all lanthanides and all actinides. To clarify, let us express two competing transfer of cations (let us name them: 1 and 2) as the double difference in Gibbs transfer energies, as in reference [START_REF] Zemb | Recycling metals by controlled transfer of ionic species between complex fluids: en route to "ienaics[END_REF]. The selectivity factor S 1-2 can be then written as

S 1-2 = ∆ tr G -• 1 -∆ tr G -• 2 .
Note that written in this way, the selectivity factor is expressed through Gibbs potential, but the link with statistical thermodynamics is through Helmholtz potential. Still, the Legendre transform does not need to be performed because resulting P V terms are automatically canceled in the ratio, which further simplifies the calculation. S 1-2 can be written as

S 1-2 = Z system 1 Z system 2 (5)
where Z system 1 and Z system 2 are system partition functions of the transfer of cation 1 and 2, respectively. After factorization of Z system 1 and Z system 2 and cancellation of similar terms, only delicate, first-coordination sphere effects give the selectivity.

Current paradigm of coordination numbers and issues concerning QM simulations in solvent extraction

A majority of QM-based works used truncated analogues of extractant molecules (shortening of the extractant side chains) in calculations, where the focus is on the description of only short-range interactions between central cation and coordinating electronegative atoms of the extractant head groups. Even if all atoms of extractants are taken as an input (corresponding to a full description of side chains), running simulations in vacuum or at continuum solvent level, yields counter-empirical idealization of many collective effects [START_REF] Roberts | Treatment and recycling of spent nuclear fuel Actinide partitioning -Application to waste management[END_REF][START_REF] Rey | Synergy in Extraction System Chemistry: Combining Configurational Entropy, Film Bending, and Perturbation of Complexation[END_REF].

A major ambiguity in separation science of lanthanides and actinides is the stoichiometry of the transfer involved with the extraction. Apart from the fact that complexation and aggregation numbers differ substantially, it is also unclear whether the coordination number and the geometry of the complex itself are preserved in the solvent phase [START_REF] Chen | Stability of reverse micelles in rare-earth separation: a chemical model based on a molecular approach[END_REF]. A generally understood mechanism is the following: the high coordination hydrated lanthanides (lanthanide is coordinated by 9 water molecules) is lost in the solvent phase, where less sterically flexible extractant molecules preferentially constitute a complex with octahedral symmetry. The entropy loss is thus compensated by more favourable enthalpic contribution, preferentially electron-donor in nature.

Due to the structural diversity of extractants which are designed for lanthanides or actinides extraction, it is difficult to state a general rule of preferred geometry of complexes, especially when branching of extractants is accounted for (increased branching generates additional steric hindrance). While it is clear that the complexation free energy of coordination chemistry is the leading force for the transfer of ion to the organic phase, it is not simple to claim whether if the control is enthalpic or entropic. Just like in coordination chemistry of aqueous solutions, the enthalpy or entropy control is the consequence of the type of the ligand (electron donor and acceptor atoms, steric effects, etc.). In that sense, it can happen that complexation and ion transfer to the organic phase is under strong enthalpic control, as in the case of HDEHP extractant [START_REF] Jing | Extraction behaviors of heavy rare earths with organophosphoric extractants: The contribution of extractant dimer dissociation, acid ionization, and complexation. a quantum chemistry study[END_REF], or the control can be entropic, as was demonstrated in the case of fluorinated aryldithiophosphinates [START_REF] Pu | sweeping" ortho substituents drive desolvation and overwhelm electronic effects in nd3+ chelation: A case of three aryldithiophosphinates[END_REF].

Addition of the colloidal aspect to solubilized ion: an emergence of weak self-assembled aggregates

In order to include higher-order effects into the estimation of the extraction efficiency one must place a complex optimized by means of QM into a much larger system. Classical simulations such as MD or Monte Carlo are obvious tools for this action, and this is exactly the next step in the multi-scale modeling approach, as demonstrated in the Classical simulations square in Figure 2. Classical simulations complemented by scattering techniques can be used to probe the solvent phase structuring at nanoscale level [START_REF] Duvail | Uo22+ structure in solvent extraction phases resolved at molecular and supramolecular scales: a combined molecular dynamics, exafs and swaxs approach[END_REF][START_REF] Liang | Nanoscale view of assisted ion transport across the liquid-liquid interface[END_REF]. While for some cases distorted hydrophilic nanodomains filled with ions are formed [START_REF] Qiao | How hydrogen bonds affect the growth of reverse micelles around coordinating metal ions[END_REF][START_REF] Duvail | Uo22+ structure in solvent extraction phases resolved at molecular and supramolecular scales: a combined molecular dynamics, exafs and swaxs approach[END_REF], for other cases distinct morphologies that resemble spherical [START_REF] Duvail | The role of curvature effects in liquid-liquid extraction: assessing organic phase mesoscopic properties from MD simulations[END_REF], or cylindrical [START_REF] Qiao | Liquid worm-like and proto-micelles: water solubilization in amphiphile-oil solutions[END_REF] reverse aggregates were identified. In the case of uranyl (UO 2+ 2 ) extraction by dialkylamide extractants, the complexity of aggregation was displayed in full and pseudophase diagram showed dozen of different aggregation regimes simply by changing UO 2+ 2 concentration in the aqueous, or the initial dialkylamide concentration in solvent phase [START_REF] Pleines | A minimal predictive model for better formulations of solvent phases with low viscosity[END_REF]. The situation got even more complicated when various diluent and extractant side chains were substituted in the study [START_REF] Pleines | A minimal predictive model for better formulations of solvent phases with low viscosity[END_REF]. Indeed these weak self-assembled aggregates can represent a most stable state of the system (this is not valid for the case of third phase formation), thus the inspiration for the Figure 3. which depicts the complexity of the structured solvent phase loaded with solutes. It must be noted that mentioned morphologies are observed when extraction systems are either in: a) low aggregation regime, i.e. dilute gas of non-interacting aggregates, or b) when type and concentration of extracted solutes do not favor the formation of dense networks governed by strong ion correlations e.g. the third phase formation.

Long-range interactions and hydration stabilize aggregates in the solvent phase

Rather than emphasizing particular examples within the diversity in aggregation, which on the first glance contradict each other, we will focus on highlighting what is common to all of them. In all cases, long-range interactions, electrostatic in nature, stabilize complex fluids [START_REF] French | Long range interactions in nanoscale science[END_REF][START_REF] Zemb | Weak aggregation: State of the art, expectations and open questions[END_REF][START_REF] Ivanov | Organic and third phase in hno3/tbp/n-dodecane system: No reverse micelles[END_REF]. Weak self-assembly is contrary to a huge loss in entropy which is needed to displace a large number of extractants from bulk and to incorporate them to the emerging pseudophase. Electrostatic forces in the low-dielectric constant hydrocarbon solvent can span over one nanometer in distance [START_REF] Qiao | How hydrogen bonds affect the growth of reverse micelles around coordinating metal ions[END_REF]. In the system of non-ionic DMDOHEMA extractant which complexes Eu 3+ and co-extracts nitric acid, the combination of hydrogen donor and acceptor, as well strong correlations among repeating Eu 3+ -NO - 3 contact ion pairs, can bridge reverse micelles and collapse them into larger aggregates [START_REF] Qiao | How hydrogen bonds affect the growth of reverse micelles around coordinating metal ions[END_REF]. In the case of UO 2+ 2 /TBP complexes, the two peaks in pair distribution functions are positioned at 6 Å and 12-15 Å, respectively [START_REF] Servis | Square supramolecular assemblies of uranyl complexes in organic solvents[END_REF]. While the first peak depicts strongly ordered spatial assembly, the second peak is consistent with more flexible, loose aggregation in the solvent phase. Longrange interactions affect the activity of cation/extractants complex, thus simplification of the effect cannot be reduced to sole cancellation of activity coefficients within the closure of the MAL (as is often done in chemical engineering). The system is made of two non-miscible liquids, namely organic solvent (pale yellow color), and the aqueous phase (pale blue color). Red circles represent target solutes e.g. lanthanide cation. Zoomed region depicts sketch the nanoscale 'view' of the system, where the equilibrium state is composed of a myriad of coexisting structures (right side). The issue with the choice of the standard state in the solvent phase is demonstrated as a transfer of hydrated cation from the aqueous to the polar core of self-assembled aggregate (orange arrow). a) In the case that solvated cation is understood as a complex which assumes the coordination of the cation by extractant molecules, counterions and coextracted water molecules, then ∆trG -• accounts only the first-coordination sphere interactions. b) If the solvated ion is the whole the aggregate, then beside interactions between first neighbors are accounted for. Additional collective effects include the formation of the extractant second-coordination sphere, hydrophobic extractant chains configurations, an organic solvent wetting and depletion effects of extractant film, the activity of the aggregate itself, etc.

V V V V V V V V V V V V V V V V V ∆ tr G • ═ μ • org -μ •
So far, we have identified few forces that are responsible for the transfer of salts: strong first coordination sphere complexation as the leading motor, and electrostatic coupling on the colloidal scale. Loss of entropy clearly opposes the self-assembly, but it is not the only such effect. The steric properties of extractants, their chain conformation and head groups crowding around the target cation restrict the extent of shortrange interactions, which gives rise to delocalization of charge even within the hydrophilic nanodomain itself. The idea of the balance between electrostatic and steric interactions is consistent with the theoretical study on non-ionic diglycolamide extractant/lanthanide system, where anions are positioned in the secondcoordination sphere to lower sterical repulsions, while stabilizing the complex by charge neutralization via long-range interactions [START_REF] Brigham | Trefoil-shaped outer-sphere ion clusters mediate lanthanide(iii) ion transport with diglycolamide ligands[END_REF]. Even though counterions such as NO - 3 or Cl -coextracted with lanthanide cation overall resemble the associated salt, simulations have shown that while the cation is in close contact with ligand (demonstrated on the case of diglycolamide extractants), counterions are at a larger distance (around 4.9 Å for NO - 3 and 4.86 Å for Cl -) but still keeping the aggregate electroneutral [START_REF] Brigham | Trefoil-shaped outer-sphere ion clusters mediate lanthanide(iii) ion transport with diglycolamide ligands[END_REF]. In the case of malonamides extracting UO 2+ 2 with sulfate SO 2- 4 as counterion, distances are smaller (average distance is around 3.2 Å) due to the stronger electrostatic interactions [START_REF] Sukhbaatar | Probing the existence of uranyl trisulfate structures in the amex solvent extraction process[END_REF]. Crowding by SO 2- 4 and malonamides in the first coordination sphere expelled water molecules to the second coordination sphere.

The absence of explicit water molecules in the aggregate makes the understanding somewhat inconsistent with 'water pool' continuum picture, but that is sometimes a consequence of the initial structure used in calculations (ligands were kept frozen in the DFT optimized geometry of the complex with counterions during MD simulations). The follow-up study which included MD simulations showed severe water uptake to stabilize the complex. Furthermore, it was found that the dynamics of complexed water onto the uranyl cation is on the scale of hundreds of nanosecond, which puts it out of reach for MD -but the effect exists [START_REF] Duvail | Uo22+ structure in solvent extraction phases resolved at molecular and supramolecular scales: a combined molecular dynamics, exafs and swaxs approach[END_REF]. Indeed, an abundant water molecules presence in second-coordination sphere does give hope for treating these aggregates like reverse micelles with 'water pool' inside the polar core, as was shown also in studies on non-ionic DMDOHEMA [START_REF] Chen | Stability of reverse micelles in rare-earth separation: a chemical model based on a molecular approach[END_REF] and TBP extractants [START_REF] Singh | Insight into acidity driven third phase formation of tbp in organic solutions by md simulation[END_REF]. The stability of Ln 3+ /DMDOHEMA aggregate is always affected by the number of water molecules in Ln 3+ first coordination sphere. The aggregate was represented like a molecular solute, where umbrella sampling allowed the calculation of the McMillan-Mayer (MM) potential, which is consistent with the Gibbs free energy, thus with calculations of successive 'hydration' stability constants. Addition of water molecule to the contact with Ln 3+ stabilized aggregate for around 0.5 k B T , while the optimum number of water molecules in the first coordination sphere was found 4. Interestingly, calculations with higher initial numbers of water molecules did not cause cleavage of the aggregate, only the migration of coordinated water to the second coordination sphere even at distances from 5-8 Å. The core of the aggregated showed the character of the liquid core, with radius up corresponding to the edge of the second coordination sphere. It must be noted that the calculated MM stability map versus the successive hydration of the core of the aggregate showed a wide flat area of negative values of MM potential, which indicates that considerable water uptake to the solvent phase can be explained by population of second coordination spheres around central cation. Even in the case of acidic HDEHP extractant, a rigorous study shows that there is substantial water partition in the organic phase, contrary to the opinion that dissociated HDEHP can fully screen electrostatics of trivalent cations thus stabilizing the aggregate [START_REF] Gullekson | Extraction of Water and Speciation of Trivalent Lanthanides and Americium in Organophosphorus Extractants[END_REF]. Monte Carlo (MC) studies were used to describe systems that resemble connected dense networks, or bicontinuous emulsions [START_REF] Mu | A novel microemulsion phase transition: Toward the elucidation of third-phase formation in spent nuclear fuel reprocessing[END_REF]. A difficulty concerning classical simulations is because solvent extraction systems usually demand longer time scales for calculations, along with a large number of atoms reproduce the aggregation effects. Huge inhomogeneity of the solvent phase in the actual extraction process implies that simulated clustering or condensation of a molecular solution are sometimes not representative of the entire system (see concept Figure 5). Yet, the new coupled MC and MD approach have shown that phase instabilities can be predicted. Parallel simulations in two boxes connected by the exchange of molecules hinted a phase separation where the dilute phase became acid and water-poor, while rich with a hydrocarbon solvent. The other, dense phase showed an increase in solutes concentration and organic solvent depletion. Results were in agreement with previous studies on third phase formation, where network analysis of hydrogen bond clusters was employed [START_REF] Servis | Network analysis and percolation transition in hydrogen bonded clusters: nitric acid and water extracted by tributyl phosphate[END_REF], as well as with the some of mesoscopic models of bicontinuous water/surfactant/oil microemulsions [START_REF] Duvail | Mesoscopic modelling of frustration in microemulsions[END_REF].

4.2.

Quantifying diluent effects and extractant chain conformations: from classical simulations to mesoscopic modeling Within the multiscale approach, classical simulations can be used to access parameters needed for later mesoscopic modeling. Especially important are solvent (diluent) effects and the branching of the extractant side chains. It is well known that effects both can radically change the efficiency of solutes transfer to the organic phase [START_REF] Rey | Synergy in Extraction System Chemistry: Combining Configurational Entropy, Film Bending, and Perturbation of Complexation[END_REF][START_REF] Knight | Subtle Effects of Aliphatic Alcohol Structure on Water Extraction and Solute Aggregation in Biphasic Water/n-Dodecane[END_REF]. Intermolecular dipolar forces between constituents of polar cores of different coexisting aggregates give rise to around 1 k B T in interaction energy, which is the same order of magnitude as the solvent effects on extractant's side chains [START_REF] Qiao | Molecular origins of mesoscale ordering in a metalloamphiphile phase[END_REF][START_REF] Guilbaud | Depletion of water-in-oil aggregates from poor solvents: Transition from weak aggregates towards reverse micelles[END_REF]. Both solvent depletion and the average extension of extractant chains are controlled by entropy [START_REF] Guilbaud | Depletion of water-in-oil aggregates from poor solvents: Transition from weak aggregates towards reverse micelles[END_REF].

By means of MD simulations, the conformation of the aggregated DMDOHEMA extractant lipophilic side chains was tracked as a function of the simulation time. Simulations of stable aggregate in vacuum yielded large variations of the average length of DMDOHEMA chains, whereas in explicit n-heptane solvent, chains maintained average extension of approximately 80% and fluctuations in lengths were small [START_REF] Duvail | The role of curvature effects in liquid-liquid extraction: assessing organic phase mesoscopic properties from MD simulations[END_REF]. A near-spherical aggregates were approximated as spherical ones to establish the curvature functionality. This resulted in a link between the chemical nature of extractant molecules which assemble into aggregates and parameters that describe the curvature of theses structures. By fitting the potential of mean force (PMF) as a function of the curvature of the spherical aggregate (near-spherical but approximated as spherical) with harmonic approximation yielded the so-called bending constant κ of the curved extractant film. The minimum of PMF was identified as the spontaneous packing parameter p 0 which is defined as the most stable packing of the extractant in the extractant film [START_REF] Duvail | The role of curvature effects in liquid-liquid extraction: assessing organic phase mesoscopic properties from MD simulations[END_REF]. Note that these parameters are both extractantand diluent-specific. The explicit values of these parameters were used later within mesoscopic modeling to successfully model the organic phase as an ensemble of distinct spherical aggregates, which resulted in the accurate prediction of the laboratory extraction experiment [START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF].

Current paradigms and issues in classical simulations of solvent extraction

Currently, there are two opinions on how to model the extraction systems. Is the solvent phase a complex fluid of interacting aggregates, or is it a molecular solution? In fact, both opinions have the same premises thus are more complementary than opposing! In all cases, long-range interactions lead to condensation of cation/anion species to fulfill the electroneutrality condition. Depending on the state of the system (in terms of all independent system variables!), as well as the nature of solute/extractant couple, the condensation may be in the form of weak self-assembled aggregate, but can be collapsed dense network, i.e. the third phase, as was demonstrated in work which dealt on N,N -dialkyl amides/water/uranyl extraction system. To bring another dimension into complexity, it was showed that the increasing steric hindrance of extractant molecules reduces the extent of the aggregation in solvent phase, which is especially important when metal cations are not present in the system [START_REF] Paquet | Aggregation in organic phases after solvent extraction of uranyl nitrate: X-ray scattering and molecular dynamic simulations[END_REF]. Dipole-dipole interactions mediate weak aggregation, whereas stronger first-cooridination sphere complexation can stabilize larger polymetallic aggregates similar to those found for Ln 3+ /HDEHP system [START_REF] Ellis | In the bottlebrush garden: The structural aspects of coordination polymer phases formed in lanthanide extraction with alkyl phosphoric acids[END_REF]. Thus it seems that the better question is: how to model these attractive interactions in hydrocarbon diluent and can the used theory provide reliable outcome? These attractive interactions resemble Keesom interactions between permanent dipoles but are obviously more important. Some ideas on how to model these interactions can be found in the following reference [START_REF] Parsegian | Van der Waals Forces: A Handbook for Biologists[END_REF].

This further brings us to issues concerning classical simulations of solvent phases. Practical applications in both hydrometallurgy and nuclear industry assume intensive usage of acidic effluent. Therefore, acid needs to be modeled accordingly. Still, the modeling the dissociation of acid is a common problem in MD (solvable, but at high computational cost!), which is why associated acid molecules are used in calculations [START_REF] Qiao | How hydrogen bonds affect the growth of reverse micelles around coordinating metal ions[END_REF][START_REF] Servis | Network analysis and percolation transition in hydrogen bonded clusters: nitric acid and water extracted by tributyl phosphate[END_REF][START_REF] Singh | Insight into acidity driven third phase formation of tbp in organic solutions by md simulation[END_REF]. This sometimes generates a practical issue: a high concentration of associated acid in the simulation box sometimes overestimates the extent of hydrogen bonding networks formation. Consequently the system collapses into dense phase (often identified as the third phase) via mechanism similar to percolation [START_REF] Antonio | Structural insights into the multinuclear speciation of tetravalent cerium in the tri-n-butyl phosphate-n-dodecane solvent extraction system[END_REF][START_REF] Servis | Network analysis and percolation transition in hydrogen bonded clusters: nitric acid and water extracted by tributyl phosphate[END_REF]. This problem was recently challenged by a study which compared results of MD simulations on a TBP/HNO 3 aggregation system. Two cases represented namely, the usual simulations with associated HNO 3 molecules, and the simulations with H 3 O + and NO - 3 ions, representing the dissociated HNO 3 [START_REF] Singh | Insight into acidity driven third phase formation of tbp in organic solutions by md simulation[END_REF]. Simulations with associated HNO 3 molecules resulted in the third phase formation, whereas simulations with ions resulted in the dispersion of hydrophilic nanodomains, which resemble reverse micelles [START_REF] Singh | Insight into acidity driven third phase formation of tbp in organic solutions by md simulation[END_REF]. Without further analyzing these results, the two distinct behaviors of solvent phase are again an example of the interplay between long-range electrostatic interactions and entropy at a given temperature.

Sometimes the usual water clusters that explain high dielectric constant in the bulk cannot be formed entirely, due to the small sizes of hydrophilic cavities [START_REF] Schlaich | Water dielectric effects in planar confinement[END_REF][START_REF] Kanduč | Aqueous nanoclusters govern ion partitioning in dense polymer membranes[END_REF]. Practically, in classical simulations, cut-off values of interacting potentials are fixed to the value which gives the best explanation of the aggregation phenomenon at hand. Problem with these values is that they are seldom transferable to other slightly different systems. Another issue is the tabulation of force fields. Due to polydispersity in aggregation, it is difficult to decouple the effects of specific species versus myriad of other coexisting ones, which complicates the fitting procedures needed to obtain reproducible force fields. Way to overcome these issues is systematical employment of joint DFT and MD simulations with extended X -ray absorption fine structure (EXAFS) and small-and wide-angle X -ray scattering (SWAXS) techniques [START_REF] Ferru | Elucidation of the structure of organic solutions in solvent extraction by combining molecular dynamics and X-ray scattering[END_REF][START_REF] Duvail | Uo22+ structure in solvent extraction phases resolved at molecular and supramolecular scales: a combined molecular dynamics, exafs and swaxs approach[END_REF].

Transition to phenomenological thermodynamics assisted by mesoscopic modeling

Quantifying collective effects in solvent extraction is possible if the free energy of every self-assembled aggregates, and its corresponding activity can be estimated. Accounting for every aggregate within multi-component MAL yields the correct speciation of all constitutes in the system. QM and MD cannot account for such diversity in aggregation due to the huge computational demand, but this can be achieved by further coarse-grained procedure, or the 'climb' to mesoscopic modeling scale. This intermediate level of theory (between atomistic simulations and macroscopic thermodynamics) focuses on establishing equations of state concerning the aggregation. Naturally, analytical expressions for speciation of the solvent phase can be transferred to macroscopic scale, including the chemical engineering.

Pseudophase model is the origin of the mesoscopic modeling of the solvent extraction

The origin of mesoscopic modeling in solvent extraction draws roots in colloidal self-assembly with the definition of coexisting pseudophases as stable entities [START_REF] Clarke | The hydrophobic effect: Formation of micelles and biological membranes, 2nd edition (tanford, charles)[END_REF]. The corresponding pseudophase model assumes the equilibrium between 'loose' extractant monomers and aggregated extractants [START_REF] Clarke | The hydrophobic effect: Formation of micelles and biological membranes, 2nd edition (tanford, charles)[END_REF]. Equilibrium can be written as

µ L,aggregated = µ L,monomer (6)
where µ L,aggregated and µ L,monomer are chemical potentials of aggregated and monomeric extractants in a given solvent. Underlying assumption of the model is that aggregates exist only at and above the micellization threshold. If the approach, which is first proposed for direct surfactant micelles in the aqueous phase, is generalized for organic phases as well, the notation of the micellization threshold dictates the usage of more accurate nomenclature: the critical aggregation concentration (CAC). The basic equation which describes equilibrium has the following form

µ -• L,aggregated = µ -• L,monomer + k B T ln c CAC (7)
where µ -• L,aggregated and µ -• L,monomer are standard chemical potentials of the aggregated and monomeric ('loose') extractant in the bulk, respectively. The left side of Eq. 7 does not depend on the concentration because in the pseudophase model, it corresponds to the chemical potential of an condensed phase which is constant (external potential). c CAC is the equilibrium concentration at which aggregation starts. Rearranging Eq. 7 gives µ -• L,aggregated -µ -• L,monomer = k B T ln c CAC where we can identify the aggregation free energy per extractant molecule ∆ agg G -• = µ -• L,aggregated -µ -• L,monomer . Within the pseudophase model ∆ agg G -• thus depicts the free energy for transferring a single extractant molecule from the bulk solvent phase to the existing aggregate i.e., its extractant film. The model is sketched in Mesoscopic modeling square in Figure 2, where the equilibrium between monomeric extractant in the bulk and extractant film (which represents the pseudophase) are shown.

In a way, all subsequent mesoscopic (sometimes referred simply as thermodynamic models) are an offspring of the original pseudophase model. In fact, µ -• L,aggregated accounts for local collective effects of the aggregated extractant in the film. Collective effects are often decoupled as sterical constraints, volume terms, and well as interfacial oil/water effects and extractant head groups lateral interactions [START_REF] Danov | Growth of wormlike micelles in nonionic surfactant solutions: Quantitative theory vs. experiment[END_REF].

Still, unlike for the case of direct o/w micelles, where surface complexation can be accounted by simple Poisson-Boltzmann description of charge inhomogeneities along with site population balance [START_REF] Müller | How do anions affect self-assembly and solubility of cetylpyridinium surfactants in water[END_REF][START_REF] Leodidis | Amino acids in reversed micelles. 4. amino acids as cosurfactants[END_REF], in the case of reverse w/o aggregates in organic solvent the cation/extractant complexation energy contribution, which takes hold in the interior of the aggregate cannot be easily accounted for. In the attempt to overcome this major difficulty, ienaic approach was introduced, where complexation energy was deduced from total aggregation free energy and every other relevant term in the free energy sum [START_REF] Zemb | Recycling metals by controlled transfer of ionic species between complex fluids: en route to "ienaics[END_REF].

It is very important to note that CAC, just like cation extraction efficiency, can be readily measured. Indeed, techniques such as vapor pressure osmometry (VPO), small-angle neutron (SANS) and X-ray scattering (SAXS) [START_REF] Pecheur | Synergism in a HDEHP/TOPO Liquid-Liquid Extraction System: An Intrinsic Ligands Property?[END_REF][START_REF] Déjugnat | Reverse aggregate nucleation induced by acids in liquid-liquid extraction processes[END_REF] or surface tension measurements [START_REF] Rey | Thermodynamic Description of Synergy in Solvent Extraction: I. Enthalpy of Mixing at the Origin of Synergistic Aggregation[END_REF] can be used to obtain an average macroscopic response of the extraction system e.g. CAC. Determination of c CAC at which some macroscopic property changes, yields ∆ agg G -• when Eq. 7 if used. Note that the aggregation free energy ∆ agg G -• and the free energy of transfer of solutes to solvent phase ∆ tr G -• are not same, even though they are often of the same order of magnitude [START_REF] Rey | Thermodynamic Description of Synergy in Solvent Extraction: I. Enthalpy of Mixing at the Origin of Synergistic Aggregation[END_REF]. Ions do not appear directly in Eq. 7, therefore, one must be careful when comparing these values: only qualitative comparability is possible.

From the theoretical point of view c CAC can be calculated by combining compatible models, such as the pseudophase model and models for activity in of the aqueous phase [START_REF] Déjugnat | Reverse aggregate nucleation induced by acids in liquid-liquid extraction processes[END_REF]. For the practical case of acids extraction by non-ionic DMDOHEMA extractant, Eq. 7 was used to express the system of two distinct aggregates: 'empty' aggregates (no electrolyte present in the aggregate) and acid-filled aggregates. Minimum aggregation number was set to 4, based on experimental data. Exploiting the property that the sum of the aggregation free energies is equivalent to competition reaction in multicomponent MAL, two competing transfers of extractant from the bulk to the distinct pseudophases were rationalized: electrolytefree and electrolyte-assisted aggregation. Few important outcomes emerged: the electrolyte lowers ∆ agg G -• for 1.5 k B T , and acid serve as a sort-of 'nucleator'. Also, aqueous phase activities for a series of acids were calculated and plugged into the model, ion-specificity of acid extraction to solvent phase. Calculated CAC as a function of the activity of the aqueous phase predicts two distinct regimes, with sharp transition among the two. This sharp transition resolves pre-micellization and micellization regime, and it is a trademark of pseudophase modeling in colloidal chemistry. This is somewhat contrary to experimental findings where extraction of solutes is observed even before CAC [START_REF] Meridiano | Aggregation in organic solutions of malonamides: Consequences for water extraction[END_REF].

An opposite way of using CAC is to fit the experimentally determined c CAC with a thermodynamic model, as was done for the case of synergistic mixtures of non-ionic TOPO and acidic HDEHP extractants [START_REF] Rey | Thermodynamic Description of Synergy in Solvent Extraction: I. Enthalpy of Mixing at the Origin of Synergistic Aggregation[END_REF]. The novelty in that study was to write µ -• L,aggregated as an equation of state, where the generalized 2D regular solution theory was employed to represent lateral interaction of extractant in the extractant film. χ interaction parameter was deduced by fitting the ∆ agg G -• (which was calculated from measured c CAC ). This was the first theoretically obtained argument that the composition of the mixed aggregated extractant film (expressed by extractant mole fraction in the aggregate) does not need to be identical to extractant mole fraction in the bulk. This was a cornerstone for the development of models which can include much-needed polydispersity in aggregation.

A recent study employed equilibrium approach between distinct pseudophases (via Monte-Carlo-like algorithm) to obtain nanoscale picture of solvent phase in the case of dialkylamides/UO 2+ 2 system and to further connect it with viscosity properties, which are of crucial importance for actual plant design [START_REF] Pleines | A minimal predictive model for better formulations of solvent phases with low viscosity[END_REF]. The authors demonstrated that the larger curvature effect favors smaller aggregates formations, while reducing the probability of the formation of dense aggregates networks. As a consequence the global viscosity of the solvent phase decreases [START_REF] Pleines | A minimal predictive model for better formulations of solvent phases with low viscosity[END_REF]. This opens a path towards the generalization of co-solvents and co-surfactants influence and applicability in certain industrial processes.

Before proceeding to more complicated models, it is worth to note that Langmuir adsorption models are consistent with pseudophase model because of the assumption that above CAC the aggregate needs to exist in the solvent phase. With this assumption the aggregates can be treated like complexation sites [START_REF] Zemb | Reverse Aggregates as Adaptive Self-Assembled Systems for Selective Liquid-Liquid Cation Extraction[END_REF][START_REF] Dourdain | Liquid-Liquid Extraction of Acids by a Malonamide: II-Anion Specific Effects in the Aggregate-Enhanced Extraction Isotherms[END_REF]. ∆ tr G -• is then equal to Langmuir binding energy. The approach was used to calculate the acid extraction [START_REF] Dourdain | Liquid-Liquid Extraction of Acids by a Malonamide: II-Anion Specific Effects in the Aggregate-Enhanced Extraction Isotherms[END_REF], and Eu 3+ extraction isotherm, which asymptotically approaches the saturation threshold of the solvent phase by cations [START_REF] Testard | Liquid-liquid extraction: An adsorption isotherm at divided interface?[END_REF]. It must be noted that only single type of aggregate is assumed to exist (in order to have reasonable fitting experimental data) which is contrary to polydispersity of actual systems. Also, it remains unclear whether complexation site corresponds to the entire aggregate (and if, what is the aggregation number!?), or the total aggregation number can be divided into multiple sites?

All pseudophase models suffer from intrinsic default: the aggregate (or the pseudophase) must exist, with pre-defined aggregation number. Accounting for different aggregates usually requires additional parameterization. Since the concentration of aggregates in the solvent phase cannot be obtained explicitly, and furthermore, there are issues concerning the inclusion of cation-specific complexation energy terms, pseudophase model was subject to upgrades.

Going beyond the pseudophase-like models

In order to calculate ∆ tr G -• rather than ∆ agg G -• as in usual pseudophase model, the first adjustment is to include the target solute directory into both sides of MAL, which then yields the partition (distribution) coefficient (see Eq. 2 and 3) [START_REF] Rey | Synergistic Extraction of Rare Earth Elements from Phosphoric Acid Medium using a Mixture of Surfactant AOT and DEHCNPB[END_REF]. There are several models which dealt with the transfer of solutes from the aqueous phase to reversed aggregates in organic solvent, where aggregates are assumed to be spherical on average. The authors would like to emphasize the importance of older, but not nearly outdated work by E. Leontidis and A. Hatton on the extraction of amino acids in bis(2-ethylhexyl) sulfosuccinate (AOT) reverse micelles in organic solvents, which was published on series of papers [START_REF] Leodidis | Amino acids in aot reversed micelles. 1. determination of interfacial partition coefficients using the phase-transfer method[END_REF][START_REF] Leodidis | Amino acids in aot reversed micelles. 2. the hydrophobic effect and hydrogen bonding as driving forces for interfacial solubilization[END_REF][START_REF] Leodidis | Amino acids in reversed micelles. 3. dependence of the interfacial partition coefficient on excess phase salinity and interfacial curvature[END_REF][START_REF] Leodidis | Amino acids in reversed micelles. 4. amino acids as cosurfactants[END_REF]. The series of papers gives a detailed description of the colloidal soft-matter theory in terms of competing forces that drive or quench the transfer of amino acids into the solvent phase, as well as highlights issues concerning the mesoscopic modeling of amino acid transfer. Questions were raised concerning the definition of standard states of solutes partitioning in w/o interface. Furthermore, in the attempt was made to rationalize the local environment of extracted amino acid. The partition ratio between 'water pool' of AOT reversed micelles, or in the interfacial film of AOT surfactants, was established for a series of distinct amino acids. Calculations showed that the decrease of the radius (an increase of the curvature) of reverse micelle depletes amino acids from the interfacial region into the bulk-like 'water pool'. Furthermore, amino acids susceptibility to hydrogen bonding at AOT interface is also one of limiting criteria i.e. more hydrogen donor/acceptor atoms, higher the partition at the interfacial film. A turn in mesoscopic modeling was made when the free energy of the aggregate F agg was calculated [START_REF] Dufrêche | Effect of long-range interactions on ion equilibria in liquid-liquid extraction[END_REF]. The practical case studied was Eu 3+ extraction by DMDOHEMA. Note that this is opposite to pseudophase model where the free energy of aggregated extractant is calculated. By definition, at infinite dilution F agg is equal to the standard chemical potential of the spherical aggregate in particular solvent µ org,-• Aagg [START_REF] Dufrêche | Effect of long-range interactions on ion equilibria in liquid-liquid extraction[END_REF]. As a result, the transfer of solutes was finally represented by equilibria as given by Eq. 4. This first model used colloidal soft-matter theory to decouple µ org,-• Aagg into contributions by namely, the free energy of the extractant film µ chain agg and the free energy of the polar core µ polar core agg . Former was identified µ chain agg = N L κ (p -p 0 ) 2 /2, where N L is the aggregation number, and κ , p, p 0 are again bending constant, packing and spontaneous packing parameters. µ polar core agg was identified as the Gibbs energy of electrolyte solution confined inside the aggregate, µ bulk agg . Additional terms were added to represent long-range interaction, such as dispersion force by the calculation of image charge and polarization, which was found to be small effect. When µ bulk agg is collected with Gibbs energy of solutes in the aqueous phase (left side of MAL), the free energy of the transfer of solutes to polar cores of aggregates is obtained, ∆µ bulk agg . The term is usually the energy penalty associated with transfer of solute to highly concentrated medium. Within ienaic approach, the energy contribution due to the coordination chemistry complexation was added in the calculation. The model was fitted to the experimental value of ∆ tr G -• = 2k B T , and complexation energy term was estimated to -32 k B T , thus clearly a leading motor for the extraction.

Since the concept included the formation of single spherical aggregate, the next step was to account for polydispersity in aggregation numbers. In the follow-up study, the benefits of multicomponent MAL were exploited, where aggregates different by the aggregation number were included [START_REF] Karmakar | Combined supramolecular and mesoscale modelling of liquid-liquid extraction of rare earth salts[END_REF]. This was possible since energy terms in the sum have functional forms that depend on the composition, while parameters are fixed and irrespective of the compositions [START_REF] Bley | A predictive model of reverse micelles solubilizing water for solvent[END_REF]. The model had slight issues with overestimation of water coextraction but agreed fairly well with experiments of cation extraction. Still, it must be noted that at the low radius of spherical aggregate (high curvature) the validity of dispersion forces questionable, therefore the model is probably more suited for larger aggregates such as microemulsions (larger w/o droplets) or vesicles.

The next improvement of the approach came with redefining the complexation energy term F complexation as a sum of individual internal energy contributions. Dispersion forces and polarization were omitted. The sketch of reverse spherical aggregate is given in Figure 4. Note that the notation of the energy terms was changed from µ agg to F agg . Contributions approximate cation/extractant bonds as additive pair-potentials, and in the model are fitted to an independent set of experimental extraction data. The analytical form of entropic part was obtained with approximating the curved o/w interface as 2D lattice (see the right side of Figure 4). N complex depicts the number of arrangements that cation can have on 2D lattice, consequently, the entropy is equal S = -k B lnN complex . Approximating the curved interface like this enables also the extractant exchange effect among the first and second coordination sphere by permutation plugged within the analytical expressions for the entropy [START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF]. This represents the effect that the aggregation number is higher than the complexation number i.e. not all aggregated extractants are involved in the direct complexation. This purely colloidal aspect was so far out of the reach by means of QM and MD but comes out naturally within mesoscopic modeling. Preferential cation binding (complexation) is still weighted by specific interactions, but the permutation effect exists and can be more than 1 k B T in magnitude, thus important. Model was applied on Eu(NO 3 ) 3 , Fe(NO 3 ) 3 and HNO 3 competing extractions by DMDOHEMA in n-heptane. F chain was calculated within multi-scale modeling, where results of MD simulations to access an average chain extension and geometrical parameters [START_REF] Duvail | The role of curvature effects in liquid-liquid extraction: assessing organic phase mesoscopic properties from MD simulations[END_REF]. Calculations in the semi-grand canonical ensemble, where thermodynamic reservoir included only water molecules i.e. fixed chemical potential of water, yielded quantitative speciation of all constituents in both phases, as well estimation of CAC and the extraction efficiency for various states of the system. Polydispersity in aggregation number was moderate, ranging from 4-5, but the number of coextracted water molecules varied severely [START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF]. Eu(NO 3 ) 3 or Fe(NO 3 ) 3 salts required more water molecules to dilute the aggregate core, thus the stabilization effect of water was quantified for the first time.

Further work included generalization acidic extractans such as HDEHP. The effect of deprotonation was captured on the basis that dissociated form of extractant is always bonded (complexed) to cation. Deduced functionality had form of -N L -ln(10)(pK • a -pH), where N L -is the number of depronotonated HDEHP which directly binds to central cation (smaller than the total aggregation number!), K • a is the effective depotonation constant. The functionality automatically satisfied the pH dependence [START_REF] Špadina | A colloidal model of prediction of the extraction of rare earths assisted by the acidic extractant[END_REF]. Deprotonated of acidic extractant generates a distinct head group in the extractant film, further enhancing N complex . The effect was taken into account as a generalized 2D Hildebrand regular solution theory. It must be noted that because of coupling of the complexation and deprotonation phenomena (deprotonation is necessary to complex multivalent cation!) have similar entropic parts, why is why double-counting in N complex needs to be prevented. Key novelties provided by this work were rationalization of complexation energy to -20.2 k B T per cation/extractant bond, which was stronger than in the case of non-ionic extractants (7.7 k B T for DMDOHEMA), extraction of Eu 3+ in pre-micellization system conditions, water depletion from the core of aggregates because of stronger complexation, proton exchange, as well as the introduction of the 'extraction landscape maps'. As a function of any independent system variable within the framework of the semi-grand canonical ensemble, it became possible to have a full speciation and extraction efficiency. CAC maps as a function of acidity and Eu(NO 3 ) 3 molality in the aqueous phase showed different regimes in aggregation. It was shown that optimized industrial formulations correspond to a 'hill' or a region of higher CAC (lower aggregation). CAC and thus the self-assembly is selectively controlled by either, Eu(NO 3 ) 3 or acid extraction [START_REF] Špadina | A colloidal model of prediction of the extraction of rare earths assisted by the acidic extractant[END_REF].

If the macroscopic molality of acid in the aqueous phase is changed by an order of magnitude, approximately 6 k B T additional complexation energy per cation/extractant bond is required to maintain the experimental extraction efficiency (expressed as ∆ tr G -• ). If we now consider that the ionization free energy of acidic extractants (which is calculated by more accurate QM at continuum solvent level of theory) is around 16 k B T , we realize that changing the macroscopic acidity of the system can serve as an extraction efficiency switch on/off.

A direct link between system independent variables and molecular nature of the cation/extractant couple is thus rationalized, further allows chemical engineering more efficient design of the extraction process. 5.3. Occurrence of the synergy in solvent extraction in the case of mixtures of extractants: a natural aspect of colloidal chemistry Some mixtures of extractants, such as tri-n-octyl phosphine oxide (TOPO)/HDEHP or DMDOHEMA/HDEHP, at certain conditions of system, show the occurrence of the synergistic extraction. Synergistic extraction is defined as higher efficiency than the sum of efficiency of both individual extractants [START_REF] Rydberg | Solvent Extraction Principles and Practice, Revised and Expanded[END_REF][START_REF] Muller | Understanding the synergistic effect on lanthanides(III) solvent extraction by systems combining a malonamide and a dialkyl phosphoric acid[END_REF][START_REF] Maangar | A microfluidic study of synergic liquid-liquid extraction of rare earth elements[END_REF]. Non-linearity of extraction efficiency versus the extractant mole fraction is ever the stepping stone in process modeling by chemical engineering [START_REF] Maangar | A microfluidic study of synergic liquid-liquid extraction of rare earth elements[END_REF]. During decades issue was not solved but rather the extra parametrization was employed which completely detached used models from chemistry i.e. non-linearity was treated as a numerical problem [START_REF] Rydberg | Solvent Extraction Principles and Practice, Revised and Expanded[END_REF]. A recent study on synergistic extraction compared some of the models based on chemical equilibria approach (formation of various complexes described by many fitted constants) and a simple reverse micelle theory [START_REF] Jackson | Synergism and aggregation in multi-extractant solvent extraction systems[END_REF]. The conclusion was drawn that in most cases the prediction fails completely. Instead of the synergistic peak, models predict linear extraction behavior.

In the attempt to resolve the phenomenon at hand, the reverse micelle theory described in previous section [START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF][START_REF] Špadina | A colloidal model of prediction of the extraction of rare earths assisted by the acidic extractant[END_REF], was generalized to account for formations of mixed aggregates as well. Main terms in the generalized expression for the free energy of particular aggregate x (F agg,x ) were written as

F agg,x = F droplet + F chain + F head group interactions + F complexation (8)
where F droplet is again the free energy of electrolyte solution confined inside the aggregate (µ bulk agg to previous studies), F chain is the free energy of mixed extractant film F head group interactions accounts for lateral head groups interactions, F complexation is generalized coordination chemistry complexation energy. Generalized expressions can be found in the reference [START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. Calculations were made in semi-grand canonical ensemble, where aqueous phase was considered as a infinite reservoir of solutes and water molecules with fixed chemical potential. The two extractant concentrations were subject to minimization procedure which resulted in finding their values at equilibrium. The approach included all possible aggregates compositions, practically tens of thousands [START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. The connectivity of cation and extractants were taken based on QM modeling of mixed aggregates [START_REF] Muller | Understanding the synergistic effect on lanthanides(III) solvent extraction by systems combining a malonamide and a dialkyl phosphoric acid[END_REF], within a multiscale approach. This robust analysis provided an insight into macroscopic control of the complex aggregation behavior at nanoscale. The complexity is sketched in Figure 5 where the extraction efficiency (the cation loading in solvent phase) is plotted as a function of acidity in the aqueous and extractant mole fraction in the solvent phase. Due to the speciation which is obtained by the model, decoupling of the extraction efficiency regimes into contributions by different aggregates is possible. For example, at low acidity and a lower concentration of non-ionic extractant, deprotonation of acidic extractant is responsible for high extraction. Consequently, only aggregates with acidic extractants form reverse aggregates, while non-ionic extractants remain dispersed in the solvent phase. Increasing acidity starts to quench the extraction, but if the extractant mole fraction is high enough (higher concentration of the second extractant in the mixture) the synergistic behavior takes hold, where the aggregation is predominantly controlled by mixed aggregates.

Furthermore, the synergy in HDEHP/DMDOHEMA extraction system was attributed to the higher number of possibilities for complexing the cation (see the right side of Figure 4). This specific gain in entropy contribution lowers the overall configurational entropy loss of bringing a large number of species into selfassembled aggregate i.e. reduces the system free energy (see MAL in Eq. 3). In coordination chemistry, this effect of huge polydispersity is not tractable, but in colloidal chemistry, it comes out naturally by a simple exchange of extractant between the first and the second coordination sphere. The most relevant conclusion was that many aggregates, similar in F agg,x but completely different in their compositions, can coexist at equilibrium.

Even at this point, a model be can used to decipher the mechanism of the phenomenon underlying the unexpected 'salting-in' effect of the mineral acid in non-ionic solvating extractants system, where after its [START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. Distinct aggregation behavior of mixtures, sketched by different types of reverse micelles, is under control of all independent system variables. The response of the system i.e. its extraction efficiency is under the strong influence of the extractant mole fraction in the solvent phase, but not only! Differences in aggregation, as well as the extraction efficiency, become transparent as can be seen, for example, that synergistic peak attributed to mixed aggregates formation is not nearly 'efficient' as the extraction at low acidity where only dissociated acidic extractants aggregate and extract heavily [START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. This example shows how mesoscopic modeling can capture the description at nanoscale and bridge it with a macroscopic domain of extraction efficiency studied by chemical engineering.

certain bulk concentration, the cation extraction efficiency decreases [START_REF] Marcus | The law of mass-action versus non-ideal behaviour in distribution equilibria[END_REF]. The colloidal theory stresses that for very large chemical potentials of the acid in the aqueous phase, the self-assembly is induced by acid/water reverse aggregates formation. Consequently, a cation extraction competition reaction is quenched within global MAL. The systems relaxes by high acid/water uptake to the solvent phase.

Current issues and challenges in mesoscopic modeling of solvent extraction

There are few issues in mesoscopic modeling of solvent extraction, namely on intra-and inter-aggregate scales.

On the intra-aggregate scale, the modeling of aggregated extractant chains in an explicit solvent is based on the concept of packing parameter and its deviation from spontaneous packing. The approach is somewhat inefficient and requires several adjusted parameters, as well as approximations of the average conformations of chains. Furthermore, the approach is a subject of a constant change of the definitions, thus values of κ , p, p 0 from different studies can be seldom compared quantitatively. Another difficulty arises since most utilized concepts do not have temperature dependency. The ones that do have it rely on extra parameter and functionality is glued by hand. Absence of temperature dependency is quite a severe issue since it is well known that weak self-assembly can be annihilated by raising the temperature and the extraction efficiency drops [START_REF] Rydberg | Solvent Extraction Principles and Practice, Revised and Expanded[END_REF]. In modeling, this property is still out of reach and presents an open challenge. Current mesoscopic models calculate the free energy of the polar core of the aggregate, which is always assumed to be liquid [START_REF] Karmakar | Combined supramolecular and mesoscale modelling of liquid-liquid extraction of rare earth salts[END_REF][START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF][START_REF] Rey | Synergistic Extraction of Rare Earth Elements from Phosphoric Acid Medium using a Mixture of Surfactant AOT and DEHCNPB[END_REF]. Since concentrations of ions in small cavities are high, often higher than solubility limit in the bulk solution, most models consider ideal solutions. In this way, the energy penalty for transferring the solute to a concentrated solution is only entropic in nature. Non-ideality is then attributed to the complexation energy [START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. With that in mind, the level of theory employed for calculation of solutes properties needs to be same in both sides of MAL (see Eq. 3). For calculation of solution properties such inside polar core of aggregates (called bulk term µ bulk agg in previous, or F droplet in later studies) as activity coefficients, or well as partial molar volumes of ions should correspond to that close to solubility limit, where electrostriction is no longer valid, rather a bare ion volume terms become more important [START_REF] Marcus | Electrostriction in electrolyte solutions[END_REF]. The generalized expression of mixtures of solutes needs to be derived, especially for systems which include lanthanide salts and mineral acids.

On inter-aggregate scales, the solvent depletion effects need to be accounted as well, since quantifying this effect can lead is necessary for the calculation of inter-aggregates activity coefficients [START_REF] Sapir | Macromolecular compaction by mixed solutions: Bridging versus depletion attraction[END_REF]. Another major challenge is to obtain an analytical expression for the free energy of aggregates characterized by different morphologies, such as ellipsoidal, cylindrical, or lamellar. Restraining to spherical morphology is severe drawback [START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. Coarse-grained MD or MC simulations within multi-scale modeling could provide potentials which describe inter-aggregate interactions, which can be recalculated into activity coefficients.

Connecting mesoscale and chemical engineering

Calculations which utilize mesoscopic models are generally fast and in performance similar to the evaluation of extraction isotherms in chemical engineering (see algorithms in references [START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF] and [START_REF] Moeyaert | Experimental and modelling study of ruthenium extraction with tri-n-butylphosphate in the purex process[END_REF]). In most cases, the dichotomy-based minimization methods are utilized to find equilibrium state of the system, which is in most cases trivial since concentration profiles are continuous monotonic functions. The output of the most recent mesoscopic models is the speciation of all species in the system for any set of independent system variables [START_REF] Špadina | Multicomponent Model for the Prediction of Nuclear Waste/Rare-Earth Extraction Processes[END_REF][START_REF] Špadina | A colloidal model of prediction of the extraction of rare earths assisted by the acidic extractant[END_REF][START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. The calculated cation loading of the solvent phase can be expressed by McCabe-Thiele isotherms for a given phase volume ratio and the thermodynamic temperature [START_REF] Huang | Selective recovery of rare earth elements from ion-adsorption rare earth element ores by stepwise extraction with HEH(EHP) and HDEHP[END_REF]. Consequently, calculations can be expanded to predict extraction/desextraction efficiency (for example in the counter-current process), as a function of the number of stages by the mere addition of another minimization loop (which represents various stages).

Conclusion

In this review, we highlighted the current paradigm in modeling of the solvent extraction of lanthanides and actinides by non-ionic and acidic extractants. Since the solvent phase is a complex fluid, a multi-scale modeling should be used to elucidate effects which occur at different length scales.

Even though the solvent phase is globally apolar, ions are still readily transferred from the aqueous phase. On nanoscale hydrophilic cavities assisted by extractant molecules are formed and ions are solubilized. These nanodomains sometimes resemble self-assembled reverse aggregates. It was found that local inhomogeneities of charge, as well as hydration influence interactions within and between self-assembled aggregates. Therefore, supramolecular and colloidal effects are important. While different theoretical approaches have been utilized, there is still a challenge to assist the chemical engineering of separation industry. Due to the huge computational demand of quantum and classical simulations, many of the collective effects cannot be transferred to the phenomenological thermodynamics. It must be noted that the calculated extraction efficiency depends explicitly on the choice of the reference state of solubilized cation in both phases.

A shift to phenomenological thermodynamics can be achieved with mesoscopic modeling which serves as a bridge between nanoscale and chemical engineering. By approximating aggregates as well-shaped colloids, whose geometric properties can be expressed in the functional form (where the composition is variable), the free energy of the transfer of ion ∆ tr G -• can be calculated. This enables an explicit closure of MAL and yields full speciation of the extraction system, in both aqueous and solvent phases. The speciation is of direct use to chemical engineers.

A prospect, or a focus, of multi-scale modeling should be directed towards the prediction of actual extraction efficiency. A concept of the multi-scale modeling scheme can include: Calculations of coordination numbers, geometry and the heats of formation of cation/extractant complexes via Born-Haber cycles from QM. Results can be compared with spectroscopic data.

Simulations of the self-assembly, phase stability, as well as explicit diluent and effects concerning branching of extractant chains from MD and MC. Results can be compared with scattering techniques.

At this point, based on two previous steps in multi-scale modeling, the type of an adequate mesoscopic model should be chosen to provide the best description of the solvent phase (e.g., reverse micelle theories or microemulsion theories). Then the total free energy of aggregates with distinct morphologies can be calculated and imported into MAL. The closure of MAL can yield a good estimate of ∆ tr G -• , obtained with a minimal number of adjusted parameters. Results can be compared with batch experiments of extraction efficiency. Calculated speciations expressed in the form of McCabe-Thiele isotherms, can be implemented into engineering flowsheets.

The coupled study of solvent extraction at different scales can help chemical engineering in the design of novel processes, with respect to the green chemistry guidelines.

Once again, it is worth emphasizing that the ienaic approach, which assumes systematic estimation of aggregation and transfer free energies of extraction systems, is coherent with transfer free energies calculations of ions into collapsed polymers matrix, hydrogels, or polymer membranes [START_REF] Kanduč | Transfer free energies and partitioning of small molecules in collapsed pnipam polymers[END_REF][START_REF] Kanduč | Aqueous nanoclusters govern ion partitioning in dense polymer membranes[END_REF]. In fact, the mechanism of ions stabilization in globally apolar medium, which assumes the emergence of polar waterfilled nanodomains that nest charged species, is not surprising at all and will be soon identified in many soft-matter-based industrial applications.

Figure 1 :

 1 Figure 1: Schematic phase behavior of malonamide extraction system as a function of the solute and extractant concentrations in the organic phase. While solute concentration depicted on y-axis starts from 0, the extractant concentration depicted on x-axis starts from certain concentration c ex,0 > 0 [p.d.u.]. The figure was inspired by similar one presented in the reference[10]. Photographs of the third phase, popularly called 'the red sphere' (left), as well as the normal biphasic solution system (right), were kindly provided by Dr. Mark Antonio from Argonne National Laboratory.

Figure 2 :

 2 Figure 2: Multi-scale modeling: time and distance scales of the phenomena and corresponding theoretical methods used to represent them. The figure was inspired by a similar one presented in the reference[START_REF] Guillaumont | Modeling selectivity in liquid/liquid extraction[END_REF]. Quantum mechanics square: the snapshot of the of density functional theory (DFT) calculations yielding the optimized structure of Nd 3+ coordinated with 6 dimethylphosphoric acid ligands (a truncated version of HDEHP with short alkyl chains) in a vacuum is reproduced from the reference[START_REF] Muller | Understanding the synergistic effect on lanthanides(III) solvent extraction by systems combining a malonamide and a dialkyl phosphoric acid[END_REF] with permission from Elsevier. Classical simulations square: the snapshot of MD simulations of stable Ln(NO 3 ) 3filled aggregate in n-heptane solvent is reprinted from the reference[START_REF] Chen | Stability of reverse micelles in rare-earth separation: a chemical model based on a molecular approach[END_REF] with permission from the Royal Society of Chemistry. The neutral polar core (cation + anions) which contains multiple H 2 O molecules is enclosed by 4 DMDOHEMA non-ionic extractant molecules. Mesoscopic modeling square: idea of the pseudophase-like models where monomeric extractants are in the equilibrium with already formed extractant film. Chemical engineering square: the concept of the chemical engineering multiple equilibrium approach. Every possible equilibria are established and constants are fitted to the experiments[START_REF] Dartiguelongue | Modelling of uranium(VI) extraction by D2EHPA/TOPO from phosphoric acid within a wide range of concentrations[END_REF].

Figure 3 :

 3 Figure3: Schematic representation of a laboratory-scale solvent extraction experiment (left side). The system is made of two non-miscible liquids, namely organic solvent (pale yellow color), and the aqueous phase (pale blue color). Red circles represent target solutes e.g. lanthanide cation. Zoomed region depicts sketch the nanoscale 'view' of the system, where the equilibrium state is composed of a myriad of coexisting structures (right side). The issue with the choice of the standard state in the solvent phase is demonstrated as a transfer of hydrated cation from the aqueous to the polar core of self-assembled aggregate (orange arrow). a) In the case that solvated cation is understood as a complex which assumes the coordination of the cation by extractant molecules, counterions and coextracted water molecules, then ∆trG -• accounts only the first-coordination sphere interactions. b) If the solvated ion is the whole the aggregate, then beside interactions between first neighbors are accounted for. Additional collective effects include the formation of the extractant second-coordination sphere, hydrophobic extractant chains configurations, an organic solvent wetting and depletion effects of extractant film, the activity of the aggregate itself, etc.

FFigure 4 :

 4 Figure4: Schematic representation of the spherical reverse micelle and its corresponding free energy terms. F polar core depicts the free energy of the electrolyte solution confined within the core of the reverse micelle. F head group interactions describes free energy associated with lateral, pairwise interactions among extractant head groups within the extractant film. The term is calculated as a generalized two-dimensional Hildebrand regular solution model. F chain accounts for the free energy change due to the confining extractant side chains within a fixed volume, in a given solvent. F complexation represents coordination chemistry complexation, as a sum of individual internal energy contributions. Analytical expressions for associated entropy are obtained when the curved interface is approximated as a 2D mixed-site lattice.

Figure 5 :

 5 Figure5: Schematic representation of the cation loading in the solvent phase (concentration of the extracted cation in the solvent phase) as a function of acidity in the aqueous, and the extractant mole fraction in the solvent phases. Note that actual numbers presented in the plot reflect qualitatively the results published in the reference[START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. Distinct aggregation behavior of mixtures, sketched by different types of reverse micelles, is under control of all independent system variables. The response of the system i.e. its extraction efficiency is under the strong influence of the extractant mole fraction in the solvent phase, but not only! Differences in aggregation, as well as the extraction efficiency, become transparent as can be seen, for example, that synergistic peak attributed to mixed aggregates formation is not nearly 'efficient' as the extraction at low acidity where only dissociated acidic extractants aggregate and extract heavily[START_REF] Špadina | Synergistic solvent extraction is driven by entropy[END_REF]. This example shows how mesoscopic modeling can capture the description at nanoscale and bridge it with a macroscopic domain of extraction efficiency studied by chemical engineering.

Acknowledgments

The authors acknowledge the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. [320915] "REE-CYCLE": Rare Earth Element reCYCling with Low harmful Emissions. Research Agency for support through grant BIFR/CEA/16-18-002 and the Slovenian Research Agency for support through program P3-0388 are acknowledged. We would like to thank M. Duvail and J.-F. Dufrêche for providing us with snapshots of MD simulations. Furthermore, the author would like to thank M. Antonio for providing us figures of the third phase formation which was studied in his laboratory, as well as for delightful discussions. The authors would like to thank S. Gourdin-Bertin, S. Pellet-Rostaing, S. Dourdain, M. Duvail, J.-F. Dufrêche and T. Zemb for useful discussions and for the encouragement to challenge the current paradigms of separations processes relevant in nuclear industry for the development of closed fuel cycles, and recycling of rare earth elements for sustainable chemical industry.