Pseudopotential-based discrete unified gas kinetic scheme for modeling multiphase fluid flows

To directly incorporate the intermolecular interaction effects into the discrete unified gas-kinetic scheme (DUGKS) for simulations of multiphase fluid flow, we developed a pseudopotential-based DUGKS by coupling the pseudopotential model that mimics the intermolecular interaction into DUGKS. Due to the flux reconstruction procedure, additional terms that break the isotropic requirements of the pseudopotential model will be introduced. To eliminate the influences of nonisotropic terms, the expression of equilibrium distribution functions is reformulated in a moment-based form. With the isotropy-preserving parameter appropriately tuned, the nonisotropic effects can be properly canceled out. The fundamental capabilities are validated by the flat interface test and the quiescent droplet test. It has been proved that the proposed pseudopotential-based DUGKS managed to produce and maintain isotropic interfaces. The isotropy-preserving property of pseudopotential-based DUGKS in transient conditions is further confirmed by the spinodal decomposition. Stability superiority of the pseudopotential-based DUGKS over the lattice Boltzmann method is also demonstrated by predicting the coexistence densities complying with the van der Waals equation of state. By directly incorporating the intermolecular interactions, the pseudopotential-based DUGKS offers a mesoscopic perspective of understanding multiphase behaviors, which could help gain fresh insights into multiphase fluid flow.

tremendous advances in computing capability, mesoscopic approaches developed upon the kinetic theory offer a penetrating perspective to comprehend the multiphase interactions. By exploring the multiphase behaviors at the mesoscopic level, the mesoscopic approaches fill the gap between the macroscopic descriptions of the multiphase dynamics and microscopic intermolecular actions [5].
Among plenty of kinetic-based mesoscopic approaches, the lattice Boltzmann (LB) method has emerged as an efficient and powerful tool for simulating a wide range of multiphase fluid flows [6][7][8][9][10]. The multiphase LB models developed in the past three decades can be generally classified into four categories: the color-gradient model [11], the phase-field model [12], the free-energy model [13], and the pseudopotential (PP) model [14]. Both the color-gradient model and the phase-field model take two sets of distribution functions, one for the interfacial property and the other for the hydrodynamic property, to depict the multiphase fluid flow. The free-energy model and the pseudopotential model, which mimic the effects of phase interactions by an additional volumetric force, employ a single set of distribution functions to describe the multiphase fluid flow. With such a treatment, the complexity of the computing program gets roughly halved compared to the program implementing the color-gradient or the phase-field model. Moreover, the mass and momentum transport process in the simulations employing the free-energy or the pseudopotential model is accomplished through the migration of identical particles described by the single set of distribution functions, which tends to be more consistent than the transport process exhibited in the simulations utilizing the color-gradient or the phase-field model, where the mass and momentum transport corresponds respectively to the migration of different particles depicted by two individual sets of distribution functions. Owing to the succinct implementation of the pseudopotential model, the pseudopotential LB method has experienced continued prosperity in a wide range of multiphase fluid flows [15][16][17][18]. Nevertheless, theoretical foundations of the pseudopotential model have remained the subject of debate since its birth. A major debating issue lies in the thermodynamic inconsistency. He and Doolean [19] first addressed the problem of the thermodynamic inconsistency and provided the simplified form of the pressure tensor induced by the pseduopotential model. Benzi et al. [20] identified the complete form of the pressure tensor at the continuum level. Later, Sbragaglia et al. [21] discovered that the continuum form of the pressure tensor does not ensure uniqueness due to the arbitrary gauge. Shan [22] further emphasized that the continuum form of the interaction force does not guarantee exact mechanical balance. The pressure tensor constructed on the discrete level should be employed to accurately predict the thermodynamic behaviors of the pseudopotential model. Kupershtokh et al. [23] managed to achieve the thermodynamic consistency by introducing a tunable interaction force. To uncover the underlying thermodynamic background, Sbragaglia and Shan [24] derived the free energy functional in terms of the pseudopotential model and established the specific expression of the interaction potential. It turns out that an implementation of the equation of state in the thermodynamic theory would inevitably result in the thermodynamic inconsistency. Thereafter, Li et al. [25] explored the mechanical stability solutions in varying conditions and figured out the appropriate parameter value for approximately fulfilling the thermodynamic inconsistency requirement. With a third-order Chapman-Enskog analysis, Huang et al. [26] proved that the classical discrete pressure tensor proposed by Shan [22] remains only conditionally correct. On the contrary, the continuum pressure tensor constructed by considering the third-order isotropic term could accurately predict the thermodynamic behaviors of the pseudopotential model. Despite the tremendous progress achieved by the pseudopotential model in multiphase fluid flow [6,27], all of the research has been confined within the framework of LB method. Considering its straightforward implementation and excellent capacity, it is tempting to develop pseudopotential-based kinetic schemes that are not restricted by the uniform Cartesian grid for efficient multiphase flow simulation.
Rooted in the kinetic theory, the discrete unified gas kinetic scheme (DUGKS) developed within the finite volume framework suffers from no restriction on grid types [28]. By considering the local Knudsen information in the construction of kinetic flux, DUGKS could accurately depict extensive fluid flows ranging from the continuum regime to the free molecular regime [29]. Over the past decade, DUGKS has demonstrated its excellent capability in modeling compressible flows [30][31][32], turbulent flows [33][34][35], solid-fluid flows [36][37][38], multicomponent gas flows [39,40], microscale gas flows [41,42], radiative heat transfer [43,44], and so forth. For the widespread application of DUGKS, readers are recommended to refer to the review literature provided by Guo and Xu [45]. Despite the tremendous progress made by DUGKS, studies centered on the multiphase fluid flows are rather limited [46,47]. Moreover, the multiphase model employed in DUGKS is generally limited to the phase-field model [48], where two sets of distribution functions are utilized to separately describe the interface and fluid dynamics. Such a treatment isolates the mass transport from the momentum transport and could induce unphysical phenomena. To avoid this undesirable feature, DUGKS using a single set of distribution functions for multiphase flow simulations should be developed. Very recently, Zeng et al. [49] proposed a well-balanced DUGKS for twophase flows by absorbing the free-energy model [50] into DUGKS. Only a single set of distribution functions has been utilized in their work for the concurrent transport of mass and momentum. Numerical results demonstrated the superior stability and accuracy of DUGKS over that of LB method. Nevertheless, the free-energy model considers the phase interactions through the chemical potential field, which typically belongs to a macroscopic description. To describe the phase interactions at the mesoscopic level, the pseudopotential model that directly mimics the intermolecular interactions is a distinct preference. Considering the excellent performance of DUGKS proved in previous studies [51,52], as well as the mesoscopic feature of the pseudopotential model, we developed a pseudopotential-based DUGKS for two-phase fluid flows by coupling the pseudopotential model with the Strang-splitting DUGKS. To simulate a realistic two-phase system, the van der Waals equation of state (vdW-EOS) is implemented for the evaluation of bulk pressure. The rest of this paper is organized as follows: Section 2 introduces the Strang-splitting DUGKS and the pseudopotential model. Section 3 presents the numerical results as well as discussions. Section 4 concludes the findings.

