Continuous demixing transition of binary liquids: finite-size scaling from the analysis of sub-systems
CContinuous demixing transition of binary liquids: finite-size scaling from theanalysis of sub-systems
Y. Pathania, D. Chakraborty, a) and F. Höfling
2, 3, b) Indian Institute of Science Education and Research Mohali, Knowledge City, Sector 81, S. A. S. Nagar, Manauli-140306,India Freie Universität Berlin, Fachbereich Mathematik und Informatik, Arnimallee 6, 14195 Berlin,Germany Zuse Institute Berlin, Takustr. 7, 14195 Berlin, Germany (Dated: 1 October 2020)
A binary liquid near its consolute point exhibits critical fluctuations of the local composition; the diverging correlationlength has always challenged simulations. The method of choice for the calculation of critical points in the phase diagramis a scaling analysis of finite-size corrections, based on a sequence of widely different system sizes. Here, we discussan alternative using cubic sub-systems of one large simulation as facilitated by modern, massively parallel hardware.We exemplify the method for a symmetric binary liquid at critical composition and compare different routes to thecritical temperature: (1) fitting critical divergences encoded in the composition structure factor of the whole system,(2) testing data collapse and scaling of moments of the composition fluctuations in sub-volumes, and (3) applying thecumulant intersection criterion to the sub-systems. For the last route, two difficulties arise: sub-volumes are open systemswith free boundary conditions, for which no precise estimate of the critical Binder cumulant U c is available. Second,the periodic boundaries of the simulation box interfere with the sub-volumes, which we resolve by a two-parameterfinite-size scaling. The implied modification to the data analysis restores the common intersection point, and we estimate U c = . ± . I. INTRODUCTION
The calculation of phase diagrams of fluid substances is ofinterest for manifold applications. A particularly challengingtask is finding the loci of critical points, with the liquid–vapourcritical point as a prominent example. Binary liquid mixtures,in addition, exhibit phase separation and the coexistence ofdifferently composed liquid phases, which leads to a line ofcritical points in the phase diagram . In the vicinity of thesecritical points, the fluid is characterised by diverging lengthand time scales , which puts considerable difficulties on theirsimulation. The growth of the correlation length interferes withthe finite size of the simulation box, with the consequencesthat the observed phase transition is rounded and shifted andcritical divergences are capped.In the past decade, the seemingly never ending growth incomputing power was, among other factors, driven by a shiftto massive parallelisation, making molecular simulations of un-precedented system sizes and run lengths broadly available .The use of huge systems mitigates artifacts due to a finitesimulation box and allows, in principle, probing the criticaldivergences directly as one would do in an experiment. Yet,an elegant and conceptually superior alternative exploits thescale invariance of the critical fluid and turns the limitation ofa finite system size into an advantage by explicitly followingthe divergences of certain fluid properties as a function of thesystem size . The latter approach requires that the simula-tions at each thermodynamic state point are repeated for a widerange of system sizes, including comparably small systems.However, on prevalent massively parallel hardware such as a) Electronic mail: [email protected] b) Electronic mail: f.hoefl[email protected] high-performance graphics processing units (GPUs) it is notpossible to run a set of independent simulations concurrently,and so the conventional finite-size scaling does not participatein substantial advances of computing hardware.Both aspects, a systematic finite-size scaling and the efficientsimulation of huge systems can be combined in the finite-sizescaling analysis of sub-systems. This old idea was carried outsuccessfully for two- (2D) and three-dimensional (3D) Isinglattice models , but its application to 2D Lennard-Jones(LJ) fluids was limited by the computing resources available 30years ago and suffered from the insufficient separation of thetwo length scales . The approach was revived only recentlyto determine critical points in the phase diagrams of 3D activesuspensions ; in these studies, additional countermeasureswere needed to avoid biased sampling by sub-volumes thatcontained an interface.Compared to the conventional finite-size scaling, there aretwo important differences: boundary conditions are not pe-riodic, and the system geometry is inherently described bytwo length scales: the simulation box and the sub-system size.We will address the latter complication explicitly in a two-parameter scaling ansatz given at the end of the article. Inthe limit of small sub-systems relative to the simulation vol-ume, the sub-systems realise open systems in the sense ofstatistical mechanics; in particular, they can exchange heatand mass with a large reservoir. Open sub-systems are re-ceiving increasing interest in the simulation community, mo-tivated by studies of small-system thermodynamics butalso liquid–vapour interfaces and by the finite-size scaling ofKirkwood–Buff integrals to calculate chemical potentials .On the other hand, methodological advances in adaptive res-olution techniques permit the direct simulation of an opensystem coupled to one or many reservoirs acting as athermodynamic mean field. The free boundary conditions, as a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p realised by open sub-systems, modify the fluctuations of den-sity and composition . In the context of critical phenomena,the critical amplitudes are known to be different for free andperiodic boundary conditions , which applies also to thecritical Binder cumulant used to locate the critical temperature.In this work, we revisit the sub-system scaling analysis in thevicinity of a critical point and compare different approaches toestimate the critical temperature. We demonstrate the methodfor a symmetric binary liquid of Lennard-Jones particles, whichhas pair interactions that are symmetric with respect to the twomolecular species A and B. Thus, the critical composition isa 1:1 mixture that exhibits phase separation into symmetricA-rich and B-rich phases, occurring in close analogy to thecontinuous symmetry breaking of the Ising model. This modelsystem (and slight variants thereof) was studied extensivelynear consolute points using combinations of semi-grand canon-ical Monte Carlo (SGMC) methods and molecular dynamics(MD) simulations . In particular, cuts of the phase diagramin the temperature–composition and in the temperature–densityplanes are available, and the critical fluctuations of the localconcentration were shown to scale as expected for the 3D Isinguniversality class. Moreover, the coarsening kinetics and thecritical singularities of the transport coefficients were charac-terised in great detail, corroborating theoretical expectationsand being in agreement with the available experimental evi-dences.After giving technical details on the model liquid and thesimulations in Section II, we analyse in Section III the criticalbehaviour of the structure factor S cc ( k ) of local concentrationfluctuations obtained from the whole, large simulation box for arange of wavenumbers k . This is then contrasted in Section IVby a scaling analysis of the concentration susceptibility χ (cid:96) calculated from sub-volumes of different edge lengths (cid:96) . Bothapproaches yield already estimates of the critical temperature T c . Binder’s cumulant method based on sub-systems is carriedout in Section V to improve these values of T c and we showthe necessity for a scaling ansatz that accounts for two lengths,the size of the sub-system and that of the simulation box. II. MODEL SYSTEM AND SIMULATION DETAILS
The model of a symmetric binary mixture consideredhere employs pair interactions given by the truncated andforce-shifted Lennard–Jones potential u αβ ( r ) = u LJ; αβ ( r ) − u LJ; αβ ( r c ) − ( r − r c ) u ′ LJ; αβ ( r c ) (1)for r (cid:54) r c with u LJ; αβ ( r ) = (cid:15) αβ [︁ ( r /σ αβ ) − − ( r /σ αβ ) − ]︁ asusual for particle species’ α, β ∈ { A , B } . All particles have thesame diameter, σ αβ = : σ , and mass m , but different interactionstrengths, (cid:15) AA = (cid:15) BB = (cid:15) AB = : (cid:15) . The number density ρ ofthe fluid is fixed at ρσ − = N A = N B ), this binary fluid exhibits a continuous demixingtransition at the critical temperature T c ≈ . (cid:15)/ k B .A typical MD workflow consists of thermalisation and equi-libration runs, followed by one or several production runs. For temperature control, we used a Nosé–Hoover thermostat asdescribed in Ref. 39 for the initial equilibration in the canoni-cal ensemble ( N A , N B , V , T fixed). Then, the thermostat wasswitched off and the system evolved microcanonically, i.e., atfixed total energy, with a time step of δ t = . τ in unitsof τ = √︀ m σ /(cid:15) . The subsequent production run was per-formed microcanonically as well and particle configurationsand observables such as Fourier modes of the density fieldwere recorded every 10 integration steps. A single simula-tion run comprised of an initial thermalisation for 20 , τ ,followed by microcanonical equilibration and production runs(at fixed total energy E ) over 50 , τ and 30 , τ , respec-tively. Further, for each temperature, the data were averagedover two independent runs. The microcanonical ensemble wasused because we have also calculated transport coefficientsand time correlation functions (not discussed here) from thesame set of simulations; the intermediate equilibration step isneeded for proper temperature control in the NVE ensemble.For the cubic simulation box, periodic boundary conditionsand an edge length of L ≈ . σ were chosen, correspondingto a total of 87 ,
808 particles. The home-grown code takesadvantage of the massively parallel architecture of high-endgraphics processing units (GPU) , and simulations wererun on a GPU of type Tesla K20Xm (Nvidia Corp.) Dimen-sionless quantities are indicated by an asterisk and are formedwith σ, τ, (cid:15) as units of length, time, and energy; for example, ρ * : = ρσ and T * : = k B T /(cid:15) . III. CORRELATION LENGTH AND CONCENTRATIONSUSCEPTIBILITY
The structural properties of the system arise from two fluc-tuating fields: the density field δρ ( r ) = δρ A ( r ) + δρ B ( r ) and theconcentration field δ c ( r ) = x B δρ A ( r ) − x A δρ B ( r ), where δρ α ( r )denotes the fluctuations of the microscopic, partial density ofspecies α , x α = N α / N is the mole fraction of the species, and N = N A + N B . The relevant order parameter is ϕ = ⟨| x A − x c |⟩ ,where the critical composition is x c = / S cc ( k ) = N ⟨ | δ c ( k ) | ⟩ with δ c ( k ) : = ∫︀ V e i k · r δ c ( r ) d r , which far away from the criticalpoint is of the Ornstein–Zernike form for small wavenumber k , S cc ( k ) ≃ ρ k B T χ [︁ + ( k ξ ) ]︁ − ; k σ ≪ . (2)This asymptotic expression serves as definition for the corre-lation length ξ ( T ) and the concentration susceptibility χ ( T ).Upon approaching the consolute point, the two quantities ex-hibit the familiar algebraic divergences ξ ( T ) ≃ ξ ε − ν and χ ( T ) ≃ χ ε − γ ; T ↓ T c , (3)as functions of the reduced temperature ε : = | T − T c | / T c . Thesecritical laws are governed by the short-ranged Ising universal-ity universality class for three dimensions, which fixes theexponents to their precisely known values, ν ≈ .
630 and (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242) (cid:230) (cid:224) (cid:242) (cid:244) k (cid:45) (cid:43)Η k Σ S cc (cid:72) k (cid:76) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) k Σ S (cid:45) cc (cid:72) k (cid:76) (cid:230) (cid:230) (cid:230)(cid:230)(cid:230) (cid:230) Ξ Ε (cid:45)Ν (cid:45) (cid:200) T (cid:45) T c (cid:200)(cid:144) T c Ξ (cid:72) T (cid:76) (cid:144) Σ (cid:230) (cid:230) (cid:230) (cid:230)(cid:230) (cid:230) Χ Ε (cid:45)Γ (cid:45) (cid:200) T (cid:45) T c (cid:200)(cid:144) T c Χ (cid:72) T (cid:76) (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) Χ (cid:72) Ξ (cid:144) Ξ (cid:76) Γ (cid:144) Ν Ξ (cid:72) T (cid:76)(cid:144) Σ Χ (cid:72) T (cid:76) (a) (b)(c) (d) FIG. 1. 246.0pt(a) Static structure factor S cc ( k ) of the binary fluid atcritical composition studied for selected temperatures approaching theconsolute point from above. The reduced temperatures T * = k B T /(cid:15) are given in the legend. The data are shown on double-logarithmicscales, the inset displays a rectification of the data according to Eq. (2).Solid lines are fits to Eqs. (2) and (5), respectively (see text), and thedashed line indicates the critical divergence, S cc ( k → ∼ k − + η at T = T c . (b–d) Critical divergences of the correlation length ξ ( T ) andthe concentration susceptibility χ ( T ), the data were obtained from fitsto S cc ( k ) as shown in panel (a). Solid lines are power-law fits usingthe known critical exponents ν and γ (see main text) to estimate thecritical temperature T c and the non-universal amplitudes ξ and χ ,respectively. γ ≈ . ; the non-universal amplitudes ξ and χ are system-specific. Moreover, the static structure factorassumes the scaling form S cc ( k ) = k − + η s ( k ξ ) ; k σ ≪ , ξ/σ ≫ , (4)where the anomalous dimension η = − γ/ν ≈ .
036 charac-terises the spatial decay of the OP correlations and s ( x ) is auniversal scaling function, which interpolates from the criticalpower-law, s ( x ≫ = const , to the convergence of S cc ( k ) atsmall wavenumber, s ( x → = x − η ; the latter encodes thescaling of the susceptibility. The anomalous scaling of S cc ( k )requires a modification of Eq. (2) to the form, S cc ( k ) ≃ ρ k B T χ [︁ + ( k ξ ) ]︁ − + η/ . (5)Our simulation data for S cc ( k ) are well described by relationsEqs. (2) and (5), respectively, for temperatures ranging from T * = .
70 down to 1 .
43, see Fig. 1(a). Linear regression of S cc ( k ) − as function of ( k σ ) for k σ (cid:46) ξ and χ (see inset); Eq. (5) was used for T * (cid:54) .
60. The datanicely follow the expected critical singularities over more thanone decade in magnitude and for ε (cid:46) . T c as well asthe amplitudes ξ and χ . Using the correlation length data,we obtained T * c = . ± .
001 and ξ = (0 . ± . σ ,while the data for the concentration susceptibility yielded T * c = . ± .
001 and χ * = . ± . T c and ξ are in excellent agreement with known results . Merely ourvalue for χ * is somewhat larger than the previously reportedvalue of χ * = . ± . IV. SCALING ANALYSIS OF SUB-SYSTEM FLUCTUATIONSA. Concentration susceptibility
In a typical semi-grand canonical Monte Carlo simulation,the fluctuations in the concentration are sampled by switchingthe particle identities at a predefined rate while keeping thetotal particle number fixed. Then by linear response theory,the susceptibility χ in the mixed phase is proportional to thevariance of the fluctuating mole fraction x A = N A / N of, e.g.,species A : ρ k B T χ = N ⟨ ( x A − x c ) ⟩ gc , T > T c , (6)provided that the symmetric mixture is at critical composition, ⟨ x A ⟩ gc = ⟨ x B ⟩ gc = x c = /
2; in practice, the averages are takenin the semi-grand canonical ensemble . The standardprocedure for the finite-size scaling analysis is then based ona sequence of simulations for different boxes of linear extent L to determine χ ( T ; L ) for each value of L and to study thebehaviour as L → ∞ .Above, we extrapolated the composition structure factorby virtue of its specific wavenumber dependence, Eq. (5), toaccess the susceptibility, S cc ( k → ≃ ρ k B T χ . Such a sys-tematic change of the wavenumber has a close analogy to thedescribed finite-size scaling analysis. A finite system of size L does not support fluctuations on length scales larger than L and, likewise, the correlations contained in S cc ( k ) for given k stem essentially from fluctuations on length scales (cid:46) π/ k .Thus, the value S cc ( k ) /ρ k B T may be interpreted as the effec-tive susceptibility of a finite system of size 2 π/ k . Observingthat a single simulation yields data for different wavenumberssuggests to perform the conventional scaling analysis but for asequence of finite sub-systems .The idea is easily linked to Eq. (6) by noting that sub-volumes of a large system represent open systems coupledto a (finite) particle reservoir. Within each sub-volume, particlenumber and concentration exhibit similar fluctuations as inthe grand canonical ensemble, although they are still locallyconserved quantities. . Specifically, we partitioned the sim-ulation box of linear extent L into m cubic sub-systems ofedge length (cid:96) = L / m , and recorded the fluctuating particlenumbers N ( (cid:96) ) A and N ( (cid:96) ) B for each sub-system; the total number ofparticles is denoted by N (cid:96) : = N ( (cid:96) ) A + N ( (cid:96) ) B . The susceptibility χ (cid:96) is then calculated from Eq. (6) by substituting x A = N ( (cid:96) ) A / N (cid:96) forthe instantaneous concentration and replacing N by the averageparticle number ⟨ N (cid:96) ⟩ = N / m in each sub-system. (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) (cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242) (cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244) k B T c (cid:144) Ε(cid:61) (cid:72) a (cid:76) (cid:230) (cid:224) (cid:236) (cid:242) (cid:244) y (cid:61) (cid:123) (cid:144) Ξ (cid:182) Γ Χ (cid:42) (cid:123) (cid:72) T (cid:76) (cid:180) (cid:45) (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) (cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242) (cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244) k B T c (cid:144) Ε(cid:61) (cid:72) b (cid:76) (cid:230) (cid:224) (cid:236) (cid:242) (cid:244) y (cid:61) (cid:123) (cid:144) Ξ (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) (cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242) (cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244) k B T c (cid:144) Ε(cid:61) (cid:72) c (cid:76) (cid:230) (cid:224) (cid:236) (cid:242) (cid:244) y (cid:61) (cid:123) (cid:144) Ξ (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) (cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242)(cid:242) (cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244)(cid:244) k B T c (cid:144) Ε(cid:61) (cid:72) d (cid:76) (cid:230) (cid:224) (cid:236) (cid:242) (cid:244) y (cid:61) (cid:123) (cid:144) Ξ FIG. 2. Scaling plots of the reduced susceptibility χ * ( T ) : = χ ( T ) (cid:15)/σ according to Eq. (7) with the correlation length ξ = ξ ( T ) obtained fromthe static structure factor as in Fig. 1. The data collapse onto the scaling function χ Z ( y ) is tested for four different choices of the criticaltemperature T c as written in each panel, which enters the reduced temperature ε : = ( T − T c ) / T c . Black solid lines represent the critical asymptoteas y →
0, amended by the leading universal correction [Eq. (9)], which employs the amplitudes z , z from the fit to f in Fig. 3. The blackdashed line in panel (c) is a fit to the large- y approximation [Eq. (8)] of χ Z ( y ) and indicates the approach to the Gaussian fixed point, y → ∞ .The amplitude χ (horizontal black line) estimated from this fit and its uncertainty (gray shaded region) are indicated. The standard finite-size scaling of the susceptibility sug-gests that χ (cid:96) for T > T c asymptotically obeys the scalingansatz χ (cid:96) ( T ) ≃ χ ε − γ Z ( (cid:96)/ξ ) , (cid:96), ξ ≫ σ , (7)where ε = ε ( T ) and ξ = ξ ( T ) as before. The function Z ( · )is the appropriate scaling function that interpolates betweenthe thermodynamic critical divergence ( (cid:96) → ∞ , ξ fixed) andthe finite-size scaling at criticality ( ξ → ∞ , (cid:96) fixed). Inthe first limit, the scaling variable y : = (cid:96)/ξ → ∞ and thescaling function approaches unity exponentially fast , and asatisfactory description of the data for y (cid:38) Z ( y → ∞ ) ≃ (︀ − z ∞ e − a y )︀ , (8)including the leading correction to scaling with z ∞ and a beingsystem-specific parameters (see Figs. 2 and 3).For the opposite limit y →
0, we recall that for a finitesystem, all physical observables are analytic in the temperature,also at T c . For the susceptibility, it implies χ (cid:96) ( T ) = χ (cid:96) ( T c ) + O ( ε ) as ε → Z ( y → ≃ z y γ/ν (︁ − z y /ν )︁ , (9)which introduces the universal amplitudes z and z . Combin-ing with ξ ≃ ξ ε − ν and specialising to T = T c , one recovers thefinite-size scaling form well-known for periodic boundaries, χ (cid:96) ( T c ) ≃ z χ (︃ (cid:96)ξ )︃ γ/ν , (cid:96) ≫ σ. (10)More generally, one finds for the critical region near T c a lineardependence on T : χ (cid:96) ( T ) ≃ ˜ χ (cid:96) γ/ν [︁ − ˜ z (cid:96) /ν ( T − T c ) ]︁ (11)with coefficients ˜ χ : = z χ ξ − γ/ν and ˜ z : = z ξ − /ν T − c for theregime σ ≪ (cid:96) ≪ ξ ( T ). The scaling form Eq. (7) can be used to determine the criticaltemperature T c , the amplitude χ , and the scaling function Z ( y ).Plotting the data for ε γ χ (cid:96) ( T ) for a range of temperatures andsub-system sizes as a function of y = (cid:96)/ξ should yield datacollapse onto the function χ Z ( y ). Figure 2 shows a sequenceof such plots for different estimates of T c , which are informedby our previous result from Section III. The quality of the datacollapse for the different choices of T c gives us a fairly goodestimate of the critical temperature. Comparing panels (b) and(c) of Fig. 2, in particular close to y = (cid:96)/ξ ≈ .
8, suggests thatthe critical temperature is close to T * c = .
422 [Fig. 2(c)].A more quantitative estimate is obtained by using the formof Z ( y ) as y →
0. However, this requires knowledge of theparameters z and z , which can independently be obtainedfrom the sub-system scaling of the second moment of the orderparameter φ : = x A − x c as discussed in Section IV B. Specif-ically, the quantity a (cid:96) β/ν ⟨ φ ⟩ (cid:96) in the limit of y = (cid:96)/ξ → z (1 − z y /ν ); for the definition of a seeEq. (13) below. Fitting this functional form to the data for ⟨ φ ⟩ (cid:96) , scaled by a (cid:96) β/ν , within the range 0 . (cid:54) y (cid:46) . z = . ± .
004 and z = . ± .
05 (Fig. 3).Equipped with these values, we plot the scaling function Z ( y )for y → a contains the critical temperature T c , why our estimate T * c = .
421 from the susceptibility data given in Section IIIwas used to fix a . It is worth noting that from the data collapseof χ (cid:96) , the value of T c obtained here agrees very well with thatobtained from the critical divergences of the correlation lengthand the susceptibility, given the measurement uncertainty ofthe values of ξ entering the scaling plots in Fig. 2.For this value of T c , the rescaled data represent the func-tion χ Z ( y ) (Fig. 2c). For values of y = (cid:96)/ξ (cid:38)
20, the dataindeed converge and, in principle, one could read off χ di-rectly. However, a more robust approach should allow forcorrections to scaling, Eq. (8), which suggests to fit this formto the data with χ , z ∞ and a as parameters. Even though therange of values for large y is severely limited, we performed anasymptotic fit as y → ∞ to estimate the amplitude χ , yielding χ = . ± .
002 and thereby improving our earlier estimateof χ = . ± . a from the best fit comesout to be a = . ± .
02 and that of z ∞ = . ± . B. Order parameter distribution
The concentration susceptibility χ is essentially the varianceof the fluctuating order parameter φ : = x A − x c , see Eq. (6),and thus merely one characteristic of the statistical distributionof φ . A more general description is in terms of the probabilitydensity P (cid:96) ( φ ; T ) of φ for sub-system size (cid:96) and temperature T ,which also admits a finite-size scaling hypothesis : P (cid:96) ( φ ; T ) ≃ a (cid:96) β/ν P (︁ a (cid:96) β/ν φ, (cid:96)/ξ ( T ) )︁ , (cid:96), ξ ( T ) ≫ σ, (12)with a universal scaling function P (ˆ y , y ), where ˆ y : = a (cid:96) β/ν φ and y : = (cid:96)/ξ ( T ), and a scale factor a such that a − = k B T c χ ξ − γ/ν . (13)By normalisation of P (cid:96) , it holds ∫︀ ∞−∞ P (ˆ y , y ) dˆ y = y . At critical composition of the mixture, ⟨ x A ⟩ (cid:96) = x c , the firstmoment of P (cid:96) vanishes due to symmetry, ⟨ φ ⟩ (cid:96) =
0. The secondmoment encodes the scaling of the susceptibility, χ (cid:96) ( T ) = (cid:96) d k B T ∫︁ d φ φ P (cid:96) ( φ ; T ) (14) ≃ (cid:96) d − β/ν a k B T ∫︁ dˆ y ˆ y P (ˆ y , (cid:96)/ξ ) , (15)by virtue of the definition of χ (cid:96) analogous to Eq. (6), with ⟨ N (cid:96) ⟩ = ρ(cid:96) d , and using Eq. (12) in the second step. For everyvalue of y = (cid:96)/ξ , the ˆ y -integral yields a constant f ( y ). Invokingthe hyper-scaling relations γ/ν = − η = d − β/ν , we recoverEq. (7): χ (cid:96) ( T ) ≃ (cid:96) γ/ν a k B T f ( (cid:96)/ξ ) (16) ≃ χ ε − γ + ε Z ( (cid:96)/ξ ) , (cid:96), ξ ≫ σ (17)upon identifying the scaling function as Z ( y ) = y γ/ν f ( y ), sub-stituting ξ = ξ ε − ν and k B T = k B T c (1 + ε ), and due to ourchoice of a .Similarly, one readily obtains for the k th moment of the orderparameter the scaling form ⟨︀ φ k ⟩︀ (cid:96) : = ∫︁ d φ φ k P (cid:96) ( φ ; T ) (18) ≃ a − k (cid:96) − k β/ν ∫︁ ∞−∞ dˆ y ˆ y k P (ˆ y , (cid:96)/ξ ) . (19) = : a − k (cid:96) − k β/ν f k ( (cid:96)/ξ ) , (20)which defines universal scaling functions f k for k = , , . . . If P (cid:96) ( φ ) is symmetric, ⟨︀ φ k ⟩︀ (cid:96) = k odd. In Fig. 3, we test this (cid:230)(cid:230) (cid:230)(cid:224)(cid:224)(cid:224) (cid:224) (cid:224) (cid:224)(cid:236)(cid:236)(cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242)(cid:242)(cid:242) (cid:242) (cid:242)(cid:244) (cid:244)(cid:244)(cid:244) (cid:244) (cid:244)(cid:231) (cid:231)(cid:231)(cid:231) (cid:231) (cid:231)(cid:225) (cid:225) (cid:225)(cid:225)(cid:225) (cid:225) (cid:225)(cid:237) (cid:237) (cid:237)(cid:237)(cid:237) (cid:237)(cid:243) (cid:243) (cid:243)(cid:243)(cid:243) (cid:230)(cid:230) (cid:230)(cid:224)(cid:224)(cid:224) (cid:224) (cid:224) (cid:224)(cid:236)(cid:236)(cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242)(cid:242)(cid:242) (cid:242) (cid:242)(cid:244) (cid:244)(cid:244)(cid:244) (cid:244) (cid:244)(cid:231) (cid:231)(cid:231)(cid:231) (cid:231) (cid:231)(cid:225) (cid:225) (cid:225)(cid:225)(cid:225) (cid:225) (cid:225)(cid:237) (cid:237) (cid:237)(cid:237)(cid:237) (cid:237)(cid:243) (cid:243) (cid:243)(cid:243)(cid:243) k (cid:61) (cid:61) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) a k (cid:123) k (cid:32) Β (cid:144) Ν (cid:88) Φ k (cid:92) (cid:123) y (cid:61) (cid:123) (cid:144) Ξ (cid:230) (cid:224) (cid:236) (cid:242) (cid:244) (cid:231) (cid:225) (cid:237) (cid:243) (cid:123) (cid:144) Σ FIG. 3. Scaling of the 2 nd and 4 th moments of the concentration fluc-tuations φ = x A − x c within the sub-volumes. The data collapse ontothe scaling functions f k ( y ) is tested by plotting a k (cid:96) k β/ν ⟨ φ k ⟩ (cid:96) against y = (cid:96)/ξ ( T ) for k = k =
4, respectively, with the critical expo-nents β and ν given by the Ising universality class (see text), the cor-relation length ξ = ξ ( T ) obtained from S cc ( k ) [Fig. 1], and a being aconstant factor [Eq. (13)]. The different symbols refer to data from thesub-volume sizes (cid:96)/σ specified in the legend. Lines show the small-and large- y asymptotes of the scaling functions f ( y ) = y − γ/ν Z ( y )and, within the Gaussian approximation, of f ( y ) ≃ f ( y ) ; the latterbecomes exact as y → ∞ . Equation (8) was used for the large- y be-haviour of Z ( y ) (solid lines), whereas the critical asymptote Z ( y → scaling prediction for k = a k (cid:96) k β/ν ⟨ φ k ⟩ (cid:96) against (cid:96)/ξ thedata collapse nicely onto the functions f k over the full range0 . (cid:46) (cid:96)/ξ (cid:46) y ≫
1, i.e., for large, near-critical sub-systems,the distribution P (cid:96) ( φ ) is Gaussian, which is inherited to thescaling function P ( · , y ) being Gaussian in its first argumentwith zero mean and variance f ( y ). In this case, the momentsand thus the functions f k are related to each other since allcumulants except the first two vanish. Specifically, ⟨ φ ⟩ (cid:96) ≃ ⟨ φ ⟩ (cid:96) for (cid:96) ≫ ξ ≫ σ , which implies f ( y → ∞ ) ≃ f ( y ) = y − γ/ν Z ( y ) . (21)For comparison, the functions f and f as predicted from thesmall- and large- y approximations to Z ( y ) (Fig. 2) have beenincluded in Fig. 3. Both functions describe the 2 nd and 4 th moments very well for the whole range of y -values available.Merely for y (cid:46)
1, the data for k = f ( y ),indicating a non-Gaussian distribution as expected close tocriticality. V. BINDER CUMULANT FOR SUB-SYSTEMS
The data collapsing approach to determine T c as describedin Section IV A has some subjective component. A superiormethod to locate a continuous phase transition was establishedby Binder and is based on the 4 th normalised cumulant U (cid:96) ( T ) of the order parameter distribution; it is closely relatedto the kurtosis used in descriptive statistics. Close to criticality,the composition fluctuations φ = x A − x c are symmetric undersign change and one defines the dimensionless ratio U (cid:96) ( T ) = − ⟨ φ ⟩ (cid:96) ⟨︀ φ ⟩︀ (cid:96) , (22)where the subscript (cid:96) , as before, denotes the linear extent ofthe (sub-)system. At high temperature, the fluctuations are ofGaussian nature, and thus U (cid:96) ( T ≫ T c ) →
0. Far below thecritical temperature, the distribution has two sharp peaks atthe coexisting compositions, and the Binder cumulant tends to U (cid:96) ( T ≪ T c ) → /
3. As the critical temperature is approachedfrom either side of T c , the correlation length in the systemdiverges and the critical divergences of the moments cancel[cf. Eq. (20)] so that U (cid:96) ( T c ) = : U c remains finite and becomesindependent of the sub-system size (cid:96) . The critical Bindercumulant U c = − f (0) / f (0) is a universal property of thecritical renormalisation group fixed point and as such dependsonly on the boundary conditions and the geometric shape ofthe sub-system .The analysed open sub-systems mimick a grand canonicalensemble in the limit of an infinitely large reservoir, L → ∞ ,which is relaxed in practice to the condition that the sub-systemsize (cid:96) does not compete with the finite extent L of the wholesimulation box, (cid:96) ≪ L ; see ref. [23] for a quantitative estimate.We consider the idealised case L → ∞ first, corrections due tothe finite simulation box are studied in Section V D. A. Common intersection point
The basis of our discussion of the 4 th cumulant is the generalscaling ansatz U (cid:96) ( T ) = U (︀ (cid:96) /ν ε )︀ , ε = ( T − T c ) / T c , (23)employing a scaling function U ( · ) that is analytic since U (cid:96) ( T )describes finite (sub-)volumes and depends smoothly on tem-perature. Scaling is expected to hold when all length scales aresufficiently large, in particular, when (cid:96) ≫ σ ; the ratio (cid:96)/ξ ( T )is controlled by the scaling variable x : = (cid:96) /ν ε . Moreover, U ( x ) fulfills U (0) = U c and interpolates between the limits U ( x → ∞ ) = U ( x → −∞ ) = /
3. Expanding Eq. (23)for small argument shows that, near criticality, the quantity U (cid:96) ( T ) varies linearly with temperature around U c : U (cid:96) ( T ) ≃ U c + u (cid:96) /ν ( T − T c ) , T → T c , (24)where u : = U ′ (0) T − c is a non-universal constant.Accordingly, a family of curves U (cid:96) ( T ) for different sub-system sizes σ ≪ (cid:96) ≪ ξ ( T ) has a common intersection pointat ( T c , U c ), in an asymptotic sense, which suggests a procedureto locate the critical temperature T c . It requires the simula-tion of one large system for a number of temperatures suffi-ciently close to T c , which in itself is challenging due to critical slowing down. Our simulation results for sub-system sizes (cid:96) = L / m with m = , , ,
12 and L ≈ . σ are shown inFig. 4(a,b). The data exhibit a common intersection point asanticipated from Eq. (24), and we read off T * c = . U c = . ± .
01, close to theearlier reported value in Ref. 14 for Ising spins on a lattice.Note that U c adopts a value that is different from the 3D Isinguniversal value . ± . . The open sub-systems studiedhere are characterised by free boundary conditions for smallsub-systems ( (cid:96) ≪ L ), and the quoted value for U c is expectedto be universal for finite-size scaling based on the analysis ofsub-systems. B. Scaling function
A more physics-adapted way of writing the scaling formEq. (23) for T > T c is U (cid:96) ( T ) = ˆ U + ( (cid:96)/ξ ) with ˆ U + ( y (cid:62)
0) : = U (︁ ( ξ y ) /ν )︁ , noting that ε ≃ ( ξ/ξ ) − /ν in the critical region. In-deed, plotting our results for U (cid:96) ( T ) against y = (cid:96)/ξ ( T ) for dif-ferent values of (cid:96) yields data collapse onto the scaling functionˆ U + ( y ), see Fig. 5. For small arguments, i.e., taking T → T c for (cid:96) fixed, the data converge to ˆ U + ( y → = U c . Near (cid:96)/ξ ≈
1, a crossover occurs to the high-temperature regime(which implies small ξ , i.e., (cid:96)/ξ → ∞ ), where the fluctuationsare Gaussian and thus ˆ U + ( y → ∞ ) → U + ( y → ∞ )? An intuitiveargument assumes that the sub-volume of linear extent (cid:96) canbe divided into independent “correlation blobs” of size ξ andinvokes a standard proof of the central limit theorem: considerthe sum Y = ∑︀ ni = X i of n independent random variables X i that are identically distributed according to some characteristicfunction ϕ X ( · ) with variance σ X < ∞ . Then Y has the character-istic function ϕ Y ( · ) = ϕ X ( · ) n and its cumulants µ k are generatedby n log ϕ X ( · ), which shows that µ k ∝ n for all k = , , . . . Thenormalised variable Y /σ Y , with σ Y : = n σ X , becomes Gaus-sian as n → ∞ since its cumulants follow µ k σ − kY ∝ n − k / andvanish for k (cid:62) Y / n instead of Y since the concentration x A in a sub-volume of size (cid:96) is given bythe arithmetic mean of the values of x A in each of the n = ( (cid:96)/ξ ) d correlation blobs. The cumulants of Y / n are proportional to n − k , such that the normalised cumulant [Eq. (22)] vanishesas U (cid:96) ∝ n − ⧸︀(︁ n − )︁ = n − . Thus, we expect that ˆ U + ( y ) ∼ y − d for y = (cid:96)/ξ → ∞ , here d =
3. Despite the limited availabilityof data for U (cid:96) in this regime, there is numerical evidence forthe scaling U (cid:96) ∼ y − for y (cid:38) U L ∼ L − d , see Eq. (14) of Ref. [14].Eventually, using x = ( (cid:96)/ξ ) /ν as the scaling variable recti-fies the critical power-law ˆ U + ( y ) − U c ∼ y /ν for y → U ( x ) − U c ∼ x [Eq. (24)], which is supported bythe inset of Fig. 5. (cid:230) (cid:230) (cid:230) (cid:230)(cid:230)(cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224) (cid:224)(cid:224)(cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:236) (cid:236) (cid:236) (cid:236)(cid:236)(cid:236)(cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242) (cid:242) (cid:242)(cid:242)(cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:72) a (cid:76) U (cid:123) (cid:72) T (cid:76) k B T (cid:144) Ε (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:72) b (cid:76) U (cid:123) (cid:72) T (cid:76) k B T (cid:144) Ε (cid:231) (cid:231) (cid:231) (cid:231) (cid:231)(cid:225) (cid:225) (cid:225) (cid:225) (cid:225) (cid:237) (cid:237) (cid:237) (cid:237) (cid:237) (cid:243) (cid:243) (cid:243) (cid:243) (cid:243) (cid:72) c (cid:76) U (cid:123) (cid:144) U (cid:123) ' k B T (cid:144) Ε (cid:230) (cid:224) (cid:236) (cid:242) (cid:123) (cid:144) Σ FIG. 4. (a) Simulation results for the Binder cumulant U (cid:96) ( T ) for different sub-system sizes (cid:96) as indicated in the legend, with the same edgelength L ≈ . σ of the overall simulation box. The data are based on the moments shown in Fig. 3, and solid lines are fits of a tanh( · )-shapeserving as a guide to the eye. (b) Close-up of the critical region of the same data as in panel (a). Black lines indicate the common intersectionpoint at the critical values ( T c , U c ). (c) Cumulant ratios U (cid:96) / U (cid:96) ′ as a function of temperature for different pairs of sub-system sizes: ( (cid:96), (cid:96) ′ )denoted by the symbols ∘ (7 . σ, . σ ), (cid:3) (7 . σ, . σ ), ◇ (7 . σ, . σ ) and △ (7 . σ, . σ ). Coloured lines are linear fits according toEq. (25), and the horizontal black line marks the fixed point U (cid:96) = U (cid:96) ′ . (cid:230) (cid:230) (cid:230)(cid:230) (cid:230) (cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:236) (cid:236) (cid:236)(cid:236) (cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242)(cid:244) (cid:244) (cid:244) (cid:244) (cid:244) (cid:244) (cid:244) (cid:244)(cid:231) (cid:231) (cid:231)(cid:231) (cid:231) (cid:231) (cid:231) (cid:231) y (cid:45) (cid:230) (cid:224) (cid:236) (cid:242) (cid:244) (cid:231) (cid:123) (cid:144) Σ y (cid:61) (cid:123) (cid:144) Ξ U (cid:123) (cid:180) (cid:45) (cid:230)(cid:224) (cid:224) (cid:224) (cid:224)(cid:236) (cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242) (cid:242) (cid:242) (cid:242)(cid:244) (cid:244) (cid:244) (cid:244) (cid:244) (cid:244)(cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) U (cid:123) (cid:72) (cid:123) (cid:144) Ξ (cid:76) (cid:144) Ν FIG. 5. Test of the scaling form Eq. (23) of the Binder cumlant U (cid:96) ( T ) plotted against y = (cid:96)/ξ ( T ) for the data in Fig. 4 but only fortemperatures above T c . The dashed line indicates an algebraic decay U (cid:96) ( T ) ∼ [ (cid:96)/ξ ( T )] − . Inset: Data collapse onto the scaling function U ( x ) is obtained by plotting the same data with x = ( (cid:96)/ξ ) /ν ε asthe abscissa [Eq. (23)]. U ( x ) is analytic in x =
0, and the linearextrapolation of the data [solid line, Eq. (24)] intersects the verticalline x = U c : = U (cid:96) ( T c ) ≈ . C. Cumulant ratios
An alternative method for estimating T c from U (cid:96) ( T ) thatdoes not depend on the value of U c is to consider the ratio U (cid:96) / U (cid:96) ′ as function of temperature, where (cid:96) > (cid:96) ′ are two dif-ferent sub-system sizes. For any choice of (cid:96) and (cid:96) ′ the ratio U (cid:96) / U (cid:96) ′ = T = T c , since U (cid:96) ( T ) is independent of (cid:96) at thecritical point, see Fig. 4(c). The behaviour for T → T c follows (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) (cid:72) b (cid:76) U (cid:123) ' U (cid:123) (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) (cid:72) a (cid:76) U (cid:123) ' U (cid:123) FIG. 6. (a) Graph of U (cid:96) ′ ( T ) vs. U (cid:96) ( T ) in the vicinity of the fixedpoint ( U c , U c ) for the investigated temperatures and for different pairsof the sub-system sizes ( (cid:96), (cid:96) ′ ) with (cid:96) = . σ fixed and (cid:96) ′ = . σ ( ∙ ), 4 . σ ( (cid:4) ), and 3 . σ ( (cid:4) ); the size of the simulation box was L ≈ . σ . Broken lines are fits of Eq. (26) to the data pointsfor the same ( (cid:96), (cid:96) ′ ) with U c as only free parameter; the solid lineindicates U (cid:96) = U (cid:96) ′ . (b) Same representation as in panel (a) with solidlines showing linear fits to data with the slope and U c as parameters.The inset provides a close-up of the intersection of these lines withthe diagonal U (cid:96) = U (cid:96) ′ ; arrows indicate the range of an anticipatedcommon intersection point. from Eq. (24) U (cid:96) ( T ) U (cid:96) ′ ( T ) ≃ + u U c (cid:96) /ν [︁ − ( (cid:96)/(cid:96) ′ ) − /ν ]︁ ( T − T c ) , (25)which allows for a linear regression of the data for U (cid:96) / U (cid:96) ′ near T c to find the intersection with unity. A larger slope isachieved for a larger sub-system size (cid:96) and a larger ratio (cid:96)/(cid:96) ′ ,whereas (cid:96) ′ ≫ σ must not be chosen too small. From this, weinferred k B T c /(cid:15) = . ± . T c given above.The determination of the critical Binder cumulant U c followsa similar approach. From the foregoing discussion it is clearthat the graph of U (cid:96) ( T ) vs. U (cid:96) ′ ( T ) for a given choice of (cid:96) and (cid:96) ′ displays a fixed point at U c . Thus, U c can be determinedfrom the intersection of this graph with the diagonal, U (cid:96) = U (cid:96) ′ .Close to criticality, solving Eq. (24) for T and substituting into U (cid:96) ′ ( T ) yields a linear relationship between U (cid:96) ′ and U (cid:96) (at thesame temperature): U (cid:96) ′ = U c + ( (cid:96)/(cid:96) ′ ) − /ν ( U (cid:96) − U c ) ; U (cid:96) → U c . (26)Thus, U c will be the only free parameter in a linear regres-sion of the data for ( U (cid:96) , U (cid:96) ′ ). The procedure is illustratedin Fig. 6 with three different choices of ( (cid:96), (cid:96) ′ ). The data for U (cid:96) ( T ) plotted against U (cid:96) ′ ( T ) for a range of temperatures T andfixed ( (cid:96), (cid:96) ′ ) do indeed fall on straight lines as inferred fromEq. (24). However, the slopes do not match with the prediction( (cid:96)/(cid:96) ′ ) − /ν of Eq. (26), which points at a deficiency of the ansatzEq. (23). Nevertheless, permitting both U c and the slope as fitparameters yields approximately a common intersection pointof the straight lines at U c ≈ .
22 [Fig. 6(b)]. A close-up of thisintersection region reveals appreciable differences between theintersection points of two of the coloured curves (for differentchoices of (cid:96) ′ ) and their intersection with the diagonal. Thefailure of Eq. (26) and this observation led us to revisit ourscaling ansatz and to rederive Eq. (24) in the next section. Theerror stems from the coefficient u , which was obtained as u : = U ′ (0) T − c , that is, as constant with respect to (cid:96) . However,taking into consideration the finite size of the simulation boxshows that in fact u depends on (cid:96)/ L . D. Finite-size corrections
So far, we have ignored the finiteness of the total simulationvolume, which can taint the estimates of the critical point, in-cluding both T c and U c . In particular, there is a competitionbetween the sub-system size (cid:96) and the length L of the wholesimulation box, leading to a kind of effective boundary con-ditions on the sub-system as (cid:96) grows. (For example, consider (cid:96) = L /
2, which implies boundary conditions that are neitherfree nor periodic.) Further, the correlation length ξ enters as athird length scale, and it may be necessary to distinguish theregimes (cid:96) ≪ ξ ≪ L and (cid:96) ≪ L ≪ ξ , in addition to ξ ≪ (cid:96), L .In the following, we will assess the importance of these cor-rections for our results combining theoretical arguments andsimulation data for a range of box sizes L .Conventional finite-size scaling is based on the fact that near T c , the correlation length exceeds the system size by far. Forthe sub-system analysis, the corresponding regime is (cid:96), L ≪ ξ and we expect that the predominant corrections due to finite L are controlled by the aspect ratio α : = (cid:96)/ L . This suggests toamend the scaling ansatz (23) by the variable (cid:96)/ L to U (cid:96) ( T ; L ) = ˜ U (︁ (cid:96) /ν ε, (cid:96)/ L )︁ , ε = ( T − T c ) / T c , (27)which is expected to hold when all length scales are suf-ficiently large, i.e., for (cid:96), L , ξ ( T ) ≫ σ . (An alternativeto this ansatz is discussed briefly in Appendix A.) Taking (cid:230) (cid:230) (cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242) (cid:242) (cid:242) (cid:242)(cid:244) (cid:244) (cid:244) (cid:244) (cid:244) (cid:244) (cid:72) a (cid:76) U Α L k B T (cid:144) Ε (cid:230) (cid:230) (cid:230) (cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224) (cid:224) (cid:224) (cid:224)(cid:236) (cid:236) (cid:236) (cid:236) (cid:236) (cid:236)(cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:242) (cid:72) b (cid:76) U Α L k B T (cid:144) Ε (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:72) c (cid:76) U c Α (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:72) d (cid:76) T c Α (cid:230) (cid:224) (cid:236) (cid:242) (cid:244) L (cid:144) Σ Α (cid:45) (cid:61) (cid:230) (cid:224) (cid:236) (cid:242) L (cid:144) Σ Α (cid:45) (cid:61) FIG. 7. (a),(b) Binder cumulant U (cid:96) ( T ; L ) as a function of temperaturefor fixed ratios α = (cid:96)/ L between the sizes of the sub-system ( (cid:96) ) andthe simulation box ( L ) as indicated in each panel. Linear fits (solidlines) in the critical region exhibit a common intersection point at (︀ ˜ U c ( α ) , T c ( α ) )︀ , which determines effective, α -dependent values for U c and T c according to Eq. (28). (c),(d) Results for U c ( α ) and T c ( α ), respectively, obtained from aspect ratios of α − = , , . . . , L ≈ . σ to 50 . σ . Solid lines aresmooth interpolations of the data to guide the eye. L → ∞ yields Eq. (23), and in terms of the scaling func-tions: U ( x ) = ˜ U ( x , α → U ( x , α ) is ana-lytic in x = α (cid:62)
0, which we infer fromthe fact, used before, that in a finite system thermodynamicobservables depend smoothly on temperature. The depen-dence on α , on the other hand, is not known albeit an exponen-tially fast approach to the thermodynamic limit is not unlikely:˜ U ( x , α ) − ˜ U ( x , = O (︀ e − /α )︀ as α → x = U (cid:96) ( T ; L ) ≃ ˜ U c ( (cid:96)/ L ) + ˜ u ( (cid:96)/ L ) (cid:96) /ν ( T − T c ) , (28)as T → T c , where we introduced an effective critical cumulantas ˜ U c ( α ) : = ˜ U (0 , α ) and the geometry-dependent coefficient˜ u ( α ) : = T − c ∂ x ˜ U ( x , α ) ⃒⃒⃒ x = . It becomes evident from Eq. (28)that the curves U (cid:96) ( T ; L ) vs. T for different (cid:96) , but the samesystem size L , do not have a common intersection point. Sucha point emerges only asymptotically for (cid:96)/ L sufficiently smallsuch that ˜ U c and ˜ u do not depend appreciably on (cid:96)/ L , i.e., inthe limit L → ∞ .Nevertheless, a common intersection point at ( T c , ˜ U c ( (cid:96)/ L ))is achieved for fixed ratios α = (cid:96)/ L , since then the (cid:96) -dependence enters only the term proportional to T − T c inEq. (28). In practice, the reciprocal m = /α is an integernumber counting the sub-volumes (along each axis) that fitinto the whole system. Thus, the refined analysis procedureaccounting for finite system sizes L < ∞ would be as follows:From a set of simulations for a few system sizes L and many T values in the critical regime, compute the sub-system statis-tics and U (cid:96) ( T ; L ) in particular for selected values of α . Datafor different α are analysed separately. For given α , plotting U α L ( T ; L ) / U α L ′ ( T ; L ′ ) as function of T the value of T c ( (cid:96)/ L )can be read off from the intersection with unity; more precisely,it follows from the linear regression according to [cf. Eq. (25)] U α L ( T ; L ) U α L ′ ( T ; L ′ ) ≃ + ˜ u ( α ) U c ( α ) ( α L ) /ν [︁ − ( L ′ / L ) /ν ]︁ ( T − T c ) (29)in the limit T → T c .Since Eq. (29) follows directly from Eq. (28), we test thisfinite size scaling analysis on the measured data for the Bindercumulant. To this end, we consider U α L ( T ) as a function oftemperature for different combinations of (cid:96) and L but with α fixed. We have carried out additional simulations with differentbox lengths varying from L ≈ . σ (32 ,
000 particles) to L ≈ . σ (131 ,
072 particles), and the results of the analysisare shown in Fig. 7.Since T c ( α ) and U c ( α ) are a priori not known, we fittedthe data U α L near the critical temperature to a linear function[Eq. (28)]. For each value of α , this lines exhibit a well-definedcommon intersection point, from which we have read off T c ( α )and U c ( α ) [Fig. 7(a,b)]. We note that for α − (cid:62)
10, no commonintersection point could be detected due to almost equal slopes,why we restricted the analysis to α − = , , . . . ,
9. The valuesobtained for U c ( α ) and T c ( α ) show a considerable dependenceon α , yet they converge as α → L ≫ (cid:96) , yields our final estimates of (i) the universal value for the critical Binder cu-mulant U c (0) = . ± .
001 and (ii) the critical temperature T c = . ± .
001 specific to the investigated binary liquid.
VI. SUMMARY AND CONCLUSIONS
We have given a comprehensive analysis of the local orderparameter fluctuations in open sub-systems and compared dif-ferent approaches to locate the critical temperature T c . Theapplicability of the procedures was demonstrated for a sym-metric binary liquid with known phase diagram, fully basedon molecular dynamics simulations of one large system. Thus,such simulations provide an alternative to (semi-)grand canon-ical Monte Carlo schemes, and arguably have the advantageof simultaneously probing the dynamic properties of the sys-tem, e.g., transport coefficients (which we have not discussedhere). Complementary, the study of the static structure factorcalculated for a large system size yields the critical divergencesof the correlation length and the susceptibility (Fig. 1), whichwere found to be compatible with the 3D Ising universalityclass as expected and from which we got a first estimate of T c .For the composition fluctuations, obtained from the particlenumber statistics in cubic sub-volumes, we have shown thatthe standard finite-size scaling procedures, as for a sequenceof simulations with periodic boundaries, are successful if theedge length L of the simulation box is replaced by that of theof sub-volume ( (cid:96) ). In particular, the data for the susceptibility χ (cid:96) and the 4 th moment ⟨ φ ⟩ (cid:96) collapse onto master curves afterappropriate rescaling [see Eqs. (7) and (20)] and if T c is chosenproperly (Figs. 2 and 3).Further, we have found that the Binder cumulant U (cid:96) ( T ) ofthe sub-systems, plotted against temperature T , yields onlyan apparent common intersection point (Fig. 4). Nevertheless,it yields T c to an accuracy of about 0.2% in our example,where we partitioned the simulation box into m cubes for m = , . . . ,
12 with (cid:96) = . σ for the smallest sub-systems.Here, it was favourable to consider the cumulant ratios U (cid:96) / U (cid:96) ′ ,which does not require knowledge of the critical value U c . Dueto the free boundary conditions, the latter adopts a universalvalue U c ≈ . U per c ≈ . U (cid:96) ( T )by the aspect ratio α = (cid:96)/ L as a second scaling variable, wehave shown that a true common intersection point is predictedasymptotically and observed in our simulation data, providedthat lines of constant (cid:96)/ L are considered [Eq. (28) and Fig. 7].A disadvantage of this more correct approach is that it requiresagain a sequence of simulations for a number of large systemsizes, so that (cid:96) can be varied while keeping (cid:96)/ L fixed. Despitethe existence of a common intersection point for aspect ratioseven close to unity (e.g., α = / U c , T c ) can sensitively depend on α , but converged for α (cid:46) . . From the extrapolation α →
0, weimproved previous estimates of the critical Binder cumulant0to U c = . ± . T c = . ± . .In conclusion, we have shown that the analysis of opensub-systems offers a reliable method to locate critical points,thereby taking advantage of large-scale simulations facilitatedby massively parallel computing hardware. A complicationarises due to the interference of the sub-system size (cid:96) withthe size L of the simulation box, which requires a finite-sizescaling procedure with (cid:96)/ L fixed; the latter defeats the idea ofsticking to a single, large value of L . Only for (cid:96)/ L (cid:46) /
10 oreven less, free boundary conditions are effectively realised onthe sub-volume surfaces and the dependence on L drops out. Inmany practical situations, this allows resorting to the simplifiedanalysis with L fixed [Fig. 4], i.e., based on one or few simula-tion runs of one large system. On the other end of the scale, wefound that sub-volume sizes as small as (cid:96) ≈ − σ still permitscaling, so that the choice L ≈ − σ yields sufficient roomfor varying the sub-system size (cid:96) . A spin-off from a large ratio m = L /(cid:96) is that the simulation data at a single instance in timepermit a statistical average over m sub-systems, similarly ascontributions at small wavenumbers k to the structure factor S cc ( k ) are self-averaging. Eventually, the presented sub-systemanalysis combined with the two-parameter scaling should alsobe applicable to the more challenging study of asymmetric,e.g., colloid–polymer mixtures . ACKNOWLEDGMENTS
We thank Siegfried Dietrich (MPI-IS Stuttgart) and Sura-jit Sengupta (TIRF Hyderabad) for useful discussions. DCacknowledges funding from BRNS vide grant number37(3)/14/10/2018-BRNS/370132 and YP acknowledges fund-ing from DST Women Scientist Scheme via grant numberSR/WOS-A/PM-36/2017(G).
Appendix A: Alternative scaling ansatz for U (cid:96) ( T ; L ) It is tempting to propose as a natural extension of the scalingansatz Eq. (23) that U (cid:96) ( T ; L ) = ≈ U (︁ (cid:96) /ν ε, L /ν ε )︁ , ε = ( T − T c ) / T c , (A1)which is expected to hold when all length scales are sufficientlylarge, i.e., for (cid:96), L , ξ ( T ) ≫ σ . Taking L → ∞ yields Eq. (23),and in terms of the scaling functions, U ( x ) = ≈ U ( x , X → ∞ )with X : = L /ν ε , so that lim x → ≈ U ( x , X → ∞ ) = U c ≈ . (cid:96) = L reproduces the conven-tional finite-size scaling with periodic boundaries , and solim x → ≈ U ( x , x ) = U per c ≈ . U c are obtained from carrying out the dataanalysis described in Section V A for various finite L . Wethus define U c ( α ) : = lim x → ≈ U ( x , α − /ν x ) where α : = (cid:96)/ L fixesthe aspect ratio of the geometry. With this, the limit x → α and interpolates between U c for α = U per c for α = ≈ U ( x , X ) is not continuous atthe origin, x = X =
0, not even speaking of analyticity. Yet,the function is analytic in x alone for any fixed L (cid:54) ∞ since U (cid:96) ( T ; L ) depends smoothly on temperature in a finite (sub-)system. The standard treatment to derive asymptotic scalingbehaviour, however, fails as it relies on a Taylor expansion ofthe scaling function around the critical point ( ε = S. Dietrich and M. Schick, “Wetting at a solid-liquid-liquid-vapor tetra point,”Surf. Sci. , 178–181 (1997). J. Köfinger, N. B. Wilding, and G. Kahl, “Phase behavior of a symmetricalbinary fluid mixture,” J. Chem. Phys. , 234503 (2006). R. Folk and G. Moser, “Critical dynamics: a field-theoretical approach,” J.Phys. A , R207 (2006). H. E. Stanley, “Scaling, universality, and renormalization: Three pillars ofmodern critical phenomena,” Rev. Mod. Phys. , S358–S366 (1999). D. Beysens, A. Bourgou, and P. Calmettes, “Experimental determinationsof universal amplitude combinations for binary fluids. I. Statics,” Phys. Rev.A , 3589–3609 (1982). H. C. Burstyn and J. V. Sengers, “Dynamical scaling and critical-pointuniversality of fluids,” Phys. Rev. Lett. , 259–262 (1980). M. Heinen and J. Vrabec, “Evaporation sampled by stationary moleculardynamics simulation,” J. Chem. Phys. , 044704 (2019). A. V. Straube, B. G. Kowalik, R. R. Netz, and F. Höfling, “Rapid onset ofmolecular friction in liquids bridging between the atomistic and hydrody-namic pictures,” Commun. Phys. , 126 (2020). F. Höfling and S. Dietrich, “Enhanced wavelength-dependent surface tensionof liquid-vapour interfaces,” EPL (Europhys. Lett.) , 46002 (2015). T. S. Ingebrigtsen, J. C. Dyre, T. B. Schrøder, and C. P. Royall, “Crys-tallization instability in glass-forming mixtures,” Phys. Rev. X , 031016(2019). P. Chaudhuri and J. Horbach, “Structural inhomogeneities in glasses viacavitation,” Phys. Rev. B , 094203 (2016). J. Roth, H.-R. Trebin, A. Kiselev, and D.-M. Rapp, “Laser ablation of Al–Nialloys and multilayers,” Appl. Phys. A , 500 (2016). J. Glaser, T. D. Nguyen, J. A. Anderson, P. Lui, F. Spiga, J. A. Millan, D. C.Morse, and S. C. Glotzer, “Strong scaling of general-purpose moleculardynamics simulations on GPUs,” Comput. Phys. Commun. , 97–107(2015). K. Binder, “Finite size scaling analysis of Ising model block distributionfunctions,” Z. Phys. B: Condens. Matter , 119–140 (1981). K. Binder, “Critical properties from Monte Carlo coarse graining and renor-malization,” Phys. Rev. Lett. , 693–696 (1981). M. Rovere, D. W. Hermann, and K. Binder, “Block density distributionfunction analysis of two-dimensional Lennard-Jones fluids,” Europhys. Lett.(EPL) , 585–590 (1988). M. Rovere, D. W. Heermann, and K. Binder, “The gas-liquid transitionof the two-dimensional Lennard-Jones fluid,” J. Phys.: Condens. Matter ,7009–7032 (1990). M. Rovere, P. Nielaba, and K. Binder, “Simulation studies of gas-liquidtransitions in two dimensions via a subsystem-block-density distributionanalysis,” Z. Phys. B: Condens. Matter , 215–228 (1993). B. Trefz, J. T. Siebert, T. Speck, K. Binder, and P. Virnau, “Estimation of thecritical behavior in an active colloidal system with Vicsek-like interactions,”J. Chem. Phys. , 074901 (2017). J. T. Siebert, F. Dittrich, F. Schmid, K. Binder, T. Speck, and P. Virnau,“Critical behavior of active Brownian particles,” Phys. Rev. E , 030601(2018). S. K. Schnell, T. J. Vlugt, J.-M. Simon, D. Bedeaux, and S. Kjelstrup,“Thermodynamics of a small system in a µ T reservoir,” Chem. Phys. Lett. , 199–201 (2011). S. Kjelstrup, S. K. Schnell, T. J. H. Vlugt, J.-M. Simon, A. Bardow, D. Be-deaux, and T. Trinh, “Bridging scales with thermodynamics: from nano tomacro,” Adv. Nat. Sci.: Nanosci. Nanotechnol. , 023002 (2014). F. Höfling and S. Dietrich, “Finite-size corrections for the static structurefactor of a liquid slab with open boundaries,” J. Chem. Phys. , 054119(2020). R. Cortes-Huerto, K. Kremer, and R. Potestio, “Communication: Kirkwood-buff integrals in the thermodynamic limit from small-sized molecular dy-namics simulations,” J. Chem. Phys. , 141103 (2016). M. Heidari, K. Kremer, R. Potestio, and R. Cortes-Huerto, “Fluctuations,finite-size effects and the thermodynamic limit in computer simulations:Revisiting the spatial block analysis method,” Entropy , 222 (2018). S. Fritsch, S. Poblete, C. Junghans, G. Ciccotti, L. Delle Site, and K. Kremer,“Adaptive resolution molecular dynamics simulation through coupling to aninternal particle reservoir,” Phys. Rev. Lett. , 170602 (2012). L. Delle Site and M. Praprotnik, “Molecular systems with open boundaries:Theory and simulation,” Phys. Rep. , 1–56 (2017). L. Delle Site, C. Krekeler, J. Whittaker, A. Agarwal, R. Klein, and F. Höfling,“Communication: Molecular dynamics of open systems: Construction of amean-field particle reservoir,” Adv. Theory Simul. , 1900014 (2019). R. Ebrahimi Viand, F. Höfling, R. Klein, and L. Delle Site, “Communication:Theory and simulation of open systems out of equilibrium,” J. Chem. Phys. , 101102 (2020). V. Privman and M. E. Fisher, “Universal critical amplitudes in finite-sizescaling,” Phys. Rev. B , 322–327 (1984). E. Brézin and J. Zinn-Justin, “Finite size effects in phase transitions,” NuclearPhys. B , 867–893 (1985). S. K. Das, J. Horbach, and K. Binder, “Transport phenomena and micro-scopic structure in partially miscible binary fluids: A simulation study of thesymmetrical Lennard-Jones mixture,” J. Chem. Phys. , 1547 (2003). S. K. Das, J. Horbach, K. Binder, M. E. Fisher, and J. V. Sengers, “Staticand dynamic critical behavior of a symmetrical binary fluid: a computersimulation.” J. Chem. Phys. , 24506 (2006). S. Das, M. Fisher, J. Sengers, J. Horbach, and K. Binder, “Critical dynamicsin a binary fluid: Simulations and finite-size scaling,” Phys. Rev. Lett. ,025702 (2006). S. Ahmad, S. K. Das, and S. Puri, “Crossover in growth laws for phase-separating binary fluids: Molecular dynamics simulations,” Phys. Rev. E ,031140–9 (2012). S. Roy and S. K. Das, “Transport phenomena in fluids: Finite-size scalingfor critical behavior,” EPL (Europhys. Lett.) , 36001–7 (2011). S. K. Das, S. Roy, S. Majumder, and S. Ahmad, “Finite-size effects indynamics: Critical vs. coarsening phenomena,” EPL (Europhys. Lett.) ,66006–7 (2012). S. Roy, S. Dietrich, and F. Höfling, “Structure and dynamics of binary liquidmixtures near their continuous demixing transitions,” J. Chem. Phys. ,134505 (2016). S. Melchionna, G. Ciccotti, and B. Lee Holian, “Hoover NPT dynamics forsystems varying in shape and size,” Mol. Phys. , 533–544 (1993). D. Chakraborty, M. V. Gnann, D. Rings, J. Glaser, F. Otto, F. Cichos, andK. Kroy, “Generalised Einstein relation for hot Brownian motion,” EPL(Europhys. Lett.) , 60009 (2011). D. Chakraborty, “Velocity autocorrelation function of a Brownian particle,”Eur. Phys. J. B , 375–380 (2011). D. Rings, D. Chakraborty, and K. Kroy, “Rotational hot Brownian motion,”New J. Phys. , 053012 (2012). H. Stanley,
Introduction to Phase Transitions and Critical Phenomena ,International Series of Monogr (Oxford University Press, 1971). A. Pelissetto and E. Vicari, “Critical phenomena and renormalization-grouptheory,” Phys. Rep. , 549–727 (2002). Note that depending on whether T c is approached from above or below,different critical amplitudes apply, e.g., ξ + and ξ − . Our analysis of ξ ( T ) and χ ( T ) is restricted to T > T c , why we refrain from indicating the + sign atthe amplitudes to keep the notation simple. M. E. Fisher, “Correlation functions and the critical region of simple fluids,”J. Math. Phys. , 944–962 (1964). To simplify the discussion here, we have tacitly assumed an arbitrarily largesimulation box that supports the limit (cid:96) → ∞ . See Section V D for analternative. B. Kastening, “Anisotropy and universality in finite-size scaling: CriticalBinder cumulant of a two-dimensional Ising model,” Phys. Rev. E , 1–4(2013). A. Malakis, N. G. Fytas, and G. Gülpinar, “Critical Binder cumulant anduniversality: Fortuin–Kasteleyn clusters and order-parameter fluctuations,”Phys. Rev. E , 042103–8 (2014). H. W. J. Blöte, E. Luijten, and J. R. Heringa, “Ising universality in threedimensions: a Monte Carlo study,” J. Phys. A , 6289–6313 (1995). H. W. J. Blöte, L. N. Shchur, and A. L. Talapov, “The cluster processor:new results,” Int. J. Mod. Phys. C , 1137–1148 (1999). Y. C. Kim and M. E. Fisher, “Asymmetric fluid criticality. II. Finite-sizescaling for simulations,” Phys. Rev. E , 041506 (2003). J. Liu, N. B. Wilding, and E. Luijten, “Simulation of phase transitions inhighly asymmetric fluid mixtures,” Phys. Rev. Lett. , 115705 (2006). J. Zausch, P. Virnau, K. Binder, J. Horbach, and R. L. Vink, “Statics anddynamics of colloid-polymer mixtures near their critical point of phaseseparation: A computer simulation study of a continuous Asakura–Oosawamodel,” J. Chem. Phys.130