Testing the isomorph invariance of the bridge functions of Yukawa one-component plasmas. II. Short range

It has been conjectured that bridge functions remain nearly invariant along phase diagram lines of constant excess entropy for the class of R-simple liquids. In the companion paper, this hypothesis has been confirmed for Yukawa bridge functions outside the correlation void. In order to complete the testing of the invariance ansatz, Yukawa bridge functions are here computed inside the correlation void with the cavity distribution method and input from ultra-long molecular dynamics simulations featuring a tagged particle pair. A general methodology is developed for the design of the tagged pair interaction potential that leads to the acquisition of uniform statistics. An extrapolation technique is developed to determine the bridge function value at zero separation. The effect of different sources of errors is quantified. Yukawa bridge functions are demonstrated to be nearly isomorph invariant also in the short range. Generalization to all R-simple systems and practical implications are discussed.


I. INTRODUCTION
One of the fundamental problems in the statistical mechanics of liquids concerns the accurate computation of static pair correlations with the knowledge of the pair interaction potential alone and without resorting to computer simulations.The integral equation theory of liquids features two formally exact equations for this problem that contain three unknown functions; an additional equation for the so-called bridge function is missing [1,2].The unknown bridge function incorporates all the elements that make a many body problem of infinite degrees of freedom unsolvable.Therefore, it is not surprising that its diagrammatic expansion is very slowly converging and its high-order terms quickly become overly complicated to calculate [3].As a consequence, numerous approximation schemes have been developed for the bridge function, whose effectiveness varies depending on the potential and can only be reliably evaluated a posteriori through a comparison with "exact" simulation results [1,4,5].It is fortuitous, though, that static correlations in liquids exhibit relatively weak dependence on the bridge function [1].
Contrary to the radial distribution function, the bridge function neither possesses microscopic representation (in terms of δ−functions and the instantaneous particle positions) nor a conditional probability interpretation.Thus, the bridge function cannot be directly extracted from computer simulations.Nevertheless, it can be computed with input from computer simulations.In particular, extraction of radial distribution functions leads to the closure of the system of integral equation theory and allows one to solve for the unknown bridge function.There are two main caveats with such an inversion method.First, the weak dependence of static correlations on the bridge function implies a strong sensitivity of the bridge function on the radial distribution function, in other words what is an asset in the direct problem becomes an obstacle in the inverse problem, which necessitates long simulations with large particle number.Second, irrespective of the achieved size of the statistical sample, the inversion method is doomed to fail at very short distances where the probability of encountering another particle is ultra low (for finite thermodynamically stable potentials) or is even zero (for unbounded potentials).Cavity simulations featuring a special tagged particle pair and utilizing umbrella sampling techniques have been conceived for the computation of the bridge function at such distances [6].
Despite the undeniable progress during half a century of investigations, few exact or approximate properties of bridge functions have so far been discovered.Recently, a novel integral equation theory approach has been formulated that is based on the conjecture that bridge functions remain invariant along phase diagram lines of constant excess entropy for a broad class of liquids known as R-simple [26].It has been coined as isomorph-based empirically modified hypernetted-chain (IEMHNC) and has been applied to Yukawa and bi-Yukawa liquids resulting in a remarkable agreement with simulations [26][27][28].
The main objective of the present study is to test the validity of the underlying ansatz of bridge function invariance for Yukawa one-component plasmas.In the present article, the intermediate and long range of Yukawa bridge functions will be computed along isentropic lines with ultra-accurate MD simulations.In the accompanying article, the short range of Yukawa bridge functions will be computed with specially designed MD simulations 1 .
The paper is organized as follows.Section II features an introduction to Yukawa one-component plasmas, isomorph theory, R-simple systems and discusses arguments in favor of the isomorph invariance of bridge functions.In section III, 16 Yukawa state points are identified with isomorph tracing techniques distributed amongst four isomorphs.In section IV, the inversion method is presented, MD simulation parameters are specified and bridge functions are computed for all state points.In section V, the bridge function sensitivity to artificial periodic & aperiodic multiplicative perturbations of radial distribution functions is studied.In section VI, the effect of statistical, grid, finite size, tail and isomorphic errors in the bridge function is quantified.In section VII, corrected bridge functions with error bars due to propagating uncertainties are presented and the degree of isomorph invariance is discussed.In section VIII, the results are summarized.

II. BACKGROUND
For the present article to be self-contained, following a brief introduction to the standard nomenclature of the Yukawa one-component plasmas, a concise primer focusing on the isomorph theory of R-simple systems and the isomorph-based empirically modified hypernetted-chain approach is given in the present section.The reader is addressed to the references cited below for further details.

A. Yukawa one-component plasmas
Yukawa one-component plasma (YOCP) systems are model systems whose constituents are charged point particles which are immersed in a neutralizing background and interact via the screened Coulomb (Yukawa) pair potential u(r) = (Q 2 /r) exp(−r/λ) where Q is the particle charge and λ the screening length determined by the polarizable background.It is convenient to specify the thermodynamic state points of the YOCP in terms of two independent dimensionless variables, the coupling parameter Γ and the screening parameter κ defined by [30][31][32][33] In the above; β = 1/(k B T ) with k B the Boltzmann constant, T the temperature and d = (4πn/3) −1/3 for the Wigner-Seitz radius with n the particle (number) density.
In the limit of a rigid background λ → ∞ or κ → 0, the Yukawa potential collapses to the unscreened Coulomb potential and the resulting model system is then known as one-component plasma (OCP).The YOCP enables the exploration of the full range of potential softness from the long range Coulomb interactions of the OCP for κ = 0 to ultra-short range hard-sphere interactions for κ → ∞.Due to its variable softness and its relevance to strongly coupled laboratory systems such as complex plasmas and colloidal suspensions [34,35], the YOCP is still being actively investigated in statistical mechanics studies.
It is worth noting that the distances are typically normalized by the Wigner-Seitz radius d = (4πn/3) −1/3 in the non-ideal plasma literature, while the distances are typically normalized by the mean-cubic inter-particle distance ∆ = n −1/3 in the liquid state and isomorph theory literature.Both normalizations will be used in the present work, but mainly the plasma normalization to remain consistent with the screening parameter definition.

B. Isomorph theory and R-simple systems
Isomorphic lines or simply isomorphs are phase diagram curves of constant excess entropy, along which a large set of structural and dynamic properties are approximately invariant when expressed in properly reduced units [36][37][38].In case of Newtonian dynamics, the length is normalized to the mean-cubic inter-particle distance ∆ = n −1/3 , the energy is normalized to the thermal energy k B T and the time is normalized to the time required for a particle that is free streaming with its thermal velocity to traverse an inter-particle distance τ = n −1/3 m/(k B T ) [38].All systems have isentropic curves in their thermodynamic phase diagram, but these are termed isomorphs only for the so-called Roskilde-simple or R-simple systems.
R-simple systems are rigorously defined as many body systems that possess the property that the ordering of the total potential energies of two configurations consistent with the same density is maintained when these two configurations are uniformly scaled to a different density [39].
for positive µ, where U (R) is the total potential energy, R is the particle configuration that is given by the collective N-particle position vector (r 1 , ..., r N ) and with R a , R b denoting two equal density configurations [39].The hidden scale invariance property is exact only for systems that are characterized by Euler-homogeneous interactions (plus a constant), such as inverse power law (IPL) systems.For other R-simple systems, hidden scale invariance should be understood to be valid for most of the physically relevant configurations reflecting the approximate nature of isomorph theory.
R-simple systems are practically identified as systems possessing strong correlations between their virial (W ) & potential energy (U ) constant-volume thermal equilibrium fluctuations [40].The degree of W − U correlations is quantified by the standard Pearson coefficient given by where ... NVT denotes canonical ensemble averaging and ∆A = A− A NVT denotes statistical fluctuations around the canonical mean.Strong W − U correlations are empirically delimited by the practical condition R W U > ∼ 0.9 that allows for straightforward identification of R-simple systems with canonical (NVT) computer simulations.
A recent computational investigation revealed that the YOCP is an R-simple system that exhibits exceptionally strong W − U correlations (R W U > 0.99) for an extended part of the fluid phase covering the entire dense liquid region of the phase diagram [41].This rationalizes a number of previous observations such as the fact that the YOCP excess internal energies conform to the Rosenfeld-Tarazona decomposition [42,43] as well as the fact that the YOCP reduced transport coefficients strongly abide to Rosenfeld's excess entropy scaling [44,45].

