A Telescoping View of Solute Architectures in a Complex Fluid System

Short- and long-range correlations between solutes in solvents can influence the macroscopic chemistry and physical properties of solutions in ways that are not fully understood. The class of liquids known as complex (structured) fluids—containing multiscale aggregates resulting from weak self-assembly—are especially important in energy-relevant systems employed for a variety of chemical- and biological-based purification, separation, and catalytic processes. In these, solute (mass) transfer across liquid–liquid (water, oil) phase boundaries is the core function. Oftentimes the operational success of phase transfer chemistry is dependent upon the bulk fluid structures for which a common functional motif and an archetype aggregate is the micelle. In particular, there is an emerging consensus that mass transfer and bulk organic phase behaviors—notably the critical phenomenon of phase splitting—are impacted by the effects of micellar-like aggregates in water-in-oil microemulsions. In this study, we elucidate the microscopic structures and mesoscopic architectures of metal-, water-, and acid-loaded organic phases using a combination of X-ray and neutron experimentation as well as density functional theory and molecular dynamics simulations. The key conclusion is that the transfer of metal ions between an aqueous phase and an organic one involves the formation of small mononuclear clusters typical of metal–ligand coordination chemistry, at one extreme, in the organic phase, and their aggregation to multinuclear primary clusters that self-assemble to form even larger superclusters typical of supramolecular chemistry, at the other. Our metrical results add an orthogonal perspective to the energetics-based view of phase splitting in chemical separations known as the micellar model—founded upon the interpretation of small-angle neutron scattering data—with respect to a more general phase-space (gas–liquid) model of soft matter self-assembly and particle growth. The structure hierarchy observed in the aggregation of our quinary (zirconium nitrate–nitric acid–water–tri-n-butyl phosphate–n-octane) system is relevant to understanding solution phase transitions, in general, and the function of engineered fluids with metalloamphiphiles, in particular, for mass transfer applications, such as demixing in separation and synthesis in catalysis science.