Numerical methodology
In this section, the macroscopic governing equations are first briefly introduced. Then we offer a detailed explanation of the Strang-splitting discrete unified gas kinetic scheme. The pseudopotential model for DUGKS will be introduced in the final part.

Macroscopic equations
The macroscopic governing equations recovered by the kinetic equation through the Chapman-Enskog theory read where t represents time, ρ indicates the fluid density, u denotes the flow velocity, p is the pressure, and µ is the dynamic viscosity. F s stands for the volumetric force that mimics the interaction effects between/among different phases whereas G indicates the gravitational or buoyant force.

Discrete unified gas kinetic scheme
In present research, the flow field is directly governed by the Boltzmann-BGK equation, which takes the form of where f = f (x, ξ , t) is the distribution function (DF) accounting for the particles residing at position x with a velocity of ξ at time t, τ is the relaxation time, f E is the equilibrium distribution function approached by f within each collision. The moments of the distribution function yield the conservative flow variables via A necessary prerequisite for the numerical evaluation of the moments is the discretization of the velocity space. In present work, the three-point Gauss-Hermite quadrature is employed to determine the discrete particle velocities along each single dimension. In two dimension the discrete velocities can be derived from the tensor product of the single dimensional velocities, which reads where ξ i is the ith discrete velocity and c s = 1/ √ 3 is the model speed of sound. To fulfill the relation of Eq. (3) at the discrete level, the equilibrium DF f E in present research is evaluated by represents the column vector constituted by the discrete equilibrium DFs, M is the transformation matrix defined as and m E signifies the macroscopic equilibria: Here α is a free parameter used to eliminate the nonisotropic effects of the scheme [26]. The relations between the conservative variables and the discrete DFs become With the discretization of velocity space, the discrete velocity Boltzmann-BGK equation takes the following form: To numerically solve Eq. (6), we first subdivide the spatial domain into a set of grid cells and integrate this equation over a certain cell, which yields where V c represents the integral cell centered at position x c , ∂V c indicates the surface bounding the cell, and n represents the unit vector normal to the surface. Integrating Eq. (7) over a time step of length t = t n+1 − t n yields where |V c | measures the volume of cell V c , f n i and n i approximate the cell averages of V c in such a way that F n+1/2 i measures the kinetic flux at the mid time t n + �t/2 by Note that the midpoint rule is utilized to compute the time integral of the kinetic flux and trapezoidal rule is employed to evaluate the time integral of the collision term in Eq. (8). To obtain a fully explicit evolution equation, two auxiliary distribution functions are introduced: Substituting Eq. (11) into Eq. (8), we have which turns out to be fully explicit.
To evaluate the kinetic flux F n+1/2 i , the primitive distribution function f i (x f , t n+1/2 ) on cell interface is needed. To this end, we integrate Eq. (6) along the characteristic line over a time step length of δt = �t/2: Here the trapezoidal rule is once again employed to evaluate the time integral of collision term. To realize the explicit treatment of Eq. (13), another two auxiliary distribution functions are introduced as follows: Eq. (13) then can be rearranged as The cell-centered auxiliary distribution function f + can be constructed according to its definition: The value of f + (x f − ξ i δt, t n ) can be interpolated from the corresponding cell-centered distribution function [53]. For the face-based reconstruction scheme (FRS), For the cell-based reconstruction scheme (CRS), f + (x f − ξ i δt, t n ) can be computed by where x L and x R correspond respectively to the center positions of the two cells adjacent to the interface located at x f . Once the value of f + is known, the primitive DF at cell interface can be updated by (11) Thereafter the kinetic flux F n+1/2 can be evaluated according to Eq. (10) and the auxiliary distribution function at time t n+1 can be updated via Eq. (12). The cell-averaged primitive DF then can be obtained according to To obtain the primitive DFs in Eq. (19) and (20), the value of equilibrium DFs that depend on the conservative variables should be first determined via Eq. (4). With the help of Eq. (5), the corresponding conservative variables can be evaluated by on cell interfaces and by at cell centers.
To date, the evolution process of DUGKS without considering force effects has been basically explained. To incorporate the influence of external forces, another discrete distribution function f S i accounting for the force effects should be introduced: To correctly recover the macroscopic equations, the moments of discrete force DF should obey where F = F s + G is the external force in total. In present research the force DF f S is evaluated by represents the column vector constituted by the discrete force DFs, M is the identical transformation matrix appeared in Eq. (4), and m S signifies the macroscopic force terms expressed as To circumvent the force effects on the interface flux, the Strang-splitting scheme is employed to evaluate the force influences. With such a treatment, the force effects are incorporated before and after the DUGKS procedure in a way that (19) As Eq. (26b) remains identical to Eq. (6), it can be solved by the DUGKS procedure addressed previously. Eq. (26a) and (26c) can be numerically solved by the Euler forward method: The conservative variables should be accordingly updated via