C. The isomorph-based empirically modified hypernetted chain approximation
The isomorph-based empirically modified hypernetted chain (IEMHNC) approximation is an integral equation theory approach that is based on the assumption of isomorph invariance of bridge functions when expressed in reduced distance units [26].The invariance ansatz closes the non-linear non-local equation system that arises in integral equation theory provided that two external inputs are also available: a closed-form expression for the dependence of the isomorphic curves on the thermodynamic state points and a closed-form bridge function expression that is valid along any phase diagram line that possesses a unique intersection point with any isomorphic curve [26].
With such input, the isentropic correspondence maps the bridge function from the initial phase diagram line to the entire phase diagram.
The IEMHNC approach has been successfully applied to dense Yukawa and bi-Yukawa liquids [26][27][28] taking advantage of an established parameterization of the OCP bridge function through the reduced distance and coupling parameter [17].Comprehensive benchmarking with computer simulations has revealed that the IEMHNC approach possesses a remarkable accuracy with predictions of structural properties within 2% inside the first coordination cell and predictions of thermodynamic properties within 0.5% in the entire dense liquid regime [26,27].In addition, systematic comparison with different advanced integral equation theory approaches has demonstrated that the performance of the IEMHNC approach is comparable to that of the variational modified hypernettedchain approach (VMHNC) [46] but with 10−80 times less computational cost depending on the state point [28].

D. Theoretical arguments in favor of the isomorph invariance of bridge functions
The excellent performance of the IEMHNC approach for YOCP systems [26] and for biYOCP systems [27] clearly suggests that the underlying conjecture of the isomorph invariance of bridge functions must hold to a high degree.This is also indicated by the fact that the IEMHNC approach preserves its OCP level of accuracy regardless of the value of the YOCP or biYOCP screening parameter.
Moreover, the VMHNC bridge function has been revealed to be implicitly isomorph invariant for the YOCP, since the effective packing fraction acquired by minimizing the respective free energy functional has been demonstrated to remain nearly constant along any isomorphic curve within the dense liquid regime [28].Finally, the output of the classic hypernetted-chain (HNC) approach, that completely neglects all the bridge diagrams, leads to approximately invariant static correlations for the YOCP [26], which implies that the addition of an isomorph invariant bridge function would be beneficial for this isomorph invariance to persist.Further arguments in support of the bridge function isomorph invariance are connected to the notion of bridge function quasi-universality [47] that forms the backbone of the powerful modified hypernetted-chain (MHNC) and reference hypernetted-chain (RHNC) approaches.This quasi-universality notion can be summarized in the statement that, in their short range, the bridge functions constitute the same universal family of curves irrespective of the interaction potential and it was based on the fact that bridge functions can be expressed as densely connected diagrams containing total correlation function bonds [47].Considering the isomorph invariance of the total correlation functions, the same reasoning can be extended to the notion of isomorph invariance.The isomorph theory has already rationalized a number of well-established quasiuniversalities of simple liquids, since the excess entropy always turned out to be the controlling parameter [38].In addition, the isomorph invariance of bridge functions is consistent with the zero-separation bridge function freezing criterion of Rosenfeld, which states that the value of the bridge function at the origin r = 0, when calculated along the liquid-solid phase transition line, is nearly constant and even independent of the pair potential [48].This criterion is known to be satisfied for the YOCP [49].
Finally, let us discuss an apparent incompatibility of bridge function isomorph invariance with the condition of unique functionality.This condition assumes that the exact functional relation between the bridge function and indirect correlation function approximately reduces to a unique function [50].This condition is implicitly invoked in most fundamental bridge function closures of integral equation theory [50].However, given the isomorph variance of the indirect correlation function (to be revealed in the following sections), it also implies that the bridge function cannot be isomorph invariant.Conversely, the approximate bridge function properties of isomorph invariance and unique functionality are incompatible.This does not constitute a contradiction, since it has been revealed (with the use of Duh-Haymet plots) that the above formulation of the unique functionality condition does not hold [51].In fact, the optimized unique functionality conditions that are invoked in more modern approaches feature a re-normalized indirect correlation function, typically stemming from a partition of the interaction potential [52][53][54][55].Such formulations are not incompatible with the isomorph invariance ansatz.

III. ISOMORPH TRACING AND STATE POINTS OF INTEREST
Different methods are available for the tracing of the isomorph curves of R-simple systems with or even without the use of computer simulations.In the present investigation, we shall employ the direct isomorph check, the small step method and the analytical method.The physical basis and the numerical implementation of these methods are briefly described below.
The direct isomorph check is based on an approximate relation valid for any state point that is a fundamental characteristic of R-simple systems and reads as [39] where U (R) is the instantaneous potential energy which depends on the configuration R consistent with any state point (n, T ), S ex ( R) is the instantaneous excess entropy function that depends on the reduced configuration R = n 1/3 R and U (n, S ex ) is the thermodynamic (ensemble averaged) potential energy.Let us suppose a (n 1 , T 1 ) reference state together with its isomorphic (n 2 , T 2 ) state of re-scaled density n 2 = (1/µ 3 )n 1 but unknown temperature T 2 .Let us also consider the configurations R 1 , R 2 of these state points that are identical in reduced units, n Application of the above relation for both the state points, first-order Taylor expansion with respect to S ex ( R) around the thermodynamic excess entropy S ex , utilization of the identity (∂U/∂S ex ) n = T as well as use of the identical reduced entropies and reduced configurations, leads to the approximate expression [39] U The above expression constitutes the basis of the direct isomorph check that in practice works as follows [56,57]; The potential energy U (R 1 ) is extracted from a (n 1 , T 1 ) simulation, the configuration is rescaled to R 2 = µR 1 and the potential energy U (R 2 ) is extracted.Repetition of this procedure for numerous R 1 configurations leads to a scatter plot between U (R 1 ) and U (R 2 ) that is well approximated by a straight line and whose linear regression slope T 1 /T 2 allows the determination of the unknown T 2 .
In the present application of the direct isomorph check, the algorithm is formulated in terms of (Γ, κ) and a fixed screening parameter jump ∆κ/κ = 3.1% is considered which translates to a |∆n|/n = 9.8% density variation between successive isomorphic state points.In the NVT MD simulations that are necessary for the slope extraction, reduced units are employed by setting the temperature and density equal to unity and controlling the length and energy parameters of the potential.The interaction potential is truncated at r cut = 10∆ with the shiftedforce cutoff method, the time step is ∆t/τ = 2.5 × 10 −3 , the equilibration time is 2 20 ∆t, the statistics duration is 2 20 ∆t, the saving period is 2 10 ∆t, and the number of particles is 8192 (20∆ for the simulation box length).
The small step method combines the thermodynamic definition of the so-called density-scaling exponent with an exact alternative expression that originates from thermodynamic fluctuation theory.The density-scaling exponent γ(n, T ) is defined in log-log density-temperature phase diagrams as the local slope of the isentropic line traversing the state point (n, T ).Thus, we have [58,59] γ(n, T ) = ∂ ln T ∂ ln n Sex , with S ex denoting the excess entropy.The density-scaling exponent is also acquired by the linear regression slope of the scatter plot between the virial and potential energy canonical fluctuations, since we also have [58,59] γ(n, T ) = ∆U ∆W NVT (∆U ) 2 NVT .
The fluctuation theory expression allows the evaluation of γ(n, T ) at any state point from canonical simulations, whereas the thermodynamic definition constitutes an explicit non-linear first-order differential equation with respect to the state points that can be solved with any numerical scheme in order to trace the respective isomorph.This procedure has been coined as small step method, because its typical applications utilize first-order numerical schemes for the solution of the differential equation that necessitate small steps in the density [60,61].However, implementation of higher-order schemes, such as the classical fourth-order Runge-Kutta method (RK4), allows for larger density increments and leads to equally accurate isomorph tracing but with far less computational cost.
In the present application of the small step method, the RK4 algorithm is formulated in terms of (ln n, ln T ) and a fixed logarithmic density step is considered which translates to a |∆n|/n = 8.8% density variation between successive isomorphic state points.In the NVT MD simulations that are necessary for γ extraction, natural units (n, T ) are employed, the Yukawa pair potential is truncated at r cut = 10d with the shifted-force cutoff method, the time step is ∆t/τ = 2.5 × 10 −3 , the equilibration time is 2 17 ∆t, the statistics duration is 2 17 ∆t, the saving period is 2 7 ∆t, and the number of particles is 8192 (32d for the simulation box length).
The analytical method exploits a number of exact properties of inverse power law potentials in order to define an approximate distance dependent IPL-like exponent for arbitrary pair potentials [62].Complemented with a realistic estimate for the effective nearest-neighbor distance, such a definition leads to an approximate relation for the density-scaling exponent and a closed-form expression for any family of isomorphic curves [41,62].The application of the analytical method to YOCP systems results in [41] Γ ISO (κ)e  denotes the ratio between the mean-cubic inter-particle distance and the Wigner-Seitz radius.The simple choice Λ = 1 has proven to be very accurate for the YOCP [41], The above expression is identical to the well-known semiempirical description of the YOCP melting line [63,64] Γ m (κ)e which accurately follows the near-exact data obtained by MD simulations [65,66].In the above, Γ OCP m = 171.8 is the OCP coupling parameter at melting [65].It is evident that all isomorph lines are nearly parallel to the melting line, an observation that is true for any R-simple system to the first order [58,67].
In the present work, the analytical method will only be invoked in order to specify the OCP members, Γ OCP ISO , of YOCP isomorphs.The mapping should be very accurate, since Eqs.(1,2) are nearly exact for κ < ∼ 1.5.In addition, the equivalence of Eqs.(1,2) results in Γ/Γ m = const.≤ 1 along any distinct isomorph.For brevity, in what follows, the constant approximate values of Γ OCP ISO or Γ/Γ m will be utilized in order to uniquely identify the isomorphs.
YOCP bridge functions will be determined for 16 state points that are equally spread amongst four different isomorphic curves.The normalized screening parameters κ The four isomorphic curves as obtained from the direct isomorph check (discrete symbols) are compared to those obtained from the approximate analytical expression (solid lines).The four numerical isomorphs begin to overshoot the respective analytical isomorphs roughly above κ ≃ 1.5, see also Table I.
of interest are κ = (1.0,1.5, 2.0, 2.5), since for κ < ∼ 1 the YOCP behavior becomes nearly OCP like while for κ > ∼ 3 the YOCP behavior becomes nearly hard sphere like.Such values are typically realized in complex plasma microgravity experiments [68,69].The OCP coupling parameters of interest are Γ OCP ISO = (160, 120, 80, 40) and span the whole dense YOCP liquid regime, since they correspond to Γ/Γ m = (0.93, 0.70, 0.47, 0.23), respectively.The κ = 1 members of these four isomorphs are then calculated with the analytical method, Eq.( 1), that leads to Γ κ=1 ISO = (205.061,153.796, 102.531, 51.265) which are the starting state points of the direct isomorph check and the small step method.The isomorphic coupling parameters are determined exactly at the remaining screening parameters (κ = 1.5, 2.0, 2.5) by a targeted jump from the closest κ point that emerges from the algorithms depending on their assumed density variations.
The YOCP state points of interest are listed in Table I.In all cases, the relative deviations between the results of the direct isomorph check and the small step method are always less than 0.53%.Since the two MD implementations are characterized by the same number of particles (2 13 ) and statistically useful configurations (2 10 ), this should be a consequence of the very high virial -potential energy fluctuations R W U ≥ 0.988.This near-unity correlations imply that the starting equation of the direct iso-morph check is exact, similar to the starting equation of the small step method.The YOCP state points that stem from the direct isomorph check were selected for bridge function computation.On the other hand, the absolute relative deviations between the results of the analytical method and the direct isomorph check can reach 4.17%.A consistent overestimation is observed for the analytical method above κ ≃ 1.5, as better illustrated in figure 1.This is rather expected because the analytical expression for the melting line, see Eq.( 2), also overshoots the MD results above κ ≃ 1.5.Notice that the density-scaling exponent γ, constant for inverse power law systems, varies from 0.5 (κ = 1.0) to 1.2 (κ = 2.5) within each isomorph.
Finally, it is worth pointing out that the above isomorph tracing methods manage to identify the isentropic points of R-simple systems without ever calculating the excess entropy.Accurate theoretical determination of the excess entropy can be formidable due to the need for either thermodynamic integration or high-order correlation inclusion (see the Nettleton-Green expansion) [70].The same applies for the computational determination owing to the inefficiency of Widom test particle insertion methods at high densities [45].For completeness, an estimate of the reduced excess entropy of each isomorph line has been attempted based on the equation of state suggested by Hamaguchi et al. [65,66].This led to (Γ OCP ISO , −s ex ) = {(40, 2.02), (80, 2.87), (120, 3.49), (160, 3.99)}.