■ INTRODUCTION
The literature of science and engineering of colloids is filled with ternary phase diagrams comprising two immiscible liquids (e.g., water and oil) and a soluble amphiphilic surfactant. 1,2 Upon mixing, the initially biphasic system can equilibrate as one of four main types of microemulsions described by Winsor. 3 In each of these oil-in-water and water-in-oil phase combinations, the micelle 4 and the reverse micelle, 5 respectively, are the fundamental building blocks in the formation of supramolecular solution structures, including cylindrical, close-packed, and lamellar mesophases. 6,7 In the vernacular of soft matter sciences, such complex (structured) fluids are built upon multiscale correlations resulting from physicochemically soft and weak interactions. The organization of solutes (i.e., neutrals, cations, and anions) into aggregated structures through weak interactions, with energies of 1 k B T and less, 8,9 drives the assembly of aggregate architectures in which a delicate balance between short-range attractive, longrange repulsive (SALR) interactions regulates pattern formation 10 and phase behaviors. Whereas the weak long-range forces are small compared the strong short-range forces in metalloamphiphile chemical bonding (approximately 300 k B T), they influence complex physical responses, often resulting in complicated fluid behaviors, including the critical phenomenon of phase splitting.
In the limit of a vanishing critical temperature and weak attraction, the splitting of SALR fluids into colloid-poor (gas) and colloid-rich (liquid) phasesknown as the gas−liquid phase separation model 11,12 is suppressed 10 in favor of percolation. 13, 14 In contrast, strong attraction energy landscapes between nanoscale solute architectures leading to gelation in colloidal SALR systems are a direct consequence of equilibrium gas−liquid phase separation, not percolation. 15 Recently, self-assembly phenomena in SALR colloids at macroscopic length scales were attributed to a previously unrecognized mechanism, dubbed ionic sphericity. 16 In previous research on complex fluids that function in ternary (water−oil−surfactant) systems for chemical separations by liquid−liquid (biphasic) extraction, the splitting of the oil phaseinto two oil phaseshas been attributed to the gas− liquid model of phase transition. 17 This interpretation, originally proposed in 1998 by Erlinger et al. 18 on the basis of small-angle X-ray scattering (SAXS), is the genesis of the micellar model of phase splitting. 19 In this, a complex fluid of sticky spheres (e.g., reverse micelles with surface adhesion) with a short-range interaction potential that is approximated by a hard repulsive core and a rectangular attractive well (as defined by Baxter) 20 was shown to be consistent with the experimental SAXS data. Despite the system approximations and simplifications required to make the data analyses tractable, the insights obtained from the use of this model with regard to the energetics of the phase transition provide a zeroth-order benchmark for all subsequent phenomenological comparisons, including those reported here.
To address aspects of the complexities with regard to organic phase architectures, in general, and the origins of phase splitting, in particular, we provide new insights into the soft matter chemistry and physics of separation science. In particular, lacking in our understanding of these systems is a robust description of the weak interactions involved in aggregation and how they scale up (and down) as manifest in the type of solution structures that are formed, and how, in combination, they influence technically relevant material behaviors. For example, self-organized, complex fluid systems are exemplified by a wide variety of familiar substances, including polymers, gels, soaps, and inks. Recent advances with X-ray and neutron scattering probes focused on fluid systems have significantly expanded the list of liquids exhibiting aggregation behaviors in a diverse range of applications. Some of these include commercial-scale processes for enhanced oil and gas recovery; production of pharmaceuticals and commodity chemicals; transportation and sequestration of natural gas and CO 2 ; phase transfer and micellar catalysis; and even nuclear fuel reprocessing operations like PUREX (plutonium uranium reduction extraction). 21 In this latter system and many others like it, there is an emerging consensus that solute aggregation can have a profound impact on the chemical behavior of the organic solution, notably on both the efficiency and selectivity of the process. 22−33 In fact, reverse micellar solvents with the macroscopic appearance of a homogeneous solution are routinely employed in organic-and biological-based separation, purification, and conversion processes with little understanding of the morphology and arrangement of the micelle domains. 34,35 These contemporary applications are revolutionizing fundamental separations concepts by shifting the focus away from the classical enthalpic drivers associated with local, molecular-scale interactions into multiscale-ordering, which is hypothesized to include large entropic contributions that may dominate the energy landscape of bulk behaviors 36 that are hallmarks of coacervation 35 and microemulsions. 37 However, the relationships between the microscopic solute structures and phase behaviors remains largely unresolved. The present study investigates the macroscopic splitting of the organic phase during liquid−liquid (water−oil) extraction by examining the microscopic structures of metal coordination chemistry and the mesoscopic architectures of supramolecular chemistry. Our findings provide a fundamental, telescoping structure perspective into processes that are pivotal to the international energy economies vis-a-vis the refining and purification of platinum group and base metals as well as rare earth and actinide elements.
In these biphasic hydrometallurgical operations, a high metal loading of the organic phase is advantageous to obtain costeffective separations. However, in practice, it is not possible to achieve this because when the concentration of the extracted metal ion exceeds the critical concentration, a ruinous phase transition occurs. In this, the organic phase splits into a light phase and a heavy one; the heavy organic phase is called the third phase. 21,38 In line with recent research on the subject, 39 we choose to cast organic phase splitting in liquid−liquid extraction as a critical phenomenon. Understanding the mechanism of third-phase formation offers the prospects for the development of advanced, high-performance separation science and technology. With specific regard to the PUREX processon a molecular scaleit is generally known that the amphiphile tri-n-butyl phosphate (TBP; see the Supporting Information showing the chemical structure) coordinates with tetravalent plutonium as the nitrate complex to form the charge-neutral tetranitrato-disolvate, Pu(NO 3 ) 4 (TBP) 2 , 38 in the organic phase. On a supramolecular scale, small-angle neutron scattering (SANS) has been modeled in the Baxter sticky-sphere framework to reveal interactions between small spherical reverse micelles containing three to five TBP molecules in n-dodecane. 40 The critical phenomenon of phase splitting was attributed to the strength of the attractive forces (with a potential energy of up to −2.6 k B T) between the polar cores of the reverse micelles. This is dubbed the micellar interaction model. 41 More recent reports involving computational studies have called into question the validity of the

