Bistable director alignments of nematic liquid crystals confined in frustrated substrates
Takeaki Araki Jumpei Nagura









    ∂ϵ≳≲ε∇δ̣
        
        ^1Department of Physics, Kyoto University, Kyoto 606-8502, Japan
        ^2CREST, Japan Science and Technology Agency, Japan
        
        We studied in-plane bistable alignments of nematic liquid crystals 
        confined by two frustrated surfaces by means of Monte Carlo simulations 
        of the Lebwohl-Lasher spin model. 
        The surfaces are prepared with orientational checkerboard patterns, 
        on which the director field is 
        locally anchored to be planar yet orthogonal between the 
        neighboring blocks. 
        We found the director field in the bulk tends to be aligned along the 
        diagonal axes of the checkerboard pattern, 
        as reported experimentally  
        [J.-H. Kim  et al., Appl. Phys. Lett.  78, 3055 (2001)]. 
        The energy barrier between the two stable orientations 
        is increased, when the system is brought to the isotropic-nematic 
        transition temperature.  
        Based on an elastic theory, we found that the bistability 
        is attributed to the spatial modulation of the director field 
        near the frustrated surfaces. 
        As the block size is increased and/or the elastic modulus is reduced, 
        the degree of the director inhomogeneity is increased, enlarging the 
        energy barrier. 
        We also found that the switching rate between the stable states is 
        decreased when the block size is comparable to the cell thickness. 
        
        
        
        
        64.70.mf,
        61.30.Hn, 
        64.60.De, 
        42.79.Kr 
        Bistable director alignments of nematic liquid crystals confined 
        in frustrated substrates
            Takeaki Araki^1,2 and Jumpei Nagura^1
            December 30, 2023
        ==========================================================================================
        
        
        
        § INTRODUCTION
        
        
        
        Liquid crystals have been utilized in many applications. 
        In particular, they are widely used 
        in optical devices such as flat panel displays 
        <cit.>. 
        Because of the softness of the liquid crystal, 
        its director field is deformed by 
        relatively weak external fields <cit.>. 
        To sustain the deformed state, the external field 
        has to be constantly applied to the liquid crystal substance. 
        In order to reduce power consumption, 
        a variety of liquid crystal systems showing 
        multistable director configurations or storage effects 
        have been developed <cit.>. 
        In such systems, 
        a pulsed external field can induce permanent changes of the 
        director configurations. 
        Liquid crystals of lower symmetries, 
        such as 
        cholesteric, ferroelectric and flexoelectric phases, 
        are known to show the storage effects 
        <cit.>. 
        
        
        
        A nematic liquid crystal 
        in a simple geometry,  e.g. that  
        sandwiched 
        between two parallel plates with homeotropic anchoring, 
        shows a unique stable director configuration if external 
        fields are not imposed. 
        By introducing elastic frustrations, 
        the nematic liquid crystals can have 
        different director configurations 
        of equal or nearly equal elastic energy 
        <cit.>. 
        For instance, 
        either of horizontal or vertical director orientation 
        is possibly formed 
        in nematic liquid crystals confined 
        between two flat surfaces of uniformly tilted but oppositely directed 
        anchoring alignments <cit.>.  
        Also, it was shown that the nematic liquid crystal confined in 
        porous media shows a memory effect <cit.>. 
        The disclination lines of the director field 
        can adopt a large number of trajectories running through the 
        channels of the porous medium 
        <cit.>. 
        The prohibition of spontaneous changes of the defect 
        pattern among the possible trajectories leads to the 
        memory effect. 
        
        
        
        Recent evolutions of micro- and nano-technologies 
        enable us to tailor substrates of inhomogeneous 
        anchoring conditions, the length scale of which 
        can be tuned less than the wavelength of visible light. 
        With them, many types of structured surfaces for the liquid crystals 
        and the resulting 
        director alignments have been reported in the past 
        few decades 
        <cit.>. 
        For example, a striped surface, in which the 
        homeotropic and planar anchorings appear alternatively, 
        was used to control the polar angle of the 
        director field in the bulk <cit.>. 
        
        Kim  et al. demonstrated 
        in-plane bistable alignments 
        by using a nano-rubbing technique with an atomic 
        force microscope <cit.>. 
        They prepared 
        surfaces of orientational checkerboard patterns. 
        The director field in contact to the surfaces 
        is imposed to be parallel to the 
        surface yet orthogonal between the neighboring domains. 
        They found that the director field far from the surface tends to be 
        aligned along either of the two diagonal axes of the 
        checkerboard pattern. 
        More complicated patterns are also possible to prepare <cit.>. 
        
        In this paper, we consider the mechanism of 
        the bistable orientations 
        of the nematic liquid crystals 
        confined in two flat surfaces of the checkerboard 
        anchoring patterns. 
        We carried out Monte Carlo simulations 
        of the Lebwohl-Lasher spin model <cit.>
        and argued their results 
        with a coarse-grained elastic theory. 
        In particular, the dependences of the stability of the 
        director patterns on the temperature, and the domain size 
        of the checkerboard patterns are studied. 
        Switching dynamics between the stable configurations 
        are also considered. 
        
        
        
        
        § SIMULATION MODEL
        
        
        We carry out lattice-based Monte Carlo simulations 
        of nematic liquid crystals confined by two parallel plates 
        <cit.>. 
        The confined space is composed of three-dimensional 
        lattice sites (L× L × H) 
        and it is denoted by ℬ. 
        Each lattice site i has a unit 
        spin vector u_i (|u_i|=1), 
        and the spins are mutually interacting with those at the 
        adjacent sites. 
        At z=0 and z=H+1, 
        we place substrates, 
        composed of two-dimensional lattices.
        We put unit vectors d_j on the site j on 𝒮, 
        where 𝒮 represents 
        the ensemble of the substrate lattice sites. 
        We employ the following Hamiltonian for u_i, 
        ℋ   =   
        -∑_i,j (∈ℬ)
        P_2(u_i·u_j)
              
        -∑_i∈ℬP_2(u_i·e)
        -w∑_i,j,i∈ℬ, j∈𝒮
        P_2(u_i·d_j), 
        
        where P_2(x)=3(x^2-1/3)/2 is the second-order Legendre function and
        ∑_i,j means 
        the summation over the nearest neighbor site pairs. 
        We have employed the same Hamiltonian to 
        study the nematic liquid crystal confined in 
        porous media <cit.>. 
        
        The first term of the right hand side of Eq. (<ref>) 
        is the Lebwohl-Lasher potential, 
        which describes the isotropic-nematic 
        transition <cit.>. 
        In Fig. <ref>, we plot the temperature dependences of 
        (a) the scalar nematic order parameter S_ b 
        and (b) the elastic modulus K 
        in a bulk system. 
        The numerical schemes for measuring them 
        are described in Appendix A. 
        We note that 
        a cubic lattice with periodic boundary conditions (L^3=128^3) 
        is used for obtaining S_ b and K in Fig. <ref>.
        As the temperature is increased, 
        both the scalar order parameter and 
        the elastic modulus are decreased and show abrupt 
        drops at the transition temperature T=T_ IN, 
        which is estimated as 
        k_ BT_ IN/≅ 1.12<cit.>. 
        
        
        
        
        
        
        
        
        The second term of Eq. (<ref>) is the coupling between the 
        spins in ℬ and an in-plane external field 
        e. 
        The last term represents the 
        interactions between the bulk spins and the surface directors, 
        that is, the 
        Rapini-Papoular type anchoring effect <cit.>. 
        w is the strength of the anchoring interaction. 
        If d_j is parallel to the substrates and w>0, 
        the planar anchoring conditions are imposed to the 
        spins at the ℬ-sites contacting to 𝒮. 
        This term not only gives the angle dependence of the anchoring effect 
        in the nematic phase, but also enhances the nematic order near the surface. 
        In Fig. <ref>(a), we also plot the scalar nematic order parameter 
        on a homogenuous surface of w=ϵ. 
        The definition of S_ w is described in Appendix A. 
        The nematic order on the surface is larger than that in the bulk S_ b
        and is decreased continuously with T. 
        Even at and above T_ IN, 
        S_ w does not vanish to zero. 
        When the temperature is far below T_ IN, 
        on the other hand, it is close to that in the bulk S_ b.
        
        
        
        
        In this study, we prepare two types of anchoring cells. 
        In type I cells, we set hybrid substrates. 
        At the bottom surface (z=0), the preferred direction d_j 
        is heterogeneously patterned like a checkerboard as given 
        by 
        d_j(x,y)
        = {[                     (0,1,0) if ([x/D]+[y/D])   is  even;                     (1,0,0)  if ([x/D]+[y/D])   is  odd;                             ].,
        
        where [X] stands for the largest integer smaller 
        than a real number X.
        D is the unit block size of the checkerboard pattern. 
        At the top surface (z=H+1), on the other hand, 
        the preferred direction is homogeneously set to 
        d_j=d_ t≡(cosϕ_ tsinθ_ t,
        sinϕ_ tsinθ_ t,cosθ_ t). 
        θ_ t and ϕ_ t are the polar and azimuthal 
        angles of the preferred direction at the top surface. 
        In type II cells, 
        both substrates are patterned like the checkerboard, 
        according to Eq. (<ref>). 
        
        We perform Monte Carlo simulations with 
        heat bath samplings. 
        A trial rotation of the i-th spin 
        is accepted, considering the local 
        configurations of neighboring spins, 
        with the probability 
        p(Δℋ)=1/(1+e^Δℋ/k_ BT), 
        where Δℋ is the difference of the 
        Hamiltonian between before and after the trial rotation. 
        The physical meaning of the temporal evolution of Monte 
        Carlo simulations is sometimes a matter of debate. 
        However, we note that the method is known to be 
        very powerful and useful for studying glassy 
        systems with slow relaxations, 
        such as a spin glass <cit.>, 
        the dynamics of which is dominated by activation processes 
        overcoming an energy barrier. 
        
        In this study 
        we fix the anchoring strengths at both 
        the surfaces to w=, for simplicity. 
        The lateral system size is L=512 and the 
        thickness H is changed. 
        For the lateral x- and y- directions, 
        the periodic boundary conditions 
        are employed. 
        
        
        
        
        
        
        
        
        § RESULTS AND DISCUSSIONS
        
        
         §.§ Bistable alignments
        
        
        First, 
        we consider nematic liquid crystals confined 
        in cells with the hybrid surfaces (type I). 
        Figure <ref>(a) 
        plots the energies stored in the cell 
        with respect to the azimuthal anchoring angle 
        ϕ_ t. 
        Here the polar anchoring angle is fixed to θ_ t=π/2. 
        The energy per unit area 
        ℰ is calculated as 
        ℰ(θ_ t,ϕ_ t)
        =ℋ(θ_ t,ϕ_ t)/L^2-
        ℰ_ min, 
        where X means the spatial average of a variable X. 
        ℰ_ min is the lowest energy defined as 
        ℰ_ min=min_θ_ t,ϕ_ tℋ/L^2 
        at each temperature (see below).   
        ℰ is obtained after 5× 10^4 Monte Carlo steps 
        (MCS) in the absence of external fields. 
        The cell thickness is H=16 and the block size is D=8. 
        The temperature is changed. 
        
        Figures <ref>(a) indicates the 
        energy has two minima at ϕ_ t=±π /4, 
        while it is maximized at 
        ϕ_ t=0 and ±π/2. 
        To see the dependence on 
        the polar angle, we plot ℰ against θ_ t 
        with fixing ϕ_ t=π/4 in Fig. <ref>(b). 
        It is shown that ℰ 
        is minimized at θ_ t=π/2 for ϕ_ t=π/4. 
        Hence we conclude that 
        the stored energy is globally lowest 
        at (θ_ t,ϕ_ t)=(π/2,±π/4), 
        so that we set 
        ℰ_ min
        =ℋ(θ_ t=π/2,ϕ_ t=π/4)/L^2 
        in Fig. <ref>. 
        This global minimum indicates that 
        the parallel, yet bistable configurations of the director field 
        are energetically preferred in this cell. 
        This simulated bistability is in accordance with the 
        experimental observations reported by Kim  et al.<cit.>. 
        When a semi-infinite cell is used, 
        the bistable alignments of the director field would be realized. 
        Hereafter, we express these two stable directions with 
        n̂_+ and n̂_-. 
        That is, n̂_±=(1/√(2),± 1/√(2),0). 
        
        
        
        
        
        Figure <ref> also indicates the temperature dependencies 
        of the stored energies. 
        When the temperature 
        is much lower than the transition temperature T_ IN, 
        the curves of ℰ are 
        are rather flat. 
        As the temperature is increased, 
        the dependence becomes more remarkable. 
        Figure <ref>(a) 
        plots the energy difference between the 
        maximum and minimum of ℰ 
        for fixed θ_ t=π/2 
        as functions of T. 
        It is defined by 
        the in-plane rotation of d_ t 
        as Δℰ=
        ℰ(θ_ t=π/2,ϕ_ t=π/4)
        -ℰ(π/2,0).  
        We plot them for several block sizes D, while 
        the cell thickness is fixed to H=16. 
        
        In Fig. <ref>(a), we observe non-monotonic dependences of 
        the energy difference on the temperature. 
        Δℰ is almost independent of T 
        when T/T_ IN<0.6. 
        In the range of 
        0.6 ≲ T/T_ IN < 0.9, 
        it is increased with increasing T. 
        When T/T_ IN≳ 0.9, 
        it decreases with T and it almost disappears if T>T_ IN. 
        When T>T_ IN, the system is in the isotropic state, 
        and it does not have the long-range order. 
        Thus, it is reasonable that Δℰ vanishes 
        when T>T_ IN. 
        When T<T_ IN, 
        on the other hand, 
        it is rather striking that the energy difference 
        shows the non-monotonic dependences on T, 
        in spite of that the long-range order and the resultant 
        elasticity are reduced monotonically with increasing T 
        (see Fig. <ref>). 
        
        We plot the energy difference Δℰ 
        as a function of D in Fig. <ref>(b), 
        where the temperature is T/T_ IN=0.89. 
        The cell thickness is changed. 
        It is shown that the energy difference Δℰ is increased 
        proportionally 
        to the block size D when D is small. 
        When the cell thickness is large, on the other hand, 
        the energy difference is almost 
        saturated. 
        The saturated value becomes smaller when the liquid crystal is 
        confined in the thicker cell. 
        
        
        
        
        In order to clarify the mechanisms of the bistable alignments, 
        we calculate the spatial distribution of the nematic order parameter. 
        In Fig. <ref>, 
        we show snapshots of xx- and xy components of a 
        tensorial order parameter at several planes parallel to the 
        substrates. 
        Using u_i(t'), the tensorial order parameter 
        Q_μν is calculated  
        by averaging 3/2(u_μ u_ν/2-δ_μν/3) 
        in  a certain period δ t as, 
        
        Q_i,μν(t)=1/δ t∑_t'=t^
        t+δ t-13/2{u_i,μ(t')u_i,ν(t')-1/3δ_μν}, 
        
        where t' means the Monte Carlo cycle, and 
        μ and ν stand for x, y and z. 
        In this study, we set δ t=10^2, 
        which is chosen so that the system is well thermalized. 
        The block size is 
        D=8 in (a) and D=64 in (b), 
        and the cell thickness is fixed to H=8. 
        The temperature is set to T/T_ IN=0.89. 
        The anchoring direction at the top surface is along n̂_+, 
        and we started the simulation with an initial condition, 
        in which the director field is 
        along 
        n̂_+, 
        so that the director field is likely 
        to be parallel to the surface 
        and along the azimuthal angle ϕ=π/4 in average.  
        
        Q_xx near the bottom surface 
        shows the checkerboard pattern 
        as like as that of the imposed anchoring directions d_j. 
        Q_xy inside the block domains is small 
        and it is enlarged at the edges between the blocks. 
        With departing from the bottom surface, 
        the inhomogeneity is reduced and the director pattern 
        becomes homogeneous along n̂_+. 
        The inhomogeneities in Q_xx and Q_xy 
        are more remarkable for the larger D than 
        those for the smaller D. 
        
        
        
        
        In Fig. <ref>, 
        we plot the corresponding 
        profiles of the spatial modulations of 
        the order parameter with respect to z. 
        The degree of the inhomogeneity of Q_μν is defined by 
        
        I(z)   =   1/L^2S_ b^2∫ dx dy {Q_μν(x,y,z)-Q̅_μν(z)}^2,
        
        where Q̅_μν(z) is the spatial average of Q_μν 
        in the z-plane and it is given by 
        Q̅_μν(z)   =   L^-2∫ dx dy Q_μν(x,y,z).
        S_ b is the scalar nematic parameter obtained in the bulk 
        [see Fig. <ref>(a)]. 
        Since Q_μν∝ S_ b in the bulk,
        the profiles are scaled by S^2_ b in Eq. (<ref>),
        in order to see the pure degree of the inhomogeneity 
        of the director field.
        In Fig. <ref>(a), we changed the block size D and the 
        temperature is fixed to T/T_ IN=0.89. 
        It is shown that the degree of the inhomogeneity 
        decays with z, and it is larger for 
        larger D as shown in Fig. <ref>. 
        Figure <ref>(a) also shows that the 
        decaying length is also increased with the block size D. 
        Roughly it agrees with D. 
        In Fig.<ref>(b), we plot the profiles of I(z) 
        for different temperatures with fixing D=16. 
        It is shown that 
        the spatial modulation is increased as the temperature is increased. 
        This is because 
        the nematic phase becomes softer as the temperature 
        is increased (see Fig. <ref>(b)). 
        When the elastic modulus is small, the director field 
        is distorted by the anchoring surface more largely. 
        
        
        
        
        
        
        
        
        
        
        Based on these numerical results, 
        we consider the bistable alignments with a 
        continuum elasticity theory. 
        The details of the continuum theory is described in 
        Appendix B. 
        In our theoretical argument, 
        the spatial modulation of the director field due to the 
        heterogeneous anchoring 
        plays a crucial role in inducing the bistable 
        alignments along the diagonal directions. 
        After some calculations, 
        we obtained an effective anchoring energy for D≪ H as 
        
        g(ϕ_0)=-cW^2D/Ksin^22ϕ_0,
        
        instead of the Rapini-Papoular anchoring energy, -Wcos^2ϕ_0/2. 
        Here ϕ_0 is the average azimuthal angle 
        of the director field on the patterned surface.
        K is the elastic modulus of the director field 
        in the one-constant approximation of the elastic theory, 
        and W represents the anchoring strength in the 
        continuum description. 
        c is a numerical factor, which is estimated as 
        c≅ 0.085 when H/D is large. 
        g(ϕ_0) has a fourfold symmetry and is lowered 
        at ϕ_0=±π/4 and ± 3π/4. 
        The resulting energy difference per unit area is given by  
        Δℰ_ th=
        π^2K/32H{1+K^2/(8cW^2DH)}.
         
        
        First we discuss the dependence of the energy difference 
        on the block size D.
        Equation (<ref>) indicates that the energy difference 
        behaves as Δℰ_ th≈π^2 cW^2D/(4K), 
        which is increased linearly with D,  
        when D is sufficiently small. 
        If D is large enough, on the other hand, 
        the energy difference converges to Δℰ_ th≈π^2 K/(32H). 
        The latter energy difference agrees with the deformation 
        energy of the 
        director field, which twists along the z axis by ±π/4. 
        It is independent of D, but 
        is proportional to H^-1. 
        The asymptote behaviors for small and large D 
        are consistent with the numerical results 
        shown in Fig. <ref>(b). 
        
        
        
        
        
        
        
        
        Next we consider the dependence of Δℰ on the 
        temperature. 
        Equation (<ref>) also suggests 
        Δℰ_ th is proportional to W^2D/K 
        when W^2DH/K^2 is small. 
        We have speculated the anchoring strength 
        is simply proportional to the nematic order as 
        W∝ S_ b. 
        If so, the energy difference is expected to be 
        independent of S_ b 
        as Δℰ_ th∝ W^2/K 
        ∝ S^0_ b, 
        since K is roughly proportional to S_ b^2. 
        This expectation is inconsistent with the dependence of 
        the numerical results of Δℰ in Fig. <ref>(a). 
        A possible candidate mechanism in explaining this discrepancy 
        is that 
        we should use the nematic order on the surface S_ w, 
        instead of S_ b, for estimating W. 
        Since S_ w is dependent on T more weakly 
        than S_ b near the transition temperature [see Fig. <ref>(a)], 
        W^2/K can be increased with T. 
        The curve of S_ w^2/K is drawn in Fig. <ref>(b). 
        Thus, 
        the director field is more largely deformed near T_ IN as shown 
        in Fig. <ref>(b), 
        so that the resulting energy difference shows the 
        increase with T. 
        
        Also, Fig. <ref>(a) shows 
        Δℰ turns to decrease to zero, 
        when we approach to T_ IN more closely. 
        In the vicinity of T_ IN, 
        K is so small that W^2/K becomes large. 
        Then Eq. (<ref>) behaves as 
        Δℰ_ th∝ K/H. 
        It is decreased to zero as K with approaching to T_ IN. 
        In Fig. <ref>(a), we draw the theoretical 
        curve of Eq. (<ref>) with taking into account 
        the dependences of W and K on the 
        temperature. 
        Here we assume W=W_0 S_ w with W_0 being a 
        constant. 
        The theoretical curve reproduces the non-monotonic 
        behavior of the energy difference qualitatively. 
        After the plateau of Δℰ_ th 
        in the lower temperature region, 
        it is increased with T. 
        Then it turns to decrease to zero when the temperature 
        is close to the transition temperature. 
        Here we use W_0=0.3, which is chosen to adjust 
        the theoretical curve to the numerical result. 
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ Switching dynamics
        
        
        
        Next we confine the nematic liquid crystals in the 
        type II cells, both the surfaces of which 
        are patterned as checkerboard. 
        As indicated by Eq. (<ref>), 
        each checkerboard surface gives rise to the effective 
        anchoring effect with the fourfold symmetry. 
        Hence, 
        the director field is expected to 
        show the in-plane bistable alignments along 
        n_+ or n_- also in the type II cells. 
        
        
        Figure <ref>(a)
        plots 
        the spatial average of the 
        xy component of Q_μν at equilibrium 
        with respect to the block size. 
        The equilibrium value of Q_xy is 
        estimated as Q_xy^∞=Q_xy|_t=5× 10^4 
        in the simulations with no external field. 
        As the initial condition, we employ the director field 
        homogeneously aligned along n̂_+, 
        so that Q_xy^∞ is likely to be positive. 
        In Fig. <ref>(a), 
        we also draw a line of 3S_ b/4, 
        which corresponds to the bulk nematic order when the 
        director field is along n_+. 
        It is shown that 
        Q_xy^∞ is roughly constant and is 
        close to 3S_ b/4 for D≪ H. 
        It is reasonable since the inhomogeneity of the director field 
        is localized within D from the surfaces. 
        When D>H, on the other hand, 
        Q_xy^∞ is decreased with D. 
        When D≫ H, the type II cell can be considered as a 
        collection of square domains each carrying the uniform 
        anchoring direction. 
        Thus, the director field 
        tends to be parallel to the local anchoring direction 
        d_j, and then, Q_xy inside each unit block 
        becomes small locally. 
        Only on the edges of the block domains, the director fields 
        are distorted and adopts either of the distorted 
        states as schematically shown in Fig. <ref>(b). 
        With scaling D by H, the plots of Q_xy^∞ 
        collapse onto a single curve. 
        
        
        
        
        
        
        
        
        
        
        
        Then we consider the switching dynamics of the director field 
        between the two stable alignments with imposing in-plane 
        external fields e in the type II cells. 
        In Fig. <ref>, 
        we plot the spatial average of the xy component of the 
        order parameter Q_xy 
        in the processes 
        of the director switching. 
        The cell size is H=16, the block size is D=16 and 
        the temperature is T/T_ IN=0.89. 
        At t=0, we start the Monte Carlo simulation with the same initial condition, 
        in which the director field is homogeneously aligned along 
        n̂_+, 
        in the absence of the external field. 
        As shown in Fig. <ref>, 
        the nematic order is relaxed to 
        a certain positive value, 
        which agrees with Q_xy^∞ in Fig. <ref>(a). 
        From t_1=10^4, 
        we then impose an in-plane external field along 
        n̂_-, and turn it off at t_2=2× 10^4. 
        After the system is thermalized during t=2× 10^4 and 
        t_3=3× 10^4 with no external field, 
        we apply the second external field along n̂_+ 
        from t_3=3× 10^4 
        until t_4=4× 10^4. 
        We change the strength of the external field e. 
        
        
        When the external field is weak (e^2≤ 0.03), 
        the averaged orientational order is slightly reduced by the external field, 
        but it recovers the original state after the field is removed. 
        After a strong field (e^2≥ 0.04) is applied and 
        is removed off, on the other hand, 
        Q_xy 
        is relaxed to another steady state value, which is close to 
        -Q_xy^∞. 
        This new state of the negative Q_xy 
        corresponds to the other bistable alignment along 
        n̂_-.
        After the second field along 
        n̂_+ is applied, 
        the averaged orientational order Q_xy 
        comes back to the positive original value, +Q_xy^∞.
        
        
        
        
        
        In Fig. <ref>(a), 
        we show the detailed relaxation behaviors of 
        Q_xy in the first switching after t_1. 
        Δ t means the elapsed time in the first switching, 
        that is Δ t=t-t_1. 
        Here 
        we change the block size D, while 
        we fix the external field at e^2=0.03 and 
        the cell thickness H=16 (type II). 
        We note that Q_xy at Δ t=0 
        depends on D as indicated in Fig. <ref>(a). 
        In Fig. <ref>(a), it is shown that the switching rate depends 
        also on the block size D. 
        Notably, the dependence of the switching behavior is not 
        monotonic against D. 
        
        
        
        In Fig. <ref>(b), we plot the characteristic switching time 
        τ with respect to the block size D in the cells of 
        H=8, 16 and 32. 
        The temperature and the field strength are the same those for 
        Fig. <ref>(a).
        The characteristic time τ is defined such that the average orientational 
        order is equal to zero at τ, Q_xy(Δ t=τ)=0. 
        Figure <ref>(b) shows the characteristic time is maximized 
        when the 
        block size is comparable to the cell thickness. 
        When D<H, the switching process is slowed down as the block size 
        is increased. 
        On the other hand, it is speeded up with D when D>H. 
        In Fig. <ref>(b), it is suggested that 
        the dependence of τ on D becomes less significant 
        as H is increased. 
        
        
        
        
        
        Figure <ref> depicts 
        snapshots of 
        Q_xy(t)
        at the midplane (z=H/2) during the first switching process. 
        The parameters are the same as those in Fig. <ref>(a), 
        so that the pattern evolutions 
        correspond to the curves of Q_xy in Fig. <ref>(a). 
        Figure <ref> shows that the switching behavior  
        is slowed down when D is comparable to H, 
        in accordance with Fig. <ref>(b). 
        When D<H, the snapshots implies the switching proceeds 
        via nucleation and growth mechanism. 
        From the sea of the positive Q_xy, where the director is 
        aligned along n̂_+, 
        the droplets of the negative Q_xy are nucleated. 
        They grow with 
        time and cover the whole area eventually. 
        Under the external field along n̂_-, 
        the alignment of the director field along n̂_- 
        is more preferred than that along n̂_+. 
        Because of the energy barrier between these bistable alignments, 
        the director field cannot change its 
        orientation to n̂_- smoothly under 
        a weak external field. 
        From Eq. (<ref>), the energy barrier 
        for the local swiching of the director field 
        between the two stable states is given by 
        Δℱ=8D^2 {g(ϕ_0=0)-g(π/4)}=8cW^2D^3/K, 
        when D<H. 
        Thus, the slowing down of the switching process 
        with D is considered to be 
        attributed to the enhancement of the energy barrier. 
        Here we note that 
        a critical field strength 
        for the thermally activated switching cannot be defined unambiguously.
        Since the new alignment is energetically preferred over the 
        original one even under a weak field, 
        the director configuration will change its orientation 
        if the system is annealed for a sufficiently long period. 
        When the field strength is moderate (e^2≅ 0.035), 
        the averaged order 
        goes to an intermediate value, 
        neither of Q_xy^∞ or -Q_xy^∞ in Fig. <ref>. 
        Such intermediate values of Q_xy reflect large scale 
        inhomogeneities of the bistable alignments (see Fig. <ref>). 
        At each block, the director field adopts either of the 
        two stable orientations. 
        The pattern of the intermediate Q_xy 
        depends not only on the field strength, but also 
        on the annealed time. 
        Under large external fields, 
        on the other hand, 
        the energy barrier between the two 
        states can be easily overcome, so that  
        the switching 
        occurs without arrested at the initial orientation 
        (not shown here). 
        
        Regarding the local director field, which adopts 
        either of the two stable orientations (n̂_+ and 
        n̂_-), as a binarized spin 
        at the corresponding block unit, we found a similarity 
        of the domain growth in our system and that in a two-dimensional 
        Ising model subject to an external magnetic field. 
        If the switching of the director field occurs locally 
        only at each block unit, 
        there is no correlations between the director fields 
        in the adjacent block units. 
        Therefore, the nucleation and growth switching behavior 
        implies the director field at a block unit prefers to be aligned 
        along the same orientation as those at the adjacent block units. 
        
        We observed  string-like patterns 
        as shown at Δ t=4000 for D=1 in 
        Fig. <ref>. 
        Here we note that they are not 
        disclinations of the director field. 
        They represent domain walls perpendicular to the substrates. 
        In the type II cells, we have not observed 
        any topological defects, 
        although topological defects are sometimes stabilized 
        in the frustrated cell <cit.>. 
        The string-like patterns remain rather stable transiently. 
        On the other hand, 
        such string-like patterns are not observed in the switching process 
        in the Ising model. 
        This indicates  the binarized spin description of the bistable director 
        alignments may be not adequate. 
        Under the external field along n̂_-, 
        the director field rotates to the new orientation 
        clockwise or counter-clockwise. 
        New domains, which appear via the clockwise rotations, 
        have some mismatches against those through the counter-clockwise 
        rotations. 
        The resulting boundaries between the incommensurate domains 
        are formed and 
        tend to suppress the coagulations of them more and less, 
        although the corresponding energy barriers are not so large. 
        
        
        
        When D>H, 
        the switching occurs in a different way. 
        The director rotations are localized around the 
        edges of the blocks 
        as indicated in Fig. <ref>(b). 
        As D is increased, the amount of the 
        director field that reacts to the field 
        is reduced. 
        Although 
        the director fields 
        around the centers of the blocks
        do not show any switching behaviors 
        before and after the field application, 
        they are distorted to orient slightly 
        toward the field. 
        It is considered that the distortion of the 
        director field inside the blocks effectively reduces 
        the energy barrier against the external field.
        We have not succeeded in explaining 
        the mechanism of the reduction of the 
        switching time with D. 
        When D>H, the inhomogeneous director field 
        contains higher Fourier modes of the distortion. 
        The energy barrier for 
        each Fourier mode becomes lower for the higher 
        Fourier modes [see Eq. (<ref>). 
        Thus, such higher Fourier modes are more active 
        against the external field and they would 
        behave as a trigger of the switching process. 
        
        
        
        
        
        
        
        
        
        
        
        § CONCLUSION
        
        
        In this article, we studied nematic liquid crystals confined 
        by two parallel checkerboard substrates by means of 
        Monte Carlo simulation of the Lebwohl-Lasher model. 
        As observed experimentally by Kim  et al.,
        we found the director field in the bulk shows the bistable alignments, 
        which are along 
        either of the two diagonal axes. 
        We attribute the bistability of the alignments to the 
        spatial modulation of the director field near the substrates. 
        Based on the elastic theory, we derived 
        an effective 
        anchoring energy with the fourfold symmetry (Eq. (<ref>)). 
        Its anchoring strength is expected to behave as 
        W^2D/K, when the block size D is smaller than 
        the cell thickness. 
        As the temperature is increased to the isotropic-nematic transition 
        temperature, the elastic modulus K of the nematic phase is reduced 
        so that the director field is deformed near the substrates more largely. 
        With this effective anchoring effect, we can explain 
        the non-monotonic dependence of the energy stored in this cell 
        qualitatively. 
        
        
        We also studied the switching dynamics of the director configuration 
        with imposing in-plane external fields. 
        Usually, the switching is considered to be associated with the 
        actual breaking of the anchoring condition. 
        Thus, the energy barrier for the switching is expected to be 
        proportional to W<cit.>. 
        In this article, we propose another possible mechanism of the 
        switching, in which the anchoring condition is not necessarily broken. 
        Since the energy barrier is increased with the block size, 
        the switching dynamics notably becomes slower when the block size is 
        comparable to the cell thickness. 
        
        By solving Δℱ=Δϵ E^2/2× (4D^2H), 
        we obtain a characteristic strength of the electric field E 
        as 
        E_ c≅{8cW^2D/(Δϵ KH)}^1/2, 
        where Δϵ is the anisotropy of the dielectric 
        constant [see Eq. (<ref>)]. 
        If we apply an in-plane external field larger than E_ c, 
        the switching occurs rather homogeneously 
        without showing the nucleation and growth processes. 
        This characteristic strength is decreased with decreasing D, 
        so that the checkerboard pattern of smaller D is preferred 
        to reduce the field strength. 
        With smaller D, however, the stability of the two preferred 
        orientations is reduced. 
        If the effective anchoring energy is lower than the thermal energy, 
        the 
        bistable alignment will be destroyed by the thermal fluctuation. 
        In this sense, the block size D should be larger 
        than D_ c≈ (Kk_ BT/8cW^2)^1/3, where we 
        assumed that the switching occurs locally in each block, 
        that is Δℱ≈ k_ BT.
        For a typical nematic liquid crystal with K=1  pN and 
        W=10^-5  J/m^2 at room temperature T=300  K, 
        it is estimated as D_ c≅ 34  nm. 
        
        
        
        In our theoretical argument, 
        we assumed the one-constant approximation of the elastic modulus. 
        However, the director field cannot be described 
        by a single deformation mode in the above cells. 
        The in-plane splay and bend deformations are localized 
        within the layer of D near the surface. 
        On the other hand, 
        the twist deformation is induced by the external field 
        along the cell thickness direction. 
        If the elastic moduli for the three deformation modes 
        are largely different from each other, 
        our theoretical argument would be invalid. 
        We need to improve both the theoretical and numerical 
        schemes to consider such dependences more correctly. 
        Also, we considered only the checkerboard substrates. 
        But, it is interesting and important 
        to design other types of patterned surfaces <cit.> 
        to append more preferred functions, 
        such as faster responses against the external field, 
        to liquid crystal devices. 
        We hope to report a series of such studies in the near future. 
        
        
        
        
        § ACKNOWLEDGEMENTS
        
        
        We acknowledge valuable discussions with J. Yamamoto, H. Kikuchi, 
        I. Nishiyama, K. Minoura and T. Shimada. 
        This work was supported by KAKENHI 
        (No. 25000002 and No. 24540433). Computation was done using the
        facilities of the Supercomputer Center, the Institute for Solid
        State Physics, the University of Tokyo. 
        
        
        
        
        § ESTIMATIONS OF THE NEMATIC ORDER AND THE ELASTIC MODULS
        
        
        In this appendix, 
        we estimate the scalar nematic order parameter 
        and the elastic modulus in Fig. <ref> 
        from the Monte Carlo simulations with Eq. (<ref>) 
        <cit.>. 
        First we consider the bulk behaviors of nematic liquid crystals, 
        which are described by the Lebwohl-Lasher potential. 
        Here we remove the surface sites 𝒮 and employ 
        the periodic boundary conditions for all the axes (x,y, and z). 
        We use the initial condition along u_i=(1,0,0) and thermalize the 
        system with the heat bath sampling. 
        The simulation box size is L^3 with L=128. 
        
        It is well known that this Lebwohl-Lasher spin model describes 
        the first-order transition between isotropic and nematic phases. 
        In Fig. <ref>, we plot the xx component of the tensorial order 
        parameter after the thermalization (t≤ 5× 10^4) 
        as a function of T. 
        Since the initial condition is along the x axis, 
        the director field is likely to be aligned along the x axis.
        We here regard Q_xx 
        as the scalar nematic order parameter S_ b.
        We see an abrupt change of S_ b around k_ BT_ IN≈ 
        1.12ϵ, which is consistent 
        with previous studies <cit.>. 
        Above T_ IN, the nematic order almost vanishes, while 
        it is increased with decreasing T when T<T_ IN. 
        
        In the nematic phase (T<T_ IN), 
        the director field n can be defined. 
        Because of the thermal noise, 
        the local director field is fluctuating around the 
        average director field, reflecting the elastic modulus. 
        The elastic modulus of the director field is obtained 
        by calculating the scattering function of the 
        tensorial order parameter as <cit.>,
        |Q̃_xμ(q)|^2_T=
        k_ BT/A+4L_Qsin^2 |q|a/2,
        
        for μ=y and z. 
        Q̃_xμ(q) is the Fourier component 
        of Q_xμ at a wave vector q. 
        a is the lattice constant, and ⋯_T means 
        the thermal average. 
        A and L_Q are the coefficients appearing in the 
        free energy functional for Q_μν. 
        
        In the case of T<T_ IN, the scattering function goes to 
        zero for |q|a≅ 0. 
        Then, we obtain the coefficient L_Q by fitting 
        |Q̃_xμ|^2^-1_T 
        with 4(k_ BT)^-1L_Qsin^2 |q|a/2. 
        L_Q is proportional to the elastic modulus K of the director field 
        n as L_Q=KS^2_ b. 
        In Fig. <ref>(b), 
        the elastic modulus K is plotted with respect to T. 
        It is decreased with increasing T, if T>T_ IN. 
        This indicates the softening of the nematic phase near the transition 
        temperature. 
        
        Next we consider the effect of the surface term. 
        The surface effect not only induces the angle dependence 
        of the anchoring effect in the nematic phase, but also leads to the wetting 
        effect of the nematic phase to the surface in the isotropic phase 
        <cit.>. 
        We set homogeneous surfaces of w=ϵ 
        at z=0 and z=H+1 as in the main text. 
        The anchoring direction is d_j=(1,0,0). 
        The periodic boundary conditions are imposed for the x- and y-directions
        and the initial condition is along u_i=(1,0,0). 
        The profile of Q̅_xx (not shown here) indicates 
        Q̅_xx at z=1 becomes larger than that in the bulk, S_ b. 
        This value is the surface order S_ w, which is also plotted in 
        Fig. <ref>(a) with red open squaress. 
        Notably, S_ w remains a finite value 
        even when T>T_ IN. 
        In Fig. <ref>(a), we cannot see any drastic change of S_ w, 
        which is continuously decreased with T. 
        
        
        
        
        § ANALYSIS WITH FRANK ELASTICITY THEORY
        
        
        Here, we consider the nematic liquid crystal 
        confined in the checkerboard substrate 
        on the basis of the Frank elasticity theory. 
        The checkerboard substrate is placed at z=0, 
        while we fix the director field at the top surface 
        like the type I cells employed in the simulations. 
        The free energy of the nematic liquid crystal 
        is given by 
        ℱ   =   K/2∫ dr
        (n)^2-Δϵ/2∫ dr(n·E)^2
        
              -W∫_z=0dxdy(n·d)^2,
        
        where n is the director field. 
        The first term in the right hand side of Eq. (<ref>) 
        is the elastic energy. 
        Here we employ the one-constant approximation 
        with the elastic modulus K. 
        E and Δϵ are external electric field 
        and the anisotropy of the dielectric constant. 
        Here
        we do not consider the effect of the electric field. 
        The third term in Eq. (<ref>) 
        represents the anchoring energy in 
        the Rapini-Papoular form. 
        W is the anchoring strength and d 
        is the preferred direction on the surface at z=0. 
        For the checkerboard substrates, 
        we set d according to Eq. (<ref>). 
        
        At the top surface, we fix the director field 
        as n(z=H)=d_ t(cosϕ_ t,sinϕ_ t,0), 
        and the bottom surface also prefers the planar anchoring. 
        From the symmetry, therefore, 
        we assume that the director field in the bulk lies parallel 
        to the substrates everywhere. 
        Then, we can write it only with the azimuthal angle ϕ as 
        n=(cosϕ,sinϕ,0). 
        
        Also, we assume that the director field 
        is periodic for x and y directions, 
        so that we only have to consider the 
        free energy in the unit block (0≤ x,  y ≤ 2D). 
        With these assumptions, the free energy per unit area
        is written as 
              ℰ
        =K/8D^2∫_0^2Ddx∫_0^2Ddy∫_0^Hdz (ϕ)^2
        
              
        -W/2D^2.∫_0^Ddx
        {∫_0^Ddy sin^2ϕ
        
        +∫_D^2Ddycos^2ϕ}|_z=0.
        
        
        In the equilibrium state, 
        the free energy is minimized with respect to ϕ(x,y,z). 
        Inside the cell (0<z<H), 
        the functional derivative of ℰ gives 
        the Laplace equation of ϕ as
        δℰ/δϕ=-K^2ϕ=0. 
        
        From the symmetry argument, 
        we have its solution as 
        ϕ(x,y,z)   =   ϕ_0+(ϕ_ t-ϕ_0)z/H+
        Δ(x,y,z),
        
        Δ(x,y,z)   =   ∑_m,n=0^∞Δ_mnsin(2m+1)π x/Dsin(2n+1)π y/D
              ×sinh(πγ_mn(H-z)/D), 
        γ_mn   =   √((2m+1)^2+(2n+1)^2)
        where 
        ϕ_0 and Δ_mn are determined later. 
        
        It is not easy to calculate the second term 
        in Eq. (<ref>) analytically. 
        Assuming |Δ|≪ 1, we approximate sin^2ϕ as 
        sin^2(ϕ_0+Δ) 
        ≈sin^2ϕ_0+Δsin 2ϕ_0 +Δ^2cos 2ϕ_0,
        
        Then, we obtain the free energy per unit area as 
        ℰ   =   K/2H(ϕ_ t-ϕ_0)^2-W/2
              +∑_m,n[
        π K Δ_mn^2sinh(2πγ_mnH/D)/16D.
              .
        -4WΔ_mnsin2ϕ_0 sinh(πγ_mnH/D)/(2m+1)(2n+1)π^2]. 
        
        
        
        First we minimize the free energy with respect to Δ_mn by solving  
        ∂ℰ/∂Δ_mn=0. 
        Then, we have 
        Δ_mn=16WDsech
        (πγ_mnH/D)sin2ϕ_0/(2m+1)(2n+1)γ_mnπ^3 K,
        
        and
        ℰ   ≈   K/2H(ϕ_ t-ϕ_0)^2-W/2
        -cW^2D/Ksin^22ϕ_0,
        
        
        c   =   ∑_mn32tanh(πγ_mnH/D)/(2m+1)^2(2n+1)^2π^5γ_mn.
        
        In the limit of H≫ D, c converges to c≈ 0.085, 
        while it behaves as c≈ 0.5H/D if H≪ D. 
        Since c is positive, 
        the last term in the right-hand side of Eq. (<ref>) 
        represents an effective anchoring condition  
        [Eq. (<ref>)] in the main text. 
        It indicates 
        the director field tends to be along the diagonal axes of 
        the checkerboard surface, ϕ_0=±π/4. 
        Then, we minimize ℰ with respect to ϕ_0 
        and obtain 
        ℰ=-W/2-cW^2D/K+K/2H(ϕ_ t∓π/4)^2/1+K^2/(8cW^2DH).
        
        It corresponds to the plots in Fig. <ref>(a). 
        Here we assumed |ϕ_0∓π/4|≪ 1, 
        so that sin^22ϕ_0≅ 1-4(ϕ_0±π/4)^2. 
        The resulting energy difference is obtained as
        Δℰ=π^2K/32H{1+K^2/(8cW^2DH)}.
        
        
        In the strong anchoring limit, 
        we can obtain ϕ rigorously as 
              ϕ(x,y,z)=(ϕ_ t∓π/4)z/H ±π/4
              +4/π∑_m,n1/(2m+1)(2n+1)sin(2m+1)π x/Dsin(2n+1)π y/D
              ×sinh(πγ_mn(H-z)/D)
        
        Its energy difference is then given by Δℰ=π^2K/(32H). 
        It is consistent with Eq. (<ref>) 
        in the limit of W→∞. 
        
        
        apsrev4-1
