Dielectric constant of the polarizable dipolar hard sphere fluid studied by Monte Carlo simulation and theories

A systematic Monte Carlo (MC) simulation and perturbation theoretical (PT) study is reported for the dielectric constant of the polarizable dipolar hard sphere (PDHS) fluid. We take the polarizability of the molecules into account in two different ways. In a continuum approach we place the permanent dipole of the molecule into a sphere of dielectric constant ∞ in the spirit of Onsager. The high frequency dielectric constant ∞ is calculated from the Clausius-Mosotti relation, while the dielectric constant of the polarizable fluid is obtained from the Kirkwood-Fröhlich equation. In the molecular approach, the polarizability is built into the model on the molecular level, which makes the interactions nonpairwise additive. Here we use Wertheim’s renormalized PT method to calculate the induced dipole moment, while the dielectric constant is calculated from our recently introduced formula [22]. We also apply a series expansion for the dielectric constant both in the continuum and the molecular approach. These series expansions ensure a better agreement with simulation results. The agreement between our MC data and the PT results in the molecular approach is excellent for low to moderate dipole moments and polarizabilities. At stronger dipolar interactions ergodicity problems and anizotropic behaviour appear where simulation results become uncertain and the theoretical approach becomes invalid.


Introduction
The classical theories of the dielectric constant founded by Kirkwood [1], Onsager [2], and Debye [3] use a continuum approach: they place the molecule in a cavity surrounded by the material treated as a continuum.The Clausius-Mosotti (CM) equation is valid for apolar molecules, while the Debye equation [3] holds good approximately to gases and dilute solutions of molecules carrying a permanent dipole.The Debye equation for a spherical shaped sample is: where α is the polarizability of the molecule, ρ = N/V is the number density, and y 0 is the dimensionless dipole strength function, y 0 = 4πρµ 2 0 /9kT , where µ 0 is the permanent dipole moment, k is the Boltzmann factor and T is the temperature (for µ 0 = 0, equation (1) yields the CM equation).Onsager [2] placed a point dipole in the centre of a cavity of dielectric constant ∞ and the effect of the surrounding dielectric is measured by the dielectric response of the polarization charges induced on the wall of the cavity.The resulting equation is where ∞ is the high frequency dielectric constant, which is related to the molecular polarizability via the CM equation The Onsager equation works quite well for liquids when the dipole moment is not too high.
When the molecules have a large permanent dipole moment, the correlation between a central dipole and the surrounding dipoles cannot be replaced by the response of a continuum dielectric.Kirkwood has introduced the so-called Kirkwood g-factor to take into account these short range correlations.The g-factor is obtained from the fluctuation of the total dipole moment of the system M, where • • • denotes ensemble average.For a liquid consisting of non-polarizable molecules, the Kirkwood-equation reads as To take into account the molecular polarizability, we can choose the continuum approach used by Onsager.In this case we obtain the Kirkwood-Fröhlich (KF) equation [4] ( It is seen that the difference between this equation and the Osager equation, (2), is that the Kirkwood factor is present on the right hand side representing the correlation between the permanent dipoles of the molecules.Nevertheless, this is still a continuum theory regarding polarizability, and y 0 is calculated by using the permanent dipole moment µ 0 (that is why the subscript 0 is used).When a molecular polarizability is added to the permanent dipole of a dipolar molecule, we take into account the polarizability on the molecular level: we call this the molecular approach.The equation that applies for this case is In this equation y is calculated by using the total (permanent plus induced) dipole moment of the molecules (subsript 0 is absent) and g K contains the correlation between permanent and induced dipoles too.
To calculate g K , a molecular model and a statistical mechanical method are needed to study it.In the case of the Kirkwood, (5) and the KF, (6), equations, the calculation of g K is relatively simple because the pair-wise additive dipole-dipole interaction is considered and the calculation of pair correlation functions is sufficient.Nevertheless, when the molecular polarizability is included in the model, the intermolecular potential is no longer pairwise additive and many-body correlations have to be taken into account.This is the characteristic problem that greatly complicates the statistical mechanical treatment of the systems containing strongly interacting polarizable dipoles.
Most theories rather avoid the treatment of these many-body correlations by mimicing the polarizable fluid with a system characterized by a pairwise additive interaction applying an effective permanent dipole moment.This approach was used by the graph-theoretical treatment of Wertheim [5], its extension to dipolarquadrupolar fluids and mixtures [6][7][8], and a self consistent mean field (SCMF) approach of Carnie and Patey [9].In these approaches, a theory is needed to calculate the properties of the system defined by the effective pair potential.For this purpose, one can use the hypernetted chain theories [9,10] or the thermodynamic perturbation theory (PT) which was developed by many authors during the second half of the twentieth century [11][12][13][14][15][16][17][18][19][20][21].A major breakthrough was the theory of  which was the first useful PT.The PT was used in our previous paper [22] where we proposed an equation for the dielectric constant of polarizable polar fluids.We use this equation to study the dielectric constant of the polarizable dipolar hard sphere (PDHS) fluid in this work.
In computer simulations, the treatment of non-additive potentials is straightforward since the pioneering works by Vesely [23,24].The problem here lies in the time consuming iteration procedure calculating the induced dipoles in every simulation step.After a few early works [10,[23][24][25][26][27], more and more computer simulation studies of polarizable fluids have been appearing in the literature due to the continously increasing speed of computers [28][29][30][31][32][33][34][35].While the earlier works dealt with simpler models such as PDHS or polarizable Stockmayer (PSTM) fluids, the latter works concentrated on more complex molecular liquid models aimed at representing real substances, first of all, water.Lately, the authors seem to have lost interest in simple potential models.This is probably due to the increasing power of computers.
Thus computer simulation became an everyday tool and investigators are capable of studying more and more sophisticated models for a specific material.Nevertheless, this specification results in a loss of generality meaning that simple models can capture a few important, characteristic features of a whole class of systems, while detailed models are restricted to the molecule for which they have been designed.
Notable exceptions in this trend of specification are the works of Kriebel and Winkelmann [36,37] who gave a systematic study of the thermodynamic properties of the PSTM [36] and the polarizable dipolar two-centred Lennard-Jones (LJ) fluid [37].Nevertheless, they did not have interest in the dielectric properties of these systems.In an earlier work [38], we studied the dielectric constant of this system using Monte Carlo (MC) simulaton and our renormalized PT equation proposed previously [22].
The PDHS model, which we consider in this work, is probably the simplest molecular model for a polarizable dipolar fluid.The core is a hard sphere (HS) which gives the molecules finite size, so contributions from excluded volume effects can be studied.Using this simple model for the core, we can concentrate on the effect of electrostatic interactions.In this study, we use both the continuum approach (equation ( 6)) and the molecular approach (equation ( 7)) to calculate the dielectric constant of the PDHS fluid over a wide range of density, dipole moment, and polarizability.The results obtained from the theories are compared to our own MC simulation data that are obtained using the pair approximation of polarization interaction (PAPI) method of Předota et al. [35] in order to accelerate calculations.
It is known since the work of Tani et al. [39] for the dipolar hard sphere (DHS) fluid that equation (5) gives poor results in the liquid phase for the dielectric constant of the DHS fluid using PT to calculate g K compared to simulation data.Nevertheless, if we expand (y 0 ) into a Taylor series as a function of the dipole strength function y 0 , we obtain much more reasonable results.We detail the nature and importance of this series expansion in subsequent sections and we develop such series expanions for the case of polarizable fluids using both the continuum and the molecular approach.We will show that these series expansions give results in better agreement with MC data than the original equations, ( 6) and (7) do.

Models
First, we introduce the PDHS model, of which the DHS fluid is a special case with α = 0. Since the potential in the PDHS system is not pair-wise additive, we have to define the total energy of the system for a complete definition of the Hamiltonian.In the case of the DHS fluid, of course, the system can be defined merely by giving the pair potential; the total energy is the sum of the pair potentials.
The total energy of the N-particle system of the PDHS fluid can be divided into three parts: The first part, the non-electrostatic spherical reference potential is a sum of pair-wise additive potentials which is the HS potential in this study: where σ is the hard sphere diameter and r ij is the magnitude of the vector r ij = r i − r j .This term merely forbids the overlapping of the spheres.The electrostatic terms of the energy consist of the interaction between the permanent dipoles and the induction energy [31] where (µ 0 ) i and (µ ind ) i are the permanent and the induced dipole moments on the i-th molecule, respectively, and (E 0 ) i is the electric field raised at the place of the i-th molecule by the permanent dipoles of all other molecules: with being the electrostatic tensor.Introducing the electric field raised at the place of molecule i by the induced dipoles of all other molecules the induced dipole moment of the i-th molecule can be written as In general, α is a tensor, but in this study, we restrict ourselves to a scalar polarizability.For the non-polarizable case U ind = 0, and U el reduces to the sum of the dipole-dipole pair potentials where (ij) is a notation for (r ij , n i , n j ) and n i = µ i /µ 0 is the unit vector in the direction of the i-th dipole.

Monte Carlo simulations
Computer simulations provide an efficient method to study many-particle thermodynamic systems.For a well defined system, the results of simulations are accepted as a gold standard to which the results of various theories, which usually contain some approximation, can be compared.By well defined system we mean that the interaction potential acting between the particles together with a few thermodynamic parameters that define the thermodynamic state are given.The details of both MC and molecular dynamics simulations are given in excellent books [40][41][42], here we only describe some aspects regarding the calculation of induced dipoles and the dielectric constant.
Normally, in an MC step only one molecule is displaced.In the case of the nonadditive polarizable dipole potential a self consistent iteration procedure based on equations ( 14) and ( 15) is needed to calculate (µ ind ) i for a given configuration.The values of the induced dipoles are needed to calculate the energy from equations (10) and (11).The computer time cost of the iteration procedure is proportional to N 2 .Therefore, MC simulation seems to be a less economic method for non-additive potentials if single particle moves and full iteration are used.
To decrease computational burden of MC simulations, Předota et al. [35] have suggested a procedure called "pair approximation for polarization interaction" (PA-PI).The basis of their idea is that after a small displacement of a single particle, a full update of the induced dipoles is necessary only in the close neighbourhood of the displaced molecule i, namely, close to the source of the disturbance that causes the reordering of the induced dipoles.The neighbourhood for which full iteration is performed is defined as two spheres around the old and the new position of molecule i with a radius R iter .The change of the induced dipole moment in an iteration is calculated from where the change of the induced electric field in the iteration is given by This means that the interactions between the displaced molecule i and an arbitrary molecule j are always taken into account in the calculation of ∆(E ind ) j .Furthermore, all interactions between pairs of molecules j, j that are within the cutoff spheres are also taken into account, while interactions between pairs where one of the two molecules is outside the spheres are ignored.The change of the electric field produced by the permanent dipoles, ∆(E 0 ) j , is nonzero only in the first iteration.
The dielectric constant can be determined from the following formula [26,31] ( This equation assumes that a spherical sample of polar molecules is surrounded by a uniform dielectric of dielectric constant RF .The notation RF refers to the reaction field produced by the polarization charges induced on the boundary of the sample and the surrounding dielectric.The high frequency dielectric constant, ∞ , can be calculated from the CM relationship (equation ( 3)).Nevertheless, replacing the first term on the right hand side by 4παρ/3 practically yields the same result for .The above formulation for the dielectric constant has to be consistent with the energy calculation in the simulations, namely, the term has to be added to the dielectric tensor in equation ( 13), where R c is the cutoff radius of the dipolar interactions (taken to be equal to the half of the width of the simulation box, L/2).Without using any summation technique to estimate contributions of periodic images to the potential (such as the Ewald technique), this procedure is the so-called reaction field method to treat long range corrections [43].In this work, the tin-foil or conducting boundary condition ( RF → ∞) was used.
The number of particles was 256 for most cases, 108 particles were used in a few cases for low density (ρ * = ρσ 3 = 0.2), and a few simulations have been performed with 512 and 1000 particles for comparison.The lengths of the simulations varied between 100 and 500 thousand MC cycles depending on the dipolar strengths and density (longer runs were used for higher densities and higher dipole moments).In an MC cycle N attempts were made to move a particle.In an MC movement, the molecule was attempted to be displaced and rotated randomly with respect to its old coordinates.The maximum displacements had to be kept at low values in order to maintain a fast convergence in the calculation of the induced dipoles in the PAPI procedure.The radius R iter of the PAPI technique was set to 2.5σ for higher densities (ρ * = 0.6 and 0.8) and 3σ for lower densities (ρ * = 0.2 and 0.4).A full iteration was made in every tenth MC cycle.The statistical errors were estimated from the block average method.

Theories
There are two basic problems that we expect from our theories to provide solution for.The first problem is the calculation of the Kirkwood factor for a dipolar fluid.The thermodynamic PT that determines it for the non-polarizable DHS fluid is described in the next subsection.The other problem is the consideration of the molecular polarizability.A continuum and a molecular approach, already loosely desribed in the Introduction, are presented in the subsequent subsections.

Non-polarizable case
To calculate the dielectric constant of the DHS fluid, the starting point is the Kirkwood equation, (5).The Kirkwood g-factor (equation ( 4)) takes care of the short range correlations between the dipoles.It can be related to the ensemble average of the function ∆(1j) = n 1 • n j by where • • • denotes ensemble average which can be determined using the paircorrelation function g (12).The notation • • • ω 1 ω 2 means averaging over orientations.
Perturbation theories [11][12][13][14][15][16][17][18][19][20][21] divide the intermolecular potential into an isotropic reference state and a perturbation part (u HS is the reference potential and u DD is the perturbation term).For the pair correlation function of the DHS system the perturbation series expansion read as where g 0 (12) and g 0 (123) are the two-and the three-particle correlation functions of the reference (HS) system, respectively, and β = 1/kT .The three-particle distribution function has been calculated by the Kirkwood superposition approximation: g 0 (123) = g 0 (12)g 0 (23)g 0 (13) as suggested by Barker et al. [15][16][17].Performing the angular integrals and using orthogonality relations the Kirkwood factor can be expressed as [39,[44][45][46][47]] where I dd∆ denotes a triple integral that can be given as where θ 3 is the interior angle of the (123) triangle.For the DHS fluid, the values of the integral depend only on the dimensionless reduced density ρ * = ρσ 3 through g 0 (123).The ρ * -dependence can be given by the following polynomial: This expansion was obtained by a fit to values that we calculated with the Fourier transform convolution theorem [47] and it provides results that are virtually identical to those given by the Padè type expression reported by Tani et al. (equation ( 17) in [39]).Substituting equation ( 23) into equation ( 5), the Kirkwood equation can be rewritten as The results calculated from equation (26) are not in good agreement with the simulation data, so Tani et al. [39] expanded (y 0 ) as a Taylor series on the basis of the above equation and obtained that Tani et al. have shown that the results from the above series expansion of the Kirkwood equation is in better agreement with the simulation data than those calculated from the Kirkwood equation itself, (26).Goldman [47] has found the same for the Stockmayer fluid.These series expansions, apart from the fact that they give better agreement with simulation results, have a strong theoretical basis.Jepsen [48,49] and Rushbrook [45,50] have shown for the DHS fluid that the next exact term in the y 0 -expansion of the Debye equation is from which the following series can be derived for the dielectric constant: This equation is exact in the limit of low densities.The Onsager [2] and the van Vleck [51] theories are valid up to the second order in y 0 .The y 0 -expansion of the dielectric constant given by Wertheim's MSA method [52] yields this equation too.
Note that in the low density limit I dd∆ = 17π 2 /9 and equation ( 27) becomes identical to equation (29).These considerations support that such a series expansion is not merely a computational trick but it is an equation that is independent of the boundary conditions and the shape of the sample.Equations written in a form of the Kirkwood, the Onsager, or the CM equations explicitly include the information about the boundary conditions (left hand side of the equations).Accordingly, the Kirkwood factor on the right hand side of the equations will also depend on the boundary conditions because the resulting dielectric constant should be independent of them.Nevertheless, in theory the bulk system is considered as infinite and the I dd∆ integral includes information only on the short range correlations of a central dipole and its local environment and knows nothing about boundary condition of the macroscopic sample.Therefore, it is reasonable to express the dielectric constant in a form that is independent of boundary conditions.
It is important to note that equation ( 27) can be obtained using a different route.Szalai et al. [53] studied the DHS fluid in the presence of an external electric field.Assuming that the shape of the sample is an infinitely prolate ellipsoid (which means the field inside the dielectric equals the external field), they have expressed the field dependent Helmholtz free energy of the system, F (E 0 ).The relation between the polarization and the free energy is from which the dielectric constant follows.The resulting equation for is the same as equation (27).The approach based on the response to an external field can be applied within the framework of the MSA [54,55] and for the calculation of the magnetic susceptibility of magnetic fluids using simulations [56][57][58][59].

The continuum approach
In the above section the polarizability α of the molecules was not taken into account.In our previous works [60,61] we have proposed a simple continuum approach based on the KF equation, (6).If we express the Kirkwood factor using the above outlined PT approach, we obtain the following equation The high frequency dielectric constant can be obtained from the CM equation ( 3).Following the reasoning given in the previous section, we can assume that the series expansion of the above equation might also give a more general equation for the dielectric constant.Performing the same series expansion for the KF equation as Tani et al. [39] did for the Kirkwood equation, we obtain [60] that Note that with ∞ = 1 equation (32) yields equation (27).

The molecular approach
If we want to treat the molecular polarizability in a self consistent way on the molecular level, we have to use a theory that is capable of treating the non-additive many-body interactions.Wertheim [5] developed a graph-theory in which a renormalized effective dipole moment µ is defined as a sum of the permanent dipole moment µ 0 and the induced dipole moment resulting from the mean electric field E acting on the molecules: Gray et al. [6][7][8] showed that the mean electric field can be given as: where f is the one-particle Helmholtz free energy of the system.The above two equations provide a self consistent iteration route that can be solved for the effective dipole moment.For the free energy, a Padé approximant can be given with free energy perturbation terms that can be calculated if one knows the g 0 radial distribution function of the reference HS system.The exact forms of these perturbation integrals can be found elsewhere [21,36].
After the renormalization iteration procedure we obtain the effective dipole moment as a function of α and µ 0 .Equation ( 7) is used as a starting point for the calculation of the dielectric constant.In our previous work [22], we have given an expression for yg K : Introducing the parameters a = 4π 3 and equation ( 7) can be expressed as Performing the same series expansion as in the case of non-polarizable fluids [39], but now as a function of ρ, the following equation is obtained for the dielectric constant (ρ

Results and discussion
The details of our results for the PDHS fluid are shown in table 1.In the case of the PDHS system, the dipole moment is reduced with kT , namely, (µ * ) 2 = µ 2 /(kT σ 3 ) and m 2 = M 2 /(NkT σ 3 ).The density and the molecular polarizability are reduced by σ 3 : ρ * = ρσ 3 and α * = α/σ 3 .Since we have three independent parameters, ρ * , (µ * 0 ) 2 , and α * , we can give a graphical interpretation (and a comparison with the theories) in various ways.
First, we show the results as a function of density for various fixed values of α * and (µ * 0 ) 2 in figure 1. Figures 1a and 1b show the data for a low polarizability α * = 0.02.For the moderate dipole moment (µ * 0 ) 2 = 2 (figure 1a), the agreement with the PDHS2 theory is very good even for the high density ρ * = 0.8.For the higher dipole moment (µ * 0 ) 2 = 3 (figure 1b), deviations appear at high density but the agreement is still good for intermediate densities ρ * 0.6.The same can be said about higher polarizability α * = 0.04 and moderate dipole moment (µ * 0 ) 2 = 2 (figure 1c).For a high polarizability α * = 0.06 and small dipole moment (µ * 0 ) 2 = 1 (figure 1d), the agreement with the expanded version of the molecular approach (PDHS2) is excellent.These results show that the renormalization part of the PT works quite well if the dipole moment is not too large (figure 1d).The perturbation part of the potential is the dipole-dipole potential, so the PT works better if the dipole moment is smaller.As far as the other theories are concerned, all curves obtained from the non-expanded equations (DHS1, KF1, and PDHS1) overestimate the MC data.The DHS2 curves underestimate the MC results because they are calculated with α * = 0.The continuum approach (KF2) also gives too small results compared to simulation and PDHS2.For larger values of α * and (µ * 0 ) 2 that are not shown in the figure, the PDHS2 theory gives less reasonable results.
If we plot the results against the dipole moment (figure 2), we can look at the problem from a different point of view.Figure 2a (low density and low polarizability) shows that the theories using series expansion (DHS2, KF2, and PDHS2) fail at high dipole moments, while the equations that do not use series expansion (DHS1, KF1, and PDHS1) seem to work better.From this we might draw the probably wrong conclusion that the PDHS1 theory works better for low density than the PDHS2 theory.It is more probable that at high dipole moments (near (µ * 0 ) 2 = 3) the whole PT approach fails.Indeed, if we increase the polarizability to α * = 0.04 (figure 2b), even the PDHS1 theory underestimates the MC data.Nevertheless, at this low density (ρ * 0.2), the values of the dielectric constant are quite small, and there is no big difference between the continuum and the molecular approach.If we consider a higher density (ρ * 0.6), the situation is different.Due to the different scale of the y-axis, a much better agreement is shown between the PDHS2 and the MC data for low polarizability (figure 2c).Increasing α * , the agreement becomes poorer (figure 2d).
Figure 3 is the one from which the most interesting conclusions can be drawn.This figure shows the results as a function of the polarizability.It is clearly seen that the DHS1 and DHS2 theories fail to account for the increase in the dielectric constant due to the increase of the polarizability of the molecules.At α * = 0, the three theories become equivalent.
Figures 3a and 3b shows results for a very small dipole moment ((µ * 0 ) 2 = 0.5).Surprisingly, especially for the low density case (ρ * 0.2), the continuum approach (KF2) gives better results than the molecular approach (PDHS2).The question arises whether it is a coincidence or it has a theoretical basis.We argue that the latter is the case.Since the density and the dipolar strength are low, the dielectric constant has small values: the molecules are weakly polar.In this regime, it would not be surprising if the continuum approach gave reasonable results.However, since the PDHS2 approach uses a renormalization procedure, it is possible that it overestimates the effective dipole moment which results in an overestimation of the dielectric constant.Increasing the density, the molecular approach becomes better (figure 3b).
If we increase the dipole moment to the higher, but still low, value of (µ * 0 ) 2 = 1, we can see that PDHS2 theory is superior over the others even for high densities (figures 3c-d).Figures 3c-d show that the theories that do not use series expansion (DHS1, PDHS1, and KF1) overestimate the relative permittivity even if α * = 0.They describe the α * -dependence of the dielectric constant similarly as their series expanded counterparts.Increasing the dipole moment even further ((µ * 0 ) 2 = 2), the failure of the DHS1, PDHS1,and KF1 theories becomes more apparent as seen in figures 3e-f.At this higher dipole moment, the PDHS2 theory becomes less accurate as the polarizability and the density is increased.
From the figures and the data of table 1 it can be concluded that the dielectric constant is very sensitive to the polarizability.Both simulations and theory predict that is about 2-3 times larger for α * = 0.06 than for the non-polarizable fluid.This is a property of the model as was pointed out by other authors before.An insight into this behaviour can be gained by examining the contributions of the permanent dipoles, the induced dipoles, and the cross term between them to the mean square polarization ( m 2 0 , m 2 ind , and 2 m 0 m ind ) as functions of α * (see table 1).The increase of m 2 0 with α * is moderate, the increase of m 2 ind is steep but the absolute value remains small.On the contrary, the cross term increases steeply with α, and it seems to bear the main responsibility for the strong α-dependence of the dielectric constant.This behaviour implies a very strong correlation between the permanent and induced dipoles.For real matters, the α-dependence seems not to be so strong as was found for this model.Moreover, in most cases, α is not a scalar, but a tensor, the effect of which we intend to study in future works.Table 2 shows the data of an analysis on the dependence of the dielectric properties on the system size and the parameter R iter .We have chosen a PDHS point for this purpose with a relatively large dipole moment and polarizability at the highest reduced density.The results obtained for the various values of N and R iter fluctuate within an error bar that roughly corresponds to that obtained from the block average method.This implies that the N-and R iter -dependence of the dielectric constant is not too strong at this state point.N STEP /10000  A more serious problem for systems of strongly interacting dipoles is the slow convergence of the simulations.This is due to ergodicity problems: the molecules tend to fall into energy traps formed by low energy configurations (these are typically dipole chains or clusters of aligned dipoles).The dielectric constant is especially sensitive to the ergodicity-problem because it is calculated from a fluctuation formula that is usually a more poorly converging quantity than the usual ensemble averages.With single particle moves of the MC simulations, the sampling of the configuration space is not efficient; biased sampling or molecular dynamics simulations are preferable in these states.Figure 4 shows the dielectric constant obtained from the simulations tabulated in table 2 as a function of the number of MC cycles.It is seen that in most cases there is a visible drift in the curves.A rough extrapolation to large values of N STEP implies that the curves converge to the same value at least within an acceptable uncertainty.This suggests that the insufficient length of the simulations (arising from nonergodicity) is more crucial in these simulations than the value of N or R iter .
Apart from considering tensorial polarizability, another interesting continuation of this study is to calculate the dielectric constant of mixtures of polar molecules.This issue was considered in our earlier paper [60], where we extended our continuum approach to mixtures and showed that it correctly reproduces the mole fraction dependence of the relative permittivity as obtained from experiments [62].With a reasonable equation for the (x) function, we can determine the composition of polar mixtures from simple dielectric measurements.

Figure 1 .
Figure 1.The dielectric constant of PDHS fluids as a function of the reduced density at various values of the dipole moment and polarizability as obtained from Monte Carlo simulations and the different theories.The explanation of the curves and symbols can be found in the (c) inset of the figure.The various theories correspond to equations as listed: DHS1: equation (26), DHS2: equation (27), KF1: equation (31), KF2: equation (32), PDHS1: equation (38), and PDHS2: equation (39).

Figure 2 .
Figure 2. The dielectric constant of PDHS fluids as a function of the square of the reduced dipole moment at various values of density and polarizability as obtained from Monte Carlo simulations and the different theories.For explanation of the curves and symbols see the caption of figure 1.

2 =2Figure 3 .
Figure 3.The dielectric constant of PDHS fluids as a function of the reduced polarizability at various values of the dipole moment and polarizability as obtained from Monte Carlo simulations and the different theories.For explanation of the curves and symbols see caption of figure 1.

Figure 4 .
Figure 4.The dielectric constant as a function of the MC steps as obtaind from Monte Carlo simulations at various values of N and R iter for the point ρ * = 0.8, (µ * 0 ) 2 = 2, and α * = 0.06.