ACS Central Science
Research Article model for both metal-free 42−44 and metal-ion-loaded 45 systems.
As with any model, new experimental and computational data obtained with contemporary methodology oftentimes reveal nuances and subtleties that are not resolved in the original efforts. 41 For example, Baldwin et al. 45 reported that the Baxter model for interpreting SANS profiles in the PUREX system can be misleading about the strength of the interaction energies between reverse micelles. That is, under selected conditions, the model tends to predict nonattractive interactions contrary to direct observations of phase splitting due to strong attractive associations. There are at least two possible reasons for these contradictions. First, to simplify analysis with the Baxter model, the extracted solutes were assumed to occupy the cores of spherical or ellipsoidal micelles without knowledge of the true coordination chemistry and particle morphologies. Second, the incoherent scattering intensity from hydrogen atoms contained in the samples was overly simplified as a fitting parameter in SANS data analyses. 41,45−47 Since the requisite coherent scattering intensity reflecting the structural information on the solutes is typically buried under the level of incoherent scattering intensity, it is crucial to evaluate the incoherent scattering intensity in a rigorous manner. Accordingly, in this study we used neutron polarization analysis, which is the best experimental method for evaluating the incoherent scattering intensity to avoid misinterpretation of the SANS data. 48 This development of an improved method to model long-range aggregate structures beyond coordination chemistry augments our understanding of the solution state of the organic phase following extraction.
In this work, extended X-ray adsorption fine structure (EXAFS) experiments with density functional theory (DFT) calculations, molecular dynamics (MD) simulations, and SANS were carried out in a complementary style to explore the micro-and mesoscopic architectures of the solvate cluster formed by Zr(NO 3 ) 4 and TBP, Zr(NO 3 ) 4 (TBP) 2 , as a nonradioactive chemical surrogate of Pu(NO 3 ) 4 (TBP) 2 in PUREX. The combination of experimentation and simulation is of great advantage to avoid the above-mentioned uncertainties by accessing a broad spectrum of length scales from 0.1 to 100 nm in a telescoping, hierarchical model 49 that includes the extracted coordination complexes (obtained from EXAFS) and also the morphology of aggregates (by SANS). This multimodal method allows unique data treatment, enabling us to combine EXAFS and SANS with theory, 50,51 to elucidate the molecular nature of the zirconium complexes and their subsequent self-assembly into supramolecular aggregates. The results advance our knowledge about the structures and energetics responsible for the phenomena of phase transitions in SALR fluids as they approach critical point concentrations.  4 ] aq,in is the aqueous, initial concentration of Zr(NO 3 ) 4 (prior to the extraction). After centrifugation and separation of the two distinct phases, the organic phase was collected; this is designated as sample no. 1. Next, 3.0 mL of sample no. 1 and 3.0 mL of aqueous, initial phases containing [Zr(NO 3 ) 4 ] aq,in ranging from 0.010 to 0.049 M in 10.5 M nitric acid solution were shaken for 1 h at room temperature. The mixtures were centrifuged to separate the aqueous and organic phases. Aliquots of the organic phases are designated as sample nos. 2−5 ( Table 1). Sample nos. 1−5 were examined by both EXAFS and SANS. Note that the phase separation and thirdphase formation occur only in cases where [Zr(NO 3 ) 4 ] aq,in > 0.049 M. Therefore, the complex concentration in sample no. 5 is very close to the critical concentration, i.e., the highest Zr concentration in the organic phase achievable without phase separation. In this extraction, Zr(NO 3 ) 4 (TBP) 2 is the predominant extracted coordination species in the organic phase, as will be explained in conjunction with analysis of the EXAFS data. The concentrations of Zr(NO 3 ) 4 (TBP) 2 Figure 1a,b show the k 3 -weighted Zr K-edge EXAFS oscillation, k 3 χ(k), as a function of the photoelectron wavenumber, k (Å −1 ), and the corresponding radial structural function, |FT[k 3 χ(k)]|, with the imaginary part of FT[k 3 χ(k)], Im{FT[k 3 χ(k)]} obtained for the organic phases (sample nos. 2−5) as a function of the radial distance from the Zr atom, r + Δr. Here, Δr is the magnitude of the phase shift resulting from a change in the photoelectron wave while traversing the potentials of the absorbing and scattering atoms. The k 3 χ(k) in Figure 1a and |FT[k 3 χ(k)]| in Figure 1b do not depend on the increase in [Zr(NO 3 ) 4 (TBP) 2 ] org,eq in the organic phase. Importantly, this indicates that the local coordination structure of the extracted zirconium complex barely changes at [Zr(NO 3 ) 4 (TBP) 2 ] org,eq ≤ 32 mM, which is close to the highest complex concentration attainable in the organic phase without third-phase formation. Even in the case of the third phase, the EXAFS data were similar to those of the organic phases, as shown in Figure S6. This mirrors the response observed in the corresponding Ce EXAFS data obtained under comparable conditions. 53 In Figure 1b, the |FT[k 3 χ(k)]| shows one significant peak with a shoulder on the longer r + Δr side. The peak, at approximately r + Δr ≈ 0.15 nm and designated as P 1 , is indicated with a thick arrow. Three further peaks are designated as P 2 , P 3 , and P 4 at 0.2 nm < r + Δr < 0.4 nm. Peak P 1 is attributed to two kinds of Zr−O bonds, namely, directly coordinating oxygen atoms of TBPs and nitrate ions. The peaks P 2 , P 3 , and P 4 are attributed to the correlations between Zr−P with TBPs, Zr−N with nitrate ions, and/or multiple scattering of Zr−N−O with nitrate ions, respectively.
Curve fitting for quantitative analysis of the EXAFS data was conducted to obtain the structural parameters of the extracted coordination species in the organic phases. In this analysis, the atomic positions from an optimized coordination structure of Zr(NO 3 ) 4 (TBP) 2 obtained by DFT calculation were used as a structural model for initial input in simulation code FEFF version 8.4. 54 The optimized coordination structure is shown in Figure 2a. All significant scattering paths between Zr and its interrelated atoms giving rise to the distinct peaks in |FT[k 3 χ(k)]| were taken into account in the curve-fitting analysis. The amplitude reduction factor, S 0 2 , was fixed as 0.9 in the fitting procedure. All EXAFS data were well reproduced by the best-fitted simulated curves of k 3 χ(k) (black solid lines), |FT[k 3 χ(k)]| (red solid lines), and Im{FT[k 3 χ(k)]} (green solid lines) with the refined structural parameters, as shown in Figure 1.
Representative structural parameters, refined in the fitting analysis of sample no. 4 ([Zr(NO 3 ) 4 (TBP) 2 ] org,eq = 24 mM), are summarized in Table 2, where CN, r EXAFS , σ DB 2 , and ΔE 0 correspond to the coordination number, bond distance, squared Debye−Waller factor, and energy shift parameter, respectively. The oxygen CN due to TBP molecules (Zr−O TBP path) and nitrate ions (Zr−O NO 3 path), resulting in P 1 , are approximately 2 and 8, respectively. The nitrogen CN due to nitrate ions (Zr−N NO 3 path), phosphorus due to TBPs (Zr− P TBP path), and multiple scattering due to nitrate ions (Zr−N− O multiple path), resulting in P 2 , P 3 , and P 4 , are approximately 4, 2, and 4, respectively. Accordingly, the CNs of the Zr−O TBP and Zr−P TBP paths indicate the coordination of two TBP molecules with the inner sphere of zirconium. The O and N CNs of Zr−O NO 3 , Zr−N NO 3 , and the Zr−N−O multiple paths indicate the bidentate coordination of four nitrate ions with the inner sphere of zirconium. As a result, we found that the coordination species in the organic phase is always Zr-(NO 3 ) 4 (TBP) 2 . This species is therefore expected to be a fundamental building unit of a higher-order, long-range structure formed in the organic phases.
The result of the DFT calculation shows that directly coordinated oxygen atom of only one nitrate ion is slightly more distant than those of the other three nitrate ions coordinated to Zr, thereby giving rise to the shoulder on peak P 1 and to the increase of σ DW 2 of Zr−O NO 3 . The average bond distances evaluated by DFT calculation, r DFT , are shown in Table 2 and are in good agreement with r EXAFS . Therefore, the optimized coordination structure of Zr(NO 3 ) 4 (TBP) 2 shown in Figure 2a can be taken as the actual state in the organic phase.
Overall SANS Features. Figure 3a shows the observed SANS intensity distribution as a function of q, I obs (q), obtained for sample nos. 1−5, which had different concentrations of Zr(NO 3 ) 4 (TBP) 2 in their organic phases, as shown in Table 1.