Pseudopotential multiphase model
In the pseudopotential multiphase model, the interaction effects between/among different phases are mimicked by a volumetric force defined as where ψ represents the interaction potential, G indicates the interaction strength, ω stands for the weights, x ′ i denotes the nearby position that is related to x by x ′ i = x + ξ i δ ′ t , among which ξ i is the ith discrete velocity and δ ′ t is the transporting time. A utilization of nine discrete velocity points leads to the following relation: In fact, the role of expression (29) is equivalent to evaluating the gradient of ψ through an isotropic finite-difference scheme [54]. A Taylor expansion of Eq. (29) gives where δ x = ξ x δ ′ t = ξ y δ ′ t = 1 measures the grid spacing. To analytically derive the pressure tensor, Eq. (31) could be reformulated as [20] (26a) where I represents the identity matrix. However, Sbragaglia et al. [21] demonstrate that the transformation of Eq. (31) into Eq. (32) does not necessarily guarantee uniqueness. As a matter of fact, Eq. (31) can be reformulated as providing the prefactors satisfy It can be identified that Eq. (32) represents a special case of Eq. (33) at a 1 = −1, a 2 = 0, a 3 = 1/2, a 4 = 1.
The continuum pressure tensor is defined as [26] where c s stands for the model speed of sound and S represents the additional term introduced from the discretization of DUGKS. Due to the reconstruction approaches utilized, the additional term S contributed from DUGKS lacks of isotropy. To balance the anisotropic influences, a free parameter α has been introduced in Eq. (4). As the discretization approaches utilized in DUGKS appear to be complex, it is quite difficult to obtain the general expression of P . Nevertheless, the normal pressure P n in such a condition could be similarly postulated as [55] where n denotes the direction normal to the interface. The normal component of the pressure tensor P n at the equilibrium state should be equal to the bulk pressure p 0 [22]. The mechanical stability condition can then be obtained as [8] where ψ ′ = dψ/dρ , ǫ = k 1 /k 2 , p 0 = p(ρ l ) = p(ρ g ) denotes the bulk pressure and p EOS represents the non-ideal equation of state (EOS) in terms of the pseudopotential model: Providing the coexistence densities ( ρ l and ρ g ) have been estimated by DUGKS, the value of the produced parameter ǫ can then be determined numerically [25]. To consider the effects of various non-ideal equations of state, He and Doolen [19] pointed out that the pseudopotential ψ should be evaluated as where p EOS should be one of the non-ideal equations of state in the thermodynamic theory. In present research, the dimensionless van der Waals equation of state (vdW-EOS) expressed as is employed, with a = 9/392 and b = 2/21 . The critical density ρ c and temperature T c hold the value of 7/2 and 1/14 in such a condition [56]. With the incorporation of the dimensionless vdW-EOS, the pseudopotential ψ could be directly calculated through Eq. (39).