IV. BRIDGE FUNCTIONS DETERMINED BY THE ORNSTEIN-ZERNIKE INVERSION METHOD
By inspecting the building blocks of the integral equation theory of liquids, it becomes evident that the bridge function acts as an additional many-body component of the pair interaction potential [47].Hence, the computational technique utilized for the deduction of bridge functions from simulation structural data is identical to the computational method employed for the deduction of interaction potentials from experimental structural data [71,72].
We will refer to it as Ornstein-Zernike inversion method, since the aforementioned determination of the pair interaction potential is known as the inverse problem [73].In this section, the inversion method will be described, the parameters of the production runs or test simulations will be provided and the numerical results will be analyzed.
A. The computational method In the case of one-component pair-interacting isotropic systems, the integral equation theory of liquids consists of the convolution-type Ornstein-Zernike (OZ) equation [1] h combined with the following exact non-linear closure condition that is derived from cluster diagram analysis [1] g where g(r) is the radial distribution function, c(r) is the direct correlation function, h(r) = g(r) − 1 is the total correlation function and B(r) the bridge function.Other static correlation functions of relevance concern the indirect correlation function γ(r) = h(r) − c(r), the potential of mean force βw(r) = − ln [g(r)] and the screening potential βH(r) = βu(r) − βw(r).A formally exact expression for the bridge function that is required to close the system of equations is unavailable due to the highly connected nature of bridge diagrams.In addition, the exact virial-type series that define bridge functions through Mayer functions or total correlation functions converge very slowly and are notoriously hard to compute [3,74].
For this reason, integral equation theory approaches invoke approximations that either prescribe B(γ) of a function of the indirect correlation function or parameterize B(r/d) as a function of the normalized distance [1,4,5].
In integral equation theory approaches, Eqs.(3,4) are solved for [g(r), c(r)] with known u(r) and an assumption for B(r) which suggests that the equations are coupled.In pair interaction reconstruction, Eqs.(3,4) are solved for [c(r), u(r)] with known g(r) and an assumption for B(r) which suggests that multiple viable solutions can emerge.In bridge function reconstruction, Eqs.(3,4) are solved for [c(r), B(r)] with known [g(r), u(r)] which suggests that the equations are decoupled and a unique solution exists.
The computation of bridge functions from MD simulations proceeds in the following manner.(1) The radial distribution function is extracted from MD simulations with the histogram method.Constant ∆/r = 0.002d bin widths are assumed, see section VI B for a detailed justification.(2) The Fourier transform of the total correlation function is calculated.Invoking the spherical symmetry, )i] for the discrete sine transform of a space resolution equal to the bin width (to avoid interpolations).In the above, W is the histogram bin number, ∆k = π/(N ∆r), k i = ı∆k, r j = j∆r − ∆r/2, r = {r i } and k = {k i }.Fast Fourier Transform (FFT) algorithms are employed in order to reduce the computational cost.(3) The Fourier transform of the direct correlation function is computed.By Fourier transforming the OZ equation and solving for C(k), we obtain The direct correlation function is calculated from the inverse Fourier transform.Invoking spherical symmetry, we acquire the expression c Inside the correlation void, that can be loosely defined as arg r {g(r) ≪ 1} or arg r {g(r) ≃ 0}, particle encounters are extremely rare but the probability remains finite.As a consequence, the histogram method could, for instance, lead to either g(r) = 10 −8 or to g(r) = 10 −12 due to the poor statistical sampling.The mathematical structure of the OZ integral equation ensures that this enormous statistical uncertainty does not propagate up to the direct correlation function.However, because of the logarithm that is present in the OZ closure equation which becomes B(r) = ln[g(r)] + 1 + c(r) + βu(r) at very short distances, it strongly impacts the bridge function leading to a significant −8 or −12 contribution for this example.In conclusion, the unavoidable insufficient statistics within the correlation void suggest that the OZ inversion method is only effective for the computation of the bridge function at intermediate and long ranges.For a given interaction potential, the validity limit mainly depends on the state point of interest, the overall statistics (number of particles and number of uncorrelated configurations) and the desired accuracy.

B. The production runs
The production runs for the extraction of the radial distribution functions, as well as the tracing of the isomorphic curves, were carried out on graphics cards with the RUMD open-source software [75].A small number of test runs were performed with the LAMMPS package [76].
In the production runs dedicated to g(r) extraction; the (canonical) NVT MD simulations utilize the shiftedforce cut-off method with the Yukawa pair potential truncated at r cut = 10d and the time-step employed for the propagation of equations of motion is ∆t/τ = 2.5 × 10 −3 .The MD equilibration time is 2 20 ∆t, the statistics duration is 2 23 ∆t and the configuration saving period is 2 7 ∆t leading to M = 2 16 = 65536 for the statistically independent configurations.The particle number is N = 54872 leading to 60d for the cubic simulation box length.The bin width size of the histogram method is ∆r/d = 0.002d.
The configuration saving period was selected so that uncorrelated radial distribution functions are always extracted, see section VI A 1. The combination of the statistics duration and the number of particles was selected so that sufficient pair correlation sample sizes are collected even up to r = 1.25d, see section VI A. The size of the histogram bin width was selected so that grid errors are much smaller than statistical errors, see section VI B. A large number of test runs were carried out in order to choose a near-optimal cut-off method, truncation radius and MD time-step.These runs were carried out at the YOCP state point Γ OCP ISO = 160, κ = 1.0 for which the bridge function exhibits the highest sensitivity to uncertainties in the radial distribution function, see section V. Finally, some test runs were performed with the RUMD and LAMMPS softwares for the same YOCP state points and for identical simulation settings as a validation check.