ACS Central Science
Research Article vector, and λ and 2θ are the wavelength of the incident neutrons and the scattering angle, respectively. The scattering intensity in the low-q region (q < 1.0 nm −1 ) gradually increases with increasing [Zr(NO 3 ) 4 (TBP) 2 ] org,eq , whereas that in the high-q region (q > 1.0 nm −1 ) barely changes as a function of [Zr(NO 3 ) 4 (TBP) 2 ] org,eq . Even in the absence of extracted Zr(NO 3 ) 4 (TBP) 2 complex in the organic phase (sample no. 1), scattering intensity is observed in the SANS profile, particularly in the low-q region. This is attributed to aggregates formed by extracted HNO 3 and TBP, which have been reported by SANS 41 and computational studies. 43−45 Our preceding research concluded that the TBP aggregates arise from the formation of hydrogen-bonding networks consisting of HNO 3 , H 2 O, and the PO group of TBP, leading to an isotropic microemulsion in the organic phase. 44 Because the scattering from the HNO 3 −H 2 O−TBP aggregates for sample nos. 2−5 can vitiate the quantitative structural analyses of the extracted Zr(NO 3 ) 4 (TBP) 2 in the organic phase, the undesirable scattering contributions were subtracted deliber-    2 in the organic phase, determined on the basis of the Debye scattering formula for randomly orientated Zr(NO 3 ) 4 (TBP) 2 using eqs S3, S6, and S7. Solid lines in part b are the best-fit theoretical SANS profiles obtained by using eqs S10−S13 together with the characteristic parameters listed in Table 3.

