Structure and thermodynamics of two dimensional Yukawa liquids

The thermodynamic and structural properties of two dimensional dense Yukawa liquids are studied with molecular dynamics simulations. The"exact"thermodynamic properties are simultaneously employed in an advanced scheme for the determination of an equation of state that shows an unprecedented level of accuracy for the internal energy, pressure and isothermal compressibility. The"exact"structural properties are utilized to formulate a novel empirical correction to the hypernetted-chain approach that leads to a very high accuracy level in terms of static correlations and thermodynamics.


I. INTRODUCTION
The two-dimensional Yukawa one-component plasma (2D-YOCP) consists of charged point particles that are confined on a two-dimensional surface and are immersed in a polarizable neutralizing background with their interactions described by the Yukawa (screened Coulomb) pair potential φ(r) = (Q 2 /r) exp (−r/λ). Here Q is the particle charge and λ the screening length defined by the polarizable background. Thermodynamic state points of the 2D-YOCP are specified by two dimensionless variables [1, 2]: the coupling parameter Γ = βQ 2 /d and the screening parameter κ = λ/d, where β = 1/(k B T ) with k B Boltzmann's constant and where d = (πn) −1/2 is the 2D Wigner-Seitz radius with n the particle density. The coupling parameter provides a measure of the strength of the unscreened particle interactions with the strong coupled (liquid) regime characterized by Γ > ∼ 1 [3,4], while the screening parameter dictates the interaction softness that varies from infinitely long-ranged Coulomb-like in the one-component plasma (2D-OCP) limit of κ → 0 to extremely short-ranged hard sphere-like in the opposite limit of κ → ∞.
This work focuses on two issues which are still not fully resolved for 2D-YOCP liquids: (a) the acquisition of an accurate equation of state through the reduced excess internal energy that can also be employed to accurately estimate other thermodynamic properties over the entire range of screening parameters relevant to experimental realizations of Yukawa systems, (b) the development of an accurate integral equation theory approach that would allow for the reliable computation of structural properties without necessarily resorting to computer simulations.
In order to address the first issue, systematic molecu-lar dynamics simulations are carried out in the entire 2D-YOCP liquid regime that are utilized for direct extraction of the internal energy, pressure and inverse isothermal compressibility. A novel approach is then presented that simultaneously utilizes these exact thermodynamic data in order to acquire an equation of state through the internal energy that is robust with respect to thermodynamic integration and thermodynamic differentiation. In order to address the second issue, systematic long molecular dynamics simulations are performed in the entire 2D-YOCP liquid regime that are employed for the direct extraction of the radial distribution function and its characteristic functional features. In the absence of a straightforward way to adapt to the 2D-YOCP advanced ultra-accurate integral equation theory approaches that are available for the 3D-YOCP [33], the exact data for the magnitude of the global radial distribution function maximum are used in order to construct an empirical modification to the hypernetted-chain approach. In spite of its simplicity, the emerging approximation leads to very accurate predictions for the structural (and thermodynamic) properties of 2D-YOCP liquids.
The paper is organized as follows. In Section II, the molecular dynamics simulations are presented and their results are discussed. In Section III, the simulation data are employed for the determination of a new 2D-YOCP liquid equation of state which is extensively compared to other equations of state that are already available in the literature. In Section IV, the simulation data are employed for the construction of a novel integral equation theory approximation of empirical nature, whose accuracy and validity region are quantified. In Section V, the results are summarized and possible future developments are discussed.

II. MOLECULAR DYNAMICS SIMULATIONS
A. Simulation parameters Molecular dynamics (MD) simulations were performed to determine the thermodynamic and structural properties of dense 2D-YOCP liquids. MD simulations were carried out for nearly 150 state points characterized by screening parameters that are relevant to experimental 2D-YOCP realizations [12,13,17], i.e. κ = {0.5, 1, 1.5, 2, 2.5, 3}, where Γ OCP m = 131.0 denotes an approximation for the 2D-OCP melting point that is consistent with the results of both computer simulations [34] and experiments [5].
All the simulations were performed with the LAMMPS package [35] and employed 4096 particles in the canonical NVT ensemble. The dynamics was resolved with a timestep of ∆τ = 0.001d √ βm while the interaction potential was truncated at r = 20d for κ = 0.5 and at r = 10d for κ > 0.5. Two set of simulations were performed: one set of short MD simulations consisting of 2 19 time-steps for equilibration followed by 2 19 time-steps for statistics (that were employed to collect 2048 samples for the thermodynamic properties) and one set of long MD simulations consisting of 2 19 time-steps for equilibration followed by 2 24 time-steps for statistics (that were employed to collect 65536 samples for the structural properties).