C. The numerical results
The sequential output of the OZ inversion method, i.e. g(r) → c(r) → B(r), is illustrated in figure 2 for the 4 isomorphic curves and 16 YOCP state points of interest.
The radial distribution functions g(r/d) along each isomorph are illustrated in the first panel within the interval r/d ≤ 6, see subfigures 2 (a), (d), (g), (j).In addition, the potentials of mean force − ln [g(r)] are plotted within the range 1.25 ≤ r/d ≤ 2 in the respective insets.As demonstrated from earlier investigations of the YOCP [41], the radial distribution function is a strongly invariant quantity along any isomorphic curve with the exception of a narrow interval around the first maximum.This invariance holds to a very good approximation outside the correlation void but is rapidly distorted at short distances.The local variance becomes especially apparent when inspecting the potentials of mean force for r/d < ∼ 1.5.To be more concrete, the g(r) deviations between isomorphic state points reach two orders of magnitude at r = 1.25d for Γ OCP ISO = 160 and this trend is expected to get further augmented at shorter distances.This behavior does not contradict the basic property of R-simple systems which states that they possess approximate invariant structural properies in reduced r/d units, because it manifests itself in short distances where the radial distribution function can be approximated with zero and its exact values are inconsequential.In other words, this behavior concerns ultra-rare structural configurations that are physically insignificant.Actually, the observed correlation void variance is a direct consequence of the asymptotic limit of a theorem derived by Widom, which states that g(r) becomes proportional to exp [−βu(r)] as r → 0 [77].
The direct correlation functions c(r/d) along each isomorph are illustrated in the second panel for the interval r/d ≤ 5, see the subfigures 2(b),(e),(h),(k).It is evident that the direct correlation function is a strongly variant quantity everywhere.This could be expected from the reduced excess inverse isothermal compressibility relation μT = −n c(r)d 3 r and the fact that μT is variant as a thermodynamic quantity that involves second order volume derivatives [57] as well as from the exact asymptotic limit c(r → ∞) = −βu(r) [1].It is worth noting that direct correlation functions reach their asymptotic limit much faster than other static correlation functions.This takes place prior to r/d = 2, close to the foot of the c(r) curve where the slope exhibits a rapid change.This observation justifies the satisfactory performance of the soft mean spherical approximation for the YOCP [78,79].
The bridge functions B(r/d) along each isomorph are illustrated in the third panel for the interval 1.5 ≤ r/d ≤ 5, see subfigures 2 (c), (f), (i), (l).The insets feature a magnification of the oscillatory bridge function pattern.It is evident that the bridge function is a strongly invariant quantity along any isomorph curve in the whole range across which it can be accurately computed with the OZ inversion method.In contrast to the radial distribution function for which there are strong variant fea- tures within the correlation void and for which longer range weak variant features are concentrated in the first maximum vicinity, the bridge function variant features appear to be uniformly spread in the whole computation range.In the long range, the observed invariance is justified by the asymptotic behavior of the bridge function.
Substituting for c(r) ≃ −βu(r) and g(r) = h(r) + 1 in the closure equation, we have Taylor expanding the logarithm with respect to h(r) ≃ 0 and retaining up to the second order term, we end up with B(r) = −(1/2)h 2 (r).In the intermediate range, the observed invariance is in accordance with the h-bond expansion that formally defines the bridge function through the infinite series B(r) = ∞ i=2 b i (r)n i where the un-known weighting functions b i (r) are given by multiple integrals only involving the total correlation function [2].Since the integration space for the weighting function of the n i term is R 3i , the introduction of reduced units r/d or r/∆ means that all the powers of the density vanish for each term, thus leading to B(r/d) = ∞ i=2 b i (r/d).To sum up, for all YOCP isomorphic lines, a high level of bridge function invariance has been observed.However, the small deviations between the bridge functions of YOCP state points that belong to the same isomorph might be comparable to the omnipresent uncertainties in bridge function determination.Therefore, in order to accurately quantify the level of isomorph invariance of the bridge function in the intermediate and long ranges, a detailed analysis of all different uncertainty sources is required.This is pursued in sections V and VI.

V. SENSITIVITY STUDIES
Bridge function uncertainties originate from uncertainties in the simulation-extracted radial distribution functions that propagate through all the steps of the OZ inversion method.For distances far from the edge of the correlation void that can be loosely defined by g(r) ≪ 1, the induced direct correlation function uncertainties are dominant.On the other hand, in the vicinity of the correlation void, the radial distribution function uncertainties themselves are the most crucial.Only the former regime is of relevance for the present work, since the OZ inversion method fails within the correlation void.
In modern computer simulations, where large numbers of particles and long observation times are routinely feasible, errors in extracted radial distribution functions are negligible.However, as discussed earlier, it is well-known that bridge functions are highly sensitive to radial distribution function uncertainties.Rigorous quantification of the sensitivity degree is complicated due to the non-linear non-local nature of the exact functional that connects the two quantities.Useful insights can be gained by inserting small controlled perturbations in the extracted g(r) and then documenting their effect in the computed B(r).From this numerical investigation, the following conclusions can be drawn: (a) Regardless of the state point and for constant error magnitude ǫ, the bridge function is mainly sensitive to aperiodic uncertainties in the radial distribution function.(b) Regardless of state point and for constant error magnitude ǫ, as the wavelength λ ǫ of periodic uncertainties in the radial distribution function decreases, the sensitivity of the bridge function becomes progressively weaker.(c) As κ increases for constant nor-malized coupling parameters Γ/Γ m , the sensitivity of the bridge function to periodic and aperiodic uncertainties dramatically decreases.Thus, bridge function uncertainties are larger near the OCP limit.(d) As Γ/Γ m increases for a constant screening parameter κ, the sensitivity of the bridge function to periodic and aperiodic uncertainties increases.Hence, bridge function uncertainties are larger near the melting line.(e) The study of aperiodic errors for the (Γ OCP ISO = 160, κ = 1.0)YOCP state point corresponds to the worst case scenario, where ǫ < 5×10 −9 was required for bridge functions to be indistinguishable.
The artificial error studies demonstrate the well-known fact that very accurate simulation-extracted radial distribution functions should be available for reliable computation of the bridge function.It is evident though that the desired level of accuracy depends strongly on the YOCP state point.The aforementioned trends will be invoked in the discussion of statistical errors and their variations with the state point, in the analysis of grid errors and the existence of a near-optimal bin width, in the investigation of explicit finite size errors and the necessity for their correction, as well as in the rationale behind the discarding of implicit finite size errors.

VI. PROPAGATION OF UNCERTAINTIES
Five types of uncertainties are relevant to investigations of the isomorph invariance of bridge functions: The first four uncertainties refer to the propagation of radial distribution function uncertainties and are relevant for all bridge function studies.The last uncertainty refers to the propagation of isomorphic state point uncertainties and is only relevant for bridge function studies which test the degree of invariance along isomorphic curves.
In the rich literature of bridge function extraction by computer simulations, finite size errors and tail errors are usually discussed, since rigorous procedures for their correction have been developed.On the other hand, statistical errors and grid errors are rarely analyzed in depth.Notable exceptions to this norm are the detailed uncertainty analysis performed by Poll and collaborators for the OCP bridge functions [16] as well as by Kolafa and collaborators for the hard sphere bridge functions [8].In the present investigation, a detailed quantification of uncertainties is crucial in order to correctly attribute the physical origin of the small deviations which are observed between the bridge functions of state points that belong to the same YOCP isomorph.To our knowledge, the following uncertainty study is the most meticulous analysis to be reported in the literature thus far.

A. Statistical errors
Statistical errors emerge in the extraction of any thermodynamic, structural or dynamic property from computer simulations as a consequence of their finite duration [80].For ergodic systems, the thermodynamic ensemble average of physical quantities can be substituted by their long time average.However, the latter can only be practically approximated in a truncated series form, i.e.
with f N (r N , p N ) the N −particle distribution function, a an arbitrary physical quantity, ∆t the simulation time step, K the total number of time steps in the equilibrium state.Due to the inherent fluctuations of a(i∆t) at each time step and the finite value of K, statistical uncer-tainties arise in the average value of a.In the following, we shall denote thermodynamic ensemble averages with ... and simulation averages with ... K where K is the number of samples.