ACS Central Science
Research Article ately as described in the Methods section (see the Supporting Information). Figure 3b shows I sub (q) plots, which are the SANS profiles of sample nos. 2−5 from which the components of the HNO 3 − H 2 O−TBP aggregates were subtracted from I obs (q). The scattering intensity in the high-q region (q ≥ 3.0 nm −1 ) increases in proportion to the increase in [Zr-(NO 3 ) 4 (TBP) 2 ] org,eq , indicating that the coordination structure of the extracted metal ion species does not depend on the concentration and maintains a specific structure, which is consistent with the EXAFS results discussed above. In addition, the scattering intensity in the low-q region (q < 1.0 nm −1 ) increases out of proportion to [Zr(NO 3 ) 4 (TBP) 2 ] org,eq suggesting that Zr(NO 3 ) 4 (TBP) 2 does not disperse homogeneously and, rather, forms supramolecular structures.
Quantitative Analyses of SANS Intensity Distributions. We first numerically analyzed the form factor of a single Zr(NO 3 ) 4 (TBP) 2 complex in sample nos. 2−5 using the optimized coordination of each atom with the DFT calculation by means of the Debye scattering formula for randomly orientated objects; 50,51 the scattering function is described in detail in the Supporting Information (see eqs S3−S9). The four dashed lines in Figure 3b indicate the same q dependence and show the contribution of the form factors of Zr-(NO 3 ) 4 (TBP) 2 for sample nos. 2−5. The small-angle scattering intensities are proportional to the number density of the extracted Zr(NO 3 ) 4 (TBP) 2 complex, n complex , and agree well with I sub (q) for q ≥ 3.0 nm −1 , whereas those for q < 3.0 nm −1 deviate from each form factor. This deviation is attributed to an ordered (aggregate) structure well beyond the inner-sphere coordination about Zr, and thus analyzing the excess scattering components for q < 3.0 nm −1 over the form factor of the coordination species will allow elucidation of the aggregate architecture. Analyses of I sub (q) with Guinier's law 55 as shown in Figure S7 provide an initial evaluation of the size of the high-order structures formed by Zr(NO 3 ) 4 (TBP) 2 in the organic phase in terms of the radius of gyration, R g . These were determined to be 3.2, 3.3, 4.5, and 5.0 nm for sample nos. 2−5, respectively. The observed power law scattering, I sub (q) ∼ q −D for 1/R g < q (nm −1 ) < 3, shows D ≈ 3 for sample nos. 3−5, which is solid evidence of the existence of an inhomogeneous spatial distribution of Zr(NO 3 ) 4 (TBP) 2 complexes inside the high-order structures.
Given the above findings, we assumed that the hierarchical structure, which comprises primary clusters consisting of extracted Zr(NO 3 ) 4 (TBP) 2 and its aggregates, is formed in the organic phases as shown in Figure 2. A small cluster of Zr(NO 3 ) 4 (TBP) 2 complexes is considered; that is, the N complexes distribute spherically and randomly around the center complex at a distance R (see Figure 2a,b). Then, the small primary clusters consisting of N complexes with radius R S form a large aggregate (supercluster) of M primary clusters with radius R L (see Figure 2c). A similar structure hierarchy was proposed to explain the pseudoliquid phase behaviors of hetero-polyanion cluster compounds. 56 Based on this hierarchical aggregate model, 49 the small-angle scattering intensity can be exactly calculated as described in the Supporting Information, and the SANS data were compared with the model. Sample no. 2 has the lowest concentration of extracted Zr(NO 3 ) 4 (TBP) 2 in the organic phase; consequently, I sub (q) can be reproduced by the contribution of the primary cluster without considering the supercluster, and thus M was fixed as 1 (the model shown in Figure 2b). By comparison, a clustering of the primary clusters is crucial to reproducing I sub (q) when The solid lines in Figure 3b show the best-fit theoretical scattering curves for the hierarchical aggregates, which exhibit good agreement between the model and the experiments. The refined characteristic parameters are summarized in Table 3. Note that the number density of the superclusters, n super , satisfies the relation n complex = n super MN, and the distributions of R S and R L were necessary to reproduce the scattering intensity distribution precisely; thus the Schultz distribution with the corresponding standard deviations, σ L and σ S , respectively, was used. 57 Remarkably, the parameters N and R S relating to the primary cluster barely change in sample nos. 3−5, even with increasing M and R L , which relate to the superclusters. This result suggests that the size and aggregation number of the primary cluster are insensitive to [Zr(NO 3 ) 4 (TBP) 2 ] org,eq for it to be stable in the organic phase when the fundamental building unit is Zr(NO 3 ) 4 (TBP) 2 . We believe that aggregation of the hydrophilic parts of Zr(NO 3 ) 4 (TBP) 2 provides the driving force to form the primary cluster, and that this is driven by coordinated nitrate ions, which form a hydrogen-bonding network involving extracted HNO 3 and H 2 O. 43,44,58 Consequently, the surface of the primary cluster should be stably covered with the hydrophobic butyl groups of TBP in the organic phase as will be explained by the MD simulations described below. Subsequently, we speculate that the process of primary cluster aggregation is similar to the formation of concentrated reverse micelle systems, as supported by previous studies. 41,59,60 MD Simulations and the Structure of the Primary Clusters. MD simulations were carried out for each of the systems given in Table 1. Trajectories from the simulations indicated that the Zr complexes tended to aggregate to form clusters, namely, the primary clusters described above. To analyze this clustering quantitatively, we calculated the Zr−Zr radial distribution function, g(r Zr−Zr ), and the corresponding coordination number of Zr with the other Zr, CN(r Zr−Zr ), shown for system no. 5 in Figure 4a. The g(r Zr−Zr ) function consists of a large peak followed by an extended tail. This large peak, which corresponds to the closest Zr−Zr distance on average, reaches its maximum at approximately 0.88 nm. This value is close to the R S value of 0.92 ± 0.03 nm in the characteristic parameters of the hierarchical aggregates model.