B. Structural properties
Concerning structural properties, the focus lied on the radial distribution function, g(r), which was extracted from the long MD simulations with the histogram method [36] for a bin-width of ∆r = 0.002d. The narrow bin-width was selected so that the magnitude of the first g(r) maximum is determined with very high accuracy, since it constitutes the MD simulation input that is employed for the construction of our integral equation theory approximation, see section IV for details. Naturally, such a narrow bin-width necessitates longer simulations, since a large number of uncorrelated samples are necessary to obtain highly-resolved radial distribution functions which are unaffected by the omnipresent statistical noise. Some examples of the radial distribution functions obtained from the long MD simulations are illustrated in Fig.2, where it is apparent that the g(r) curves are subject to negligible statistical errors. Further support for the accuracy of the present radial distribution functions comes from the observation that some key figures of merit including the magnitude and position of the first maximum, first nonzero minimum and second maximum are consistent with results available in the literature from Langevin Dynamics simulations for κ = {0.5, 1.0} [21]. Tabulated values of the basic g(r) figures of merit are provided in the supplementary material [37] for all the state points summarized in Fig.1.

C. Thermodynamic properties
Concerning thermodynamic properties, the focus lied on the internal energy (U ), pressure (P ) and inverse isothermal compressibility [K T = −V (∂P/∂V ) T ] with V = N/n the volume of a homogeneous system with N particles and n density. In two dimensions, V is strictly the area. However, aiming to be consistent with the nomenclature developed for three-dimensional systems, we shall still refer to it as volume also for two-dimensional systems. In what follows, we shall discuss normalized (reduced) thermodynamic properties with the internal energy expressed as u = U/(N k B T ), the pressure as p = P/(nk B T ) and the inverse isothermal compressibility as µ = K T /(nk B T ). In addition, since the ideal gas contributions are all known, we shall ignore them and discuss exclusively the reduced excess thermodynamic properties that are emerging from the interaction part of the Hamiltonian.
The reduced excess internal energy is related to the ensemble averaged total potential energy per particle, u ex = U /N with U = N i N j>i βφ(r ij ) and r ij the distance between any particle pair [38]. The reduced excess pressure is related to the (two-dimensional) microscopic virial W = −(1/2) N i N j>i βw(r ij ) with w(r) = r[dφ(r)/dr], through the virial equation p ex = W /(N V ) [38]. The reduced excess inverse isothermal compressibility is obtained from the so-called hypervirial theorem from which   it follows that µ ex = [ W − (δW) 2 + X ]/N , with (δW) 2 = W 2 − W 2 the microscopic virial fluctuations, X = (1/4) N i N j>i βx(r ij ) the (two-dimensional) microscopic hypervirial and x(r) = r[dw(r)/dr] [36].
Therefore, in order to obtain u ex , p ex , µ ex from MD simulations, it is sufficient to collect samples for U, W, X at regular time-intervals in the course of the simulation, to invoke the ergodic hypothesis for the computation of the ensemble averages and to utilize the aforementioned expressions. This procedure was followed for the short MD simulations, since it was observed that relatively few samples are required for the determination of u ex , p ex , µ ex with a negligible statistical uncertainty. Fig.3 illustrates two sample collection examples for all three thermodynamic properties at two 2D-YOCP state points characterized by κ = 3.0. It is evident that all samples fluctuate around a constant average value without systematic deviations; a behavior that confirms that sufficient time was provided before the sampling procedure for the simulated system to efficiently equilibrate. In addition, the magnitude of the fluctuations is extremely small, which confirms that the properties are accurately determined. It is worth noting that compressibility samples cannot be collected directly in the course of the simulation, in contrast to energy and pressure samples. In particular, the hypervirial theorem for µ ex involves the virial fluctuations, which can be evaluated only after the simulation is completed and the average virial is known. Therefore, in order to generate a set of compressibility samples which could be used for uncertainty analysis, we adopted the bootstrap resampling technique [39] that allowed us to construct 2048 compressibility samples from the simulation data at each state point.
The resulting u ex , p ex and µ ex values for all the roughly 150 state points of interest have been tabulated in the supplementary material [37]. Both the average value and the standard deviation of each quantity are provided for the 2D-YOCP state points depicted in Fig.1. All three thermodynamic properties obtained from the MD simulations are determined with a negligible statistical uncertainty that is quantified by a relative standard deviation (defined as the ratio between the standard deviation and the average value) which never exceeds 10 −4 .
It should be pointed out that, when the radial distribution function is known, then the reduced excess internal energy and reduced excess pressure can be computed from the following integral relations [38] This indirect extraction procedure was followed in earlier MD simulation works focusing on the thermodynamics of 2D-YOCP liquids [15,25]. Our direct extraction procedure was preferred because, apart from being much faster, it allows for the quantification of statistical uncertainties in the determination of thermodynamic properties and does not suffer from tail or truncation errors.