Level of radial distribution functions
In order to calculate the statistical uncertainties in the radial distribution function, the number of time steps that is required for radial distribution functions to become uncorrelated at all distances should be first identified.This was accomplished with the application of the block averaging method proposed by Flyvbjerg and Petersen [81] at all distances and led to the conclusion that 64 MD time steps are required, regardless of the YOCP state point.Adding a safety margin, after the equilibration phase, the configurations were saved every 128 time steps for the extraction of g(r) in the production runs.It is evident that any larger period of configuration saving would unnec- essarily increase the computational cost.Given the lack of time-correlations, well-known statistical formulas can be employed for the average radial distribution function and its standard deviation, with M the total number of saved (uncorrelated) configurations (M = 65536) and where σ denotes the standard deviation of the average or standard error of the mean.The uncertainty σ[g(r)] depends on the simulation duration through the sample size M and the particle number through the inherent g(r, i∆t) fluctuation level, since for a single configuration the total g(r, i∆t) statistics are equal to the N (N − 1) number of pairs.Figure 5 illustrates these dependencies for the relative standard error σ[g(r)]/g(r).Notice that the relative standard error exhibits a divergent behavior at short distances, since close particle encounters are extremely rare.In particular, distances within the correlation void are so poorly sampled that the error in g(r) becomes comparable to the g(r) average value prior to one inter-particle distance d.Our production runs were designed in a manner that guarantees σ[g(r)]/g(r) < 0.01 up to r = 1.25d.

Level of bridge functions
In principle, knowledge of σ[g(r)] should allow for a determination of σ[B(r)] by standard error propagation anal-ysis.However, the non-locality of the OZ equation prohibits such calculations, since σ[c(r)] depends on the g(r) values at all possible distances that are naturally correlated with each other.Nevertheless, the application of error analysis to the closure equation alone can be used in order to reaffirm the inapplicability of OZ inversion at short distances.Substituting for g(r) → g(r) + σ In view of the calculation of the g(r) statistical errors, it would seem appropriate to compute the bridge function B(r, i∆t) for all extracted g(r, i∆t) and then to employ These expressions would lead to erroneous results both for the average bridge function and for its standard deviation.The reason is that the bridge function and direct correlation function lack microscopic representation, in contrast to the radial distribution function and structure factor.These static correlation functions are only defined in the thermodynamic limit and their reliable computation from OZ inversion requires a sufficiently smooth radial distribution function as input.Thus, strictly speaking, the quantity B(r, i∆t) does not make physical sense and should not be employed in the calculation of the average bridge function and its standard deviation.In fact, since the bridge function is a highly non-linear functional of the radial distribution function, there are large deviations between B[g(r, i∆t)] M and B[ g(r, i∆t) M ] with the second relation representing the correct way to calculate the average bridge function.
Based on these observations, the following procedure was devised to estimate the propagation of statistical uncertainties to bridge functions.First, the total dataset of M uncorrelated configurations was divided into N b blocks each containing N g configurations (M = N b × N g ).The N g blocking ensured the extraction of sufficiently smooth radial distribution functions and the N b blocking ensured large sample sizes for the calculation of the statistical deviations.Then, for each block, the block-averaged radial distribution function and bridge function were computed: Finally, the (N b , N g ) values should be determined.The N g value should be as large as possible for g i (r) Ng to be sufficiently smooth, otherwise B i (r) would be unphysical.
The N b value should be as large as possible for the B i (r) sample size to be sufficiently large, otherwise the estimate of σ[B(r)] would be unreliable.Within the constraint of M = N b × N g and assuming that the above requirements are of equal significance, then the near-optimal values are The above choice for the N b , N g combination was also confirmed to be near-optimal by the objective empirical analysis described below.For all 16 YOCP state points of interest, the statistical uncertainties in the bridge function were determined for 11 symmetric (N b , N g ) combinations.The statistical error was observed to strongly fluctuate for combinations that featured low values of N b , N g < ∼ 64, but remained nearly constant for the combinations (N b , N g ) = {(512, 128), (256, 256), (128, 512)}, see figure 6 for an example.Since (N b , N g ) = (256, 256) lies at the centre of the stability neighbourhood regardless of the state point, it was selected for the quantification of statistical errors.It should be emphasized that the near-optimal (N b , N g ) combination not only depends on the number of uncorrelated configurations M but also on the number of simulated particles N .Hence, the above choice does not constitute a general recommendation.
Analysis of the statistical errors for all 16 state points that belong to the 4 YOCP isomorphs led to the following conclusions, all in accordance with the sensitivity studies.

B. Grid errors
Grid errors emerge as a consequence of the finiteness of the bin width that is employed in the histogram method extraction of radial distribution functions from computer simulations.In other words, grid errors in the radial dis-tribution function are generated by the fact that only average values in intervals of finite length can be extracted from MD simulations.In mathematical terms, grid errors originate from the discretization of the microscopic representation of the radial distribution function.It is instructive to revisit the histogram method starting from the g(r) δ−function representation [16] This expression is integrated within thin spherical shells (r n , r n + ∆r) with ∆r ≪ r n .The integral of the double series is obtained from the post-processing of the MD trajectories and denoted with N (r n , r n + ∆r).The integral of ng(r) is evaluated by employing spherical coordinates, expanding the radial distribution function around the effective bin position x n and retaining up to the first order term.Introducing ∆V = 4 3 π[(r n + ∆r) 3 − r 3 n ], we get To close the system, we require that the first order term identically vanishes.This results in the conventional histogram method relations that are exact for infinitesimal bin widths ∆r → 0 and lead to g(r) errors due to the neglected high-order terms.This generates numerical errors in the calculation of the direct correlation function and thus in the computation of the bridge function.There is clear trade-off between grid errors and statistical errors; smaller bin widths reduce the grid errors but increase statistical errors (certainty in the distances but strong fluctuations in average values), while larger bin widths enhance grid errors but decrease statistical errors (weak fluctuations in the average values due to better grid statistics but uncertainty in distances).Grid errors should not be confused with numerical errors that stem from the use of discrete Fourier transforms instead of continuous Fourier transforms during the OZ inversion procedure.Within the heuristic assumption of a deterministic radial distribution function (independent of the bin width), Fast Fourier transform routines lead to negligible errors in the YOCP bridge function already for any discretization ∆r ≤ 0.1d.Grid errors stem from the fact that the histogram-extracted g(r) exhibit a nonnegligible dependence on bin width variations, i.e. it has a stochastic character.
In contrast to statistical errors, grid errors cannot be quantified.Hence, an efficient strategy should be focused on ensuring that sufficiently small bin widths are used so that grid errors are much smaller than statistical errors.The exact value should be decided by an empirical analysis of the dependence of the bridge function on the bin width [8].Our investigation revealed that as bin widths become smaller, the bridge functions start to become independent of their value.A near optimal value exists that is close to the largest bin width for which convergence has been achieved, since further reduction would only slightly alter the average bridge function but strongly increase its standard deviation (due to statistical errors).
Specifically, for each of the 16 YOCP state points of interest, the bridge functions were computed from radial distribution functions that were extracted from the same MD simulations but with varying histogram bin widths.Eight bin width values were considered, namely ∆r/d = {0.0005,0.001, 0.002, 0.004, 0.006, 0.008, 0.01, 0.04}.A characteristic example is illustrated in figure 7. Results reveal that: (a) Regardless of the YOCP state point, the bridge functions are always nearly identical for ∆r ≤ 0.002d bin widths.(b) As Γ/Γ m increases with a constant screening parameter κ, the bridge function depends more strongly on the bin width.(c) As κ increases with a constant normalized coupling parameter Γ/Γ m , the bridge function dependence on the bin width becomes weaker.(d) The near-optimal bin width depends on the YOCP state point.It lies within the interval ∆r/d = 0.002−0.01for the probed state points and acquires its smallest value when Γ OCP ISO = 160, κ = 1.0.Nevertheless, a state point independent bin width of ∆r = 0.002d was preferred.It should be emphasized that the near-optimal bin width should depend strongly on the number of statistics relevant to the extraction of the radial distribution function, i.e. on the number of particles and the number of uncorrelated time steps in the MD simulation.Thus, the above choice does not constitute a general recommendation.
The qualitative aspects of the above conclusions were anticipated from the sensitivity study.For instance, let us consider the existence of a near-optimal bin width in greater detail.The grid error can be roughly viewed as an oscillating correction to the radial distribution function with a periodicity of λ ǫ ≃ ∆r.As ∆r gets smaller, then the error ǫ in the average g(r) becomes larger.As λ ǫ gets smaller, then the bridge function becomes gradually insensitive to the error ǫ.Consequently, a near-optimal bin width exists where convergence is achieved with these two competing effects cancelling each other and the average bridge function becoming insensitive to the bin width.