Numerical tests
The capacity of pseudopotential-based DUGKS is validated by three benchmark tests. Firstly, the fundamental capability to predict coexistence densities is verified by the flat interface test. Subsequently, the nonisotropic property is investigated by the quiescent droplet test and the isotropy-preserving parameter α is tuned to cancel out the nonisotropic effects. Finally, the spinodal decomposition test is conducted to examine its isotropy-preserving property in transient conditions. For steady tests, the computing process terminates when the L 2 -norm-based velocity error E(u) falls below a critical value e: where u denotes the velocity field, t n−1000 indicates the moment 1000 time steps ahead of t n and e is given as 10 −8 if not otherwise specified. The CFL number is fixed at 0.8 across all the tests.

Flat interface
The flat interface problem serves as a fundamental benchmark for validating the basic capability of newly proposed multiphase models. Consider an infinitely long horizontal channel filled with binary fluids of different phases. The liquid resides in the middle half of the channel while the gas occupies the upper and bottom region. The computational domain is confined to a 16×256 rectangular region uniformly divided into Cartesian grids, where the grid spacing holds a fixed value of unity. Periodical boundary conditions are applied to all sides. The density field is initialized by where ρ l and ρ g correspond to the liquid and gas densities, y L and y H represent the lower and upper bound of the fluid in liquid phase, and W denotes the interface thickness. Figure 1 presents the coexisting densities produced respectively by DUGKS and LBM at τ = 0.3, α = 1.0. It can be identified that the results provided by DUGKS agree well with the theoretical solutions obtained through the Maxwell equal-area law [27] while the results offered by LBM deviate significantly from the theoretical solutions. Moreover, DUGKS operates properly in conditions of low temperatures whereas LBM fails to work when T r < 0.75 . The inferior results produced by the standard LBM could be attributed to the superfluous terms recovered in the momentum equation [56]. The stability superiority of DUGKS might result from the coupling of transport and collision process. The influences of isotropy-preserving parameter α on the coexisting densities are investigated and the corresponding results are presented in Fig. 2. Numerical results remain unchanged despite the varying parameter α . As the density distributions in the horizontal direction stay unaltered, which has no impact on the fluid behavior, it is reasonable to obtain identical results with varying α . The density profiles along the vertical direction produced by DUGKS at different temperatures are illustrated in Fig. 3. With the increasing temperature, the flat interface gets sharpened, which is due to the strong coupling of physical properties originating from the pseudopotential model [27]. The varying interface thickness also suggests that W in Eq. (42) has no concrete meaning and remains useful only during the initialization process.