Research Article
To define a cluster, we used the algorithm of Sevick, 61 with the criterion that two Zr complexes are connected if the Zr−Zr distance is less than 1.2 nm. This distance is approximately where the first peak of the radial distribution function ends. This choice is somewhat arbitrary, but the resulting analysis is in reasonable agreement with the conclusions drawn from observing many snapshots. The distribution of the primary cluster aggregation numbers is shown in Figure 4b for system no. 5. It is clear that primary clusters form, several with aggregation numbers of 5−7, in agreement with the SANS analysis. The shapes of the primary clusters, according to our simulations, are not uniform, as they vary from spherical to elongated aggregations. This observation is supported by the coordination number analysis, CN(r Zr−Zr ), given in Figure  Smaller clusters are also seen, but it is noteworthy that clusters of three or four Zr atoms are relatively uncommon. To check that this was not an equilibration issue, we ran 10 independent MD simulations, starting from random initial configurations. All the simulations used the same simulation protocols as described in the Supporting Information (Computational Methods). The results, which are shown in Figure S2 in the Supporting Information, indicate that the cluster size distribution is identical, to within statistical uncertainty, for all 10 systems.
The system sizes are too small to observe and analyze hierarchical clustering of these primary clusters, but the long broad tail in the radial distribution function in Figure 4a indicates that the primary clusters in the simulated systems are not homogeneously distributed. An examination of simulation snapshots shows the loose aggregation of a small number of primary clusters. An example is given in Figure 5a, where two primary clusters are seen in close proximity, providing an indication of incipient hierarchical structuring.
We now consider the structure of the primary clusters. Examples are shown in Figure 5b. Within a primary cluster there are hydrogen bonds, involving uncoordinated TBP, HNO 3 , and H 2 O. These play a role in binding the Zr complexes together. There are, in addition, nonspecific polar interactions, for example, involving the charged species within the Zr complexes. Hydrophobic interactions between the nonpolar groups are also in evidence. Turning now to the interactions between primary clusters, any inferences based on simulation results must be extremely tentative, due to the system size issues noted previously. That said, the indications are that there is no direct hydrogen bonding between primary clusters and no specific bridging molecules. This is exemplified by Figure 5a. The available simulation data suggest, instead, that primary clusters interact via a combination of nonspecific electrostatic interactions between the charged groups and hydrophobic interactions between the aliphatic groups. It is noteworthy that recent research by Servis et al., 62 on the structure of the aggregates formed by uranyl nitrate and TBP, also indicated that such electrostatic interactions play a key role in stabilizing the clusters.
Finally, we return to the jagged nature of the cluster distribution (Figure 4b). Our analysis shows that this is not related to the choice of the distance criterion for defining a cluster. The results are insensitive to this. As discussed previously, the fact that 10 independent runs give the same distribution leads us to believe this is an equilibrium distribution. To gain insight into what factors might lead to this distribution, we carried out a hydrogen-bonding analysis, counting the number of hydrogen bonds between NO 3 − , HNO 3 , and H 2 O. The criterion for determining the presence of a hydrogen bond is that the donor−acceptor distance should be no more than 0.35 nm, and the acceptor−donor−hydrogen angle should be no more than 30°. The oxygen atoms that have covalent bonds with the hydrogen atoms in the H 2 O and HNO 3 molecules were regarded as potential donors, and the electronegative atoms that possess a lone electron pair were regarded as potential acceptors. Stable primary clusters, defined for these purposes as clusters with a lifetime greater than 1 ns, were identified and indexed for hydrogen bonding. These primary clusters contained 1, 2, 3, 5, and 6 Zr complexes. Tetramers are apparently too unstable to survive for this length of time. The average number of hydrogen bonds per Zr atom is plotted against cluster size in Figure 6a.
While the number of hydrogen bonds increases when two monomers form a dimer, there is no increase in the number of hydrogen bonds on forming higher-order clusters. It would thus appear that hydrogen bonding is probably not a governing factor in determining the number of Zr complexes in a primary cluster.
Another factor that might influence the aggregation number is the degree to which the octane solvent is shielded from polar

