Effect of Mechanical Loads on Stability of Nanodomains in Ferroelectric Ultrathin Films: Towards Flexible Erasing of the Non-Volatile Memories

Intensive investigations have been drawn on nanoscale ferroelectrics for their prospective applications such as developing memory devices. In contrast with the commonly used electrical means to process (i.e., read, write or erase) the information carried by ferroelectric domains, at present, mechanisms of non-electrical processing ferroelectric domains are relatively lacking. Here we make a systematical investigation on the stability of 180° cylindrical domains in ferroelectric nanofilms subjected to macroscopic mechanical loads, and explore the possibility of mechanical erasing. Effects of domain size, film thickness, temperature and different mechanical loads, including uniform strain, cylindrical bending and wavy bending, have been revealed. It is found that the stability of a cylindrical domain depends on its radius, temperature and film thickness. More importantly, mechanical loads have great controllability on the stability of cylindrical domains, with the critical radius nonlinearly sensitive to both strain and strain gradient. This indicates that erasing cylindrical domain can be achieved by changing the strain state of nanofilm. Based on the calculated phase diagrams, we successfully simulate several mechanical erasing processes on 4 × 4 bits memory devices. Our study sheds light on prospective device applications of ferroelectrics involving mechanical loads, such as flexible memory devices and other micro-electromechanical systems.

impact on the domain stability of the nanofilm. For example, the film gradually loses switchable domains after repetitive electrical switching, known as fatigue problem, due to the happening of electric-driven processes such as defect generation 22 . For the above complicated factors, although quite an amount of relevant works can be found in literature, stability of ferroelectric domains still needs further investigations to achieve an effective and quantitive control on it. For the simplest domain pattern for memory applications, i.e., 180u nanodomains, previous researches paid attention to dependence of their stability on SPM electric field (e.g., strength and period) [4][5][6][7][8][9] . The minimal domain size that can be achieved in SPM experiments so far is limited by the radius curvature of SPM tip (typically 20-50 nm) 9 . To reach a lower spatial limit of domain stability, exploiting the effects of other factors, such as mechanical loads would be beneficial.
Due to the natural coupling between polarization and electric field, almost all reported reading/writing/erasing mechanisms of ferroelectric domains are electrical. In literature, it has been understood that film-substrate mismatch strains can be utilized to obtain novel ferroelectric phases and domain structures in ferroelectric thin films 17,23,24 . However, few attentions have been paid to the evolution of an existed domain pattern in response to mechanical loads (i.e., treating mechanical loads in a similar role of electric field). Being a phenomenon originated from the collective distortions of lattice unitcells, ferroelectric polarization is intrinsically coupled with mechanical field. Consequently, external strain can change the polarization state, known as piezoelectricity. It has been also realized that strain gradient can affect polarization in a similar way of electric field, namely flexoelectricity (see reviewed paper Ref. 25). These electromechanical effects provide the possibility of processing ferroelectric domains mechanically. Indeed, by exploiting the flexoelectric field of local strain gradient, mechanical writing ferroelectric domains has been recently demonstrated 26 . Considering the possible problems of electrical means (e.g., leakage, heat, dielectric breakdown and fatigue), and potential applications of ferroelectrics in non-electrical environments, effective strategies for processing ferroelectric domains by mechanical means would be meaningful.
In this letter, based on systematical phase field simulations, we demonstrate that stability of polar domains in ferroelectric ultrathin films can be effectively controlled by macroscopic mechanical loads. As a consequence, ultrahigh density of stable polar domains can be obtained. Moreover, erasing of domains can be realized by adjusting the strain state of the nanofilm. With carefully taken into account the effects of inhomogeneous electromechanical fields, ambient temperature and surfaces, stability of 180u cylindrical domains has been comprehensively revealed as a function of domain size, temperature, film thickness and particularly various mechanical loads, such as uniform strain, cylindrical bending and wavy bending. The revealed mechanical controllability of domain stability should provide insight into the fundamental limits of ferroelectric memories and indicate prospective applications of ferroelectrics involving mechanical loads.