Quiescent droplet
The quiescent droplet provides another fundamental benchmark for validating the model's capability. With this test, we specially investigated the isotropic property of pseudopotential-based DUGKS. Initially, a circular droplet is placed at the center of a L 0 ×L 0 square domain according to where ρ l and ρ g correspond respectively to the liquid and gas density, (x c , y c ) represents the center position of the square domain, and R d denotes the droplet radius. The parameters in this test are set as L 0 = 256 , R d = 0.24L 0 , W = 5 . The computing process terminates once the L 2 -norm-based velocity error meets the relation defined by Eq. (41). Figure 4 shows the density distribution at the initial moment. The ultimate density contours produced by DUGKS using different reconstruction schemes are illustrated in Fig. 5. The square droplet presented in Fig. 5a is obtained with the second-order face-based reconstruction scheme, i.e., the central-difference scheme [28]. The nearly circular droplet shown in Fig. 5b is produced by the second-order cell-based reconstruction scheme, also known as the upwind scheme [57]. As a matter of fact, it is Dr. Wang Peng [51] who first identified the nonisotropic property of the pseudopotential-based DUGKS. In the spring of 2018, Wang and Zhu [57] dropped in at NWPU and gave a brief presentation themed around DUGKS. During a casual conversation, we talked about the nonisotropic property of the pseudopotential-based DUGKS. Wang suggests that the upwind reconstruction approach, together with the Strang-splitting scheme, should be employed to obtain an isotropic interface. Following this idea, we conducted a few simulations with the pseudopotential-based DUGKS. It turned out that although the nonisotropic problem could be relieved to some extent, the pseudopotential-based DUGKS still fails to produce a perfectly isotropic interface, which has been demonstrated by the interface profile presented in Fig. 5b. The droplet produced with a third-order accuracy of CRS is illustrated in Fig. 5c. It can be observed that employing a high-order upwind Through a third-order Chapman-Enskog analysis, Huang et al. [26] identified the nonisotropic and isotropic terms in the pseudopotential model. Numerical results proved that the free parameter α played a key role in preserving the isotropic property. With this proof, we introduced the isotropy-preserving parameter α in the calculation of the equilibrium distribution function. By adjusting α to an appropriate value, a perfectly isotropic interface can be produced and maintained. In the condition of τ = 0.3, T r = 0.95, the corresponding value of α is 1.304. The circular droplet produced is illustrated in Fig. 5d. In practical computations, the density criterion and the velocity criterion have been compared to identify the optimal criterion for the determination of α . The density criterion compares the final density field with the initial density field and recognizes the value of α that creates the minimum density difference as the most appropriate choice. The velocity criterion considers the value of α that generates the minimum L 2 norm of velocity as the most appropriate choice. An optimal criterion should help maintain the quiescent interface as isotropic as possible. The isotropic level of an interface is assessed by comparing the horizontal density profile with the diagonal density profile produced at a final moment. Through a series of comparisons, it has been identified that the velocity criterion best determines the value of α . Figure 6 illustrates the density contour produced by DUGKS on different criteria. It can be observed that the interface generated under the density criterion suffers from a nonisotropic problem while the interface created under the velocity criterion maintains excellent isotropy. Numerical experiments also revealed that the appropriate value of α varies along with changes in the relaxation time τ or the reduced temperature T r . Through repeated calibration, we ascertained the specific values of α for preserving an isotropic interface. Figure 7 illustrates the ternary relation diagram for α , τ and T r . Table 1 presents the detailed data. With this information, it is now possible to interpolate the specific value of α for a wide range of T r and τ . Considering the situation of τ = 0.63 and T r = 0.72 , the most appropriate value of α determined by numerical simulations is 1.025 while the interpolated value of α in such a condition is 1.02494, which is pretty close to the calibrated value. Figure 8 compares the density profiles extracted from various directions. The density profile along the horizontal direction coincides completely with the profile along the vertical direction. The profile along the diagonal direction deviates slightly from them, which could be partially attributed to the additional error introduced during data extraction. Nevertheless, the nearly indistinguishable deviation demonstrates the well-maintained isotropic property of the interface. It is worth mentioning that we failed to find an analytical expression to describe the mapping relation between τ , T r , and α , which might be beneficial for the application of pseudopotential-based DUGKS. To quantitatively examine the accuracy of the current scheme, we further validated the Laplace's law by simulating a number of quiescent droplets with varying radii. Figure 9 illustrates the relation between the pressure jump P and the reciprocal of radius. Obvious linear relation could be identified from the results produced by DUGKS with different temperatures, which is in accordance with the Laplace's law: �P = σ/R d . Due to the strong coupling effects of the pseudopotential model, it is generally difficult to determine the theoretical value of surface tension coefficient  . 7 Variation of isotropy-preserving parameter α with regard to relaxation time τ and reduced temperature T r σ . For the same reason, the droplet radius obtained at the final moment could deviate slightly from the initialization value. Nevertheless, the linear relation reflected by the numerical results could surely demonstrate the fundamental capability of pseudopotential-based DUGKS.