ACS Central Science
Research Article groups in the clusters. Such interactions would be hydrophobically unfavorable. We therefore calculated the radial distribution functions for octane carbon atoms surrounding a central Zr atom and for octane carbon atoms surrounding the N atom of a nitrate ligand. The plots are shown in Figure 6b,c. These two plots show that there is considerably more unfavorable polar−nonpolar contact for the trimer, as compared to either the dimer or the hexamer. Thus, it would appear that the polar groups in the trimer are less well shielded from the diluent than is the case for either the dimer or the hexamer, and this may partly explain the relative stabilities of these clusters.
Thermodynamic Considerations. We relate the extracted parameters from the model fittings of I sub (q) in SANS to the thermodynamic properties of the system. We start with the van't Hoff equation given by 63 where Π is osmotic pressure, v is the occupied volume of the solute per molecule, k B is the Boltzmann constant, and T is the thermodynamic temperature. The osmotic pressure can be related to the forward scattering intensity, I sub (q = 0), by 64 with the difference between the scattering length densities of the Zr(NO 3 ) 4 (TBP) 2 complex and the matrix, Δρ, and ϕ Zr(NO 3 ) 4 (TBP) 2 . In the studied system, the volume v in eq 1 can be assumed to be the occupied volume per supercluster, and hence, v ≈ 1/n super . The calculated v and the osmotic pressure obtained by eq 1, Π van't-Hoff , are listed in Table 4. Because the vapor pressure of octane at 20°C is 1.33 kPa, and the matrix consists of a mixture of extracted HNO 3 , extracted H 2 O, uncoordinated TBP, and n-octane-d 18 (Table 1), the calculated Π van't-Hoff is a good estimation. The osmotic pressures can also be calculated with eq 2 by using I sub (q = 0) from the experimental data, which were estimated by Guinier plots (Table 4). The osmotic pressures from I sub (q = 0), Π I(0) , and Π van't-Hoff agree within 10%, which indicates that our structural model is correct. The van der Waals equation of state with hard spheres with radius R HS modifies eq 1 to give 63

ACS Central Science
Research Article i k j j j y where β = 4ω with ω = 4πR HS 3 /3 and α = 27k B T c β/8 with the critical temperature, T c , for phase separation. By assuming T c / T = 1, R HS for each sample can be calculated with assigning Π I(0) and v, which are shown in Table 4. The R HS values agree well with the R L values. The parameter R HS is slightly affected by T c /T around T c /T = 1, although the variation of R HS with a ±10% variation of T c /T = 1 is less than 5%. This result suggests that R L obtained by the model is reasonable.
It is worth noting, at this point, some curious features of this aggregation process. First, as shown in Table 1, the molar concentrations of zirconium in the organic phase are relatively small. The highest concentration studied was sample no. 5, where [Zr(NO 3 ) 4 (TBP) 2 ] org,eq = 0.032 M, and as noted previously, third-phase formation sets in at concentrations just a little greater than this. A normal solution at such low concentrations might be expected to be approximately idealdilute, with Henry's law applying, but this is clearly far from the case here. To account for the clustering and for the phase transition, it must be the case that the mixture is very far from ideal, implying remarkably strong interactions between the zirconium complexes and, indeed, between the primary clusters. As shown in Table 3, the aggregation number of the primary clusters is approximately 7 and does not change with increasing [Zr(NO 3 ) 4 (TBP) 2 ] org,eq . This is reminiscent of spherical micelles formed in lyotropic systems, where again there is a favored aggregation number, but it is not clear why this should apply here, where simulation suggests the primary clusters to be far from nicely organized spheres. This is certainly worthy of future study.
Even though the concentrations of primary clusters are small, the results of Table 3 show that they aggregate, indicating strong interactions between them. One may note that the primary clusters, with solvating water, nitric acid, and uncoordinated TBP molecules, are polar entities immersed in a nonpolar, low dielectric medium. This may mean that Coulomb interactions between primary clusters are strong, contributing considerably to the attractive forces between them. Once these superclusters have formed, their concentration is very small, and as shown above, the osmotic pressure may be well estimated by assuming one has an ideal suspension of superclusters.
For values of [Zr(NO 3 ) 4 (TBP) 2 ] org,eq a little greater than that of sample no. 5, the system separates into two phases. One phase is diluent-rich with little Zr, while the coexisting phase, the third phase, is rich in Zr and has a reduced concentration of diluent. Our data indicate that the aggregation number, M, of the superclusters rises rapidly with [Zr(NO 3 ) 4 (TBP) 2 ] org,eq , presaging an oncoming phase instability as more and more primary clusters condense. Unfortunately, these data alone are insufficient to cast light on the nature of this phase transition, for this would require information about the structure of the third phase. In the literature, this third phase transition has been treated both as a gas−liquid transition and as a reverse micelle to microemulsion or gel transition. Taken literally, 65,66 the gas−liquid transition model would correspond to a gas of primary clusters (or reverse micelles) condensing to a liquid of primary clusters, but with the primary clusters retaining their integrities. In the second scenario, the primary clusters/reverse micelles would lose their individual identities and merge to form a connected, bicontinuous structure. Of course, another possibility is that the third phase is a bicontinuous emulsion of polar and nonpolar regions, such as appears to be the case in metal-free systems, 44 with distinct Zr clusters embedded in the polar regions of this emulsion. Further study is needed to clarify this situation, but what we have shown is that the strong attractions between primary Zr clusters are the driving force that makes the "gas" phase unstable.