Results
The basic idea of our investigation is illustrated in Figure 1. The model systems are (001) oriented PbTiO 3 (PTO) ferroelectric ultrathin films divided into rectangular memory units. The conclusions of this investigation are believed to be applicable to other systems such as Pb(Zr x Ti 1-x )O 3 and BaTiO 3 nanofilms, which have also been predicted to exhibit 180u domains. The PTO films are initially poled into single domain state, with all memory units initially carrying bit information ''0''. Some memory units are then written with 180u inverted cylindrical domains to carry bit information ''1''. The domain size is characterized by a lateral radius r. To see the effects of both strain and strain gradient on the domain stability, and to explore effective mechanical erasing mechanisms, we apply various mechanical loads to the nanofilm written with cylindrical domains, including (i) uniform strain, (ii) cylindrical bending and (iii) wavy bending, respectively. Note that cylindrical domain and control of its size can be achieved in SPM experiments by controlling strength and period of the SPM electric field [4][5][6][7][8][9] . Moreover, from an ideal point of view, the three mechanical loading conditions can be also realized in experiments by growing the PTO thin films on compliant substrates according to their elastic property. We therefore believe that the simulation results in this study can be verified in experiment, which we would nevertheless leave in a future wok.
The critical size of stable cylindrical domain. It is necessary to first gain an insight into the characteristics of 180u cylindrical domain in the nanofilm, which is important to understand its stability behavior. We simulate the equilibrium domain pattern of a 128 nm 3 128 nm 3 8 nm simulation cell at room temperature, which is free of external strain and initially written with a cylindrical domain of size r 5 16 nm (see Supplementary Information on line for the generation of such a domain pattern). After polarization relaxation, the system is found to maintain the cylindrical domain pattern. Figure 2a and 2b depict the distributions of polarization and strain field viewed from   Figure S1 on line). The equilibrium domain pattern has nonzero in-plane polarization components P 1 and P 2 of about 0.03 C/m 2 near the domain wall. Such a polarization rotation at 180u domain wall is in consistent with theoretical calculation 27 and has been observed in experiment 28 . The domain wall is sharp with a maximum width of polarization rotation at the surface being about 2 nm, similar with the experimental measurement 28 . The strain components also have large magnitudes near the domain wall. Due to the large variation of P 3 component across the domain wall, e 33 , e 13 and e 23 have much larger magnitudes near the domain wall than the other strain components. Long-range elastic interaction between the cylindrical domains in adjacent cells is also seen from the distributions of normal strain components. It is worth noting that the distributions of polarization and strain field break the cylindrical symmetry around the central axis of the initial written pattern, due to polarization relaxation mediated by the anisotropic elastic effect. From these results, it can be seen that our simulation can well capture two important features of real 180u domain walls, i.e., polarization rotation near the domain wall and domain wall anisotropy. It is known that these two features of domain wall play important roles in its transport property 29,30 . In the following, we will show that they have great impact on the stability behavior of the nanodomain.
We continue to investigate the stability of cylindrical domains with considering the effects of film thickness and temperature. To achieve this, we write cylindrical domains with different radii into 128 nm 3 128 nm 3 h simulation cells, which are free of external strain at different temperatures from T 5 0 K to T 5 350 K, with film thickness h ranging from 6 nm to 16 nm. The systems are relaxed to equilibrium to find whether the written domains can be maintained. The radius of the written cylindrical domains ranges from 1 nm to 64 nm. The simulated results are depicted in Figure 2c-f. For each investigated film thickness or temperature condition, it can be seen that there exists a critical size of cylindrical domain (denoted as r c ), below which it is instable, shrinks and disappears in the simulation cell. Above the critical size, cylindrical domains with various sizes can maintain stable, indicating that they are actually metastable states. Obviously, this result reflects the multi-valleys morphology of the system's free energy. It is the barriers between the energy valleys that lead to a pinning force against domain wall shrinkage. We believe that this domain wall pinning effect is not an artificial effect but an intrinsic effect of real systems with periodic lattice potential. While it is difficult to be taken into account in analytic models, this intrinsic domain wall pinning effect is important for domain stability in defect-free systems.
Importantly, the stability of cylindrical domain is sensitive to film thickness. At T 5 300 K, the critical size of cylindrical domain decreases from 22 nm to 6 nm when the film thickness increases from 6 nm to 16 nm (see Figure 2c and 2e). Considering the collective nature of domain, this result is easily understood as the domain has a larger volume in thicker film to withstand the devouring force from the surrounding domain. Apparently, such a dependence of domain stability on film thickness is against the well known Kittel's law 31 , which tells that the stable domain size increases as the film thickness increases. It is due to that our systems are under electrical short-circuit condition, thus the domain energy is mainly contributed by its bulk free energy rather than the depolarization energy. The latter is assumed to be dominant in the domain energy when deriving Kittel's law. A tradeoff between the increase of nano-  For 128 nm 3 128 nm 3 8 nm simulation cells under different temperatures (see Figure 2d and 2f), it is found that at low temperatures the written cylindrical domain is rather stable with its critical size as small as 2 nm at T 5 0 K. However, the critical size increases with the increase of temperature. At T 5 350 K, r c is found to be 48 nm. It is worth noting that there are strong nonlinear dependences of the cylindrical domain stability on both film thickness and temperature. Large margin of control on domain stability can be achieved at the high slope region of the r c -thickness or r c -temperature curves. On the other hand, if we were to make use of stable cylindrical domains, such nonlinear region should be kept away. Moreover, it should be pointed out that our predicted critical domain size is not necessary related with those obtained in SPM experiments [4][5][6][7][8][9] . This is largely due to the fact that the minimal domain size achieved in SPM experiments is limited by the radius curvature of SPM tip and is likely larger than our predicted critical domain size. Improved SPM technology and film grow method (to decrease the density of defects) may be necessary for inspection of the predicted ''intrinsic'' domain stability.
Equilibrium domain patterns under different external biaxial strains. After obtaining a basic insight into the stability of cylindrical domain as a function of film thickness and temperature, we focus on 8 nm-thick nanofilms at room temperature to explore the effects of mechanical loads. Mechanical loads effects on films with other thicknesses and/or at other temperatures can be qualitatively inferred. Before writing cylindrical domain patterns, in Figure 3, we first investigated the equilibrium domain patterns of a 128 nm 3 128 nm 3 8 nm simulation cell with initial random polarization distribution. The simulation cell is under different external biaxial strains, i.e., e a 11~e a 22~{ 0:01, 20.005, 0 and 0.005. It can be seen that the domain morphology of the nanofilm is significantly dependent on the mechanical load. Specifically, it shows that in the first three loading conditions, the nanofilm is in c-domain (i.e., P 1 5 P 2 5 0, P 3 ? 0) structure, whereas in the last loading condition it is in a/bdomain (i.e., P 1 ? 0, P 2 5 P 3 5 0 or P 2 ? 0, P 1 5 P 3 5 0) structure. This is due to that c-domain structure and a/b-domain structure can respectively accommodate the compressive and tensile biaxial strains, leading to a decrease of elastic energy, as has been experimentally verified from the domain abundance in epitaxial PTO thin films 33 . Importantly, comparing the first three domain patterns, we notice that c-domain size is sensitive to the external strain. For example, as indicated by the green arrows, small domains are stabilized by compressive strain (20.01). However, they disappear when the film is in a less compressive strain state (20.005). More interestingly, the formed domain patterns are also different when the nanofilm is subjected to different strain gradient (see Supplementary Figure S2 on line). These results clearly indicate the sensitivity of 180u domain stability in nanofilms to external mechanical loads. In the following, we will make a comprehensive exploration on the 180u cylindrical domains stability as a function of different external mechanical loads. In order to summarize the regularity of uniform strain effect on cylindrical domain more clearly, in Figure 4d and 4e, we draw the ''phase diagram'' depicting the equilibrium domain pattern of the simulation cell as a function of written domain size and external strain for the two strain conditions, respectively. It is important to see that the stability of cylindrical domain strongly depends on external strain. Similar with the effects of film thickness and temperature, such dependences are highly nonlinear, indicating a large margin of controllability. Specifically, for the investigated range of strain and temperature, critical size of written cylindrical domain ranges from 6 nm to a value larger than 64 nm (at e a 11~0 :002 for both uniaxial and biaxial strain conditions). Note that instability of 64 nm cylindrical domain is tested in an extended simulation cell of size 256 nm 3 256 nm 3 8 nm to avoid domain expansion caused by the cooperation between adjacent cells. Otherwise, for 128 nm 3 128 nm 3 8 nm simulation cell as shown in the phase diagrams, the 64 nm cylindrical domain would evolve into a single domain state due to the cooperation between adjacent cells at e a 11~0 :002 for both uniaxial and biaxial strain conditions. Furthermore, for the cases of largest strain, i.e., uniaxial strain (e a 11~0 :004, e a 22~0 ) and biaxial strain (e a 11~e a 22~0 :004), instead of remaining stable, the cylindrical domain would either backswitch or evolve into domain patterns with a/b-domains. Comparing the uniaxial and biaxial strain cases, it can be also seen that the abundance of a/b-domains in those patterns is much larger in biaxial strain condition.
In Figure 4f and 4g, we calculated the average polarization in zdirection over the initial cylindrical region of the written domain, denoted as ,P 3 .. The critical size of stable cylindrical domain and its dependence on external strain are clearly seen, as indicated by the transformation from instable state to stable state with a changing sign of ,P 3 .. It can be also inferred that those stable cylindrical domains probably have slight change in shape, particularly when the domain size is small, according to the round corner section right after transformation point of the ,P 3 . curves. Moreover, significant rumplings are observed in the ,P 3 . curves for the cases of largest strain, due to the fact that a/b-domains are induced and their morphologies are sensitive to the size of the initial cylindrical domain.
Note that Figure 4d and 4e present the equilibrium phase diagrams of cylindrical domain as a function of domain size and external strain. From these phase diagrams, we can readily infer the size range of stable cylindrical domain in an initially poled nanofilm for memory applications. Importantly, this size range can be effectively controlled by external strain. Through applying compressive external strain, the stability of cylindrical domain would be significantly enhanced and thus increases the storage density and data stability. Even more important indication from these phase diagrams is that we can erase the domain by changing the strain state. As indicated by the arrow in Figure 4d, if we write relatively small domain at compressive strain condition, by changing the strain condition to be less compressive, we can erase the domain. Interestingly, further simulation shows that the induced a/b-domains by uniaxial strain (e a 11~0 :004, e a 22~0 ) are instable if the strain is removed, leading to a single (poly) c-domain state when the initial cylindrical domain size is small (large). This indicates that mechanical erasing can be achieved even though a/b-domains were induced. Nevertheless, this seems only possible when the abundance of induced a/b-domains is relatively small. For those patterns with a large abundance of induced a/b-domains at biaxial strain condition (e a 11~e a 22~0 :004), most of them would relax to poly c-domain state or maintain stable, leading to an unsuccessful erasing. Typical evolution paths of domain patterns with induced a/b-domains after the strain load is off can be found in Figure S3 in Supplementary Information on line.
Control of domain stability by cylindrical bending. Besides uniform external strain, a ferroelectric nanofilm can be subjected to mechanical bending loads, which impose strain gradient and/or strain to the nanofilm. It is of practical significance to explore an effective mechanical way to erase nanodomains. The investigation of mechanical bending effects on the domain stability of ferroelectric nanofilm is also important for applications of flexible ferroelectric devices. Therefore, in the following we will inspect another two mechanical loading conditions, i.e., cylindrical bending condition and wavy bending condition. Note that the existence of strain gradient e jk,l may lead to a flexoelectric field, i.e., E f lexo i~f ijkl e jk,l , with f ijkl being the flexocoupling coefficients 25 . This flexoelectric field brings a mixed mechanical and electrical feature to the strain gradient effect. Currently, flexoelectricity is still a puzzled phenomenon in literature despite its growing relevance. Particularly, for flexocoupling coefficients, there are orders of magnitude in discrepancy between the theoretical estimations and experimental measurements 25 . Considering this uncertainty of the flexocoupling coefficients, and to see the strain gradient effect on domain stability more clearly, in the following we take values of the flexocoupling coefficients in between those of theoretical estimations and experi- For the cylindrical bending condition, we mean that the nanofilm is under an external bending strain in form of e a 11~e top z hz e bot 1{z=h ð Þ , where e top and e bot are the magnitudes of e a 11 at the top and bottom surfaces of the film, respectively, and coordinate z is measured with respect to the bottom surface of the film. This bending strain would cause a cylindrical-shape bending of the nanofilm along the x-direction as shown in Figure 5a. Typically, we consider two loading cases, i.e., the pure bending case with zero membrane strain, i.e., e top~{ e bot , and mixed bending-strain case with nonzero membrane strain, i.e., e top ? 0 and e bot 5 0. The distribution of e a 11 as a function of e top in the x-z plane of an 8 nm-thick nanofilm for the two loading cases is depicted in Figure 5b and 5c, respectively. Now we write cylindrical domains with different radii, i.e., ranging from 1 nm to 64 nm, into 128 nm 3 128 nm 3 8 nm simulation cells under different cylindrical bending strains at room temperature. The corresponding phase diagrams depicting the equilibrium domain pattern as a function of domain size and bending strain are shown in Figure 5d and 5e, for the pure bending case and mixed bending-strain case respectively. Note that the flexoelectric field has been neglected to see the pure mechanical effect. From the phase diagrams, it can be seen that cylindrical bending can also effectively control the stability of cylindrical domain. Firstly, it shows that strain gradient decreases the stability of cylindrical domain, similar with the effect of tensile strain. Nevertheless, strain gradient can largely depress the formation of a/b-domains. Moreover, dependence of the domain stability on cylindrical bending also exhibits a strong nonlinear feature. For the pure bending case (Figure 5d), it shows that the critical size of stable cylindrical domain remains almost the same when the bending strain is small, i.e., e top 5 2e bot , 0.004. However, as the bending strain further increases, the critical size of stable cylindrical domain rapidly increases to ,30 nm at e top 5 2e bot 5 0.006, and exceeds 64 nm at e top 5 2e bot 5 0.008 (To see the instability of 64 nm cylindrical domain, the simulation cell size is extended to be 256 nm 3 256 nm 3 8 nm). This nonlinear feature is even more significant in the mixed bending-strain condition (Figure 5e), due to the superposition of strain gradient effect and membrane strain effect. For this loading case, the membrane strain is compressive when e top , 0, which cancels the effect of strain gradient, leading to small change of domain stability. Meanwhile, the membrane strain becomes tensile when e top . 0, which has similar effect of strain gradient, and leads to large change of domain stability. This indicates a larger margin of controllability on the cylindrical domain stability of this type of mechanical loading.
In Figure 5f and 5g, corresponding to the two cylindrical bending cases, we also calculated the average polarization in z-direction over the initial cylindrical region of the written domain. The critical size of stable cylindrical domain and its dependence on cylindrical bending strain are clearly seen, as indicated by the transformation from instable state to stable state with a changing sign of ,P 3 .. It should be also noted that the ,P 3 . of single domain state almost remains the same despite the large change of bending condition. Moreover, compared with the cases of uniform strain loading conditions as shown in Figure 4, the stable cylindrical domains change their shape more significantly in cylindrical bending conditions, as indicated by the round corner section right after transformation point of the ,P 3 . curves.
Control of domain stability by wavy bending. The investigated mechanical loading conditions so far are with uniform strain or uniform strain gradient. It has been recently demonstrated that a wavy bending configuration of ferroelectric nanofilm can be achieved by integrating it on prestrained elastomeric supports, which is important for developing flexible ferroelectric devices 15,16 Considering the fact that the strain gradient of wavy bending is not uniform, an exploration on its effect is not trivial to the cases of uniform strain and uniform strain gradient. For a nanofilm with a wavy bending along x-direction (Figure 6a), its strain state can be approximately described by a cosine form strain, i.e., characterizing the largest defection of the nanofilm from its neutral plane z 5 h/2, and l is the wave length along the xdirection. Note that for simplicity we fix the neutral plane at middle plane of the film to only consider the case of zero membrane strain. To see how wavy bending affects the stability of 180u cylindrical domain, simulations are conducted on 128 nm 3 128 nm 3 8 nm simulation cells under different wavy bending strains (l 5 128 nm) at room temperature. The simulation cells are initially written into cylindrical domains, with their radii ranging from 1 nm to 64 nm. Figure 6b depicts the distributions of e a 11 as a function of A b in the x-z plane of the simulation cell under a wavy bending with l 5 128 nm. Here we also investigate the effect of wavy bending in the presence of a flexoelectric field, E f lexo 3~f 12 Le 11 =Lz, with f 12 5 10 V. The distribution of E f lexo 3 as a function of A b in the x-z plane of the simulation cell is depicted in Figure 6c. Thus in Figure 6d and 6e, we respectively calculate the phase diagram of equilibrium domain patterns in the simulation cells for two wavy bending cases, i.e., with zero and nonzero flexoelectric field, respectively. The average polarization in zdirection, i.e., ,P 3 ., over the initial cylindrical region of the written domain, are also calculated and depicted in Figure 6f and 6g, corresponding to the two bending cases.
Firstly, it can be seen that wavy bending has significant impact on the stability of cylindrical domains. As shown in Figure 6d, where the flexoelectric field is switched off, the critical size of stable cylindrical domain changes from 11 nm to 30 nm when the bending amplitude A b increases from 0.2 nm to 1.0 nm. Note that although the maximum strain of a wavy bending with A b 5 10 nm is about 0.01 according to e max~2 p 2 A b l 2 h, the controllability of investigated wavy bending is less significant than the uniform in-plane stains ( Figure 4) and cylindrical bending ( Figure 5). This is due to the fact that, for a large size of domain, its stability is more likely to be affected by the overall strain and strain gradient state. Therefore, if the flexoelectric field is negligible, for the margin of controllability on cylindrical domain stability, we have uniform in-plane strain . cylindrical bending (with zero membrane strain) . wavy bending.
Nevertheless, significant difference can be found in the controllability on domain stability of the bending loads, if a flexoelectric effect is not trivial. Specifically, flexoelectric effect imposes an electric field proportional to the strain gradient, and changes sign as direction of the strain gradient changes. As a consequence, both enhancement and depression of the stability of the 180u cylindrical domain are possible, dependent on the direction of the domain polarization and the flexoelectric field. This feature is clearly seen in Figure 6e, where a flexoelectric field is included. Specifically, the critical size of stable cylindrical domain is tuned from 3 nm to 48 nm, dependent on the bending direction. Note that upper limit of 48 nm exists in the upper bending case due to the cosine form of the strain gradient. A much larger critical domain size can be reached by increasing the wave length or in case of cylindrical bending. Note also that the maximum flexoelectric field is about 2.5 3 10 7 V/m in our investigation (Figure 6b). In normal ferroelectrics, polarization switching occurs at much lower field strengths than that predicted by Landau-type phenomenology. The flexoelectric field needed to significantly affect the stability of cylindrical domains might be much less than our prediction. Due to the modulated flexoelectric field, large change in domain shape not only happens at region near domain instability but also happens when the domain size is comparable to the wave length, as indicated by the ,P 3 . curves as shown in Figure 6g.
Mechanical erasing of the information in nano-ferroelectric memories. Finally, we would like to conduct simulations on multibits memory system with 4 3 4 memory units at room temperature, to clearly show how to erase the information by exerting or removing mechanical loads. In Figure 7a, a 256 nm 3 256 nm 3 8 nm memory system at free of external strain condition is initially written with 4 3 4 bits of information, with the size of written cylindrical domains being 16 nm. The information is found stable at such a strain state, in consistent with the phase diagrams in Figures 3-6. Then we apply external strain to the memory system. Figure 7b shows that evolution of the information as the film is under uniform biaxial strain, e a 11~e a 22~0 :002. It can be seen that the information is totally erased, due to the instability of 16 nm cylindrical domain at this loading condition. The erasing of the information can be also done by applying cylindrical bending (e top 5 2e bot 5 0.006) to the memory system ( Figure 7c). Interestingly, making use of the flexoelectric effect, we can also partially erase the information by applying a wavy bending (A b 5 1.0 nm and l 5 128 nm) to the system (Figure 7d).
Moreover, the remaining information can be also erased by changing the phase of bending wavy, e.g., by p. Another mechanical erasing example exploiting the flexoelectric effect is illustrated in Figure 7e and 7f. In Figure 7e, it shows that a 256 nm 3 256 nm 3 6 nm memory system under wavy bending (A b 5 0.2 nm and l 5 128 nm) can maintain 16 nm cylindrical domains at the valleys, due to enhance of domain stability by the flexoelectric effect. After removing the wavy bending, as 16 nm cylindrical domain is no longer stable, the information is successfully erased (Figure 7f). These results clearly demonstrate the possibility of erasing the information by mechanical loads.

Discussion
We have conducted a systematical investigation to explore the effects of macroscopic mechanical loads on 180u cylindrical domain patterns in ferroelectric ultrathin films, and the possibility of mechanical erasing the ferroelectric domains. Dependence of cylindrical domain stability on domain size, film thickness and temperature has been revealed. Subsequently, control on the cylindrical domain stability of different mechanical loads, including uniform strain, cylindrical bending, and wavy bending, has been comprehensively investigated and compared. Our results show that both strain and strain gradient have great impact on the stability of 180u cylindrical domain, with its critical size sensitive to strain state of the nanofilm. Particularly, we found that if the flexoelectric field is neglected, the cylindrical domain stability is most sensitive to uniform in-plane strain, and least sensitive to wavy bending. Meanwhile, if a flexoelectric effect is not trivial, significant difference can be found in the controllability on domain stability of the bending loads. All these results indicate erasing 180u cylindrical domain can be achieved by changing the strain state of the nanofilm. According to the calculated phase diagrams depicting the cylindrical domain stability as a function of mechanical loads and domain size, we successfully simulated possible mechanical erasing processes of a memory device storing 4 3 4 bits information. We believe our finding is instructive to prospective device applications of ferroelectric ultrathin films, such as flexible memory devices and other micro-electromechanical systems based on nanodomains.
Here we would like to make a further discussion about the indications of large controllability of mechanical loads on cylindrical domain stability. Besides utilizing this controllability to increase storage density of cylindrical domains and to realize mechanical erasing of the nanodomains, some interesting points can be further drawn. Obviously, the coercive field of domain switch, and thus the operating voltage and speed, should be also strong functions of mechanical loads. In general, there is always a tradeoff in the pursuit of high density, low coercive field/operating voltage, and high operating speed. By exploring the dependence of related factors on mechanical loads, an optimized point at which the memory device is of high-density, high-speed and low-power consumption can be reached. An even more clever strategy is to maintain the nanodomain of memory device at a mechanical strain state and to process the domain at another mechanical strain state. Moreover, for memory devices under the wavy bending, our simulation indicates that the domain stability is spatial dependent. Making use of this property, we can realize addressable erasing of part of the information if we can precisely control the wavy pattern (e.g., wave length and wave phase). Moreover, information can be stored at different regions according to their operating condition (for example, information that needs frequent operation can be stored at less stable region). This would be of practical importance for applications. All of these, although there is still a long path to go through, indeed present us wonderful possibility in memory applications involving mechanical loads.

Methods
Typically, the PTO nanofilms are modeled as infinite along the in-plane x-and ydirections. We mainly restrict ourselves to simulation cells containing one inverted cylindrical domain at the center region (Figure 1), and apply periodic conditions in x-and y-directions. Simulation cells with 4 3 4 memory units to mimic multi-bits memory devices will be also investigated to show how to erase the information by adjusting mechanical loads. The nanofilms are assumed to be under in-plane strain constraint with free top-bottom surfaces and under electrical short-circuit condition 18 . Periodic conditions are applied in the in-plane directions. Meshing grids of n x Dl 3 n y Dl 3 n z Dl are used to simulate different simulation cells, with scale Dl equal to 1 nm.
The domain structure is solved by a generally accepted and powerful theory for simulating domain structures of ferroelectrics, i.e., the phase field model based on the Landau-type formalism, in which we also carefully takes into account the effects of inhomogeneous mechanical field, electric filed and surfaces. Such approach has been widely used to simulate ferroelectric domain structure, and represents well the experimental results of ferroelectric domain structures in a wide range of materials and conditions [33][34][35][36] . Specifically, in our phase field model, evolution of the polarization field is described by the Time-Dependent Ginzburg-Landau (TDGL) equations, i.e., hP i /ht 5 2MdF/dP i , where F is the free energy of the system, P i are the polarization components, t is time, and M is the kinetic coefficient related to the domain wall mobility. The evolution of elastic and electric fields are treated in an ''adiabatic'' way. At each time step of polarization evolution, these fields are calculated according to the mechanical and electrostatic equilibrium equations, using a fast Fourier transformation (FFT) technique based on Khachaturyan's microscopic elasticity theory 37 and Stroh formulism of anisotropic elasticity 38 . In literature, based on this approach, the domain structures in epitaxial ferroelectric thin films have been successfully simulated 17,18 . To guarantee the convergence of domain structure to equilibrium, the simulation time is set to be sufficient long (up to 10 6 time steps). The material parameters in calculations are listed Table 1. For PTO, a commonly used sixorder Landau-potential is adopted in this study. The potential coefficients were determined by fitting the experimental results thus the potential captures well the phase transition behavior of PTO 39