Finite-size effects at first-order isotropic-to-nematic transitions
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J un Finite-size effects at first-order isotropic-to-nematic transitions
J.M. Fish and R. L. C. Vink
Institute of Theoretical Physics, Georg-August-Universit¨at G¨ottingen,Friedrich-Hund-Platz 1, 37077 G¨ottingen, Germany (Dated: December 17, 2018)We present simulation data of first-order isotropic-to-nematic transitions in lattice models ofliquid crystals and locate the thermodynamic limit inverse transition temperature ǫ ∞ via finite-sizescaling. We observe that the inverse temperature of the specific heat maximum can be consistentlyextrapolated to ǫ ∞ assuming the usual α/L d dependence, with L the system size, d the latticedimension and proportionality constant α . We also investigate the quantity ǫ L,k , the finite-sizeinverse temperature where k is the ratio of weights of the isotropic to nematic phase. For anoptimal value k = k opt , ǫ L,k versus L converges to ǫ ∞ much faster than α/L d , providing an economicalternative to locate the transition. Moreover, we find that α ∼ ln k opt / L ∞ , with L ∞ the latentheat density. This suggests that liquid crystals at first-order IN transitions scale approximately as q -state Potts models with q ∼ k opt . PACS numbers: 05.70.Fh, 75.10.Hk, 64.60.-i, 64.70.mf
I. INTRODUCTION
The investigation of the isotropic-to-nematic (IN) tran-sition in liquid crystals via computer simulation is longestablished. Decades ago Lebwohl and Lasher (LL) in-troduced a simple lattice model, the LL model, to studythis transition [1]. At each site i of a cubic lattice theyattached a three-dimensional unit vector ~d i (spin) inter-acting with its nearest-neighbors via H = − ǫ X
0, with k B the Boltzmann constantand T the temperature. Despite its simplicity, the LLmodel captures certain aspects of liquid crystal phase be-haviour remarkably well, and has consequently receivedconsiderable attention [2, 3].A common problem in locating the IN transition viasimulation is the issue of finite system size. Phase tran-sitions are defined in the thermodynamic limit, whereassimulations always deal with finite particle numbers. Inorder to estimate the thermodynamic limit transitionpoint, it is typical to perform a number of simulationsfor different system sizes and to subsequently extrapo-late the results following some finite-size scaling (FSS)procedure. Which procedure to use depends on the typeof transition, i.e. whether it is continuous or first-order.In three spatial dimensions, the IN transition is typicallyfirst-order; in two dimensions both continuous [4, 5, 6]and first-order [7, 8] IN transitions can occur, dependingon the details of the interactions [9, 10, 11]. In this paperwe focus on the first-order case.The literature on FSS at first-order transitions is quiteextensive, for a review see [12], the majority of whichdeals exclusively with the Potts model [13]. An impor-tant result is that the “apparent” transition inverse tem-perature ǫ L,CV , obtained in a finite system of size L , is shifted from the thermodynamic limit value ǫ ∞ as [14, 15] ǫ L,CV = ǫ ∞ − α/L d + O (1 /L d ) , (2)with proportionality constant α >
0. Here ǫ L,CV is theinverse temperature where the specific heat in a finitesystem of size L attains its maximum, d is the spatial di-mension of the lattice and L denotes the linear extensionof the simulation box, generally square or cubic, withperiodic boundary conditions.We emphasize that Eq.(2) was derived for the Pottsmodel where the proportionality constant is known to be α Potts = ln q/ L ∞ , (3)with q the number of Potts states and L ∞ the latentheat density in the thermodynamic limit [14, 15]. Inter-estingly, simulations of the LL model have shown thatthe functional form of Eq.(2) also works well for IN tran-sitions [16, 17]. That is, meaningful extrapolations of ǫ L,CV can be performed, although the significance of α isnot obvious. It certainly cannot be related to the numberof spin states, i.e. conform to Eq.(3), since the LL modelis a continuous spin model, in contrast to the discrete spin variables of the Potts model [18].In any case, based on the success of Eq.(2) in describingfinite-size effects in the LL model, it could be hoped thatother scaling relations, originally derived for the Pottsmodel, also remain valid. Of particular interest is theresult of Borgs and Kotecky, who showed that for thePotts model exponentially decaying finite-size effects arealso possible [15, 19, 20]. The obvious advantage of expo-nential decay is that ǫ ∞ is approached much faster withincreasing L , compared to the power law decay of Eq.(2).This means that moderate system sizes may suffice tolocate the transition, thereby saving valuable computertime. As liquid crystal phase transitions are in any caseexpensive to simulate, such a gain in efficiency would cer-tainly be highly desirable.We will show in this paper that it is indeed possible tolocate first-order IN transitions from finite-size simula-tion data with shifts that vanish much faster than 1 /L d .This is possible by considering ǫ L,k , the inverse tempera-ture at which the “ratio-of-weights” of the isotropic andthe nematic phases is equal to a value k . This ratio-of-weights is obtained from the order parameter distri-bution P L,ǫ ( λ ), defined as the probability to observe anorder parameter λ , when simulating a system of size L at inverse temperature ǫ . In the vicinity of the IN tran-sition the distribution becomes bimodal, with one peakcorresponding to the isotropic phase and the other to thenematic phase. The ratio-of-weights is simply the ratio ofthe peak areas. Provided k is chosen optimally ǫ L,k ap-proaches ǫ ∞ extremely rapidly as L increases, yielding aneconomic alternative over Eq.(2). A prerequisite is thatthe transition must be strong enough first-order for theratio-of-weights to be meaningfully calculated [21]. Forthis reason, we do not consider the original LL model, asthe transition is extremely weak here, but a variation ofit.In this paper we firstly provide the details of the mod-ified LL model in Section II, together with a descriptionof the simulation method that was used to obtain theorder parameter distribution. Next, we measure ǫ ∞ us-ing the “standard” approach of extrapolating ǫ L,CV viaEq.(2), as well as using the “new” approach based on ǫ L,k . In particular we demonstrate how to locate the op-timal value k opt , along which finite-size effects are min-imal. As expected, both approaches are in good agree-ment, with the essential difference that ǫ L,k converges to ǫ ∞ already for very small systems. This fast convergenceproperty was observed at all transitions studied by us, ir-respective of space and spin dimension. We also considerthe finite-size scaling of the latent heat density and showthat, for IN transitions, k opt becomes the “analogue” ofthe number of Potts states q . Finally we present a sum-mary of our findings in Section IV. II. MODEL AND SIMULATION METHODA. Modified LL model
In order to study finite-size effects at phase transitions,simulation data of high statistical quality are essential.This sets a limit on the complexity of the models thatcan be handled as well as on the system size. For ourpurposes already the simple LL model is too demand-ing, the problem being that the IN transition in thismodel is extremely weak. Generally, in computer sim-ulations, first-order phase transitions are identified bymeasuring the probability distribution of the order pa-rameter [22]. At the transition, this distribution displaystwo peaks: one corresponding to the isotropic phase andthe other to the nematic phase. In the thermodynamiclimit, the peaks become sharper, and ultimately a distri-bution of two δ -functions is obtained. In finite systems,however, the peaks are broad and possibly overlapping,especially when the transition is very weak. Such be- haviour is observed in the LL model: even in simulationboxes of L = 70 lattice spacings the peaks strongly over-lap and the logarithm of the peak height, measured withrespect to the minimum in between, is less than 2 k B T [17]. Since the peaks overlap one never truly sees purephases, which complicates the analysis. In order to yieldreasonable results we require in this paper that the peaksin P L,ǫ ( λ ) be well-separated. More precisely, it must bepossible to assign a “cut-off” separating the peaks, onwhich the final results may not sensitively depend. Forthis reason we do not consider the original LL model butrather a generalization of it, where the exponent p ofEq.(1) exceeds the LL value. We expect this will lead toa much stronger first-order IN transition [9, 10, 11], sodistributions will display non-overlapping peaks alreadyin moderately sized systems. In fact by using a large ex-ponent p in Eq.(1) strong first-order IN transitions maybe realized even in purely two-dimensional systems [7, 8].Hence, the model that we consider is just the LL modelof Eq.(1) but with p >
2. Note the absolute value | · | such that the system is invariant under inversion of thespin orientation. We thus impose the symmetry of liquidcrystals although we believe that our results also applyto magnetic systems.Note that the use of a large exponent p in Eq.(1) mayalso yield a better description of experiments on con-fined liquid crystals. The latter systems are quasi two-dimensional. If one studies the LL model in two dimen-sions, i.e. with p = 2, and three-dimensional spins, a truephase transition appears to be absent [23]. In contrast,experiments clearly reveal that transitions do occur. Infact, these transitions appear to be of the IN type andare quite strong, as manifested by pronounced coexis-tence between isotropic and nematic domains [24]. Suchbehaviour cannot be reproduced easily with the standardLL model, but it can be using the modified version con-sidered in this work, with a sufficiently large exponent p . B. Transition matrix Wang-Landau sampling
Following earlier work on the LL model [16, 17] oursimulations are based on the order-parameter distribu-tion. We use the energy E of Eq.(1) as order parameterand aim to measure P L,ǫ ( E ) as accurately as possible.Recall that P L,ǫ ( E ) is the probability to observe energy E , in a system of size L , at inverse temperature ǫ . De-pending on the case of interest, the simulations are per-formed on square or cubic lattices of linear size L , usingperiodic boundary conditions.In order to obtain P L,ǫ ( E ) we use Wang-Landau (WL)sampling [25, 26] additionally optimized by recordingsome elements of the transition matrix (TM) [27, 28].The aim of WL sampling is to perform a random walkin energy space, such that all energies are visited equallyoften. To this end, we use single spin dynamics, wherebyone of the spins is chosen randomly and given a new ran-dom orientation. The new state is accepted with proba-bility p ( E I → E J ) = min (cid:20) g ( E I ) g ( E J ) , (cid:21) , (4)with E I and E J the energies of the initial and final statesrespectively and g ( E ) the density of states. The densityof states is unknown beforehand and g ( E ) is initially setso g ( E ) = 1. Upon visiting any particular energy thecorresponding density of states is multiplied by a modifi-cation factor f ≥
1. We also keep track of the histogram H ( E ), counting the number of times each energy E isvisited. Once H ( E ) contains sufficient information overthe range of energy of interest, the modification factor f is reduced and the energy histogram H ( E ) is reset tozero. These steps are repeated until f has become closeto unity, after which changes in the density of states be-come negligible. The sought order parameter distributionis then obtained from P L,ǫ ( E ) ∝ g ( E ) exp( − ǫE ).The above procedure is the standard WL algorithm,which works extremely well in many cases [1]. However,it has been noted [28, 29] that the WL algorithm in itsstandard form reaches a limiting accuracy, after whichthe statistical quality of the data no longer improves, nomatter how much additional computer time is invested.Hence, these authors also propose to measure the TMelements T ( E I → E J ). These are defined as the numberof times that, being in a state with energy E I , a statewith energy E J is proposed, irrespective of whether thenew state is accepted. From the TM elements one canestimate Ω( E I → E J ) = T ( E I → E J ) P K T ( E I → E K ) , (5)which is the probability that being in state with energy E I , a move to a state with energy E J is proposed. Thisis related to the density of states via [27] g ( E I ) g ( E J ) = Ω( E J → E I )Ω( E I → E J ) . (6)Hence, by recording TM elements the density of statescan also be constructed, the great advantage being thatrejected moves also give useful information.To combine WL sampling with the TM method wesomewhat follow [28]. At the start of the simulation thedensity of states g ( E ) is set to unity, while the energyhistogram H ( E ) and the TM elements are set to zero.We perform one WL iteration, i.e. accepting moves con-form Eq.(4), using a high modification factor ln f = 1.At each move both H ( E ) and the TM elements are up-dated. We continue to simulate until all bins in H ( E )contain at least n entries over the chosen energy range.We then use the TM elements to construct a new densityof states, which serves as the starting density of states forthe next WL iteration. For the next iteration H ( E ) is re-set to zero, the modification factor is reduced to (ln f ) /l but the TM elements remain untouched. These steps arerepeated until ln f ≈ O (10 − ), after which we store the corresponding density of states g P ( E ). This marks theend of the “prepare” stage.Next we proceed with the “collect” stage. The TM ele-ments are set to zero, whereas H ( E ) is no longer needed.During collection we sample according to Eq.(4) using g P ( E ) as estimate for the density of states. However,only the TM elements and not g P ( E ) are further up-dated. As collection proceeds the accuracy of the TM el-ements increases indefinitely, as does the accuracy of thedensity of states obtained from them. The reason to havea separate “collect” stage is because during “prepare” de-tailed balance is not strictly obeyed, due to the initiallylarge modification factor f [28]. For this reason, g P ( E )could be biased and we are reluctant to perform finite-size scaling with it.During the “prepare” stage small values n ≈
10 to-gether with large values l ≈ −
10 can be used. Thissignificantly speeds up the simulation and similar obser-vations have been made in other works [28]. Histogramswere collected by discretizing the energy in bins of reso-lution ∆ = 1 k B T . In order to avoid “boundary effects”during WL sampling states are counted as in [30]. Toreduce memory consumption, only the nearest-neighborelements T ( E I → E I ± ∆) of the TM along with thenormalization P K T ( E I → E K ) are stored. Since sin-gle spin dynamics are used, these are the dominant en-tries. Constructing the density of states using Eq.(6)and recursion is then a straightforward matter. If allTM elements were to be used, constructing the densityof states becomes more complex while not yielding sig-nificantly higher accuracy [27], so this is not attemptedhere. The required computer time depends sensitively onthe size of the system. For small systems consisting of ≈ ,
000 spins or more this can take more than oneweek. In these cases it is necessary to collect the densityof states over a number of separate energy intervals witha single processor assigned to each interval. Such a par-allelization is trivially implemented. The “collect” phasetypically lasts as long as the “prepare” phase except forvery large systems, where it is found that it takes a muchlonger time to obtain an equivalently accurate density ofstates.
III. RESULTS AND ANALYSIS
We have performed extensive simulations of Eq.(1)varying both the space and spin dimension, as well asthe exponent p . More precisely, the following scenariosare considered:1. three-dimensional lattices, three-dimensional spinswith p = 5 − p = 20 −
50 and l n P L , e ( r ) r = −E / L d Dr D F FIG. 1: Logarithm of P L,ǫ ( E ) using p = 10 in Eq.(1) for sys-tem sizes L = 10 , ,
15 (from top to bottom), cubic lattices,and three-dimensional spins. In each of the distributions ǫ was tuned so the peaks are of equal height. The barrier ∆ F ,here marked for the L = 10 system, is defined as the height ofthe peaks measured with respect to the minimum in between.The peak-to-peak distance ∆ ρ corresponds to the latent heatdensity. Note that we have plotted the distributions as afunction of the negative energy density: the left peak thuscorresponds to the isotropic phase and the right peak to thenematic phase.
3. two-dimensional lattices, two-dimensional spinswith p = 150 − p of Eq.(1) is sufficientlylarge [9, 10, 11]. Hence, in two dimensions, one generallyneeds p ≫ →
2, even greater exponents p are re-quired. For this reason, the chosen p ranges vary signifi-cantly between the three scenarios. A. Determining the order of the transition
In order to verify the presence of a first-order transitionwe use the scaling method of Lee and Kosterlitz [31, 32].Recall that the order parameter distribution becomes bi-modal in the vicinity of the IN transition; see Fig. 1 foran example. The idea of Lee and Kosterlitz is to monitorthe peak heights ∆ F of the logarithm of the order param-eter distribution, measured with respect to the minimum“in-between” the peaks. At a first-order transition ∆ F corresponds to the formation of interfaces between co-existing isotropic and nematic domains [33]. In d spatialdimensions it is therefore expected that ∆ F ∝ L d − , pro- D F L d−1 p=5p=8p=10p=15p=20p=45 FIG. 2: Variation of ∆ F versus L d − for d = 3 dimensionallattices, three-dimensional spins, and the various values of p as indicated. D F L d−1 p=20p=25p=30p=35p=40p=45p=50 FIG. 3: Similar to Fig. 2 but for two-dimensional lattices andthree-dimensional spins. D F L d−1 p=150p=1000 FIG. 4: Similar to Fig. 2 but for two-dimensional lattices andtwo-dimensional spins. e L , C V (b) p=20 d=2 3D−spins2.0252.0352.0452.0552.0650.000 0.002 0.004 0.006 0.008 0.0101/L d (c) p=150 d=2 2D−spins FIG. 5: Estimation of ǫ ∞ via extrapolation of ǫ L,CV usingEq.(2). Shown is ǫ L,CV versus 1 /L d using (a) p = 10, cu-bic lattices and three-dimensional spins, (b) p = 20, squarelattices and three-dimensional spins and (c) p = 150, squarelattices and two-dimensional spins. The open symbols aresimulation data and the lines are fits to Eq.(2). viding a straightforward recipe to identify the transitiontype: a linear increase of ∆ F versus L d − indicates that afirst-order transition is taking place, with the slope yield-ing the interfacial tension [33], whereas for a continuoustransition ∆ F becomes independent of L , or vanishesaltogether if no transition takes place at all in the ther-modynamic limit. In Fig. 2 ∆ F is plotted for the purelythree-dimensional case; the linear increase is clearly vis-ible, confirming the presence of a first-order transition.For two-dimensional lattices the results have been col-lected in Fig. 3 and Fig. 4. Once again the presence ofa first-order transition is confirmed. Note that on two-dimensional lattices the slopes of the lines correspond toline tensions. B. Extrapolation of ǫ L,CV
Next, we measure the thermodynamic limit inversetemperature ǫ ∞ by means of extrapolating ǫ L,CV viaEq.(2). Recall that ǫ L,CV is the finite-size inverse tem- perature where the specific heat C L = h E i − h E i L d (7)attains its maximum. Shown in Fig. 5(a) is ǫ L,CV versus1 /L d for the purely three-dimensional case using p =10 - results for different p are qualitatively similar andtherefore not explicitly shown. In agreement with earliersimulations of the “original” LL model [16, 17], the dataare well described by Eq.(2) and from the fit ǫ ∞ can bemeaningfully obtained. The resulting fit parameters arecollected in Table I. Repeating the same analysis for two-dimensional lattices yields similar results; some typicalplots are shown in Fig. 5(b) and (c) with the resulting fitparameters collected in Tables II and III. C. Extrapolation of ǫ L,k
We now arrive at the main result of this paper, namelythe estimation of ǫ ∞ by monitoring ǫ L,k . Recall that ǫ L,k is defined as the finite-size inverse temperature where theequality W N /W I = k (8)is obeyed, with W N and W I the areas under the nematicand isotropic peaks of the order parameter distributionrespectively. No matter what value of k is used, providedit is positive and finite, we expect that lim L →∞ ǫ L,k = ǫ ∞ . The reason is that in the thermodynamic limit abimodal order parameter distribution survives only at ǫ ∞ and not anywhere else [32]. Hence, keeping the arearatio fixed at some value of k whilst increasing L , ǫ L,k will definitely approach ǫ ∞ . The rate of the convergence,however, does depend on k . Assuming that the predictionof Borgs and Kotecky for the Potts model also holds atfirst-order IN transitions, it should be possible to locatean optimal value k opt at which the convergence to ǫ ∞ is fastest and hopefully faster than 1 /L d . Therefore, wepropose to manually inspect the convergence of ǫ L,k usingseveral values of k .A prerequisite for numerically solving Eq.(8) is thatthe transition must be sufficiently first-order in order forbimodal distributions P L,ǫ ( E ) with well-separated peaksto appear. By this we mean that the barrier ∆ F , definedin Fig. 1, is large enough. The areas of the nematic andisotropic peaks may then be calculated using W N = Z E c −∞ P L,ǫ ( E ) dE, W I = Z E c P L,ǫ ( E ) dE, (9)where we remind the reader that the energy in our modelis negative. The details of defining the “cut-off” energy E c are somewhat arbitrary, but as states around E c con-tribute exponentially little to the peak areas, the preciseform does not matter [34]. In this work E c is taken to bethe average E c = R EP L,ǫ ( E ) dE , with P L,ǫ ( E ) obtained TABLE I: Properties of the IN transition of Eq.(1) for three-dimensional lattices and three-dimensional spins versus p . Listedare the fit parameters ǫ ∞ ,CV and α of Eq.(2), the best estimate ǫ ∞ ,k obtained from the convergence of ǫ L,k along k opt , thelogarithm of k opt with uncertainty ∆ k , the latent heat density L ∞ and the ratio ln k opt / L ∞ . p ǫ ∞ ,CV α ǫ ∞ ,k ln k opt ± ∆ k L ∞ ln k opt / L ∞ . ± .
001 2 . ± . . ± .
010 5 . − .
28 1.5207 5.28 1 . ± .
001 5 . ± . . ± .
005 5 . − .
710 1.5862 4.28 1 . ± .
001 5 . ± . . ± .
005 4 . − .
015 1.7126 3.94 1 . ± .
001 6 . ± . . ± .
002 3 . − .
120 1.8063 3.76 1 . ± .
001 6 . ± . . ± .
002 3 . − .
945 2.0838 3.55 2 . ± .
001 8 . ± . . ± .
001 3 . − . p ǫ ∞ ,CV α ǫ ∞ ,k ln k opt ± ∆ k L ∞ ln k opt / L ∞
20 2.7695 5.18 2 . ± .
001 4 . ± . . ± .
001 5 . − .
325 2.8678 4.84 2 . ± .
001 4 . ± . . ± .
001 4 . − .
330 2.9517 4.69 2 . ± .
001 4 . ± . . ± .
001 4 . − .
035 3.0240 4.58 3 . ± .
001 5 . ± . . ± .
001 4 . − .
440 3.0882 4.52 3 . ± .
001 5 . ± . . ± .
001 4 . − .
745 3.1456 4.47 3 . ± .
001 5 . ± . . ± .
001 3 . − .
850 3.1976 4.50 3 . ± .
001 5 . ± . . ± .
001 4 . − . at equal-height, i.e. as in Fig. 1. Once E c has been setits value is kept fixed whilst solving Eq.(8).For the purely three-dimensional case, the behaviourof ǫ L,k is shown in Fig. 6. Using a number of exponents p in Eq.(1), we have plotted ǫ L,k versus L for several val-ues of k . The data are consistent with the expectationthat, regardless of k , ǫ L,k converges to a common value,corresponding to ǫ ∞ . Note also that ǫ ∞ is approachedfrom above for large k , and from below for low k . Hence,we can indeed identify an optimal value k opt along whichfinite-size effects are minimal. The optimum can be esti-mated by locating, for a pair of system sizes L i and L j ,the inverse temperature ǫ ij where for both system sizesthe same ratio k ij of the peak areas is observed. By con-sidering all available pairs of system sizes, the averageand root-mean-square fluctuation in ǫ ij and k ij can becalculated, which then yield ǫ ∞ and k opt with uncertain-ties, shown in Table I. Although k opt itself is not knownvery precisely, since ∆ k is quite large, very accurate esti-mates of ǫ ∞ can still be obtained as this quantity is ratherinsensitive to the precise value of k opt being used. Thismeans that the series ǫ L,k also provides a valid methodfor locating IN transitions. The corresponding estimatesof ǫ ∞ are in good agreement with those obtained via ex-trapolation of ǫ L,CV , as inspection of the various tablesindicates. The practical advantage of using ǫ L,k with k = k opt is that the L -dependence is very weak, so muchso that ǫ ∞ is captured already in small systems. Similarfindings were obtained using two-dimensional lattices, ofwhich some typical plots are provided in Fig. 7 and Fig. 8with the corresponding numerical estimates collected in Tables II and III.For non-optimal values k = k opt , we observe that theshift ǫ ∞ − ǫ L,k ∝ /L d , i.e. the shift vanishes as a powerlaw in the inverse volume, similar to ǫ L,CV . At the opti-mal value k = k opt , finite-size effects in ǫ L,k are typicallytoo small in order for a meaningful fit to be carried out.Hence, our data confirm Borgs and Kotecky in the sensethat optimal estimators can be defined which convergeonto ǫ ∞ faster than 1 /L d ; whether the optimal conver-gence is indeed exponential requires more accurate data,which is currently beyond our reach [35].An alternative, but completely equivalent, methodto investigate the convergence of ǫ L,k is presented in[34, 36, 37, 38], albeit for the Potts model. The ideais to plot the area ratio W N /W I versus ǫ for several sys-tem sizes. The resulting curves are expected to revealan intersection point at the transition inverse tempera-ture; the value of the area ratio at the intersection thenyields k opt . For completeness we have prepared one suchplot, see Fig. 9. The curves indeed intersect and giveestimates of ǫ ∞ and k opt that are fully consistent withthose reported in Table I. D. Latent heat density
It appears that scaling relations derived for the Pottsmodel also work remarkably well at IN transitions. Inagreement with earlier simulations of the LL model [16,17], the validity of Eq.(2) is confirmed additionally byus. Furthermore, our data suggest that an analogue of
TABLE III: Similar to Table I but for two-dimensional lattices and two-dimensional spins. p ǫ ∞ ,CV α ǫ ∞ ,k ln k opt ± ∆ k L ∞ ln k opt / L ∞
150 2.063 3.67 2 . ± . . ± . . ± .
005 1 . − . . ± . . ± . .
270 + 0 .
005 2 . − . e L , k (c) p=20(b) p=10(a) p=51.58401.58501.58601.58701.58801.5890 e L , k (c) p=20(b) p=10(a) p=51.38501.39001.39501.40001.40501.4100 e L , k (c) p=20(b) p=10(a) p=5 FIG. 6: Variation of ǫ L,k versus L , for three-dimensional lat-tices and three-dimensional spins, using different exponents p as indicated. The symbols are simulation data, the linesserve to guide the eye. The central curves in each plot show ǫ L,k using k = k opt along which finite-size effects are minimal;also shown is ǫ L,k using k = k opt − k (lower curves), and k = k opt +5∆ k (upper curves). The methods for locating k opt and ∆ k are explained in the text and the resulting values, aswell as the estimates of ǫ ∞ , are listed in Table I. the Borgs and Kotecky prediction, namely that finite-sizeeffects vanish faster than 1 /L d at appropriate points, canbe defined. In this case, it is needed to measure ǫ L,k using the optimal value k = k opt . In the Potts modelit holds that k opt = q , where q is the number of Pottsstates. In other words, finite-size effects in the Pottsmodel are minimized when the ratio of the peak areas inthe order parameter distribution is held fixed at q . Basedon our results, it seems reasonable to assume that scalingrelations for the Potts model also hold at IN transitions,but with q replaced by k opt .To test this assumption we consider the proportional-ity constant α from Eq.(2), which is given by Eq.(3) forthe Potts model. If q can be replaced by k opt , α shouldcorrespond to ln k opt / L ∞ , where L ∞ is the latent heatdensity. The latter can be obtained independently from C L, max = L ∞ L d /
4, where C L, max is the maximum value e L , k (c) p=50(b) p=35(a) p=202.983.003.023.043.06 e L , k (c) p=50(b) p=35(a) p=202.742.762.782.80 e L , k (c) p=50(b) p=35(a) p=20 FIG. 7: Similar to Fig. 6 but for two-dimensional latticesand three-dimensional spins. Numerical estimates are givenin Table II. e L , k (b) p=10001.98 e L , k (a) p=150 FIG. 8: Similar to Fig. 6, but for two-dimensional latticesand two-dimensional spins. Numerical estimates are given inTable III. −30−20−10 0 10 20 30 401.5850 1.5855 1.5860 1.5865 1.5870 l n W N / W I e L=10L=12L=15L=20L=25L=30
FIG. 9: Plot of ln W N /W I versus ǫ for several system sizes L . This plot was generated using p = 10 in Eq.(1), three-dimensional lattices and three-dimensional spins. Dr L,1 Dr L,2 l a t en t hea t den s i t y (b) p=20 d=2 3D−spins Dr L,1 Dr L,2 Dr L,1 Dr L,2
FIG. 10: Finite-size variation of the latent heat density es-timators ∆ ρ L, and ∆ ρ L, . Results are shown for Eq.(1)using (a) p = 10, three-dimensional lattices and three-dimensional spins, (b) p = 20, two-dimensional lattices andthree-dimensional spins and (c) p = 150, two-dimensional lat-tices and two-dimensional spins. of the specific heat in a finite system of size L [15]. Hence,we introduce the latent heat estimator∆ ρ L, = q C L, max /L d , (10)which should approach L ∞ as L → ∞ . Addition-ally, the latent heat density can be read-off directly,as the peak-to-peak distance in the energy distribution,marked ∆ ρ in Fig. 1. Numerically this is expressed by M L = 2 h| E − h E i|i /L d ; plotting M L versus ǫ gives amaximum ∆ ρ L, , which in the limit L → ∞ also ap-proaches L ∞ . Typical behaviour of ∆ ρ L,i is shown inFig. 10. As expected, both latent heat estimators con-verge to a common value, which can be read-off reason-ably accurately; the resulting estimates of L ∞ are givenin the various tables. Note also that L ∞ is approachedfrom below in three dimensions, whereas in two dimen-sions, it is approached from above. If an appropriatenumber of two-dimensional lattice layers stacked on topof each other were simulated, it is likely that a cross-overregime could be found where ∆ ρ L,i depends only weaklyon L , as these systems are effectively in-between two andthree dimensions.Having measured L ∞ , the ratio ln k opt / L ∞ is eas-ily obtained, which may then be compared to α ,see Tables I-III. The uncertainty is admittedly ratherlarge, but within numerical precision, and the relationln k opt / L ∞ ∼ α appears to hold. IV. SUMMARY
In this paper we have presented simulation data offirst-order isotropic-to-nematic transitions in lattice liq-uid crystals with continuous orientational degrees of free-dom for various space and spin dimensions. As with ear-lier simulations of this type [16, 17], we find that the ex-trapolation of the finite-size inverse temperature of thespecific heat maximum ǫ L,CV can be consistently per-formed assuming a leading α/L d dependence, exactly asin the Potts model. Inspired by this result, we have inves-tigated an alternative approach to locate the transitioninverse temperature using estimators ǫ L,k , defined as thefinite-size inverse temperature where the ratio of peak ar-eas in the energy distribution is equal to k . In agreementwith the Potts model, ǫ L,k converges to ǫ ∞ much fasterthan 1 /L d , provided an optimal value k = k opt is used.Moreover, the ratio k opt / L ∞ , with L ∞ the latent heatdensity, is remarkably consistent with the proportional-ity constant α from the scaling of ǫ L,CV . This leads usto conclude that finite-size scaling predictions originallyproposed for first-order transitions in the Potts model re-main valid at first-order IN transitions too, but with thenumber of Potts states q replaced by k opt .It is perhaps somewhat surprising that a continuousspin model at a first-order transition, such as the LLmodel, scales in the same way as the Potts model, whichis, after all, a discrete spin model. In fact, Borgs andKotecky have remarked that the derivation of their scal-ing results cannot be easily extended to continuous spinmodels [20]. Nevertheless, the LL model and its vari-ants may be more closely connected to the Potts modelthan one may initially think. Note that for large p theHamiltonian of Eq.(1) becomes increasingly Potts-like,in the sense that the pair interaction approaches a δ -function: lim p →∞ | ~d i · ~d j | p = δ ( ~d i , ~d j ). This implies thatneighbouring spins only interact when they are closelyaligned and are otherwise indifferent to each other, justas in the Potts model. It has indeed been suggested thatsuch models approximately resemble q -state Potts mod-els, with q ∼ √ p [9]. The observed trends in this work arecertainly consistent with this interpretation. For all casesconsidered the strength of the transition increases with p ,as manifested by the growing latent heats and interfacial tensions, exactly as in the Potts model with increasing q .Also the upward shift of ǫ ∞ with p is consistent with thePotts model. However, it is clear that new theoreticalapproaches are needed to fully understand finite-size ef-fects at first-order transitions in the models studied here.We hope that the present simulation results may inspiresuch efforts. Acknowledgments
This work was supported by the
Deutsche Forschungs-gemeinschaft under the Emmy Noether program(VI 483/1-1). We thank Marcus M¨uller for stimulatingdiscussions. [1] P. A. Lebwohl and G. Lasher, Physical Review A , 426+(1972).[2] U. Fabbri and C. Zannoni, Molecular Physics , 763(1986).[3] C. M. Care and D. J. Cleaver, Reports on Progress inPhysics , 2665 (2005).[4] M. A. Bates and D. Frenkel, The Journal of ChemicalPhysics , 10034 (2000).[5] D. Frenkel and R. Eppenga, Phys. Rev. A , 1776(1985).[6] H. Kunz and G. Zumbach, Physical Review B , 662+(1992).[7] R. L. C. Vink, Physical Review Letters , 217801+(2007).[8] H. H. Wensink and R. L. C. Vink, Journal of Physics:Condensed Matter , 466109+ (2007).[9] E. Domany, M. Schick, and R. H. Swendsen, Phys. Rev.Lett. , 1535 (1984).[10] A. C. D. van Enter and S. B. Shlosman, Phys. Rev. Lett. , 285702 (2002).[11] A. C. D. van Enter, S. Romano, and V. A. Zagrebnov, J.Phys. A (2006).[12] K. Binder, Reports on Progress in Physics , 487 (1997).[13] The exceptions appear to be: M. Fisher and V. Privman,J. Appl. Phys. , 3327 (1985); Phys. Rev. B , 447(1985); Commun. Math. Phys. , 527 (1986).[14] M. S. Challa, D. P. Landau, and K. Binder, PhysicalReview B , 1841+ (1986).[15] C. Borgs, R. Koteck´y, and S. Miracle-Sol´e, Journal ofStatistical Physics , 529 (1991).[16] Z. Zhang, O. G. Mouritsen, and M. J. Zuckermann, Phys.Rev. Lett. , 2803 (1992).[17] N. V. Priezjev and R. A. Pelcovits, Physical Review E , 062702+ (2001).[18] F. Y. Wu, Rev. Mod. Phys. , 235 (1982).[19] C. Borgs and R. Koteck´y, Journal of Statistical Physics , 79 (1990).[20] C. Borgs and R. Kotecky, Phys. Rev. Lett. , 1734(1992).[21] Strictly speaking, Eq.(2) also requires a strong enough first-order transition. For example, in the derivation ofRef. 14, it is assumed that the order parameter distribu-tion consists of two non-overlapping Gaussians.[22] K. Vollmayr, J. D. Reger, M. Scheucher, and K. Binder,Zeitschrift f¨ur Physik B Condensed Matter , 113(1993).[23] C. Chiccoli, P. Pasini, and C. Zannoni, Physica A: Sta-tistical and Theoretical Physics , 298 (1988).[24] R. Garcia, E. Subashi, and M. Fukuto, Physical ReviewLetters , 197801+ (2008).[25] F. Wang and D. P. Landau, Physical Review Letters ,2050+ (2001).[26] F. Wang and D. P. Landau, Physical Review E (Statisti-cal, Nonlinear, and Soft Matter Physics) (2001).[27] J.-S. Wang and R. H. Swendsen, Journal of StatisticalPhysics , 245 (2002).[28] S. M. Shell, P. G. Debenedetti, and A. Z. Panagiotopou-los, The Journal of Chemical Physics , 9406 (2003).[29] Q. Yan and J. J. de Pablo, Physical Review Letters ,035701+ (2003).[30] B. J. Schulz, K. Binder, M. M¨uller, and D. P. Landau,Physical Review E , 067102+ (2003).[31] J. Lee and J. M. Kosterlitz, Phys. Rev. Lett. , 137(1990).[32] J. Lee and J. M. Kosterlitz, Physical Review B , 3265+(1991).[33] K. Binder, Phys. Rev. A , 1699 (1982).[34] W. Janke, Physical Review B , 14757+ (1993).[35] For the Potts model, the exponential size dependence hasbeen accurately resolved [38].[36] C. Borgs and W. Janke, Physical Review Letters ,1738+ (1992).[37] W. Janke, in Computer Simulations in Condensed Mat-ter Physics VII , edited by D. P. Landau, K. K. Mon,and H. B. Schuettler, pages 29+ (Springer-Verlag, Berlin,Heidelberg, 1994).[38] A. Billoire, R. Lacaze, and A. Morel, Nuclear Physics B370