■ CONCLUSION
We have clarified the microscopic structure of the organic phases containing extracted Zr(NO 3 ) 4 (TBP) 2 using a method combining EXAFS with DFT calculations, MD simulations, and SANS observations. The coordination structure of the extracted Zr(NO 3 ) 4 (TBP) 2 in the organic phase does not depend on the concentration of Zr(NO 3 ) 4 (TBP) 2 and maintains a specific structure as a fundamental building unit to form high-order (aggregate) architectures. Moreover, quantitative analyses of SANS data using an accurate evaluation of the incoherent scattering intensity and a form factor of randomly oriented Zr(NO 3 ) 4 (TBP) 2 provide a clearer picture of the high-order structures on the basis of scattering theory. Zr(NO 3 ) 4 (TBP) 2 was found to form a hierarchical aggregate composed of small primary clusters comprising a supercluster (Figures 2 and 5). A hybrid interaction consisting of a hydrogen-bonding network for the primary clusters and an attractive interaction for the superclusters induces the formation of the hierarchical aggregate. An increase in [Zr(NO 3 ) 4 (TBP) 2 ] org,eq has little effect on the size of the primary cluster but increases the size of the superclusters. Furthermore, MD simulations provided direct evidence about the formation of primary clusters, which are assembled by the hydrogen-bonding network involving uncoordinated TBP, HNO 3 , and H 2 O. In the clustering analyses of the MD snapshots, the distribution of the primary cluster aggregation numbers of 5−7 agreed well with the SANS data analysis. Accordingly, we conclude that growth of the superclusters is due to an increase in the number of small primary clusters, causing third-phase formation. In this paper, we quantitatively applied the simple aggregation model to analyze SANS profiles. The number of the clusters and their sizes, obtained by SANS data analysis, are fully compliant with van der Waals equations of state. These findings suggest that such behaviors are general ones, extending beyond zirconium in the PUREX process to a wide variety of extraction systems (e.g., DIAMEX, TRUEX, ALSEP) 38 that show multiscale ordering. There is then the intriguing possibility that, by tailoring the interactions between these superclusters, we may afford a new entry into the design of next-generation ionic separation techniques with complex fluids typical of soft matter chemistry.
It is known that the addition of phase modifiers, such as aliphatic alcohols, into organic phases during liquid−liquid extraction allows a higher loading of metal ions (without thirdphase formation) than in unmodified systems. 67,68 The fact that the separation selectivity and efficiency are also impacted leads us to speculate that modifiers contribute to the formation of hydrogen-bonding networks inside the primary clusters, and these networks vary in accordance with the intra-and intercluster interactions. In view of the largely empirical use of modifiers in practical chemical separations, the combination of SANS, EXAFS, and computational methods used here can help to elucidate the role of the modifiers in other solvent extraction systems, contributing to the development of highperformance processes in chemical separation science.

ACS Central Science
Research Article ■ ASSOCIATED CONTENT

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acscentsci.8b00669.
Method section and additional data and figures including cluster probabilities, configurations, radial distribution function, chemical structure, EXAFS data, Guinier plots, theoretical modeling of the hierarchical aggregates for SANS data analysis, simulated SANS profiles, and coordination structure of zirconium nitrate complex in initial aqueous phase (PDF) ■ AUTHOR INFORMATION Corresponding Authors *E-mail: motokawa.ryuhei@jaea.go.jp. *E-mail: andrew.masters@manchester.ac.uk. ORCID