C. Finite size errors
Finite size errors ultimately emerge from the finite number of particles considered in computer simulations [82].In the case of static equilibrium correlations, two types have been identified in the literature [83][84][85]; the explicit (or ensemble) size errors that emerge in the passage from the canonical ensemble to the grand-canonical ensemble . Determination for 8 different histogram bin widths ∆r/d = {0.0005,0.001, 0.002, 0.004, 0.006, 0.008, 0.01, 0.04} It is evident that a near-optimal bin width emerges below which bridge functions overlap.The near-optimal bin width depends slightly on the screening parameter and obtains the rough values 0.002d (κ = 1.0), 0.002d (κ = 1.5), 0.004d (κ = 2.0), 0.006d (κ = 2.5).Unless the bin widths are very large (∆r/d > ∼ 0.04), grid errors are rather small and better discerned in the zoomed-in insets.
and the implicit (or anomalous) size errors that emerge in the passage from a finite simulated system with infinite boundary conditions to an infinite macroscopic system.

Explicit finite size errors
The explicit size effect in NVT MD simulations is a direct consequence of the suppression of particle number fluctuations in the canonical ensemble.The correspondence between the NVT and µVT radial distribution functions can be derived by expressing the grand canonical twoparticle densities via the canonical two-particle densities, Taylor expanding the canonical quantity with respect to the average particle number, retaining up to the second order term and employing the fluctuation relation for the isothermal compressibility.The final expression is the socalled Lebowitz-Percus correction and reads as [86][87][88] where g c (r; n, T ) is the corrected radial distribution function, g MD (r; n, T ) the (NVT) MD-extracted radial distribution function, χ T the reduced isothermal compressibility and N the particle number.After switching from the (n, T ) to the (Γ, κ) state variables and expanding the derivatives, the Lebowitz-Percus correction for the YOCP becomes The final term involves all possible first and second order partial derivatives of the radial distribution function with respect to the state variables.In principle, these can be evaluated with simulations by employing finite difference approximations.Nevertheless, the standard second-order central-difference scheme would require nine MD simulations for any state point which correspond to all the first neighbouring grid points of a two-dimensional computational stencil centered at the state point of interest.To avoid this formidable task, the relative importance of the final term has been assessed with the aid of IEMHNCgenerated radial distribution functions.For all 16 state points of interest, the results revealed that the final term has no impact on the computed bridge functions.Hence, a simplified version of the Lebowitz-Percus correction can be safely employed that reads as [16,88] We note that discarding of the final term is equivalent to stating that the radial distribution function is relatively insensitive to very small density variations; a reasonable assumption above the Fisher-Widom line and especially close to crystallization.The only obstacle remaining in the implementation of the correction for the explicit size effect has to do with the calculation of the isothermal compressibility.In absence of an established YOCP equation of state that remains very accurate in the entire liquid domain, two alternative paths were followed.(a) Integral equation theory.Given the well-documented successes of the IEMHNC approach for YOCP liquids, IEMHNC-produced radial distribution functions were employed.The IEMHNC approximation was solved in very dense extended grids of ∆r/d = 0.0001 and max (r/d) = 100.The virial route as well as the statistical route were considered.In the virial route to χ T , the first derivatives of the pressure with respect to (Γ, κ) were computed with the second-order central difference scheme and ∆Γ = ∆κ = 0.001 finite steps.In the statistical route to χ T , the integral relation that contains the direct correlation function was preferred and the asymptotic limit was added and then subtracted to reduce the truncation errors [26].(b) Hypervirial route.This path is only realizable in MD simulations, since the ensemble averaged hypervirial is a non-thermodynamic quantity.The hypervirial relation reads as [89,90] (δW In the above; the dimensionless virial function is defined by W = −(1/3) i j>i βw(r ij ) with the interatomic virial given by w(r) = rdu/dr, the dimensionless hypervirial function is defined by X = (1/9) i j>i βx(r ij ) with the interatomic hypervirial given by x(r) = rdw/dr, the operator ... denotes canonical ensemble averaging and δW = W − W for the canonical deviations from mean.The MD implementation only requires the recording of the virial and hypervirial quantities at each time step.For all the 16 YOCP state points of interest, the χ T values stemming from the three routes are reported in Table II.The MD-hypervirial route will be preferred, but it is worth pointing out that the IEMHNC-virial route provides very accurate results while the IEMHNC-statistical route typically leads to slight underestimations.Given its aperiodic form, see Eq.( 5), the explicit finite size error can be expected to weakly affect the computed bridge functions despite the large simulated particle number (N = 54872).The magnitude of this size error is controlled by the isothermal compressibility χ T , whose thermodynamic state dependence is opposite to that of the bridge function sensitivity.In particular: (a) As Γ/Γ m increases for constant screening parameter κ, χ T decreases but the sensitivity increases.(b) As κ increases for a constant normalized coupling parameter Γ/Γ m , χ T increases but the sensitivity decreases.In both cases, the sensitivity variations are far more dramatic than the compressibility variations, which implies that the bridge functions of the state points that lie closer to the crystallization line and closer to the OCP limit should be more susceptible to explicit finite size errors.
The numerical results confirmed our theoretical expectations.The explicit finite size errors only influence the bridge functions along the isomorph line that is closer to the melting line and mainly affect the κ = 1 member.The effect is minor and mostly manifested as a correction of the asymptotic behavior of the bridge function which now more properly converges to zero.A characteristic example is illustrated in figure 8.

Implicit finite size errors
The implicit size effect in computer simulations is a direct consequence of the spurious correlations which are introduced by imposing periodic boundary conditions [91][92][93].In simple liquids that are simulated with standard cubic boxes, the infinite system should be viewed as a primitive cubic crystal which results from periodic repetition of the complex composite unit cell in all directions [92].Such a picture demonstrates that the imposed symmetry of the periodic boundary is equivalent to a progressively weaker rigid bond between each unit cell particle and its infinite periodic images, which induces orientational order in the liquid.This leads to slightly non-isotropic radial distribution functions and to systematic errors in the spherically averaged radial distribution functions [94].
Within the grand-canonical ensemble and with the aid of cluster expansion techniques, Pratt and Haan devised a formally exact theory for the implicit finite size effect under the assumption that particles do not directly interact with any of their periodic images [92], which is valid for pair interaction potentials that are truncated within the unit cell.When neglecting a class of bridged graphs, an approximate expression was derived that connects the slightly anisotropic (explicit effect corrected) g c (r) of the simulated system with the actual g(r) of the bulk system.T , the absolute relative deviations between the hypervirial and virial results are denoted by ǫIEMHNC,v and the absolute relative deviations between the hypervirial and statistical results are denoted by ǫIEMHNC,s.The MD-hypervirial results are extremely accurate as inferred from the smallness of e MD T /µ MD T and have been preferred.The IEMHNC-virial results are also very accurate with relative deviations that are always less than 0.17%, while the IEMHNC-statistical results are less accurate with relative deviations ranging from 0.3% up to 8.2%.This is expected given the small degree of thermodynamic inconsistency of the IEMHNC approach [26,28].The expression has the superposition form and reads as where r 1 − r 2i refers to the displacement vector between the particle 1 and any of particle's 2 infinite periodic images [92].It is apparent that the above expression cannot be inverted with respect to the unknown bulk g(r).In contrast to the explicit size effects, the implicit size effects cannot be directly corrected.However, under some reasonable additional assumptions, it has been shown that implicit size errors for the spherically averaged radial distribution function g c (r) = (1/2π) g c (r)dΩ scale as .
On the other hand, see Eq.( 5), explicit size errors scale as g c (r)/g MD (r) ∼ 1 + O [1/N ].This implies that, for sufficiently short-ranged potentials ν > 3, implicit size errors should decrease much faster with the number of particles.Since in our simulations the explicit size effect was very weak, the implicit size effect should be expected to be insignificant for the κ ≥ 1 Yukawa potentials considered, which have shorter range than the ν = 3 IPL potential.
In addition, implicit size errors are not aperiodic like the explicit size errors, but display a λ ǫ ≃ 2d periodicity that is comparable to the distance between successive coordination cells [93].As demonstrated earlier, the sensitivity of the bridge function to such radial distribution function errors is much smaller.Hence, it can be safely concluded that implicit size errors have negligible effect on the computed bridge functions.

D. Tail errors
Tail errors emerge from the finite extent of the primary cell of the simulations.These errors are generated by the fact that radial distribution functions can only be reliably extracted inside a restricted domain with the maximum distance naturally equal to half the length of the cubic box.Unless MD-extracted radial distribution functions have truly reached their asymptotic limit of unity prior to r = L/2, their direct implementation will lead to errors in the direct correlation functions that will magnify while propagating to the bridge functions.Even for MD simulation domains of moderate size, tail errors can be nearly eliminated with a long-range extrapolation method.There are three common types of long-range extrapolation methods.chain closure c(r) = g(r)−1−βu(r)−ln [g(r)] [96,97] and mean spherical approximation closure c(r) = −βu(r) [98] are typically invoked.The extrapolated g(r) is obtained by numerically solving the OZ equation that is complemented with mixed closure conditions.(b) Baxter methods.These techniques are also based on integral equation theory and assume the long-range validity of a traditional closure equation, but they consider the equivalent Baxter system of equations that emerges from the Wiener-Hopf factorization of the OZ equation [99,100].(c) Asymptotic methods.Such techniques conjecture that the long-range (near-asymptotic) rh(r) can be approximated by a discrete sum of exponentially decaying oscillating functions.The solution of the OZ equation or the Baxter equations is circumvented.The unknown algebraic parameters are chosen to fit the simulation results in a subinterval of the intermediate range that extends up to r = L/2 [96,101].
The N = 54872 particles employed in the present MD simulations correspond to L/2 ≃ 30d, which was considered long enough for g(r) to be sufficiently close to unity so that the influence of tail errors on the bridge functions is completely negligible.This expectation was verified by employing the Verlet method and the asymptotic method for all the 16 YOCP state points of interest.
The application of the asymptotic method was based on the assumption that rh(r) can be approximated by a single exponentially decaying oscillatory term.The four unknown parameters {a 0 , a 1 , a 2 , a 3 } that emerge in the expression g fit (r) = 1+(1/r) exp (a 0 − a 1 r) cos (a 2 + a 3 r) were determined by least square fitting to the extracted g MD (r) in an interval r f,d < r < r f,u .The upper fitting limit was always retained at r f,u = L/2 ≃ 30d, while the lower fitting limit varied within r f,d = 15−10d depending on the state point.The negligible fitting errors supported the above functional assumption for g fit (r).The extrap-olated radial distribution function reads as where r m /d = 60 was selected in order to double the numerical extraction domain and where r x ∈ [r f,d , r f,u ] was selected so that the discontinuity between g MD (r), g fit (r) is minimized (r x /d = 13.6 − 28.3 depending on the state point).The bridge functions that result from g MD (r) up to r = L/2 ≃ 30d and also from g c (r) up to r = 60d were computed, corrected for finite size errors and compared.For all YOCP state points, the bridge functions totally overlapped, confirming that tail errors are negligible.The Verlet method was utilized in combination with the hyper-netted chain (HNC) closure.The OZ equation has been supplemented with the closure condition where we selected r m /d = 60 in order to double the numerical extraction domain and where we probed different r x /d = {5, 8, 10, 15, 20} values in order to investigate how limited the MD extraction domain should be for the tail errors to be impossible to correct by long range extrapolations.The system of the OZ integral equation and the non-linear closure condition has been solved with Picard iterations in Fourier space; the real space resolution was ∆r/d = 0.002 that is equal to the optimum bin width employed to extract g MD (r), the reciprocal space resolution was ∆q = πd/r m and at each iteration n convergence was assumed when the criterion ||c n (r)−c n−1 (r)|| < 10 −8 was satisfied for all distances.The bridge functions that result from g MD (r) up to r = L/2 ≃ 30d and also from g c (r) up to r = 60d were computed, corrected for finite size errors and compared.For all the YOCP state points, the bridge functions totally overlapped at least for r x /d ≥ 8 which demonstrates that tail errors are negligible for our simulations where r x /d = 30.Even for r x /d = 5, the deviations were very small along the Γ OCP ISO = 160 isomorph and remained negligible along the Γ OCP ISO = 40 isomorph, see figure 9 for an illustration.The results suggest that the radial distribution functions did not need to be extracted up to 30d; an extraction up to 8d combined with a reliable extrapolation method would have led to identical bridge functions with far less computational cost.

E. Isomorphic errors
Isomorphic state point uncertainties emerge due to the fact that isomorphic curves are not traced out exactly.To be specific, the direct isomorph check is based on an approximate expression and it requires input from NVT MD simulations, whereas the small step method is based on an exact differential equation that is solved numerically and it also requires input from NVT MD simulations.As isomorphic curves are being traced out, both methods involve a controlled density re-scaling.As a consequence, densities are known exactly and only temperature uncertainties exist.For the YOCP, this implies that there are only uncertainties in the coupling parameters of the isomorphic state points.The isomorphic errors emerge from the propagation of these uncertainties to the bridge function.There is a direct propagation through the presence of the pair interaction potential and an indirect propagation through the presence of the radial distribution and direct correlation functions in the closure equation.
The isomorphic coupling parameter uncertainties have been quantified from analysis of the local truncation errors and the global accumulated errors of the fourth order Runge-Kutta method that also accounted for statistical errors in MD-extracted density scaling exponents.
Isomorph tracings based on the small step method were carried out with different logarithmic density increments (|∆n|/n = 2.3%, 4.8%, 8.8%) and particle numbers (N = 8192, 17576) in order to verify the output of this analysis.The empirical deviations between the isomorphs that stem from the direct isomorph check and the small step method were also considered.The emerging upper uncertainty thresholds varied within ∆Γ/Γ < 0.5%−1.0%with the exact values depending on the screening parameter.
Isomorphic errors can be quantified by extracting the radial distribution function g(r; Γ±∆Γ, κ) from MD simulations, computing the bridge function B(r; Γ ± ∆Γ, κ) and comparing with B(r; Γ, κ).It is preferable that grid errors are minimized as well as that finite size errors and tail errors are corrected in both bridge functions prior to comparison.This was pursued for the four state points corresponding to the κ = 2.5 members of all YOCP isomorphs.These state points are subject to the largest isomorphic coupling parameter uncertainties ∆Γ/Γ < 0.01, since the sequential tracing of the isomorphs with the direct isomorph check and the small step method initiates at κ = 1.0.Comparison revealed that isomorphic errors are negligible, see figure 10.In particular, the two bridge functions were indistinguishable for Γ OCP ISO = 40, 80, 120 and very minor deviations were observed for Γ OCP ISO = 160 that were much smaller than the statistical errors.
A similar procedure was followed in an attempt to answer the closely-related question of how much isomorphic coupling parameter uncertainties would be required in order to make the bridge functions of different isomorphic state points fully invariant.This problem was pursued for the Γ OCP ISO = 160 isomorph in the following manner.The coupling parameter of the κ = 2.5 member was reduced by 2.5%, 5.0%, 7.5%, 10%, the respective radial distribution functions were extracted from MD simulations and the respective bridge functions were computed.The latter were then compared with the bridge functions of the κ = 1.0, 1.5, 2.0 members of the same isomorph.Comparison revealed that Γ uncertainties exceeding 10% are required in order to make the κ = 2.5 and κ = 1.0 bridge functions invariant, see figure 11.

VII. CORRECTED BRIDGE FUNCTIONS INCLUDING UNCERTAINTIES
As demonstrated in section VI for the MD simulation parameters of the production runs; the histogram bin width is selected in a manner which guarantees that grid errors are much smaller than statistical errors, the nature of the Yukawa pair potential together with the large number of simulated particles ensure that implicit finite size errors are negligible, the extraction domain is extended enough so that tail errors are negligible, the accuracy of isomorph tracing methods is high enough so that isomorphic errors are negligible.Hence, in order to incorporate the uncertainties, it suffices to correct for explicit finite size errors and to include error bars that stem from statistical errors.The finite-size corrected bridge functions including statistical uncertainties are computed in the following manner for all the 16 YOCP state points of interest.First, the reduced isothermal compressibility is calculated from the production runs with the hypervirial route of Eq.( 6).The simplified version of the Lebowitz-Percus correction, see Eq. (5), is now applied to the average radial distribution function g(r) M and all block-averaged radial distribution functions g i (r) Ng .Then, the average bridge function B(r) = B[ g c (r) M ] and all the block-averaged bridge functions B i (r) = B[ g c,i (r) Ng ] are computed with the OZ inversion method from the respective sizecorrected radial distribution functions.Afterwards, the standard deviation of the average bridge function σ[B(r)] is computed with the block averaging method outlined in section VI A 2 by assuming the near optimal combination of (N b , N g ) = (256, 256).Finally , the standard error of the mean σ[B(r)] is utilized for the determination of confidence intervals.In particular, the selected error bars for the statistical uncertainties correspond to 95% confidence intervals.
The finite-size corrected bridge functions including error bars are illustrated in figure 12 for the 4 isomorphic curves and the 16 YOCP state points of interest.It is evident that the errors cannot account for the small deviations that are observed between the bridge functions of the different members of the same isomorph.Therefore, the observed isomorph invariance of bridge functions in the long and intermediate range r ≥ 1.5d is only approximate.Note that the relative invariance holds to nearly the same degree regardless of the YOCP isomorph.Note also that the approximate bridge function invariance extends up to the edges of the correlation void, where the approximate radial distribution function invariance begins to break down.Finally, it is worth pointing out that, for all the YOCP state points investigated, the bridge function becomes slightly positive close to r ≃ 2d well within the statistical uncertainties.In rough accordance with the asymptotic behavior of B(r) = −(1/2)h 2 (r) and with the ∼ 1.5d periodicity of the total correlation function, additional shallower positive maxima appear with a periodicity slightly less than ∼ 0.8d.The emergence of sign switching and of a positive maximum within the first coordination cell seems to be a rather standard feature of the bridge functions of dense fluids that has also been observed for hard-sphere systems [8], Lennard-Jones liquids [10,11], IPL-12 systems [14] and OCP liquids [17].It is a salient feature of bridge functions that cannot be captured by the VMHNC approximation [46] that utilizes the non-positive analytical Percus-Yevick hard-sphere bridge function.This deficiency has been suggested as responsible for the minor structural inaccuracies of the VMHNC approach that are observed in the vicinity of the first peak of the radial distribution function [28].

VIII. SUMMARY AND CONCLUSIONS
Yukawa bridge functions were systematically computed aiming to confirm or disprove the validity of the ansatz of reduced unit bridge function invariance along isomorphs, i.e. phase diagram lines of constant excess entropy.16 state points were selected that belong to four isomorphs and cover the entire dense liquid YOCP phase diagram.The YOCP isomorphic curves were traced out with the small step method as well as the direct isomorph check.The intermediate and long bridge function ranges were made accessible after application of the Ornstein-Zernike inversion method with radial distribution function input from ultra-accurate molecular dynamics simulations that employed carefully selected parameters.
In order to accurately quantify the level of isomorph invariance of the bridge functions in the reliable extraction range, meticulous analysis of all sources of error was carried out.A detailed investigation of the sensitivity of bridge functions to aperiodic and periodic multiplicative perturbations in radial distribution functions led to important insights concerning the propagation of uncertainties.Regardless of the state point, it was consistently observed that YOCP bridge functions are far more sensitive to aperiodic perturbations.A strong phase diagram dependence was also discerned with the YOCP state points that lie in the vicinity of the melting line or near the OCP limit possessing more sensitive bridge functions.The use of these controlled artificial errors facilitated the understanding of all types of naturally emerging errors.
The statistical errors were quantified with a block averaging procedure that ultimately revealed a near-optimal combination of the block number and the sub-block configuration number after exhaustive trials.The grid errors were minimized compared to the omnipresent statistical errors due to the utilization of a near-optimal histogram bin width that emerged after comprehensive testing.The explicit finite size errors were observed to influence the bridge function asymptotes and were corrected with the simplified version of the Lebowitz-Percus expression valid for radial distribution functions that are relatively insensitive to small density variations, while the implicit finite size errors were deduced to be negligible after inspecting an approximate scaling derived from the exact Pratt and Haan rigid bond theory.The tail errors were confirmed to be negligible by the application of the Verlet and asymptotic long range extrapolation methods.The isomorphic errors that emerge from slight state point excess entropy mismatches due to inaccuracies in the isomorph tracing techniques were also shown to be negligible.
The final YOCP bridge functions, corrected for explicit finite size errors and featuring error bars stemming from the statistical uncertainties, were observed to be nearly isomorph invariant in the intermediate and the long range for all four excess entropies probed.This invariance was concluded to be approximate, since the small deviations observed between isomorphic bridge functions always exceed the quantified level of uncertainties.However, the bridge functions remain nearly isomorph invariant even at the edge of the correlation void, where the radial distribution functions and potentials of mean force already exhibit strong variance.The isomorph invariance of YOCP bridge functions well within the correlation void up to the origin is investigated in the accompanying paper [102].

FIG. 1 :
FIG.1:The 4 targeted isomorphs Γ OCP ISO = (160, 120, 80, 40) together with the approximate melting line Γ OCP ISO = 171.8for dense YOCP liquids in the log Γ − κ phase diagram.The four isomorphic curves as obtained from the direct isomorph check (discrete symbols) are compared to those obtained from the approximate analytical expression (solid lines).The four numerical isomorphs begin to overshoot the respective analytical isomorphs roughly above κ ≃ 1.5, see also TableI.
for the residue.FFT algorithms are again employed.(5) The bridge function is computed.The closure equation is solved for the bridge function leading to B(r) = ln[g(r)] − h(r) + c(r) + βu(r).

1 .
statistical errors due to the finite simulation duration, 2. grid errors due to the finite histogram bin width, 3. size errors due to the finite simulated particle number, 4. tail errors due to the finite simulation box length, 5. isomorphic errors due to excess entropy mismatches.

ϵ = 6 × 0 FIG. 3 :
FIG. 3: Sensitivity of the YOCP bridge function to artificial aperiodic uncertainties in the radial distribution function at two state points along the isomorph Γ OCP ISO = 160.(a) For κ = 1.0; error magnitudes ǫ > 2 × 10 −7 lead to large perturbations in the bridge function, error magnitudes ǫ > 10 −8 lead to small but observable perturbations and error magnitudes ǫ < 5 × 10 −9 are required for the bridge function to become insensitive.(b) For κ = 2.5; error magnitudes ǫ < 6 × 10 −7 suffice for the bridge function to become insensitive to perturbations.There is a dramatic dependence of the sensitivity on the screening parameter.

7 = 10 FIG. 4 :
FIG. 4: Sensitivity of the YOCP bridge function to artificial periodic uncertainties in the radial distribution function at the state point Γ OCP ISO = 160, κ = 1.0 for different wavelengths.(a) For λǫ = d, error magnitudes ǫ < 10 −7 lead to very small bridge function perturbations.(b) For λǫ = 0.1d, error magnitudes ǫ < 10 −5 lead to very small bridge function perturbations.There is a strong dependence of the sensitivity on the error wavelength.

FIG. 5 :
FIG. 5: Relative standard errors in the MD-extraction of radial distribution function at the YOCP state point Γ OCP ISO = 160, κ = 1.0.Results for varying particle number (N ) and uncorrelated configuration number (M ).Here, rcut = 16d and Tsave = 64∆t was employed for the potential cut-off and configuration saving period instead of the standard rcut = 10d, Tsave = 128∆t of the production runs.See the inset for the insufficient statistics within the correlation void.
[g(r)], B(r) → B(r) + σ[B(r)], c(r) → c(r) + σ[c(r)], linearizing with respect to the small uncertainties and disposing the averages via the unperturbed closure equation, we obtain σ[B(r)] = σ[g(r)] g(r) − σ[g(r)] + σ[c(r)] .The first term represents the relative standard error in the radial distribution function and stems from the logarithmic term of the closure equation.As illustrated in figure 5, this term rapidly diverges at short distances where it should overcome contributions from other terms, thus demonstrating that B(r) cannot be reliably computed inside the correlation void with the OZ inversion method.
(a) As Γ/Γ m increases for a constant screening parameter κ, the statistical errors in the bridge function increase.(b) As κ increases for a constant normalized coupling parameter Γ/Γ m , the statistical errors in the bridge function decrease.(c) Regardless of the state point, the statistical errors in the bridge function are larger than the statistical errors in the radial distribution function.The largest difference is achieved for Γ OCP ISO = 160, κ = 1.0 (∼ 50×) and the smallest difference for Γ OCP ISO = 40, κ = 2.5 (∼ 3×).

FIG. 8 :
FIG. 8: Finite size errors in the computation of the YOCP bridge function at the state point Γ OCP ISO = 160, κ = 1.0.Results with and without the Lebowitz-Percus correction for the explicit finite size error.

FIG. 9 :
FIG. 9: Tail errors in the computation of the YOCP bridge function at the state point Γ OCP ISO = 80, κ = 1.0.Results for the extended radial distribution function extraction up to 30d without compensation for the tail errors and two limited radial distribution function extractions (rx/d = 5, 8) combined with the Verlet HNC extrapolation method.Tail errors can be totally compensated for g(r) extractions at least up to 8d, as seen by the indistinguishable respective bridge functions.

5 FIG. 12 :
FIG. 12: Finite-size corrected bridge functions featuring error bars due to the statistical uncertainties (95% confidence intervals) resulting from the application of the Ornstein-Zernike inversion method for dense YOCP liquids with radial distribution function input from NVT MD simulations (see the production runs).(a, b, c, d) Intermediate and long range results for the 4 members of the Γ OCP ISO = 160, 120, 80, 40 isomorphs, respectively.The isomorphic deviations and the error bars are rather small and can be better discerned in the zoomed-in insets.The same applies for the quasi-periodic sign switching of the bridge functions.

TABLE II :
The reduced inverse isothermal compressibility µT = 1/χT for all the 16 YOCP state points of interest.Results from the MD-hypervirial route (superscript MD), the IEMHNC-virial route (superscript IEMHNC,v) and the IEMHNC-statistical route (superscript IEMHNC,s).The standard deviations of the mean values of µ MD T are denoted by e MD