III. EQUATION OF STATE
Practical equations of state specify the analytical relation between the reduced excess internal energy and the 2D-YOCP state variables, i.e. u ex (Γ, κ). Once the equation of state is determined, all other thermodynamic properties of the system follow from standard thermodynamic identities. In particular, for 2D-YOCP systems, an analytical expression for u ex (Γ, κ) allows the computation of the reduced excess Helmholtz free energy from the reduced excess pressure from and the reduced excess inverse isothermal compressibility from

A. Equations of state available in the literature
In strongly coupled liquids, the reduced excess internal energy can be conveniently decomposed into the sum of two contributions [40]: a static part u st describing the energy of the system with its constituents frozen in a regular structure (at zero temperature) and a thermal part u th accounting for the finite temperature effects that cause the particles to be displaced from such regular structure. The reduced excess internal energy decomposition reads as u ex (Γ, κ) = u st (Γ, κ) + u th (Γ, κ) for the YOCP with the static component given by For the 3D-YOCP, the Madelung constant is given by a simple closed-form expression [41] that can be obtained from the unitary packing fraction limit (also known as asymptotically high density limit) of the Percus-Yevick approximation for hard spheres [42,43] or the ion-sphere model [44,45]. In addition, Rosenfeld and Tarazona (RT) have shown that the thermal component obeys the particularly simple scaling u th (Γ, κ) ∝ [Γ/Γ m (κ)] 2/5 in the dense fluid region [41,46] where the 3D-YOCP Γ m (κ) is given by Eq.(4) of Ref. [47] and should not be confused with the 2D-YOCP Γ m (κ) that is described by Eq.(1). A particularly attractive feature of the RT scaling has to do with its validity for a variety of three dimensional systems characterized by different interactions and molecular topology [46,[48][49][50]. Furthermore, there exists a deep connection between isomorph theory and the RT scaling, with systems that follow the RT scaling often also being R-simple [49]. In fact, it has been demonstrated that 3D-YOCP liquids are R-simple in an extensive region of their phase diagram [51]. R-simple systems possess isomorph curves, i.e. lines of constant excess entropy along which a large set of thermodynamic, structural and dynamic properties are approximately invariant when expressed in properly reduced units [52,53]. The static part of the excess internal energy naturally produces no entropy, thus the excess entropy is exclusively computed from the thermal part of the excess internal energy via Hence, the RT scaling is compatible with isomorph theory only if Γ/Γ m (κ) is an accurate representation for the isomorphs. The latter is true for the 3D-YOCP [51], but in general the melting line constitutes an isomorphic line only to a first order approximation [54].
For the 2D-YOCP, the situation is more complicated. The Madelung constant does not possess an analytical expression that can be derived from purely theoretical considerations, while the existence of a RT scaling for the thermal component is still open for debate and the functional form of the scaling is unknown. Moreover, no analytical representation is available for the 2D-YOCP iso-morphs, neither is it even known whether the 2D-YOCP is R-simple. This lack of rigorous theoretical foundation for the construction of a 2D-YOCP equation of state has led, over the years, to the emergence of various functional forms for the analytical parametrization of the reduced excess internal energy.
Earlier attempts to obtain an analytical u ex (Γ, κ) expression include the equation of state proposed by Hartmann and collaborators [15], (3,4,7) of Ref. [15] as well as the equation of state presented by Vaulina [55] [55]. In spite of a satisfactory accuracy for κ ≤ 1.5, these u ex (Γ, κ) expressions have two major problems: they become very inaccurate for larger screening parameter values and do not lead to accurate thermodynamic properties via Eqs. (5,6). Thus, in what follows, we focus on two more accurate u ex (Γ, κ) expressions.
Kryuchkov and collaborators have proposed the following u ex (Γ, κ) equation of state that reads as [25] with M (κ) the Madelung constant for 2D-YOCP crystals with triangular lattice that can be fitted with [22,25] while the unknown coefficients of the thermal component are a K (κ) = 0.357 + 0.094κ, b K (κ) = 1.655 exp(−0.769κ), s K (κ) = 0.688 − 0.052κ. It is worth noting that an alternative fit for the thermal component was provided, where all the κ-dependence was absorbed in the form Γ/Γ m (κ) with Γ m (κ) given by Eq.(1) [25]. The accuracy of the latter fit, which predicts u ex within a few percent over the entire dense fluid region of the 2D-YOCP for κ ≤ 3.0 [25], supports the possibility of a modified RT scaling that is applicable to the 2D-YOCP. Finally, Feng and co-workers have proposed the following u ex (Γ, κ) equation of state that reads as [26] where the unknown coefficients for the static and thermal part are described by a F (κ) = 2(0.8394 + 0.5162κ) −6.4 and by b F (κ) = 2 exp(−1.579 − 0.3935κ). Particular care should be taken during the application of these equations of state in the OCP limit (κ = 0). In this limit, it is necessary to replace the reduced excess internal energy u ex with u ex − Γ/κ in order to explicitly take into account the diverging background contribution, Γ/κ. It is evident that the equations of state proposed by Kryuchkov, Hartmann or Vaulina can be safely applied in the OCP limit by simply removing the Γ/κ term, while the equation of state proposed by Feng and collaborators should not be applied in the OCP limit since the divergence is not removable.

B. A new equation of state
The equations of state for u ex (Γ, κ) discussed in Section III A were all obtained by fitting simulation data for the reduced excess internal energy alone. In what follows, a novel approach is presented that determines the equation of state for u ex (Γ, κ) by simultaneously fitting simulation data for the reduced excess internal energy, pressure and inverse isothermal compressibility with the aid of the thermodynamic Eqs. (4,5,6). Initially, the internal energy is expressed as u ex (Γ, κ) = M (κ)Γ + u th (Γ, κ) with M (κ) as given by Eq.(8) and the thermal component defined as This parametrization is applicable in the OCP limit by simply removing the Γ/κ term in the M (κ) expression. The first two terms in Eq.(10) were inspired from the successful equation of state proposed by Hamaguchi and collaborators for 3D-YOCP liquids [56,57] which contains terms proportional to Γ, Γ s and Γ −s with s = 1/3. After trial and error, a different s exponent was adopted and the Γ −s term had to be dropped, since it led to a nonmonotonic pre-factor with respect to κ. The linear term acts as correction to the static component M (κ)Γ. The logarithmic term acts as a residual introduced to adjust the values for small coupling parameters (Γ/Γ m (κ) ≈ 0.1) and was inspired from known low coupling expansions of the 3D-OCP internal energy in terms of Γ ln(Γ) [58]. The κ-dependent coefficients were expressed as Pade' approximants a(κ) = a n 0 + a n 1 κ 1/2 + a n 2 κ + a n 3 κ 2 + a n 4 κ 5/2 1 + a d 1 κ 1/2 + a d 2 κ + a d 3 κ 2 + a d 4 κ 5/2 , Alternative expressions for the Pade' approximants containing only powers of κ or of κ 1/4 were also tested, but proved to be less accurate than the above approximants. The OCP coefficients a n 0 , b n 0 and c n 0 were determined by fitting the OCP simulation results for the thermal component of the reduced excess internal energy that have been tabulated in Table II of Ref. [34]. The remaining coefficients were determined as follows. Starting from the reduced excess internal energy from MD simulations, u MD ex , the thermal component of the excess internal energy was computed as u MD th = u MD ex −M (κ)Γ. Then, u MD th was fitted with Eq.(10) six times, one for each value of κ i belonging to the set κ = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0} and for all corresponding Γ values depicted in Fig.1. The resulting six values of c(κ i ) were fitted with the Pade' approximant given in Eq. (13) to define the coefficients c n 1 , c n 2 , c n 3 , c d 1 and c d 2 . The resulting six values for a(κ i ), b(κ i ) were stored for later analysis. Afterwards, the leading contribution  (11,12,13) that determine the κ-dependent coefficients a(κ), b(κ), c(κ) appearing in the parametrization of the thermal part of the reduced excess internal energy, see Eq. (10). On the other hand, by applying the thermodynamic Eqs. (4,5,6) to the leading component of the fit for the excess internal energy, a(κ)Γ + b(κ)Γ 2/5 , we obtained that the leading contribution to the thermal pressure could be expressed as and that the leading contribution to the thermal inverse isothermal compressibility could be parameterized with The four κ-dependent coefficients in Eqs. (14,15) are connected to the a(κ), b(κ) coefficients in Eq.(10) via Therefore, p MD th,l and µ MD th,l were fitted with Eqs. (14,15) for six κ i belonging to κ = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0} producing twenty-four values for a p (κ i ), b p (κ i ), a µ (κ i ) and b µ (κ i ), six for each type of coefficient. Finally, the sets of coefficients {a n 1 , a n 2 , a n 3 , a n 4 } and {a d 1 , a d 2 , a d 3 , a d 4 }, were determined by simultaneously fitting the six values for a(κ i ) with Eq.(11), the six values for a p (κ i ) with Eq. (16) and the six values for a µ (κ i ) with Eq. (18). An analogous procedure was adopted to define the sets of coefficients {b n 1 , b n 2 , b n 3 , b n 4 } and {b d 1 , b d 2 , b d 3 , b d 4 } which were found by simultaneously fitting the eighteen values for b(κ i ),b p (κ i ) and b µ (κ i ) with Eqs. (12,17,19). The coefficients are summarized in Table I.

