Theoretical Estimate of the Glass Transition Line of Yukawa One-Component Plasmas

The mode coupling theory of supercooled liquids is combined with advanced closures to the integral equation theory of liquids in order to estimate the glass transition line of Yukawa one-component plasmas from the unscreened Coulomb limit up to the strong screening regime. The present predictions constitute a major improvement over the current literature predictions. The calculations confirm the validity of an existing analytical parameterization of the glass transition line. It is verified that the glass transition line is an approximate isomorphic curve and the value of the corresponding reduced excess entropy is estimated. Capitalizing on the isomorphic nature of the glass transition line, two structural vitrification indicators are identified that allow a rough estimate of the glass transition point only through simple curve metrics of the static properties of supercooled liquids. The vitrification indicators are demonstrated to be quasi-universal by an investigation of hard sphere and inverse power law supercooled liquids. The straightforward extension of the present results to bi-Yukawa systems is also discussed.


Introduction
When liquids are quenched below their melting point by cooling or compression in a manner that suppresses crystallization [1], they exhibit a dramatic slowdown in dynamics and remarkable increase in their viscosity. Since quenching is typically caused by cooling, these metastable liquids are known as supercooled and, for a sufficiently low temperature, they can undergo dynamical arrest and transform into a glass [2]. The process of liquidglass transition or more simply glass transition has been the source of various questions concerning the nature of the transition and the microscopic mechanisms driving it [3][4][5].
The physics of the glass transition have been addressed with a mix of experiments, computer simulations and theoretical approaches. In experiments, the glass transition has been investigated in colloidal systems [6][7][8], granular media [9,10] and organic compounds [11]. Simulations have led to important insight in the physics of supercooled liquids in regimes often not accessible in experiments by adopting simplified models such as the Kob-Andersen [12][13][14] and hard-sphere binary mixtures [15,16]. Theoretical approaches such as mode coupling theory [17,18], random first-order transition theory [19] and dynamic facilitation theory [20] have rationalized some experimental findings and even predicted previously unobserved features of the vitrification process [21,22].
In this paper, mode coupling theory (MCT) is employed to estimate the glass transition line of Yukawa one-component plasmas (YOCP). MCT is an entirely first-principle approach for the investigation of the dynamic processes occurring in glass-forming liquids, which allows to localize the glass transition point requiring only the knowledge of the static properties of the supercooled liquid as an input. MCT is known to lead to good predictions for the dynamics of supercooled liquids [22], which also applies for the glass transition point in spite of the fact that the MCT bifurcation (associated with the glass transition) should rather be interpreted as a cross-over point from a non-activated into an activated dynamics regime [4,22]. The YOCP comprises of equally charged point particles immersed in a neutralizing background that interact via the pair potential u(r) = (Q 2 /r) exp(−r/λ). Here Q is the particle charge and λ the screening length defined by the polarizable background. Thermodynamic YOCP states are uniquely specified by two dimensionless variables [23]: the coupling parameter Γ = βQ 2 /d and the screening parameter κ = λ/d. This allows to re-write the interaction potential as βu(x) = (Γ/x) exp(−κx), where d = (4πn/3) 1/3 is the Wigner-Seitz radius, n is the particle number density, β = 1/(k B T) and x = r/d is a normalized distance. The YOCP possesses a well-understood phase diagram in terms of the κ and Γ variables [24].
The main motivation of this work lies in the relevance of the YOCP model to the experimental realization of complex plasmas, a novel state of soft matter composed of charged particles of mesoscopic size that are immersed in a weakly ionized plasma [23]. The YOCP has been suggested as a promising tool to investigate the dynamics of glassy systems [25], but glass formation has remained experimentally elusive in three dimensional complex plasmas. In particular, complex plasmas have already been employed in order to study supercooled fluids near the vitrification point in two dimensions [26,27], but threedimensional glassy structures exhibiting dynamical arrest still remain to be observed. The accurate estimate of the glass transition line provided in this work should help to guide current [28] or future complex plasma experiments in microgravity conditions that are or will be actively searching for the glassy state of plasmas. It is worth pointing out that MCT calculations of the YOCP glass transition line are already available in the literature [29]. However, there is space for drastic improvement over the existing prediction due to the use of oversimplified structural input that should be grossly inaccurate within the supercooled liquid regime.

Theoretical Background
This section provides an overview of the theoretical background upon which the remaining part of this work is constructed. The equations and approximations which characterize mode coupling theory are presented, the basics of isomorph theory are discussed and the integral equation theory of liquids employed to compute the static structural properties is presented.

