Experimental and modelling evidence for structural crossover in supercritical CO 2
Cillian J. Cockrell, Oliver Dicks, Ling Wang, Kostya Trachenko, Alan K. Soper, Vadim V. Brazhkin, Sarantos Marinakis
EExperimental and modelling evidence for structural crossover in supercritical CO Cillian J. Cockrell, ∗ Oliver Dicks, † Ling Wang, ‡ and Kostya Trachenko § School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London, E1 4NS, UK.
Alan K. Soper ¶ ISIS Facility, STFC Rutherford Appleton Laboratory,Harwell Campus, Didcot, Oxon OX11 0QX, UK.
Vadim V. Brazhkin ∗∗ Institute for High Pressure Physics, RAS, 108840, Moscow, Russia.
Sarantos Marinakis †† School of Health, Sport and Bioscience, University of East London,Stratford Campus, Water Lane, London E15 4LZ, UK andDepartment of Chemistry and Biochemistry, School of Biological and Chemical Sciences,Queen Mary University of London, Joseph Priestley Building, Mile End Road, London E1 4NS, UK (Dated: May 1, 2020)Physics of supercritical state is understood to a much lesser degree compared to subcritical liquids.Carbon dioxide in particular has been intensely studied, yet little is known about the supercriticalpart of its phase diagram. Here, we combine neutron scattering experiments and molecular dynamicssimulations and demonstrate the structural crossover at the Frenkel line. The crossover is seen atpressures as high as 14 times the critical pressure and is evidenced by changes of the main featuresof the structure factor and pair distribution functions.
Usage:
Secondary publications and information retrieval purposes.
PACS numbers:
May be entered using the \pacs{ command.
INTRODUCTION
Supercritical fluids have unique properties that haveled to a rich variety of applications [1]. Rare gases, nitro-gen, CO and H O are among the most common super-critical fluids. CO in particular is an important green-house gas of Earth’s atmosphere, and in its supercriticalstate is the main component (97%) in the atmosphereof Venus. Supercritical CO is used in a great varietyof applications (see, e.g., applications in solubility, syn-thesis and processing of polymers [2–4], dissolving anddeposition in microdevices [5], green chemistry and sol-vation [6–12], green catalysis [9, 13–15], extraction [16],chemical reactions [17], green nanosynthesis [18] and sus-tainable development including carbon capture and stor-age [19]). It has been widely appreciated that improvingfundamental knowledge of the supercritical state is im-portant for the reliability, scale-up and widening of theseapplications (see, e.g., Refs [1, 7, 9, 12, 15, 17]).Compared to subcritical liquids, the supercritical stateis not well understood. Traditional understandingamounted to a general assertion that this state is phys-ically homogeneous, with no qualitative changes takingplace anywhere above the critical point [1]. The first chal-lenge to this view was the Widom Line (WL). Close tothe critical point, the WL characterises persisting near-critical anomalies such as the maximum in the heat ca-pacity [20], which can be used to stratify different statesin the supercritical region. A different subsequent pro- posal was based on the Frenkel line (FL) separating twodistinct states in the supercritical state with liquid-likeand gas-like dynamics. Differently from the WL, the FLextends to arbitrarily high pressures and temperatures(as long as chemical bonding is unaltered), is unrelatedto the critical point and exists in systems with no boilingline or critical point [21–23]. The FL is also of practi-cal importance because it corresponds to the solubilitymaxima in supercritical CO [24].Here, we combine neutron scattering experiments andmolecular dynamics (MD) simulations and show evidencefor the structural crossover of supercritical carbon diox-ide at the Frenkel line. The crossover extends to pressureas high as 14 times the critical pressure and is evidencedby changes of the main features of the structure factorand pair distribution functions. The neutron scatteringexperiments evidencing a crossover at highly supercriti-cal pressure are the first of its kind for CO . METHODS
We recall that particle dynamics combine solid-likeoscillations around quasi-equilibrium positions and dif-fusive jumps between different positions below the FL,the typical character of molecular motion in liquids [25].Above the line, particle dynamics lose this oscillatorycomponent and become purely diffusive. This givesa practical criterion to calculate the FL based on the a r X i v : . [ c ond - m a t . o t h e r] A p r disappearance of minima of the velocity autocorrela-tion function (VAF) [22]. This criterion coincides withthe thermodynamic criterion c v = 2 k B correspondingto the disappearance of transverse-like excitations in amonatomic system [26]. Since structure and dynamicsare related [27], the FL crossover was predicted to resultin a crossover of the supercritical structure .The pressures we consider for both experiments andMD simulations are 500 and 590 bar. The FL in CO was previously calculated using the VAF criterion [24],giving us the following two state points of the predictedcrossover: (500 bar, 297 K) and (590 bar, 302 K). We re-call that the FL extends to arbitrarily high pressure andtemperature above the critical point, but at low temper-ature it touches the boiling line at around 0.8 T c , where T c is the critical temperature [22] (note that the systemdoes not have cohesive liquid-like states at temperaturesabove approximately 0.8 T c [28], hence crossing the boil-ing line at around 0 . T c and above can be viewed as agas-gas transition [22].) The critical point of CO is (73.9bar, 304.3 K), hence our state points correspond to near-critical temperatures and pressures well above critical. Inthis regard, we note that the supercritical state is oftendefined as the state at P > P c and T > T c . This defini-tion is loose, not least because an isotherm drawn on ( P , T ) diagram above the critical point crosses the meltingline, implying that the supercritical state can be foundin the solid phase. As a result, one can meaningfullyspeak about near-critical part of the phase diagram onlywhen discussing the location of the supercritical stateon the phase diagram [29]. As far as our state points areconcerned, they correspond to temperatures much higherthan the melting temperature and pressures extending to14 times the critical pressure where near-critical anoma-lies are non-existent [29].A cylinder of carbon dioxide was obtained from BOC,CP grade, and used without further purification. Thepressure of the cylinder was around 50 bar and a SITECintensifier and a SITEC hand pump gas was used to raisethe pressure. Capillaries were used to connect intensifiermanifold system to the cell. The flat plate pressure cellwas made from an alloy of Ti and Zr in the mole ratio0.676:0.324, which contributes almost zero coherent scat-tering to the diffraction pattern [30]. The cell consistedof a flat section that was 12 mm thick and had four 6 mmdiameter holes running through it, so the occupied gasspace was 6 mm thick and the wall thickness was 3 mmeither side. The container was placed at right angle tothe neutron beam, which was approximately 30 mm x30 mm in cross section. A bottom loading closed cyclehelium refrigerator was used to control the temperaturewithin ± ∼
20 mbar toprovide temperature uniformity. The employed temper-atures and pressures are shown in Table I, where thedensities were calculated from the data available in the NIST database [31].
TABLE I. T − P − d state points for neutron scatteringmeasurements. The values of d are taken from [31]. T exp (K) P exp (bar) d (g/mL) P exp (bar) d (g/mL)250 500 1.1676 590 1.1821270 500 1.1131 590 1.1306290 500 1.0573 590 1.0784310 500 1.0003 590 1.0257330 500 0.9426 590 0.9729340 500 0.9137 - -350 500 0.8848 590 0.9204360 500 0.8560 - -370 500 0.8276 590 0.8688380 500 0.7996 590 0.8436390 500 0.7722 590 0.8188 Total neutron scattering measurements were per-formed on the NIMROD diffractometer at the ISIS pulsedneutron source [32]. Absolute values of the differentialcross sections were obtained from the raw scattering databy normalising the data to the scattering from a slab ofvanadium of known thickness, and were further correctedfor background and multiple scattering, container scat-tering and self-attenuation, using the Gudrun data anal-ysis program [33]. Finally the data were put on absolutescale of barns per atom per sr by dividing by the numberof atoms in the neutron beam (1 barn = 10 − m ). F ( Q ) = 19 b H CC ( Q ) + 49 b H OO ( Q ) + 49 b C b O H CO ( Q )(1)where b α is the neutron scattering length of atom α ,and the partial structure factor H αβ ( Q ) is the three-dimensional Fourier transform of the corresponding site-site radial distribution function: H αβ ( Q ) = 4 πρ (cid:90) ∞ r ( g αβ ( r ) −
1) sin
QrQr dr (2)and ρ is the atomic number density. Note that the H OO ( Q ) and H CO ( Q ) terms include both the intra- andinter-molecular scattering. The results are shown inFig. 1.The molecular dynamics (MD) simulation packageDL POLY [34] was used to simulate a system of 30752CO particles with periodic boundary conditions. Thepotential for CO is a rigid-body non-polarizable po-tential based on a quantum chemistry calculation, withthe partial charges derived using the distributed multi-pole analysis method [35]. The electrostatic interactionswere evaluated using the smooth particle mesh Ewaldmethod in MD simulations. The potential was derivedand tuned using a large suite of energies from ab initio -0.4-0.200.20.40.60.8 10 -0.4-0.200.20.40.60.8
250 300 350 4001.661.681.71.721.741.761.781.81.821.841.86
FIG. 1. Weighted sum of the experimental weighted partialstructure factors F ( Q ) for CO at (a) 500 bar and (b) 590 bar.The curves for higher temperatures have been shifted alongthe y-axis by 0.05 per set. (c) Position of the maximum of thefirst peak in the total (weighted sum) experimental structurefactor as a function of temperature for CO at 500 bar (blue)and 590 bar (red). The straight lines are visual guides. density functional theory calculations of different molec-ular clusters and validated against various sets of ex- perimental data including phonon dispersion curves and P V T data. These data included solid, liquid and gasstates, gas-liquid coexistence lines and extended to high-pressure and high-temperature conditions [35]. We alsoused another rigid-body non-polarizable potential devel-oped by Zhang and Duan [36] and found the same results.The MD systems were first equilibrated in the constantpressure and temperature ensemble for 500 ps. The datawere subsequently collected from production runs in theconstant energy and volume ensemble. In order to reducenoise and see the crossover clearly, data were averagedover 500,000 frames, involving production runs of further500 ps.
RESULTS AND DISCUSSION
Before analyzing the data, we recall that the FL cor-responds to the qualitative change of particle dynamics:from combined solid-like oscillatory and diffusive dynam-ics below the line to purely diffusive gas-like dynamicsabove the line. Therefore, the supercritical structure ispredicted to show the crossover between the liquid-likeand gas-like structural correlations. This is predictedto be the case for functions characterizing the structure,such as pair distribution function (PDF) and structurefactor (SF). In our analysis, we focus on meaningful fea-tures such as maxima positions of PDFs and SFs.The experimental weighted sum of the partial struc-ture factors are plotted in Fig. 1 for two pressures. Weplot the first peaks position of SFs vs temperature in Fig.1 and observe that it undergoes the crossover at temper-atures close to 320 K and around 12% larger than the FLcrossover temperature predicted from the VAF criterionmentioned earlier.The SFs were Fourier transformed to obtain the exper-imental PDF. As with previous experimental and mod-elling results on Ar [27], Ne [37], CH [38], and especiallyon H O [39] pronounced changes of first peak positionin the PDF with temperature are observed, indicatinga well-defined crossover. When a system is compressedor expanded, one expects the first nearest-neighbour dis-tance, r fnn (given by the radial position of the first peakin g(r)), and the system’s “length” ( V / ) to be propor-tional to each other unless the system undergoes a struc-tural change. In other words the system structure under-goes uniform compression. The first PDF peak positiondivided by the position ( r ) at the Frenkel temperature vsthe cube root of the volume divided by the volume ( V ) atthe same reference temperature is shown in Fig. 2. For asystem undergoing uniform compression, V /V and r/r will be equal. If there is a phase transition at higher den-sities, as there is in liquids across the melting line, thisrule cannot be extrapolated down to arbitrarily low vol-umes and hence there will be an intercept: V = αr + β ,upon which the gradient of V /V vs. r/r will depend.However as long as no structural changes occur, the gra-dient will remain constant. Specifically, in a simple cubiccrystalline solid (atomic packing fraction 0.52) the con-stant of proportionality between V /N and r is unity, in aFCC lattice (packing fraction 0.74) the constant is 0.89,and in a diamond cubic lattice (packing fraction 0.34),the constant is 1.2. In gases, the fnn distance is largelydetermined by the size, geometry, and interaction of theconstituent molecules (see, e.g., [40]) rather than the den-sity. This linear relationship has been experimentally ob-served in molten group 1 elements [41, 42] and liquid CS [43]. The fnn distance is most readily extracted from thepartial C-C PDFs. The experimental data give the totalPDF, but the peak corresponding to the fnn distance isnot profoundly changed, therefore the total PDF givesa qualitative approximation of the fnn distance. In Fig.2 we observe the crossover of the first PDF peak at atemperature within 10% of the predicted crossover tem-perature at the FL, signified by the change of gradients.This qualitative behaviour is seen more clearly in the MDresults (Fig. 4). r / r ( a ) P =
500 bar ( T = K ) F L ( V / V ) r / r ( b ) P =
590 bar ( T = K ) F L FIG. 2. First peak positions of experimental weighted sumPDF for CO at (a) 500 bar and (b) 590 bar as a function ofvolume. The vertical dashed lines show the reduced volumesat the FL, and the straight lines are visual guides. We now discuss the MD results. Examples of C-CPDFs from MD simulations are shown in Fig. 3. Weobserve a reduction in height, and corresponding broad-ening of peaks with increasing temperature as expected.The steepness of the first peak is related with the soft-ness of the effective intermolecular potential, and its re- T = = = = = = = g ( r ) ( a ) P =
500 bar T = = = = = = = ( A ) g ( r ) ( b ) P =
590 bar
FIG. 3. (Colour online) Evolution of the simulated C-C pairdistribution functions with temperature at (a) 500 bar and(b) 590 bar. duction can be quantitatively related to the reduction ofthe viscosity [44]. Fig. 4 displays the radial positionsof the first PDF peaks as a function of volume, as dis-cussed above, which shows a crossover at densities nearthe FL. Because of the reduced noise and abundance oftemperature points we can perform statistical analysisof the data to quantify the crossover. We see the samebehaviour, including a much clearer crossover, for bothpressures. The constant of proportionality between
V /N and r is ≈ .
4, implying a much more open arrange-ment than the crystal systems quoted above. This is inaccordance with the density of CO at the FL (23346mol/m ), less than half that of water (56501 mol/m ) atthe same pressure and temperature. In order to quantifythe crossover, we fitted the data to two different typesof function. The first was a single functional dependenceover the entire range. In order to avoid the extrapola-tion errors associated with high order polynomials, thetrial functions we used were quadratic, or log plus linear: f ( x ) = a + bx + c log( x ), with a , b , and c the fitting pa-rameters. The second set of functions was linear below acertain crossover volume V c , and either quadratic or logplus linear above that volume ( i.e. a piecewise function): f ( x ) = Θ( V c − V )[ a + bx ] + Θ( V − V c )[ α + βx + γ log( x )],with Θ[ V ] the Heaviside step function and a , b , β , γ , and V c the fitting parameters ( α depends on the other pa-rameters in order to ensure continuity of the function).Generally speaking, adding more parameters to a fit-ting function improves the numerical quality of the fit. A priori , one can penalise having too many parameters- this prevents the extreme situation of a perfect fit ac-quired using a piecewise function with a number of sub-domains equal to the number of data. The two closelyrelated quantitative measures of goodness of fit withpenalty terms for the number of parameters are the AIC(Akaike Information Criterion) [45] and BIC (BayesianInformation Criterion) [46]. Applied to our data, at bothpressures and with the quadratic and log plus linear vari-ants, the AIC and BIC were substantially lower than − V c ). This volumeis shown in the vertical dotted line (Fig. 4) and corre-sponds at both pressures to a temperature close to 350 K,which is within 12-15% of the predicted crossover value.Also plotted as insets in Fig. 4 are the residuals of thelow-volume linear fits which show a sharp and suddenincrease above the crossover volume, which would notbe the case if we had simply interpolated a straight linebetween non-linear data.Fig. 5 shows theoretical PDF peak heights. We notethat the PDF peak heights of a solid h = g ( r peak ) − h ∝ − log T with h = g ( r ) − [38],and the combination of neutron and Raman scattering insupercritical N [49] . Only one small-angle neutron scat-tering experiment had been used to study the FL in CO in the vicinity of the critical point only [50]. Our currentneutron scattering experiment detecting the crossover atthe FL at highly supercritical pressures is the first ofits kind and importantly widens the range of techniquesused to detect the FL. It will stimulate further neutronscattering experiments in important systems such as su-percritical H O where a pronounced crossover at the FLwas recently predicted on the basis of MD simulations r / r ( a ) P =
500 bar F L ( T = K ) F i tt ed C r o ss o v e r ( T = K ) ( V / V ) / Δ ( r / r ) ( V / V ) r / r ( b ) P =
590 bar F L F i tt ed C r o ss o v e r ( T = K ) ( V / V ) / Δ ( r / r ) FIG. 4. First peak position of simulated C-C PDF. Thestraight lines are fitted to data below the FL and serveas visual guides. The vertical dashed lines show the fittedcrossover volume and the volume at the FL. The insets showthe relative trend of the residuals of the linear fit. [39].In summary, our combined neutron scattering andmolecular dynamics simulations study has detected thestructural crossover in CO at pressures well above thecritical pressure and temperatures well in excess of melt-ing temperature. The crossover is seen in the main fea-tures of the SF and PDFs and corresponds to the pre-dicted crossover at the FL. Apart from the fundamentalimportance of understanding the supercritical state, theFL corresponds to the solubility maxima of several so-lutes in supercritical CO [24] and is therefore of practicalimportance.Neutron beam time at ISIS and project fundingwas provided by the Science and Technology FacilitiesCouncil (RB1720056). We acknowledge Chris Good-way, Thomas Headen and Damian Fornalski (Ruther-ford Appleton Laboratory) for help with the experi-ments. This research utilized Queen Mary’s MidPluscomputational facilities, supported by QMUL Research-IT, http://doi.org/10.5281/zenodo.438045. ∗ [email protected] † [email protected] ‡ [email protected] g ( r ) - ( a ) P =
500 bar F L
200 300 400 500 6000.40.60.81.01.2 T ( K ) g ( r ) - ( b ) P =
590 bar F L FIG. 5. First peak height of C-C PDF for simulated CO at(a) 500 bar and (b) 590 bar as a function of volume. Thevertical dashed lines show the temperatures of the predictedcrossover, and the solid lines show the fitted straight lines. § [email protected] ¶ [email protected] ∗∗ [email protected] †† [email protected][1] E. Kiran, P. G. Debenedetti and C. J. Peters, Super-critical Fluids: Fundamentals and Applications (NATOScience Series E: Applied Sciences vol 366) (Boston:Kluwer, 2000).[2] T. Sarbu, T. Styranec and E. J. Beckman, Nature ,165 (2000).[3] J. M. DeSimone, Z. Guan and C. S. Elsbernd, Science , 945 (1992).[4] A. I Cooper, J. Mater. Chem. , 207 (2000).[5] J. M. Blackburn, D. P. Long, A. Caba˜nas and J. J.Watkins, Science , 141 (2001).[6] J. M. DeSimone, Science , 799 (2002).[7] C. A. Eckert, B. L. Knutson and P. G. Debenedetti, Na-ture , 313 (1996).[8] C. J. Li and B. M. Trost, PNAS , 13197 (2008).[9] W. Leitner, Acc. Chem. Res. , 746 (2002).[10] P. T. Anastas and M. M. Kirchhoff, Acc. Chem. Res. ,686 (2002).[11] P. Anastas and N. Eghbali, Chem. Soc. Rev. , 301(2010).[12] E. J. Beckman, J. Super. Fluids , 121 (2004).[13] D. J. Cole-Hamilton, Science , 1702 (2003).[14] P. G. Jessop, Y. Hsiao, T. Ikarita and R. Noyori, J. Am.Chem. Soc. , 344 (1996).[15] P. G. Jessop, T. Ikariya and R. Noyori, Chem. Rev. ,475 (1999). [16] E. Reverchon, J. Super. Fluids , 1 (1997).[17] P. E. Savage, S. Gopalan, T. I. Mizan, C. J. Martino, E.E. Brock, AIChE J. , 1723 (1995).[18] J. A. Dahl, B. L. S. Maddux and J. E. Hutchison, Chem.Rev. , 2228 (2007).[19] C. Song, Catalysis Today , 2 (2006).[20] L. Xu, P. Kumar, S. V. Buldyrev, S.-H. Chen, P. H.Poole, F. Sciortino and H. E. Stanley, PNAS , 16558(2005).[21] V. V. Brazhkin, Yu. D. Fomin, A. G. Lyapin, V. N.Ryzhov and K. Trachenko, Phys. Rev. E , 031203(2012).[22] V. V. Brazhkin, Yu. D. Fomin, A. G. Lyapin, V. N.Ryzhov, E. N. Tsiok and K. Trachenko, Phys. Rev. Lett. , 145901 (2013).[23] V. V. Brazhkin and K. Trachenko, Physics Today , 68 (2012).[24] C. Yang, V. V. Brazhkin, M. T. Dove and K. Trachenko,Phys. Rev. E , 012112 (2015).[25] J. Frenkel, Kinetic Theory of Liquids (New York: Dover1955).[26] K. Trachenko and V. V. Brazhkin, Rep. Prog. Phys. ,016502 (2016).[27] L. Wang, C. Wang, M. T. Dove, Yu. D. Fomin, V. V.Brazhkin and K. Trachenko, Phys. Rev. E , 196 (1993).[29] V. V. Brazhkin, A. G. Lyapin, V. N. Ryzhov, K. Tra-chenko, Yu. D. Fomin, E. N. Tiosk, Physics Uspekhi , 1061 (2012).[30] A. K. Soper and D. T. Bowron, Chem. Phys. Lett. , 214507(2005).[37] C. Prescher, Yu. D. Fomin, V. B. Prakapenka, J. Stefan-ski, K. Trachenko, and V. V. Brazhkin, Phys. Rev. B ,134114 (2017).[38] D. Smith, M. A. Hakeem, P. Parisiades, H. E. Maynard-Casely, D. Foster, D. Eden, D. J. Bull, A. R. L. Marshall,A. M. Adawi, R. Howie, A. Sapelkin, V. V. Brazhkin, andJ. E. Proctor, Phys. Rev. E , 052113 (2017).[39] C. Cockrell, O. Dicks, V. V. Brazkhin, and K Trachenko,arXiv:1905.00747.[40] J. M. Ziman, Models of Disorder, (Cambridge: Cam-bridge University Press, 1979).[41] K. Tsuji, Y. Katayama, Y. Morimoto, and O. Shimo-mura, J. Non-Crystalline Solids , 295-298 (1996).[42] K. Tsuji and Y. Katayama, J. Phys: Condensed Matter , 6085-6103 (2003).[43] S. Yamamoto et al, J. Chem. Phys. , 144511-5 (2006).[44] J. Krausser, K. H. Samwer, and A. Zaccone, PNAS , ,712-723 (1974).[46] G. Schwarz, Ann. Statist. , 461-464 (1978).[47] A. A. Maradudin, E. W. Montroll, G. H. Weiss, and I.P. Ipatova, Theory of Lattice Dynamics in the Harmonic Approximation (New York: Academic, 1971).[48] L. Wang, C. Yang, M. T. Dove, V. V. Brazhkin, and K.Trachenko, J. Phys.: Condens. Matt. , 225401 (2019).[49] J. E. Proctor, C. G. Pruteanu, I. Morrison, I. F. Croweand J. S. Loveday, J. Phys. Chem. Lett. , 6584 (2019).[50] V. Pipich and D. Schwahn, Phys. Rev. Lett.120