Spinodal decomposition
The previous two benchmark tests belong to steady-state problems. To examine the capacity of DUGKS in dealing with transient problems, spinodal decomposition tests are conducted. Spinodal decomposition, also referred to as phase separation, occurs when a homogeneous mixture contains compositional fluctuations. The computational region is a L 0 ×L 0 square domain, with the side length L 0 set to 512. The relaxation time τ is fixed at 0.3 and the reduced temperature T r is given as 0.8. The density field is initialized according to where random(0, 1) generates the fluctuations by returning a random number between 0 and 1. Usually the evolution time should be scaled by a reference time relating to the surface tension coefficient σ . Due to the strong coupling effects of the pseudopotential model, it is difficult to determine the reference time. Here the simulation time t is directly used to indicate the time evolution. Figures 10, 11 (44) ρ(x, y) = (ρ l + ρ g )/3 + random(0, 1)/100, DUGKS at τ = 0.3 , T r = 0.8 . The complete separation process is successfully predicted and no instability phenomenon has ever been detected. In the early stages depicted by Fig. 10, the tiny fluctuations evolve into local inhomogeneities that initialize the phase  separation. As the system evolves, the inhomogeneities induce the nucleation of heavy liquid, which can be observed in Fig. 11. With the continual development, interfaces separating different phases can be clearly detected in Fig. 12. Then the small droplets  keep coalescing into large droplets. Eventually, a single quiescent droplet illustrated in Fig. 19 is formed. At the final moment, the interface produced with α = 1.035 appears to be isotropic while the interface produced with α = 1.0 suffers from a lack of isotropy. The same phenomenon can also be detected in the process of droplet coalescence depicted by Figs. 14, 15, 16, 17 and 18. This fact demonstrates the effectiveness of the isotropy-preserving property in a transient condition. The separation process produced with α = 1.0 deviates slightly from the corresponding process built with α = 1.305 , which could be partially attributed to the differences in initial fluctuations.

Conclusion
A pseudopotential-based discrete unified gas-kinetic scheme for multiphase fluid flows is proposed by coupling the pseudopotential model into the Strang-splitting DUGKS. Due to the strict requirements of the scheme isotropy by the pseudopotential model, a direct coupling of pseudopotential model into DUGKS could not maintain an isotropic interface. To cancel out the nonisotropic terms introduced during the flux reconstruction process, the equilibrium distribution functions are expressed in a moment-based form and an isotropy-preserving parameter α is introduced. By adjusting this parameter to the appropriate value, pseudopotential-based DUGKS managed to produce and maintain an isotropic interface. With a sequence of numerical experiments, we calibrate the value of α to facilitate DUGKS simulating isotropic interfaces in a wide range of conditions. Results produced in the spinodal decomposition tests further demonstrate the effectiveness of its isotropy-preserving property in the transient state. Comparative results of coexistence curves also proved the superior stability of DUGKS over that of LBM. The fundamental capacity of pseudopotential-based DUGKS has been demonstrated by the basic benchmark tests. Further investigations are needed to explore its capacity in specific fields. (56) f (x + �x, t n + �t) =f n (x) + 2�t 2τ + �t f E,n (x) −f n (x) + 2τ �t 2τ + �t f S,n (x),