Mode Coupling Theory of the Glass Transition
In order to distinguish glasses from stable and supercooled fluids, it is necessary to consider the temporal evolution of the microscopic dynamics. One common probe for such dynamics is the intermediate scattering function, F(k, t), which quantifies, at time t and for the wavenumber k, the correlation of the density fluctuations over a length scale ∼k −1 . The temporal dependence of F(k, t) exhibits three distinguished behaviors between the fluid, supercooled and glassy state [30]: for stable fluids, the intermediate scattering function relaxes exponentially in time, F(k, t) ∼ e −t/τ . On the other hand, supercooled liquids exhibit a multi-stage relaxation in which an initial exponential decay is followed by the so-called β-relaxation regime in which the particles are trapped in cages and the intermediate scattering function exhibits a plateau, i.e., F(k, t) ≈ const.r Such plateau is eventually destroyed during the final α-relaxation regime when the particles escape from their cages, ergodicity is restored and the intermediate scattering function relaxes towards zero following a stretched-exponential law, F(k, t) ∼ e −(t/τ)γ . In the above τ and γ are two constants which are both wavenumber and temperature dependent. Finally, for glasses, there is no full relaxation of the dynamics, particles escape their cages only in rare events and the intermediate scattering function is characterized by a persistent plateau which leads to a positive asymptotic limit F(k, t → ∞) > 0.
MCT provides an equation of motion for the intermediate scattering function [17,18]. In what follows, we will briefly describe the derivation of the fundamental MCT equation with the main purpose of highlighting the assumptions that are adopted in MCT, the reader is addressed to Ref. [30] for a detailed derivation of the MCT equation of motion. Let us start by considering a system with N particles of mass m that are enclosed in a volume V with average number density n = N/V and temperature T. Using r j (t) to denote the position of a particle j at a time t, we introduce the time-dependent microscopic density ρ(r, t) = ∑ N j=1 δ[r − r j (t)] and the microscopic density fluctuations δρ(r, t) = n − ρ(r, t). The intermediate scattering function is then defined as the autocorrelation function of the density fluctuations, F(k, t) = (1/N) δρ(−k, 0)δρ(k, t) , where ... denotes the ensemble average operator and where δρ( is the spatial Fourier transform of the density fluctuations (we have implicitly assumed that the system is isotropic). Employing the Zwanzig-Mori projection formalism [31,32], it is possible to obtain the following exact integro-differential equation for the intermediate scattering function [18] is the static structure factor and M(k, t) is a memory kernel which can be expressed as the ensemble average of a fluctuating force that depends on density fluctuations and pair density products of the form δρ(−k, t)δρ(k, t) [30].
Equation (1) has the structure of an oscillator equation in which Ω(k) plays the role of the characteristic frequency, while the memory kernel acts as a generalized and time-dependent friction coefficient. However, without a simplified expression for the memory function, Equation (1) cannot be solved. Within MCT, such simplified expression for M(k, t) is derived by adopting the procedure explained below. In MCT, the memory kernel is decomposed into the regular M reg (k, t) part that describes the short-time conventional liquid dynamics and the asymptotic M MCT (k, t) part that describes the long-time dynamics dominated by the interplay between caging and ergodicity-restoring effects [18], M(k, t) = M reg (k, t) + M MCT (k, t). The regular part is neglected [30] or approximated with M reg (k, t) = νδ(t) with ν a friction constant [17,33], since the glassy state is mainly concerned with the long-time behavior of the density correlation function. The second part is treated under the assumption that the most relevant contribution stems from the pair density products of the fluctuating force. Hence, M MCT (k, t) is projected onto a basis of pair density products ρ(k 1 , t)ρ(k 2 , t) with a properly defined projection operator running over all (k 1 , k 2 ) pairs relevant to the system [30]. The projection leads to the emergence of triplet correlation functions containing only static properties and of a time-dependent four particle density correlation function; the triplet correlation functions are factorized into static structure factor products within the convolution approximation [34], while the four particle density correlation function is approximated as the product of density pair correlation functions with Kawasaki's approach [35]. These lead to the MCT intermediate scattering function equation where τ(κ) = ν/[Ω(k)] 2 and the memory kernel is given by [30] In the above p = |p|, p = k − k and c(k) is the Fourier transform of the direct correlation function which is related to the static structure factor via the Ornstein-Zernike equation (see Section 2.3), S(k) = 1/[1 − nc(k)]. Provided that the static structure factor and the viscosity are known, Equation (2) is a self-consistent equation for F(k, t) which can be solved subject to the initial conditions F(k, 0) = S(k) and ∂F(k, 0)/∂t = 0. Equation (2) is sometimes reported in its over-damped form by assuming that the viscosity is so large that the inertial second order derivative term can be neglected [18,33,36]. However, while the over-damped form of Equation (2) is appropriate to study colloids, it would be inaccurate for complex plasmas in which the dilute background gas [23] makes viscous damping less predominant.
Given that glasses can be distinguished from conventional liquids from the asymptotic limit of the intermediate scattering function, it is sufficient to obtain an equation for F(k, t → ∞) in order to investigate the glass transition properties. Rewriting Equation (2) in terms of the normalized density autocorrelation function φ(k, t) = F(k, t)/S(k), Laplace transforming and employing the final value theorem to extract the asymptotic limit, leads to the form factor f (k) = lim t→∞ φ(k, t) MCT equation The solution of the MCT equation for the form factor allows to distinguish between liquids and glasses, since the former are characterized by f (k) = 0 while the latter are characterized by f (k) > 0. It should be noted that, at the glass transition point, the form factor changes discontinuously from zero to some positive critical value f c (k) > 0. This discontinuity in the form factor happens despite the fact that the static properties do not exhibit any discontinuity between supercooled liquids and glasses. This bifurcation phenomenon is a manifestation of the feedback between the force fluctuations in the memory kernel and the density fluctuations in the form factor [18,21].

Isomorph Theory
Isomorphic curves are phase diagram lines of constant excess entropy along which a substantial set of structural and dynamical properties remain approximately invariant when expressed in dimensionless units where the length is normalized to a = n 1/3 and the energy to k B T [37,38]. While it is possible to identify lines of constant excess entropy in the phase diagram of any system, only in R-simple systems the isentropic lines are also isomorphs. R-simple systems are characterized by the property that the ordering of the potential energies of two configurations corresponding to the same density is preserved when the two configurations are re-scaled uniformly to a different density [39]. In other words, denoting with U the potential energy and with R the set of N particle's position {r 1 , ...r N }, R-simple systems satisfy the relation U(R a ) < U(R b ) =⇒ U(ζR a ) < U(ζR b ), where R a and R b are two equal-density configurations and ζ a positive scaling factor.
A recent investigation has shown that the YOCP is an R-simple system whose isomorphs can be accurately described with the following analytical parameterization [40] where α = a/d = (4π/3) 1/3 is the ratio between the mean-cubic inter-particle distance employed in isomorph theory and the Wigner-Seitz radius used for the characterization of the YOCP state point, while Λ is a constant close to unity which depends weakly on the state point. It should be noted that, with the assumption Λ = 1, the isomorph parameterization described by Equation (5) is equivalent to the semi-empirical expression utilized in Ref. [41] in order to fit MD data of the YOCP melting line [24,42] Γ m (κ) = Γ OCP m e ακ 1 + ακ + (ακ) 2 where Γ OCP m = 171.8 is the coupling parameter at melting in the one-component plasma (OCP) limit (κ = 0) [42]. The fact that the isomorphs and the melting line can be approximated with the same analytical expression is a reflection of the fact that for R-simple systems the melting line constitutes an isomorphic curve to the first-order [37,43].

Integral Equation Theory of Liquids
The integral equation theory (IET) of liquids gives access to the static liquid properties without resorting to computer simulations. For one-component liquids with pair-wise isotropic interactions, IET comprises of two exact equations: the Ornstein-Zernike equation [34] h(r) = c(r) + n c(r )h(|r − r |)d 3 r and the non-linear closure condition derived from cluster diagram analysis [34] In the above g(r) denotes the radial distribution function, h(r) = g(r) − 1 the total correlation function, c(r) the direct correlation function and B(r) the bridge function. Knowledge of the total correlation function h(r) gives direct access to the static structure factor via the Fourier space relation S(k) = 1 + nh(k). Thus, the solution of the above equations provides all the input that is necessary for the calculation of the form factor via MCT. Nevertheless, the solution of the IET system of equations requires an expression for the bridge function. The latter can be formally represented as a power series of the density, but such series is known to converge slowly already at moderate densities [44,45]. Combined with the fact that the calculation of the series coefficients becomes extremely cumbersome beyond the third order [46][47][48], this calls for the adoption of approximate bridge function expressions.
Among the wide variety of bridge function approximations or closures that have been proposed over the years [49], here we shall consider only three approximations which have been previously employed to investigate the fluid properties of the YOCP: the hypernetted chain (HNC) approximation, the isomorph-based empirically-modified hypernetted chain (IEMHNC) approach and the variational modified hypernetted chain (VMHNC) approximation. Within the HNC approximation, the bridge function is neglected altogether by assuming B(r) = 0 [50]. The HNC approximation is straightforward to implement and flexible, since it is not taylored to any specific interaction potential. For systems featuring Coulomb interactions, it has been shown that the HNC approximation produces qualitatively correct results for the static and thermodynamic properties of state points far from crystallization, but that its performance rapidly degrades near the freezing line [51].
The IEMHNC approach is an advanced closure to the IET system of equations that is built upon the ansatz that reduced-unit bridge functions remain exactly invariant along isomorphic lines [52]. Systematic indirect bridge function extractions along isomorphic curves with the aid of molecular dynamics simulations have confirmed the approximate validity of the invariance conjecture for Yukawa systems [53]. This ansatz is then combined with two external inputs, a closed-form expression for the isomorphs and a closed-form bridge function expression valid along any phase diagram line that possesses a unique intersection point with any isomorph, in order to construct an expression for the bridge function valid for the whole phase diagram [52]. The IEMHNC approach has been applied to Yukawa and bi-Yukawa fluids showing a remarkable agreement with computer simulations in the entire dense fluid region up to crystallization [52,54]. Note that the IEMHNC approach is only applicable to R-simple systems for which the aforementioned external input is available. In addition, the IEMHNC bridge function expression might be accurate only in a sub-region of the phase diagram which depends on the region of validity of the external inputs. For instance, with the currently available external inputs, the IEMHNC bridge function for the YOCP is strictly applicable for coupling parameters which satisfy 5.25 ≤ Γ iso (Γ, κ) ≤ 171.8 [52].
The VMHNC approximation is an advanced closure to the IET system of equations which is based on the ansatz of bridge function quasi-universality [55]. The quasiuniversality conjecture justifies the following two-step procedure adopted to define the bridge function [56]: first, the unknown bridge function B(r) is replaced with the Percus-Yevick bridge function of the hard sphere system, B HSPY (r; η), for which an exact analytical representation in terms of the inter-particle separation r and of the packing fraction η = πnσ 3 /6 (with σ the hard sphere diameter) is available. Second, the value of the packing fraction which ensures the correct mapping between B(r) and B HSPY (r; η) is determined with a robust method based on the minimization of a properly defined free energy functional. Since it is not constructed for a specific class of systems, the VMHNC approach shares the same flexibility of the HNC approximation and can be applied without any major modification to any one-component system characterized by purely repulsive interactions. For the specific case of the YOCP, the VMHNC is known to produce results that compare exceptionally well with computer simulations [57] and are as accurate as those obtained with the IEMHNC approach over the entire fluid region of the phase diagram [58]. The main drawback of the VMHNC approach lies in its computational cost, since the minimization of the free energy functional makes the algorithm for the solution of the IET more cumbersome, eventually causing the VMHNC approach to become up to 80 times slower than approximations which do not involve minimization [58].

Computational Approach
This section describes the computational methods employed in the solution of the MCT equation for the form factor and in the solution of the IET system of equations for the static structural properties. The advantages of combining MCT with structural input obtained with advanced IET closures are discussed and the present algorithm is benchmarked against the literature results.

Combining MCT with Advanced IET Approaches
Within MCT, the static structural properties of the supercooled fluid constitute the only external to the theory input that is required for the calculation of the glass transition line and the critical form factors. In the literature, such properties have been obtained either from computer simulations [59][60][61] or from IET calculations combined with elementary closures (when simulation input was unavailable) which include the aforementioned HNC approximation for soft long range interaction potentials [29,62] or the Percus-Yevick (PY) approximation for hard sphere short range interactions [16,17,33,[63][64][65]. More advanced IET closures which enforce thermodynamic consistency through free parameters or resemble the VMHNC approximation by featuring an optimized correspondence rule have also been considered [66][67][68]. However, they have received comparatively much less attention than the elementary HNC and PY closures. The objective of this section is to demonstrate that the adoption of advanced IET approximations for the calculation of the static properties is an essential ingredient for reliable predictions of the MCT glass transition line. In what follows, the HNC approximation will be compared to the IEMHNC approach. The PY approximation will not be discussed owing to the long range interaction potential of interest, while the VMHNC approach will not be addressed here owing to its similar accuracy to the IEMHNC approach [58].
The discussion will center around stable rather than supercooled fluids. The main reason is that it has proven to be extremely challenging to simulate supercooled liquids due to the necessity of preventing crystallization (especially for one-component systems) [69][70][71][72][73] as well as due to the high computational cost necessary to obtain an equilibrated configuration and to effectively sample the phase space [74][75][76][77], while it is relatively cheap from a computational point of view to perform computer simulations of the stable fluid state in order to quantify the accuracy of the different IET approximations. The secondary reason is that, formally, the region of validity of the IEMHNC approach is the stable fluid region. Thus, any attempt to compare the HNC and IEMHNC approximations outside this region involves extrapolations which inevitably cast a doubt over the result of the comparison. However, as will shall deduce later, the IEMHNC approach can indeed be extrapolated deep into the supercooled liquid regime without losing its accuracy.
The comparison is illustrated in Figure 1. Panels (a) and (b) feature a radial distribution function comparison between the results of the HNC, IEMHNC approximations and the results of molecular dynamics simulations. It is evident that the IEMHNC approach outperforms the HNC approximation leading to radial distribution functions that are nearly indistinguishable from the results of computer simulations within the first and second coordination cells [52]. In addition, the IEMHNC approach maintains the same level of accuracy throughout the stable fluid region [58]. On the contrary, the HNC approximation exhibits strong deviations from the MD results within the first coordination cell. More important, it becomes more and more problematic upon approaching the melting point [52,58]. The latter observation is a direct manifestation of the fact that the bridge function contribution to the static structural properties gradually becomes more prominent as the liquid-solid phase transition is approached. Finally, it is reasonable to assume that this trend continues also beyond the melting point, in the supercooled portion of the phase diagram, which implies that the HNC approximation leads to grossly inaccurate estimates for state points close to the MCT glass transition. Panels (c) and (d) of Figure 1 feature a static structure factor comparison between the results of the HNC and IEMHNC approximations. These panels shed light on another property of the HNC approximation, which is here observed for the stable fluid region but also holds beyond the melting point (as we shall deduce in what follows); namely, the fact that the HNC approximation produces quantitatively correct structural properties but for state points of much stronger coupling than that of the actual state point. In particular, in panel (c), the IEMHNC and HNC approaches are compared for the same state points and, as anticipated, lead to different results. Given the high accuracy of the IEMHNC approach, the IEMHNC results can be considered as nearly exact which confirms that the HNC approximation introduces noticeable distortions in the static properties. On the other hand, panel (d) demonstrates that the two IET approximations result to nearly identical structural properties, if the HNC approximation is applied at Γ = Γ HNC and the IEMHNC approach at Γ = Γ IEMHNC with Γ HNC Γ IEMHNC . In order to rationalize this result, it is convenient to rewrite the IET diagrammatic closure, Equation (8) . This leads to the physical interpretation of the bridge function B(r) as an additional repulsion which is superimposed on the interaction potential u(r), given that it is always negative except for a series of small positive peaks in its long-range behavior [79][80][81][82]. Thus, it becomes apparent that the HNC approximation, which assumes that B(r) = 0, would produce quantitatively accurate results, only when the state point is adjusted in such a way that u * (r) ≈ u(r). For the specific case of the YOCP, this amounts to say that the HNC approximation can be expected to reproduce the structural properties of the state point (Γ, κ), only if it is applied to a different state point (Γ HNC Γ, κ), as demonstrated in panel (d).
At this point, it is worth analyzing the consequences of the results discussed above on the MCT glass transition calculations. The MCT equation for the form factor, Equation (4), does not explicitly include any information about the interaction potential whose effect appears only indirectly via the static properties in the non-linear kernel stemming from the memory function. Therefore, it can be expected that the HNC approximation, which produces nearly exact static properties if the state point is artificially re-scaled towards the stronger coupling region, produces unreliable estimates for the location of the MCT glass transition line but accurate predictions for the shape of the critical MCT form factors. Therefore, for a precise determination of the MCT glass transition line, it is necessary to adopt more advanced IET closures. Here we consider two such closures, the IEMHNC approach and the VMHNC approximation. The reason behind considering both advanced closures is that none of them has been extensively tested in the supercooled regime, so it is important to confirm that they produce consistent results for a phase diagram region that falls outside their normal range of applicability.

Numerical Implementation
The following four-step procedure was devised for the efficient localization of the YOCP MCT glass transition line: (i) For a given value of the screening parameter κ, ten values of the coupling parameter Γ that belong to the YOCP supercooled region were considered according to the prescription Γ i = Γ i−1 + ∆Γ with i = 2, 3, ...10. Here Γ 1 and ∆Γ are κ-dependent free parameters that should be chosen in a manner that allows the exploration of a sizable portion of the YOCP supercooled region. (ii) For each (κ, Γ i ) combination, Equations (7) and (8) were solved with one of the three closures described in the Section 2.3 leading to the determination of the static structural properties of the supercooled fluid. (iii) For each (κ, Γ i ) combination, the form factor was computed from Equation (4) supplemented with the structural properties obtained from the IET equations. (iv) In the case of positive form factor for a state point characterized by Γ = Γ i , the update Γ 1 → Γ i−1 was employed, ∆Γ was divided by ten and the procedure was repeated until the coupling parameter of the glass transition point was determined within four digits of accuracy. In the case of zero form factors for all Γ = Γ i state points, the update Γ 1 → Γ 10 + ∆Γ was employed and the procedure was restarted. In the case of positive form factors for all Γ = Γ i state points, the update Γ 1 → Γ 1 − 11∆Γ was employed and the procedure was restarted. For the above procedure to prove successful, robust algorithms should be available for the solution of the IET system of Equations (7) and (8) and for the solution of the MCT form factor Equation (4). The algorithms implemented are outlined in the following paragraphs.
The IET equations were solved with a well-established algorithm [52,54,58] that is based on Picard iterations in Fourier space combined (when necessary) with mixing and long-range decomposition techniques to facilitate convergence. The Fourier transforms are computed on a discretized domain extending up to R max = 80d with a real space resolution ∆r = 10 −3 d and a reciprocal space resolution ∆k IET = π/R max = 0.039/d. The convergence criterion for the Picard iterations reads as |γ m (k) − γ m−1 (k)| < 10 −5 ∀k where γ(k) is the Fourier transform of the indirect correlation function γ(r) = h(r) − c(r). When the IET equations are solved within the VMHNC closure, the effective packing fraction η which appears in the VMHNC bridge function is found with a dedicated iteration cycle that is terminated when the convergence criterion |η m − η m−1 |/η m < 10 −5 is satisfied. A more detailed description of the algorithm employed in the solution of the IET equations can be found in Ref. [58].
The MCT equation for the form factor was solved over a discretized wavenumber domain with resolution ∆k = 0.1/d containing 400 equally-spaced points distributed between ∆k and k m = 40/d . Taking advantage of the fact that the long-time limit of the form-factor obeys the maximum property [33], Equation (4) is a short-hand notation for the non-linear term which appears in the right hand-side of Equation (4). Concerning this non-linear term, it was first rewritten in a form more suitable for numerical integration by applying the bipolar convolution theorem [12,17] (9) and then it was evaluated with the adaptive Gauss-Kronrod quadrature rule as implemented in the GNU Scientific Library [83]. The discretization employed in the solution of the MCT equation has two straightforward consequences. First, the fact that ∆k > ∆k IET suggests that it is necessary to interpolate the static properties obtained from IET when they are employed in the solution of the MCT equations. This, however, should not impact the final result of the MCT calculations because there is no extrapolation involved and because such results have proven to be independent from the grid resolution if ∆k ≤ 0.2/d (see Section 3.3). Second, the integrand values within the long wavelength limit (k = 0) and within the short wavelength limit (k → ∞) are inaccessible, since both these limits fall outside of the grid employed in the solution of the MCT equations. This implies that in Equation (9)

Benchmarking and Convergence Study
Before proceeding with the presentation of the results stemming from the systematic solution of the MCT equations over the supercooled YOCP phase diagram, it is important to discuss how the present algorithm was benchmarked against results available in the literature.
The MCT benchmarking exercise initially focused on hard-sphere (HS) systems, since the glass transition of HS systems has been the focal point of numerous investigations over the years [15][16][17]33] making it relatively easy to find reference values for the critical form factor and the glass transition point. In addition, it is possible to employ analytical expressions for the static structural properties of the HS systems, which allows to test the isolated performance of the MCT solver without worrying about the performance of the IET solver. Regarding the static properties of such systems, they have either been obtained within the Percus-Yevick approximation by employing the Wertheim-Thiele analytical solution [84,85] or expressed through the accurate semi-empirical parameterization proposed by Verlet and Weis [86]. In what follows, the former case will be coined as the Hard-Sphere Percus-Yevick (HSPY) system (this is the same system employed for the definition of the bridge function in the VMHNC approximation), while the latter case will be addressed as the Hard-Sphere Verlet-Weis (HSVW) system. Panel (a) in Figure 2 illustrates how the present numerical solutions for the critical form factor of the HSPY and HSVW systems compare against results available in the literature. It is evident that the present numerical implementation is able to correctly reproduce both the glass transition points (which occur at η = 0.516 for the HSPY system and η = 0.525 for the HSVW system) and the corresponding critical form factors.  Figure 2 of Ref. [17] (discrete circles). The inset compares the form factor at the glass transition point for the HSVW system as obtained from the present numerical implementation (solid line) and from the results reported by Wu and Cao in Figure 1 of Ref. [63] (discrete squares). The HS calculations were performed with k m = 50/σ and ∆k = 0.1/σ where σ denotes the HS diameter. Panel (b) compares the form factor at the glass transition point for the OCP system as obtained from the present numerical implementation (solid line) and from the results reported by Yazdi and collaborators in Figure 2 of Ref. [29] (discrete diamonds). Note that for the present implementation we report the form factors at two values of the coupling parameter: Γ = 575 (magenta) which is the present glass transition point prediction and Γ = 590 (blue) which is the Yazdi glass transition point prediction.
The MCT benchmarking exercise was also extended to YOCP systems, where the investigation of Yazdi and collaborators [29] provided the only literature results that can be employed as a reference. In panel (b) of Figure 2, the present solution of the MCT form factor equation in the OCP limit is compared with the corresponding OCP solution (as reported in Figure 2 of Ref. [29]). In both cases, the static properties of the OCP were computed by closing the IET system of equations with the HNC approximation. Two problems can immediately be noticed: the predictions for the glass transition point do not match (the present implementation predicts a glass transition at Γ = 575, while the reference implementation predicts a glass transition at Γ = 590) and the two implementations generate different form factors for the same coupling parameter. Unfortunately, Ref. [29] contains very little information regarding the algorithms and the parameters used to solve the MCT form factor equation. Hence, in order to rule out the presence of numerical errors in the present solution, we had to resort to a systematic convergence study on the two free numerical parameters that emerge in the solution of Equation (4) with the memory kernel expressed via Equation (9), namely the wavenumber resolution ∆k and the cutoff wavenumber k m .
In the convergence study, the wavenumber resolution was varied between ∆k = 0.05/d and ∆k = 0.4/d, while the cutoff wavenumber was varied between k m = 10/d and k m = 50/d. The results of this investigation are reported in Figure 3 and reveal that a wavenumber resolution ∆k ≤ 0.2/d and a cutoff wavenumber k m ≥ 20/d are sufficient to obtain a form factor and glass transition point Γ g which are independent from the parameters employed in the numerical implementation. The results in Figure 3 refer to the OCP limit (κ = 0) and to structural input from the HNC approximation, but similar conclusions were obtained also for selected higher values of κ and structural input from the IEMHNC or VMHNC approximations. To summarize, this section demonstrated that the present numerical implementation of the MCT equations is able to reproduce results available in the literature for HS systems. Regarding YOCP systems, a small mismatch was observed between the present results and the literature results. However, while it was not possible to pin-point the origin of this mismatch due to insufficient information, the results of a thorough convergence study ensured that the present MCT calculations are accurate also for Yukawa systems.

The MCT Glass Transition Line
The glass transition line of the YOCP was computed for fifteen screening parameters which belong to the part of the YOCP phase diagram that is most relevant to experimental realizations of Yukawa systems and for which molecular dynamics based calculations of the liquid-solid phase transition are available [24,42]   The observation that suggests that the IEMHNC approximation and the VMHNC approach maintain the consistency and accuracy which characterizes them for stable YOCP liquids also for supercooled YOCP liquids. In addition, it indirectly suggests that the defacto extrapolations of the IEMHNC and VMHNC bridge functions deep into the supercooled regime are rather safe. As a consequence, both the Γ MCT-V g (κ) and Γ MCT-I g (κ) estimates can be considered to be an accurate approximation to the "exact" MCT glass transition line which would have been obtained if the MCT equations were supplied with an "exact" structural input extracted from computer simulations. On the other hand, the observation that Γ MCT−H g (κ) ≈ 2Γ MCT-I g (κ) confirms that adopting the HNC approximation would lead to inaccurate estimates for the MCT glass transition line, see the detailed discussion in Section 3.1. In fact, assuming that the MCT-I or MCT-V calculations are able to predict the "exact" MCT glass transition line within a few percent, then the MCT-H calculations lead to an overestimation of the MCT glass transition line of roughly a factor two, regardless of the value of the screening parameter. This justifies our revisiting of the YOCP glass transition results earlier reported by Yazdi and collaborators [29], since the latter were exclusively based on MCT-H calculations.
Perhaps, the most significant result of Ref. [29] was the observation that, for YOCP systems, the glass transition line is almost parallel to the melting line. As a consequence, it can be well approximated by the following analytical semi-empirical formula [29] Γ g (κ) = Γ OCP g e ακ 1 + ακ + (ακ) 2 where Γ OCP g is the glass transition point in the OCP limit. It is important to test the validity of the scaling given in Equation (10) against the accurate glass transition predictions obtained from the MCT-I and the MCT-V calculations not only because it provides a simple convenient parameterization but also because its validity is a strong indication that the glass transition line constitutes an isomorphic curve, compare Equations (10) and (5). As we shall see in the following sections, the isomorph invariance of the MCT glass transition line implies the validity of numerous state-independent symmetries that can be exploited in different ways.
In Figure 4, the analytical scaling of Equation (10) is compared against the predictions of the three MCT glass transition calculations. At this point, it should be noted that during the comparison with a type of MCT calculation, the Γ OCP g pre-factor which appears in Equation (10) is adapted to match to the corresponding value that is given in the first line of Table 1. This means that Γ OCP g = 575.5 for MCT-H, Γ OCP g = 289.8 for MCT-I and Γ OCP g = 279.7 for MCT-V. It is apparent that the semi-empirical expression performs well for the MCT-H calculations (as it was already observed in Ref. [29]) and that it performs even better for the MCT-I calculations and the MCT-V calculations. To be more quantitative, for the set of screening parameters considered in this work (κ ≤ 5.0), the relative deviations between the analytical scaling and Γ MCT-H g (κ) can become as large as 15%, while the relative deviations between the analytical scaling and Γ MCT-I g (κ) or Γ MCT-V g (κ) never exceed 5%. The high accuracy of analytical semi-empirical expression for the MCT glass transition and the equivalence of Equation (10) to Equation (5) with Λ = 1, suggest that the glass transition line can be considered to be an isomorphic curve. This observation also rationalizes why the analytical scaling agrees better with the MCT-I calculations which feature an explicitly isomorph-invariant bridge function [52] and with the MCT-V calculations which feature an implicitly isomorph-invariant bridge function [58] rather than with the MCT-H calculations which neglect the bridge function altogether. It should be mentioned that a modified version of Equation (10) was also tested, in which the analytical scaling was made fully equivalent to the expression for the YOCP isomorphs by replacing ακ with Λακ and then by determining Λ via least-square fitting. The results of such modified scaling will not be analyzed further, because for all three calculations the parameter Λ was found to be very close to unity, suggesting that the choice Λ = 1 is already nearly optimal, and that the modified scaling possesses essentially the same accuracy of the standard scaling of Equation (10).
The above results indicate that the MCT glass transition line of YOCP systems constitutes an isomorphic curve. Given that isomorphs are phase diagram lines of constant reduced excess entropy, this suggestion can be supported by computing the reduced excess entropies of the YOCP state points that lie along the glass transition line, s ex (κ, Γ g ), and then by checking to which extent they satisfy the condition s ex (κ, Γ g ) = const. The accurate computation of the excess entropy requires to perform computer simulations of supercooled fluids combined with free-energy calculation techniques [87], but the correct implementation and execution of such simulations falls outside the scope of the present work. Hence, here we computed s ex (κ, Γ g ) by utilizing two YOCP equations of state available in the literature which provide the reduced excess internal energy, u ex (κ, Γ), from which the reduced excess entropy can be straightforwardly computed by thermodynamic integration from s ex (κ, Γ) = u ex (κ, Γ) −  [24,42]. The equation of state proposed by Rosenfeld and Tarazona, within the framework of an asymptotically-high density expansion for purely repulsive potentials, fits well computer simulation results in the strongly coupled stable liquid portion of the YOCP phase diagram and leads to the reduced excess internal energy expression u RT ex (κ, Γ) = M(κ)Γ + 3.0[Γ/Γ m (κ)] 2/5 [88,89]. The κ-dependent parameters a(κ), b(κ), c(κ), d(κ) and M(κ) are given in the respective references, while the coupling parameter at the melting point, Γ m (κ), can be conveniently expressed via Equation (6) [54]. It should be emphasized that both equations of state are strictly valid only for the stable fluid phase of the YOCP. In order to confirm that they can be safely extrapolated to the supercooled regime, the reduced excess internal energies u H ex (κ, Γ) and u RT ex (κ, Γ) were compared against those computed from the IET structural properties, u ex (κ, Γ) = (3/2) ∞ 0 x 2 βu(x)g(x)dx, along the glass transition line that stems from the MCT-I and MCT-V calculations. The comparison showed that both equations of state could predict the reduced excess internal energy within 2% along the MCT-I glass transition line and within 1% along the MCT-V glass transition line leading to the conclusion that both can be employed to obtain accurate estimates for thermodynamic properties of supercooled Yukawa fluids.
The reduced excess entropies along the MCT-I and MCT-V glass transition lines are reported in Table 2. The MCT-H calculations are not included due to their poor glass transition line prediction. There are deviations in the entropy predictions of the equations of state. The Rosenfeld-Tarazona expression predicts that both glass transition lines satisfy s ex (κ, Γ g ) = −5.5 within 3%, whereas the Hamaguchi expression predicts reduced excess entropies that are not constant along the glass transition lines (see the relatively large deviations reported in the last two rows of Table 2) and change noticeably between the MCT-I and MCT-V glass transition lines. At this point, one should recall that the Hamaguchi equation of state was constructed on an empirical basis by observing the functional dependencies within the stable fluid range, while the Rosenfeld-Tarazona equation of state was derived on the basis of an asymptotically-high coupling parameter expansion that ignores the liquid-solid phase transition. Thus, the Rosenfeld-Tarazona equation of state should be more accurate for supercooled YOCP liquids. Overall, from the results of Table 2, it can be concluded that the MCT glass transition line constitutes an isentropic curve with s ex = −5.5 (within 3%) reduced excess entropy. Table 2. Excess entropy along at the MCT-I glass transition line (columns 3 and 4) and the MCT-V glass transition line (columns 6 and 7) as obtained from the Hamaguchi equation of state s H ex , and the Rosenfeld-Tarazona equation of state, s RT ex . All the excess entropies are reported in units of Nk b , i.e., these are reduced excess entropies. In the last three rows, the designation AVE denotes the average value of the reduced excess entropy along the respective glass transition line (rows 1-15), while the quantities a and m report the average and maximum relative deviations between the reduced excess entropies along the respective glass transition line and their average value reported in AVE.

The MCT form Factors
The critical form factors resulting from the MCT-H, MCT-I and MCT-V calculations are compared in Figure 5 for two screening parameters. Important observations from Figure 5 are that the MCT-H, MCT-I and MCT-V calculations lead to similar form factors at their respective glass transition points (compare the solid curves) as well as that there are enormous form factor differences between MCT-H and MCT-I or MCT-V calculations at the same state point (compare the dashed curves to the blue solid curve). These observations are in line with the discussions of Sections 2.3 and 3.1. They are both manifestations of the fact that the MCT non-linear kernel on the right hand side of Equation (4) requires only the static structure factor as external input and does not contain any explicit dependence on the interaction potential. The form-factor similarity at the glass transition point can be explained by the fact that the IEMHNC, VMHNC and HNC approximations can all be expected to produce similar structure factors at their respective glass transition points; the IEMHNC and VMHNC approximations due to their high accuracy (see Section 2.3) and the HNC for the reasons discussed in Section 3.1. On the other hand, the large discrepancies observed between the MCT-H calculations and the other two sets of calculations at the same state point are caused by the poor performance of the HNC approximation in the dense fluid region which was also discussed in Section 3.1. Figure 6 illustrates how the critical YOCP form factor changes with the screening parameter according to MCT-I and MCT-V calculations; similar results (not shown here) were also obtained from MCT-H calculations. Regardless of the IET approximation employed to compute the static structural properties, it is apparent that the critical form factor is characterized by a main peak whose magnitude and position are practically independent of the screening parameter under consideration. The main peak is preceded by a rapid decay towards a state-point-dependent long wavelength limit and is followed by slowly decaying oscillations which also appear to be state-point-independent. The results of Figure 6 are rather anticipated in light of our previous discussions. Since the MCT glass transition line constitutes an isomorphic curve, then a large number of reduced unit dynamic and structural properties should be nearly invariant while traversing the glass transition line. The static structure factor is known to be an isomorph invariant quantity, thus, given the deep connection between S(k) and f (k), it can be expected that also the form factor exhibits some degree of invariance. On the other hand, the variance of f (k) near the k = 0 limit is caused by the fact that the long wavelength limit of the MCT memory kernel is given by a relation of the form αS(0) = αχ T where the pre-factor α is described by Equation (9b) of Ref. [29] and where χ T = 1 + n h(r)d 3 r is the isothermal compressibility, which is known to strongly vary along an isomorph [90] (strictly speaking, S(0) = χ T does not hold for the OCP for which the correct infinite wavelength limit of the structure factor, S(0) = 0, must be retrieved from the small argument expansion S(q → 0) = (3Γ/q 2 + 1/χ T ) −1 with q = kd [91]). While the isomorph variance of S(0) has a negligible effect on the overall degree of invariance of the static structure factor due to the nearly incompressible nature of supercooled YOCP liquids (see Figure 7), the isomorph variance of S(0) has a strong effect on the degree of invariance of the form factor over a sizeable wavenumber interval due to the large value of the proportionality constant α which augments the small compressibility values.

The MCT Vitrification Indicators
The YOCP static structure factor, radial distribution function and direct correlation function along the MCT glass transition line (MCT-H, MCT-I, MCT-V) are illustrated in Figure 7. Let us first focus on the results for the state points pertaining to the MCT-I and MCT-V glass transition lines. In accordance with the predictions of isomorph theory, the static structure factor and radial distribution function are almost invariant. On the other hand, owing to the asymptotic limit c(r → ∞) = −βu(r) and due to the connection to the compressibility via (1/χ T ) = −n c(r)d 3 r, the direct correlation function is strongly variant. Albeit not directly evident from Figure 7, the radial distribution function and the static structure factor are similar to f (k) in the sense that their short-range values are also non-invariant; g(r) is variant due to the asymptotic limit g(r → 0) ∝ exp[−βu(r)] [92], while S(k) is variant due to its connection to the compressibility S(k → 0) = χ T . However, while the short-range variance has a noticeable effect on the critical form factor, it is inconsequential for both the radial distribution function and the static structure factor since it occurs in a region where both functions are approximately zero. Finally, we point out that the invariance of the radial distribution and the static structure factor is slightly violated in the vicinity of their first maxima, but it holds to a nearly exact degree at their subsequent minima and maxima. Regarding the static properties along the MCT-H glass transition line, it can be noticed that S(k) and g(r) both show a reduced degree of invariance, in accordance to an earlier observation which reported that assuming B(r) = 0 has a detrimental effect on the invariant properties along an isomorphic line [52]. The high degree of isomorph invariance which characterizes the static structure factor and the radial distribution function at the MCT glass transition line opens up the possibility of determining a group of empirical vitrification indicators. Such indicators are inspired by the successful application of freezing indicators which allow to locate the liquid-solid phase transition line simply by monitoring some characteristic features of the liquid state S(k) and g(r). Empirical vitrification indicators would serve as phenomenological criteria for the localization of the glass transition line based solely on the structural properties of the supercooled liquid. In fact, it was observed that two commonly used freezing indicators, namely the magnitude of the first peak of the static structure factor S max employed in the Hansen-Verlet freezing rule [93] and the amplitude ratio of the first nonzero minimum to the first maximum of the radial distribution function g R employed in the Raveché-Mountain-Streett freezing rule [94] perform reasonably well also for the prediction of the glass transition line, albeit at different values. In particular, it was revealed that the MCT-I glass transition line is characterized by S max = 4.49 ± 0.15 and by g R = 0.14 ± 0.01, while the MCT-V glass transition line is characterized by S max = 4.00 ± 0.08 and by g R = 0.13 ± 0.01.

MCT-H glass transition points
Given the fact that most of the non-invariant features of S(k) and g(r) are concentrated around their first peak and taking into account that an effective vitrification indicator should remain as constant as possible along the glass transition line, it should be possible to identify even better vitrification indicators by generating simple curve metrics that do not involve the magnitude of the first peak of these quantities. For this reason, two additional prospective vitrification indicators were considered which refer to the amplitude ratio of the first nonzero minimum to the second maximum of the radial distribution function, g R2 , and the amplitude ratio of the first nonzero minimum to the second maximum of the static structure factor, S R2 . The values of these prospective vitrification indicators along the glass transition lines stemming from the MCT-I and the MCT-V calculations are reported in Table 3. Table 3. Prospective vitrification indicators along the MCT-I glass transition line (columns 3 and 4) and the MCT-V glass transition line (columns 6 and 7). S APP R2 denotes the amplitude ratio of the first nonzero minimum to the second maximum of the static structure factor obtained from approximation APP, while g APP R2 denotes the amplitude ratio of the first nonzero minimum to the second maximum of the radial distribution function stemming from approximation APP. In the last three rows, the designation AVE denotes the average value of the vitrification indicator along the glass transition line (rows 1-15), while the quantities a and m report the average and maximum relative deviations between the vitrification indicators along the respective glass transition line and their average value reported in AVE. Both prospective vitrification indicators remain almost constant along the MCT glass transition line, with minor deviations from their average value which do not exceed 1% for g R2 and 3% for S R2 . These vitrification indicators possess another desirable property, since they exhibit respectable variations within the supercooled regime prior to and post the glass transition line (see Figure 8). This characteristic highlights their potential practical use for the localization of the MCT glass transition point. Overall, the indicators g R2 and S R2 perform better than S max and g R , since they exhibit slightly smaller variations along the MCT glass transition line and result to more consistent predictions between the two IET approaches employed in the computation of the MCT glass transition line. In this regard, the VMHNC approximation and the IEMHNC approach produce really close but not identical values for g R2 and S R2 at the glass transition point. We argue that the results obtained with the VMHNC approximation should be preferred over the ones obtained with the IEMHNC approach, since the former approximation is known to produce more accurate predictions for the second coordination cell [58]. Combining the above, it can be concluded that the MCT glass transition point is characterized by g R2 = 0.30 or equivalently by S R2 = 0.35 and that these two conditions can be employed to obtain an accurate guess for the MCT glass transition line of the YOCP without having to solve the MCT equation. On a side note, it is worth pointing out that the vitrification indicator g R2 could also be employed as a freezing indicator, since the state points along the YOCP melting line obtained via computer simulations are all characterized by g R2 = 0.4 ± 0.01 (see panel b in Figure 8).    for different screening parameters as a function of the coupling parameter normalized by its MCT glass transition value. Results for stable liquids (Γ/Γ g (κ) 0.6), supercooled liquids prior to the glass transition (0.6 Γ/Γ g (κ) < 1.0) and supercooled liquids post the glass transition (Γ/Γ g (κ) > 1.0). For each indicator, the superscript denotes the IET approximation employed, i.e., the VMHNC approach or IEMHNC approach. The panels report the S R2 or g R2 values for six screening parameters κ = {0, 1, 2, 3, 4, 5} and two MCT glass transition lines, namely those stemming from MCT-V calculations where Γ g (κ) ≡ Γ MCT-V g (κ) (main plot) and from MCT-I calculations where Γ g (κ) ≡ Γ MCT-I g (κ) (inset). The numerical values of Γ MCT−V g (κ) and Γ MCT-I g (κ) have been reported in Table 1. The full symbols represent the values of S R2 and g R2 at the YOCP melting point predicted by computer simulations [24].

Summary of the Results
In the present work, the glass transition line of Yukawa one-component plasmas was computed by combining the mode coupling theory of the glass transition with highly accurate structural input obtained from two advanced closures to the integral equation theory of liquids, namely the isomorph-based empirically modified hypernetted chain approach and the variational modified hypernetted chain approach. It was observed that both closures lead to consistent values for the YOCP glass transition line and it was concluded that the present results offer a greatly improved estimate compared to earlier estimates that are available in the literature. Besides the improvement upon existing results, the highly accurate structural input adopted in the present calculations allowed the identification of two vitrification indicators which can be employed to obtain an accurate guess for the YOCP glass transition line without necessitating a determination of the bifurcation point. The existence of vitrification indicators is important from a theoretical and practical standpoint; from a theoretical perspective, the possibility to identify reliable vitrification indicators that are solely based on the structural properties of the supercooled fluid is a direct manifestation of the fact that the glass transition line is an isomorph. From an experimental perspective, the vitrification indicators can be used to guide experiments aimed at reaching the glassy state since they can be readily estimated from the radial distribution function or the static structure factor, two quantities which are often easily measured in the course of an experiment either by direct camera observation in the case of soft matter [95] or by neutron diffraction in the case of atomic or molecular systems [96].

Quasi-Universality Aspects
A brief investigation of the glass transition point for the HSPY system and for three inverse power law (IPL) systems with exponents m equal to 4, 9 and 12 revealed that the structural vitrification indicators S R2 and g R2 have a quasi-universal character, i.e., the conditions S R2 = 0.35 and g R2 = 0.30 produce an accurate estimate for the glass transition point regardless of the system under consideration. In the investigation of the HSPY and IPL-m glass transition point, we proceeded as follows: (a) The static properties of the HSPY system were computed via the Wertheim-Thiele analytical solution [84,85], while the static structure factors for the IPL-4, IPL-9, IPL-12 systems were computed with the VMHNC approach. (b) MCT was employed for the IPL systems leading to the glass transition points n = 8.103 for IPL-4,ñ = 1.648 for IPL-9 andñ = 1.322 for IPL-12 whereñ = (β ) 3/m nσ IPL is a temperature-scaled density which fully specifies the state point of an IPL−m system with interaction potential u(r) = (σ IPL /r) m . On the other hand, it was not necessary to solve the MCT for the HSPY system, since it is known that it vitrifies when the packing fraction becomes η = 0.516, see Section 3.3. (c) The vitrification indicators were computed and it was revealed that the four systems satisfied S R2 = 0.35 and g R2 = 0.30 within 3% at their respective glass transition points.
In addition, the quasi-universal character of the thermodynamic vitrification indicator s ex = −5.5, obtained in Section 4.1, was tested. The reduced excess entropy of the HSPY system was computed from the Wertheim-Thiele equation of state [85] and the reduced excess entropy of the IPL systems was computed with the Rosenfeld equation of state [88].