C. Level of accuracy of different equations of state
The predictions of the new equation of state, see Eq.(10), and the two recent literature equations of state, see Eqs. (7,9), have been extensively compared against the MD thermodynamic property results presented in Section II C and tabulated in the supplementary material [37]. The summary of this comparison is reported in Table II.
Concerning the reduced excess internal energy, our new equation of state and the Kryuchkov equation of state are visibly more accurate than the Feng equation of state. In particular, the Feng equation leads to exceptionally large errors for κ > 2.0, while the other two equations are both able to predict u ex within 1% over the entire κ ≤ 3.0 range, with the Kryuchkov equation of state having a slight edge. Concerning the reduced excess pressure, the situation is somewhat similar with the important difference that the present equation of state is accurate within 0.5% for any value of κ, whereas the performance of the Kryuchkov equation of state abruptly degrades when κ ≥ 2.5. Concerning the reduced excess inverse isothermal compressibility, the new equation of state remains superior being accurate within 0.3% for any value of κ, while the Kryuchkov equation of state has a high accuracy up to κ = 2.5, but it becomes accurate only within ∼ 20% at κ = 3.0.
Overall, it is concluded that the new equation of state proposed in Section III B leads to improvements over all the other dense 2D-YOCP liquid equations of state currently available in the literature, especially when it comes to predictions of thermodynamic properties for κ > 2.0. Nevertheless, it must be noted that the equation of state proposed by Kryuchkov and collaborators exhibits an excellent agreement with MD simulations, despite possessing a simple κ-dependence for the thermal component. In light of such good agreement, it would seem reasonable to replace Eq.(10) with u th (Γ, κ) =ã(κ) ln[1 +b(κ)s (κ) ]. This possibility was tested but eventually discarded because the coefficientsã(κ),b(κ) ands(κ) showed un unfavorable dependence over κ characterized by changes of sign and stationary points.

