{'fraction_non_alphanumeric': 0.034794794794794796, 'fraction_numerical': 0.007247247247247247, 'mean_word_length': 4.5305580159433125, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical Mechanics of Splay Flexoelectricity in Nematic Liquid Crystals

We develop a lattice model for the splay flexoelectric effect in nematic liquid crystals. In this model, each lattice site has a spin representing the local molecular orientation, and the interaction between neighboring spins represents pear-shaped molecules with shape polarity. We perform Monte Carlo simulations and mean-field calculations to find the behavior as a function of interaction parameters, temperature, and applied electric field. The resulting phase diagram has three phases: isotropic, nematic, and polar. In the nematic phase, there is a large splay flexoelectric effect, which diverges as the system approaches the transition to the polar phase. These results show that flexoelectricity is a statistical phenomenon associated with the onset of polar order.


I. INTRODUCTION
Flexoelectricity is a coupling between elastic deformation and electrostatic polarization in a liquid crystalline medium. In general, a splay or bend deformation of the nematic director leads to an electrostatic polarization, which can be observed as a macroscopic dipole moment of the liquid crystal. Conversely, an applied electric field induces an electrostatic polarization, which leads to a combination of splay and bend distortions in the nematic director. Since its discovery in 1969 by Meyer [1], the flexoelectric effect has drawn great interest because of its possible applications [2,3] in strain gauges, transducers, actuators, micro power generator and electro-optical devices.
There have been many experimental and theoretical studies to determine the flexoelectric coefficients of nematic liquid crystals [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18], using a range of different approaches. For typical calamitic (rod-shaped) liquid crystals, the splay and bend flexoelectric coefficients are in the range of 3-20 pC/m. However, in recent experiments, Harden et al. [2,3] found that bent-core liquid crystals have a surprisingly large bend flexoelectric coefficient, up to 35 nC/m, roughly three orders of magnitude larger than the typical value. With this large bend flexoelectric coefficient, bent-core liquid crystals may be practical materials for converting mechanical into electrical energy.
For theoretical physics, a key question is how to explain the large flexoelectric effect found in bent-core nematic liquid crystals, so that it can be exploited for technological applications. Our conjecture is that the large flexoelectric effect is a statistical phenomenon associated with nearby polar phase. Near a polar phase, a nematic liquid crystal is on the verge of developing spontaneous polar order, and hence any deformation of the director should induce a large polar response. To test this conjecture, we would like to build a model with nematic and polar phases, and determine the behavior of flexoelectric effect as a function of temperature above the nematic- * Electronic address: jselinge@kent.edu polar transition. In this paper, we begin the study by investigating the splay flexoelectric effect in a system of uniaxial pear-shaped molecules. In a subsequent paper, we will investigate the more complex case of bend flexoelectricity in bent-core liquid crystals, as in the experiments.
To study the splay flexoelectric effect, we generalize the Lebwohl-Lasher lattice model of nematic liquid crystals [19]. In the original Lebwohl-Lasher model, each lattice site i has a spin n i , which represents the local nematic director, with the symmetry n i → −n i . In our generalization, the spins represent the orientations of pear-shaped molecules, which do not have that symmetry. For that reason, the interaction between neighboring spins includes three terms-one term favoring nematic order, another term favoring polar order, and a third term that couples polar order with splay of the nematic director. With this interaction, we find a phase diagram with isotropic, nematic, and polar phases, as illustrated in the snapshots of Fig. 1. The nematic phase has a flexoelectric effect, which increases as the system approaches the polar phase. Thus, this calculation demonstrates explicitly that the flexoelectric effect is a collective, statistical phenomenon, which is strongest near the transition to a phase with spontaneous polar order.
The plan of this paper is as follows. In Sec. II, we set up the theoretical framework leading to the model interaction and discuss the relevant order parameters. In Sec. III, we present the Monte Carlo simulation methods and results, for both the phase diagram and the flexoelectric effect. In Sec. IV we present a mean-field theory for this model, and compare the simulation results with the mean-field approximation. Finally, in Sec. V we discuss and summarize the conclusions of this study.

II. MODEL
In this study, our goal is to simulate the splay flexoelectric effect in a system of uniaxial pear-shaped molecules. For these simulations, we construct a lattice model that can represent both nematic and polar order. In this model, the local molecular orientation at lattice site i is represented by a unit vectorn i . If the system has ne- matic order, then the molecular orientations tend to be aligned along a preferred axis; i.e. there is a nonzero order parameter P 2 (n i ·d) , whered is the overall director and P 2 is the second Legendre polynomial. If the system has polar order, then the molecular orientations tend to point in a particular direction; i.e. there is a nonzero order parameter P 1 (n i ·d) , where P 1 is the first Legendre polynomial. Note that the system can have nematic order without polar order, but it cannot have polar order without nematic order.
The lattice Hamiltonian must have four terms: one term that favors nematic order, one term that favors polar order, one term that gives a coupling between polar order and an applied electric field, and a final term that gives a coupling between polar order and splay of the nematic director. The term favoring nematic order can be written simply as −A(n i ·n j ) 2 , summed over all pairs of neighboring sites i and j, as in the Lebwohl-Lasher model [19]. The term favoring polar order can be written even more simply as −B(n i ·n j ), again summed over all pairs of neighboring sites i and j, as in the Heisenberg model of magnetism. The coupling between polar order and an applied electric field can be written as −E ·n i , summed over i.
The coupling between polar order and nematic splay is somewhat more subtle. For this coupling we need a lattice expression for the local splay between neighboring sites i and j. Our expression for the local splay should depend only on the nematic director, and hence it should be invariant under the transformationn → −n. We cannot describe splay by the scalar ∇ ·n, because it is not invariant under that transformation. Rather, we must describe splay by the vectorn(∇ ·n), which has the correct symmetry. In the following calculation, we let Latin letters refer to lattice sites and Greek letters refer to directions.
On a continuum basis, the splay vectorn(∇·n) can be written in terms of the local nematic order tensor Q αβ (r), or equivalently in terms of the dyad n α (r)n β (r), as Hence, a lattice approximation to the splay vector between sites i and j can be written as + n iα n iγ + n jα n jγ 2 r ijβ (n jβ n jγ − n iβ n iγ ) − n iβ n iγ + n jβ n jγ 2 r ijβ (n jα n jγ − n iα n iγ ) , wherer ij = r j − r i is the unit vector from site i to j on the lattice. After some algebra, this expression simplifies to [n(∇ ·n)] ij = 1 2 n j (r ij ·n j ) −n i (r ij ·n i ) Note that this expression is invariant under the transformationsn i → −n i ,n j → −n j , and i ↔ j. Now that we have found an expression for the local splay vector, we can couple it with the local polar order. The coupling term in the lattice Hamiltonian can be written as the dot product of the splay between sites i and j with the average polar order on these sites, Simplifying with the use of Eq. (3), the coupling term between splay and polar order in the Hamiltonian is Combining all these terms, our final expression for the lattice Hamiltonian is At this point, we want to use this lattice Hamiltonian to calculate the nematic order parameter P 2 , the polar order parameter P 1 , and the average splay vector n(∇ ·n) as functions of the parameters A, B, and C and the electric field E. In the following sections, we will do this calculation through Monte Carlo simulations and mean-field theory.

III. MONTE CARLO SIMULATION
As a first step in exploring this model, we carry out Monte Carlo simulations of a system of pear-like molecules interacting with the lattice Hamiltonian of Eq. (6). In these simulations, we use a simple cubic lattice of size 20 × 20 × 20. When an electric field is applied, it is in the z direction, so that the molecules tend to align along z, with splay in the x and y directions. The lattice has periodic boundary conditions in z, but free boundaries in x and y, so that it can form splay in those directions.
The usual Metropolis algorithm was used for lattice updates. In each Monte Carlo step, a lattice site is chosen randomly, its orientation is changed slightly, and the change in energy ∆E is calculated. If ∆E < 0 the move is accepted, and if ∆E > 0 the move is accepted with probability exp(−∆E/k B T ). Starting from the hightemperature isotropic phase, the system is cooled down slowly with temperature steps of ∆T = 0.02. The final configuration at each temperature is taken as the initial configuration for the next lower temperature. Typical runs take about 10 5 steps to come to equilibrium, while runs near phase transitions take about 6×10 5 steps. The nematic and polar order parameters and the splay vector are calculated and time-averaged during the production cycle.
The nematic order parameter P 2 is calculated by the usual method using the 3D nematic order tensor where α and β = x, y, z, and N is the total number of lattice sites. The largest eigenvalue of this order tensor corresponds to P 2 .
To calculate the polar order parameter P 1 , we assume that polar order is oriented along the same axis as nematic order, as expected for uniaxial molecules. The eigenvector corresponding to the largest eigenvalue of the nematic order tensor Q αβ is the instantaneous director d. Hence, the polar order parameter is calculated as the average dot product of the director with the molecular orientation, The splay vector is calculated from Eq. (3), averaged over the four bonds in the (x, y) plane. The magnitude of this vector gives the average angle between the molecular orientations on neighboring lattice sites. For that reason, we report this magnitude as ∆θ . Figure 2 shows plots of the order parameters P 2 , P 1 , and ∆θ as functions of temperature for several values of the interaction parameters. In Fig. 2(a), for a small polar coupling B and no applied electric field, we see an isotropic-nematic transition at high temperature followed by a nematic-polar transition at low temperature. At the isotropic-nematic transition, the nematic order parameter goes from zero to a nonzero value. Here the transition is rounded by finite-size effects; we would expect a sharp first-order transition for an infinite system. Throughout the nematic temperature range, the polar order parameter and splay are both zero. At the nematic-polar transition, the polar order parameter becomes nonzero, and this polar order induces an accompanying splay. The nematic order parameter decreases as the system moves into the polar phase, because the splayed molecular orientation partially averages out the alignment, as shown in the snapshot of Fig. 1(c). In Fig. 2(b), for a larger polar coupling B, we see a direct transition from the isotropic to the polar phase, with no intervening nematic phase. In this case, the nematic and polar order parameters both become nonzero at the same transition temperature. Once again, the polar order induces a splay, which inhibits the growth of the nematic order parameter.
The simulation results for the phase diagram are shown in Fig. 2(c). In this phase diagram, the vertical axis shows temperature while the horizontal axis shows the polar coupling B for a constant nematic coupling A. For small B the phase diagram shows isotropic, nematic, and polar phases, with a nematic range that decreases as B increases. At a sufficiently large value of B, the nematic phase disappears and there is a direct transition from isotropic to polar. Note that this phase diagram is quite similar to the phase diagram found in recent work on the 2D isotropic, tetratic, and nematic phases [21].
The question is now: What happens to the nematic phase when an electric field is applied? To answer this question, Fig. 2(d) shows the simulation results for the same parameters as Fig. 2(a), but in the presence of a small electric field. In the high-temperature isotropic phase, the field induces some polar and nematic order, but this effect is very small. However, in the nematic phase, the field induces a more substantial polar order, and that polar order induces a splay in the nematic di- rector, i.e. a converse flexoelectric effect. Both the polar order and the splay are quite temperature-dependent, increasing as the system approaches the nematic-polar transition temperature, as would be expected for a divergent susceptibility above a second-order transition. The nematic-polar transition is rounded off by the applied field, and the polar order parameter and splay saturate in the low-temperature polar phase.
To provide further insight into the effect of an applied electric field, Fig. 3 shows the splay as a function of temperature for several values of the field. Within the nematic phase, the splay increases as the electric field increases, as expected for the converse flexoelectric effect. This trend is reasonable because an increasing electric field enhances polar order. For small field the splay is quite sensitive to temperature, but for large field it becomes less temperature-dependent, as the induced polar order grows larger and approaches saturation. In the low-temperature polar phase, the splay shows the opposite trend with electric field; it now decreases as the field increases. This trend is reasonable because an in-creasing electric field cannot enhance the polar order, which is already saturated; it only aligns the direction of polar order. This alignment reduces the induced splay, since splay necessarily involves some misalignment of the molecular orientation.
For comparison, Fig. 4 presents a plot of splay as a function of temperature for several values of the interaction coefficient C. This graph shows that the splay increases as C increases, over the full temperature range, in all phases. This result is reasonable because the coefficient C represents the flexoelectric coupling between polar order and splay.
As a final point, note that the behavior presented here can only occur in the limit of small splay ∆θ < ∼ π/N , where N is the system size. In the opposite limit ∆θ > ∼ π/N , the system is too large for one single splay from side to side. Instead, it must break up into modulated structures consisting of regions of splay separated by domain walls. These modulated structures might be splay stripes or even more complex two-or threedimensional arrangements of splay cells [22]. We have observed such modulated structures in our simulations, but we have not explored them in detail because they are not likely to occur in experiments, where the magnitude of splay is generally small.

IV. MEAN-FIELD CALCULATION
In this section we discuss two approximate analytic approaches to solve the problem. First, we map the interaction onto an Ising model, and use this Ising model to calculate the splay and polar order as functions of temperature and electric field. Second, we present a more general mean-field calculation with full rotational degrees of freedom, and use it to calculate the full phase diagram with isotropic, nematic, and polar phases.

A. Ising model
For a simple Ising-type model of the splay flexoelectric effect, we suppose that the system has well-defined nematic order, with variable amounts of splay and polar order. Consider a particular site i surrounded by six nearest neighbors on a cubic lattice. We suppose that site i has its director along the z-axis, as do the two neighbors above and below, while the four nearest neighbors in the xy-plane have directors that are splayed outward by a small angle ∆θ. The polar order at any site j is represented by an Ising spin variable σ j = ±1, which indicates whether the molecular orientation is pointing up or down along the local director. Thus, the central site i has the molecular orientationn i = σ i (0, 0, 1), while the six neighbors have the orientationsn +x = σ +x (sin ∆θ, 0, cos ∆θ),n −x = σ −x (− sin ∆θ, 0, cos ∆θ),n +y = σ +y (0, sin ∆θ, cos ∆θ), n −y = σ −y (0, − sin ∆θ, cos ∆θ),n +z = σ +z (0, 0, 1), and n −z = σ −z (0, 0, 1).
We now substitute these expressions for the molecular orientations into the lattice Hamiltonian of Eq. (6). As usual in mean-field theory, we assume that all the neighbors of site i have polar order given by σ j = M . (This quantity is called P 1 in the other sections; here we use the symbol M to emphasize the analogy with the Ising magnetization.) The mean potential experienced by site i, expanded to second order in the small splay ∆θ, is then Hence, the effective field acting on the Ising spin σ i is As a result, the polar order parameter must satisfy the self-consistency equation = tanh Furthermore, minimization of the mean potential over the splay ∆θ gives Solving Eqs. (11) and (12) simultaneously gives the equilibrium values of the splay ∆θ and polar order M , as functions of electric field E, temperature T , and energetic parameters A, B, and C.
To calculate the response to an electric field in the nematic phase, we assume that E, M , and ∆θ are all small, and expand Eqs. (11) and (12) to linear order in this quantities. These equations imply Note that Eq. (13) gives the polar order parameter induced by an applied electric field, while Eq. (14) gives the converse flexoelectric effect induced by the field. Both of these responses increase as the temperature decreases toward the second-order nematic-polar transition at the temperature At this transition, they diverge as (T − T N P ) −γ , with critical exponent γ = 1, as expected for the susceptibility to an applied field, in mean-field theory for the Ising model.
For a more precise calculation, we solve Eqs. (11) and (12) numerically as functions of temperature and field. The numerical results for splay ∆θ and polar order M are shown in Fig. 5. As in the approximate analytic calculation above, we see a second-order nematic-polar transition. The high-temperature nematic phase has no polar order or splay without a field, but an applied field induces both of these quantities. By contrast, the lowtemperature polar phase has both spontaneous polar order and spontaneous splay, and they both increase moderately when a field is applied.
Although the Ising model is successful in explaining some features of our Monte Carlo simulations, it is incomplete because it assumes perfect nematic order-the molecules can have only two possible orientations, up and down. It cannot describe the behavior of the nematic order parameter as a function of temperature. For that reason, we proceed to a more general mean-field theory, in which each molecule has full rotational degrees of freedom.

B. General Mean-Field Calculation
In mean-field theory, the free energy can be written as averaged over the single-particle distribution function ρ. Thus, our goal is to express the single-particle distribution function in terms of some variational parameters, calculate the energetic and entropic terms in the free energy, and then minimize the free energy over those variational parameters. For this mean-field calculation, we write the molecular orientation at each lattice site in terms of the polar angle θ i and azimuthal angle φ i with respect to the local director. We assume the distribution function depends only on the polar angle θ i , and hence write where v 1 and v 2 are variational parameters. The order parameters are then P 1 = π 0 P 1 (cos θ)ρ(θ)dθ and P 2 = π 0 P 2 (cos θ)ρ(θ)dθ, and the partition function is Z = π 0 exp[v 1 P 1 (cos θ) + v 2 P 2 (cos θ)]dθ. The entropic contribution to the free energy is therefore (18) per lattice site.
As in the previous section, we suppose that site i has its director along the z-axis, as do the two neighbors above and below, while the four nearest neighbors in the xy-plane have directors that are splayed outward by a small angle ∆θ. To calculate the average energy, we combine our distribution function of Eq. (17) with the lattice Hamiltonian of Eq. (6). After averaging over all the angles, neglecting terms involving the third-order Legendre polynomials, and expanding to second order in the small splay ∆θ, we obtain per lattice site. We now have an expression for the total free energy per lattice site. In this expression, note that P 1 , P 2 , and Z are all determined by the parameters v 1 and v 2 in the distribution function. Hence, the free energy is a function of just three variational parameters: v 1 , v 2 , and ∆θ. Thus, in the mean-field calculation, we must minimize the free energy numerically with respect to those three parameters. After this minimization, we can calculate the order parameters P 1 and P 2 , and hence determine whether the system is in an isotropic, nematic, or polar phase. Figure 6 shows the numerical mean-field results for the order parameters P 1 and P 2 and splay ∆θ, as well as a complete phase diagram as a function of temperature T . For a small polar interaction B = 0.09 there are two transitions, first from the high-temperature isotropic phase ( P 1 = 0, P 2 = 0) to the intermediate nematic phase (( P 1 = 0, P 2 = 0), and then from the nematic phase to the low-temperature polar phase ( P 1 = 0, P 2 = 0). The isotropic-nematic transition is first-order, while the nematic-polar transition is second-order. On increasing the polar interaction strength B, the polar phase becomes stable even at higher temperature. For B = 0.36, these two transitions merge into a single first-order transition directly from the high-temperature isotropic phase to the low-temperature polar phase. If E = 0, the polarization and splay are nonzero even in the nematic phase and scale with the magnitude of the field. For nonzero field, the magnitude of the splay increases on reducing temperature and is enhanced greatly near the transition to the polar phase.
Note that the numerical mean-field results of Fig. 6 are very similar to the Monte Carlo simulation results of Fig. 2, both in the overall phase diagram and in the splay and polar response to an electric field. This similarity demonstrates that the mean-field theory captures the essential physics of this model.

V. DISCUSSION
In conclusion, we have developed a lattice model for splay flexoelectricity in a system of uniaxial pear-shaped molecules. This model predicts a phase diagram showing isotropic, nematic, and polar phases, and it further predicts a converse flexoelectric effect in the nematic phase. The converse flexoelectric effect is proportional to the applied electric field, and it increases dramatically as the temperature decreases toward the nematic-polar transition. Indeed, we can regard this effect as a susceptibility to an applied field, which diverges at the second-order transition to a polar phase. Thus, flexoelectricity is not just a molecular effect arising from the microscopic inter-action of liquid crystals with a field. Rather, it is a statistical effect associated with the response of correlated volumes of molecules, which increases as one approaches a polar phase.
The recent experiments of Harden et al. [2,3] have found an anomalously large bend (rather than splay) flexoelectric effect in systems of bent-core liquid crystals. We speculate that the same considerations discussed in this paper can explain the large bend flexoelectric coefficient in those experiments. The bent-core liquid crystal should be close to a polar phase, with order in the transverse dipole moments of the molecules. As a result, there should be large correlated volumes of molecules, leading to a high susceptibility to an applied field, which induces both polar order and bend. The detailed theoretical model for bend flexoelectricity will necessarily be more complex than the model for splay flexoelectricity presented here, because the bent-core molecules are not uniaxial and hence their orientations must be characterized by two vectors or three Euler angles. This model for bend flexoelectricity will be the subject of a future paper.