A. Method
For an isotropic pair-interacting one-component system, the integral equation theory of liquids enables the computation of two-particle equilibrium correlation functions by combining the Ornstein-Zernike integral equation [38] h(r) = c(r) + n c(r ′ )h(|r − r ′ |)d 2 r ′ (20)  (9) (superscript F) and the results of the MD simulations discussed in Section II C. The average thermodynamic quantity deviation for each screening parameter was computer over all the corresponding coupling parameters, see the summary of Fig.1. Deviations for the reduced excess internal energy are reported in columns 2-4 (ǫu), deviations for the reduced excess pressure are reported in columns 5-7 (ǫp) and deviations for the reduced excess inverse isothermal compressibility are provided in columns 8-10 (ǫµ).
In the above, h(r) = g(r)−1 is the total correlation function and c(r) is the direct correlation function. An expression for the bridge function B(r) is necessary to complete the theory. It is generally prescribed by approximations that are constructed on theoretical grounds [38,59] or that are aided by computer simulations [60,61]. Popular approximations include the hypernetted-chain (HNC) approach which assumes that B(r) = 0 and has proven to be successful for systems with soft interaction potentials [18,19,62,63] or the Percus Yevick approach which assumes that B(r) = ln[1+γ(r)]−γ(r) and has proven to be successful for hard-sphere-like systems [64,65]. Here γ(r) = h(r) − c(r) is the indirect correlation function. When combined with a B(r) assumption, Eqs.(20,21) form a system of equations to be solved for g(r). For its numerical solution, a validated algorithm was followed that was previously applied to three-dimensional systems [33,61]. It is based on Picard iterations in Fourier space combined with a standard mixing technique [62] and a long-range decomposition method [66] in the OCP limit. The convergence criterion was formulated in terms of the Fourier transform of the indirect correlation function, γ(k), and chosen to be ||γ m (k) − γ m−1 (k)|| < 10 −5 ∀k. The two-dimensional Fourier transforms were first expressed as one-dimensional Hankel transforms and were then computed with the Quasi-Fast Hankel Transform algorithm [67] over a discrete grid of N p points. The chosen algorithm allowed to circumvent the unfavorable O(N 2 p ) scaling of direct Hankel transform calculations and did not feature the long-wavelength deficiencies which were observed for similar algorithms that were adopted in earlier works [68][69][70]. Nevertheless, it required the introduction of a computational grid with a constant logarithmic spacing both in real and Fourier space. We employed a grid of N p = 32768 points that extended from 3.5×10 −6 d up to 50d in real space and from 3.5 × 10 −6 /d up to 50/d in Fourier space, with both grids featuring the same logarithmic spacing, i.e. log(r i /r i−1 ) = log(k i /k i−1 ) = 5 × 10 −4 . The algorithm was successfully benchmarked against HNC results for the 2D-YOCP that are available in the literature [19,66].

B. Scaled HNC approach
The HNC approach possesses a reasonable accuracy for the 3D-YOCP, being able to reproduce the main features of the radial distribution function with an accuracy of ∼ 20% and the thermodynamic properties within 5% [61]. However, the 2D-YOCP is known to be richer in structure than its three-dimensional counterpart, which translates to a decline in the accuracy of the HNC predictions. This is demonstrated in Fig.4, where the MD-extracted and HNC-generated radial distribution functions are illustrated for constant (κ, Γ/Γ m ) pairs in the 3D and 2D case. It is evident that the maxima and minima of the radial distribution function become more pronounced and that the deviations between MD results and HNC predictions become larger as the dimensionality decreases.
In spite of these deficiencies of the HNC approach, no advanced integral equation theory approximations have been developed that would lead to more accurate structural predictions for 2D-YOCP liquids. This is in stark contrast to 3D-YOCP liquids, for which two very accurate approximations are available: namely the IEMHNC approach based on the isomorph invariance property of the bridge functions of R-simple systems [61,71] and the VMHNC approach based on the notion of bridge function quasi-universality [72,73]. Such advanced approaches are characterized by an accuracy of < 2% within the first coordination cell of the radial distribution function [33]. Unfortunately, the IEMHNC and VMHNC approach are not directly applicable to the 2D-YOCP: the IEMHNC approach because of the lack of a parameterized 2D-OCP bridge function to be used for the construction of the 2D-YOCP bridge function with the aid of the isomorph mapping [61] and the VMHNC approach because of the lack (a) 3D, κ = 3.0 Results for two state points characterized by κ = 3.0 and two values of the normalized coupling parameter Γ/Γm, namely Γ/Γm = 0.6 (green) and Γ/Γm = 0.9 (magenta). For the 3D-YOCP, Γm is given in Eq.(4) of Ref. [47] and d = (4πn/3) −1/3 . For the 2D-YOCP, Γm is given by Eq.(1) and d = (πn) −1/2 . of a reference system with a known analytical solution to be used for the construction of the VMHNC free energy functional [73]. It is worth to mention the crossover approach of Ref. [74], which constitutes a notable attempt to improve the HNC accuracy for the 2D-OCP.
In light of the above, we opted to take advantage of the vast corpus of MD results presented in Section II C in order to construct an empirical correction to the HNC approach which improves its predictions without requiring any additional input other than that already available. The sought-for correction was constructed by considering that: (a) the main shortcoming of the HNC approach within the first coordination cell refers to the large underestimation of the magnitude of the global g(r) maximum, see the lower panel of Fig.4, (b) our earlier 3D-YOCP work has demonstrated that the HNC approach produces highly accurate structural properties provided that the state point is rescaled towards the stronger coupling region [33]. Given the above, our scaled hypernetted-chain approach (SHNC) was based on retaining the assumption of a vanishing bridge function provided that the interaction strength is up-scaled in a manner that reproduces the exact first peak of g(r). In other words, a mapping is employed from the actual state point (Γ, κ) to another state point (Γ SHNC ≥ Γ, κ) with the unknown Γ SHNC determined by the condition that the HNC approach at (Γ SHNC , κ) leads to the MD-extracted first peak of g(r) at (Γ, κ). The idea of interaction strength rescaling within the HNC approach dates back to the seminal work of Ng for the 3D-OCP [62], whereas an oversimplified version of this idea is encountered in the recent T/2-HNC approach proposed for supercooled dipolar binary mixtures [75].
The SHNC approach is based on the implicit assumption that tempering with the HNC interaction strength in order to ensure an exact global maximum magnitude does not have a detrimental effect on other features of the radial distribution function, especially within the first coordination cell. Since there is no rigorous way of justifying such assumption, the discussion on its validity is postponed to Section IV C, where the predictions of the SHNC approach are compared with "exact" MD results.
In integral equation theory, the bridge function and dimensionless interaction potential appear only in the closure equation combined as βu(r) − B(r). Therefore, our tempering of the interaction potential within the HNC approach is equivalent to approximating the bridge func- . The unknown function Γ SHNC (Γ, κ) is obtained in the following manner: (a) for each (Γ, κ) state point, the coupling parameter is gradually up-scaled and the HNC approach is numerically solved until the first g(r) peak coincides with the respective MD result, (b) this procedure is repeated for all the (Γ, κ) state points considered in the computer simulations reported in Section II A as well as for all the OCP state points simulated in Ref. [21] and a dataset for (Γ, κ, Γ SHNC ) is generated, (c) a closed-form expression for Γ SHNC (Γ, κ) is acquired by sequential least-square fitting with respect to Γ, κ. The SHNC mapping reads as where It is important to emphasize that, because the mapping of Eq.(22) was obtained by fitting, the SHNC approach should not be extrapolated beyond the original fitting region of Γ/Γ m (κ) ∈ [0.1, 1.0] and κ ∈ [0.0, 3.0]. However, since the fit was constructed in such a way that the SHNC approach reduces to the HNC approach for Γ → 0, weak coupling extrapolations are permissible. Therefore, it can be concluded that the 2D-YOCP phase diagram region of validity of the SHNC approach is Γ/Γ m (κ) ≤ 1.0 and κ ≤ 3.0. In other words, provided that the screening parameter is not large, the SHNC can be employed in the entire stable fluid region but not for metastable states.
It is worth noting that we explored the possibility to construct the SHNC mapping by taking advantage of the effective coupling parameter introduced in Ref. [15], Γ * (Γ, κ) = Γf (κ), with f (κ) = 1 − 0.388κ 2 + 0.138κ 3 − 0.0138κ 4 . This effective coupling parameter is related to the melting line parametrization of Eq. (1) which, in fact, can be also expressed via Γ m (κ) = Γ * (Γ OCP m , κ)/f 2 (κ). Since the magnitude of the first maximum of g(r) is approximately constant for state points with the same Γ * (Γ, κ) [15], it should have been possible to construct an SHNC approach where Γ * (Γ, κ) is employed to map any YOCP state point to an effective OCP state point (Γ * , κ = 0) which is then rescaled so that the HNC result for the first g(r) peak coincides with the results of MD simulations. Such a mapping was tested, but was eventually discarded because it showed pronounced deviations from the simulation results in the region characterized by Γ/Γ m (κ) ≤ 0.2 and κ = 3.0 that were traced back to inaccuracies in the parametrization of the effective coupling parameter Γ * (Γ, κ) [15].

C. Structural properties
The HNC and SHNC approaches were numerically solved for all the 2D-YOCP state points illustrated in Fig.1 and all the 2D-OCP state points simulated in Ref. [21]. The computed radial distribution functions have been compared to the ones extracted from computer simulations.
As illustrated in Fig.5, multiple advantages are gained by adopting the SHNC in place of the HNC approach. In fact, modification of the HNC interaction strength to ensure an exact first peak magnitude has a positive effect on all features of the radial distribution function within the first and second coordination cells. To be more specific, apart from the expected enormous improvement concerning the magnitude of the first peak, there is also a strong improvement in the correlation void, the magnitude of the first trough and the magnitude of the second peak as well as a slight improvement in the positions of all peaks and troughs.
The superior performance of the SHNC approach compared to the HNC approach is already evident at small coupling, but it becomes much more pronounced in the vicinity of the melting line. As demonstrated in Fig.6, the SHNC approach is capable of producing accurate predictions for the radial distribution function throughout the whole dense fluid region of the 3D-YOCP phase diagram.
For a more quantitative assessment of the SHNC and HNC accuracy, we report the relative deviations in key g(r) figures of merit: the location of the edge of the correlation void (assumed to be given by the first location where g(r/d) = 0.5) as well as the magnitudes and positions of the first maximum, first non-zero minimum and second maximum. The values of these quantities for "exact" radial distribution functions are available in Table  1 1 of Ref. [21] for κ = 0.0 and in the supplementary material [37] for 0.5 ≤ κ ≤ 3.0. The most noticeable SHNC improvements take place in the location of the correlation void and in the magnitude of the first maximum. The former is predicted within 1% from the SHNC approach and only within 5% from the HNC approach. The latter is predicted within 0.5% from the SHNC approach, whereas it is consistently strongly underestimated in the HNC approach with the relative deviations even exceeding 25% close to the melting line. Concerning the locations of the first and the second maximum as well as of the first non-zero minimum, the SHNC slightly improves the HNC predictions which are anyways well within 5%. Furthermore, both approaches produce rather poor estimates for the magnitude of the first non-zero minimum, but the HNC approach is much more accurate than the SHNC approach with 20% vs 40% mean relative deviations from the MD results. Finally, concerning the magnitude of the second maximum, the SHNC is accurate within ∼ 2% and the HNC only within ∼ 10%.

D. Thermodynamic properties
The performance of the HNC and the SHNC integral equation theory approximations was also evaluated at the level of the thermodynamic properties. For this purpose, Eqs.(2,3,6) were employed to compute the excess internal energy, excess pressure and excess inverse isothermal compressibility from the radial distribution functions obtained with the two approaches. The results were then compared with the thermodynamic properties that were extracted from MD simulations and tabulated in the supplementary material [37].
Near the OCP limit, both approaches are able to reproduce all the three thermodynamic properties within 1%, namely the SHNC within ∼ 0.2% and the HNC within ∼ 0.6%. However, the accuracy of the SHNC approach remains almost constant with the screening parameter, whereas the performance of HNC approach promptly degrades as the screening parameter increases. Concerning the excess internal energy, the mean deviations reach 2% for the SHNC and 10% for the HNC for κ = 3. Concerning the excess pressure, the mean deviations reach 2% for the SHNC and 8% for the HNC for κ = 3. Concerning the excess inverse isothermal compressibility, the mean deviations remain < 1% for the SHNC and reach 5% for the HNC for κ = 3. To sum up, it can be concluded that the SHNC approach reproduces the excess internal energy, pressure and isothermal compressibility within 2% over the entire dense fluid region, while the HNC approach produces estimates which are accurate within 10%.
It is important to point out that the virial route was followed for the computation of the excess inverse isothermal compressibility. The utilization of the so-called statistical route, i.e. µ ex = −n c(r)d 2 r [38], would result to large deviations from the exact results for both approximations and especially for the SHNC approach. The SHNC bridge function does not obey the correct asymptotic limit B(r) → −h 2 (r)/2 but decays as B(r) ∝ βu(r), which suggests that the exact asymptotic limit of the direct correlation function c(r) → −βu(r) should also be violated. In fact, this can be rigorously proven from the asymptotics of the non-linear closure condition, see Eq. (21). It is known that the asymptotic range provides large contributions to the above µ ex expression, which explains why the statistical route should be avoided. Naturally, this brings forth an inherent problem of the SHNC approach: its thermodynamic inconsistency.

E. Comments on the T/2-HNC approach and the metastable states
The T/2-HNC approach utilizes the HNC approximation at a reduced half temperature. For the YOCP, this approach corresponds to a simplified version of the SHNC approach for which Γ SHNC (Γ, κ) = 2Γ. This empirical approach was recently applied to determine the glass transition line of the 2D-YOCP with the aid of mode coupling theory [18]. This was carried out without any discussion concerning the validity of the T/2-HNC approximation for metastable or even for stable 2D-YOCP liquids. Such an analysis will be performed in what follows. The T/2-HNC approach was initially proposed for twodimensional binary mixtures of point-dipoles interacting via a ∝ r −3 pair potential. For such systems it was observed empirically that, if the state point temperature is rescaled from T to T /2, then the HNC approach can be employed for fairly accurate estimates of the structural properties [75]. Given the dependence of bridge functions on the softness (see the successes of the HNC and of the Percus-Yevick approaches) and the empirical nature of the re-scaling, it is evident that the T/2-HNC approach should not be applied to other systems without prior verification of its accuracy. Some work in this direction was performed in Ref. [18], where some qualitative agreement between Monte Carlo simulations and the T/2-HNC approximation was reported for the 2D YOCP. However, (a) these simulations were performed in the stable fluid regime, thus the extension to the supercooled regime involves an unjustified extrapolation, (b) the documented accuracy of the SHNC approach within the stable fluid region and the strong deviations of the SHNC mapping from the T/2-HNC mapping (see Fig.7) prove that the T/2-HNC approach does not lead to accurate 2D-YOCP structural properties even in the stable fluid region. Generally speaking, any approximate integral equation theory closure that is derived from computer simulations should only be used within its range of validity (determined by the simulation input employed to construct it). Extrapolations outside the original range of validity are sometimes possible [76], but should always be performed with great care. Hence, considering that the range of validity of the SHNC approach is Γ/Γ m (κ) < 1.0 (see the dashed gray line in Fig.7), that the T/2-HNC approach performs poorly even in the stable fluid region and that the HNC approach is expected to perform poorly in the supercooled regime [76], it can be concluded that, at the moment, there is no integral equation theory approximation which can be employed to accurately predict the structural properties of supercooled 2D-YOCP liquids.

V. SUMMARY AND FUTURE WORK
The structural and thermodynamic properties of dense two-dimensional Yukawa liquids were extensively investigated with molecular dynamics simulations. The "exact" thermodynamic properties were employed in order to construct a new equation of state for the excess internal energy valid in the parameter regime most relevant for contemporary experiments. Our equation of state exhibited excellent agreement with the simulation results and, contrary to most 2D YOCP equations of state available in the literature, proved to be robust with respect to thermodynamic integration and differentiation. The "exact" structural properties were employed to formulate the scaled hypernetted-chain approach that is constructed by up-scaling the interaction strength until the bare HNC recovers the exact magnitude of the first peak of the radial distribution function. The SHNC was demonstrated to significantly improve the HNC structural predictions and to achieve a 2% accuracy in thermodynamic quantities.
For future improvement of the present results, it would be important to confirm that the 2D-YOCP is R-simple, to numerically trace out multiple isomorphic lines and to determine an accurate analytical parameterization of these isomorphs. This would allow for the construction of more accurate equations of state in the spirit of the Rosenfeld-Tarazona scaling and would also allow for the development of the isomorph-based empirically modified hypernetted chain approach for the 2D-YOCP. The latter integral equation theory approximation should lead to unprecedented levels of accuracy superior to that of the scaled hypernetted chain approach, but it requires an analytical parametrization for the 2D-OCP bridge function whose extraction from simulations is a formidable task.