Anna N Morozovska
email: anna.n.morozovska@gmail.com

Riccardo Hertel
email: riccardo.hertel@ipcms.unistra.fr

Salia Cherifi-Hertel

Victor Yu Reshetnyak

Eugene A Eliseev

Dean R Evans
email: dean.evans@afresearchlab.com

Chiral Polarization Textures Induced by the Flexoelectric Effect in Ferroelectric Nanocylinders

Polar chiral structures have recently attracted much interest within the scientific community, as they pave the way towards innovative device concepts similar to the developments achieved in nanomagnetism.

Despite the growing interest, many fundamental questions related to the mechanisms controlling the appearance and stability of ferroelectric topological structures remain open. In this context, ferroelectric nanoparticles provide a flexible playground for such investigations. Here, we present a theoretical study of ferroelectric polar textures in a cylindrical core-shell nanoparticle. The calculations reveal a chiral polarization structure containing two oppositely oriented diffuse axial domains located near the cylinder ends, separated by a region with a zero-axial polarization. We name this polarization configuration "flexon" to underline the flexoelectric nature of its axial polarization. Analytical calculations and numerical simulation results show that the flexon's chirality can be switched by reversing the sign of the flexoelectric coefficient. Furthermore, the anisotropy of the flexoelectric coupling is found to critically influence the polarization texture and domain morphology. The flexon rounded shape, combined with its distinct chiral properties and the localization nature near the surface, are reminiscent of Chiral Bobber structures in magnetism. In the azimuthal plane, the flexon displays the polarization state of a vortex with an axially polarized core region, i.e., a meron. The flexoelectric effect, which couples the electric

polarization and elastic strain gradients, plays a determining role in the stabilization of these chiral states.

We discuss similarities between this interaction and the recently predicted ferroelectric Dyzaloshinskii-Moriya interaction leading to chiral polarization states.

I. INTRODUCTION

Research on ferroelectric materials has received growing interest over the past years, driven in part by the potential of these material systems for low-power technological applications in a broad spectrum of domains [1,2], ranging from high-density data storage to optical nanodevices. A central aspect of this field of research is the formation of ferroelectric domain structures [3], and more generally the micro-and nanoscale structure of the polarization field [4].

Traditionally, research on ferroelectrics is centered on the study of bulk materials and thin films [5,6,7], but recently ferroelectric nanoparticles have also attracted increasing interest [8,9,10,11,12,13,14]. In ferroelectric thin films and nanoparticles, the polarization structure is strongly affected by electrostatic (depolarizing) fields [15,16,17,18], as well as by strain and strain gradients [19,20,21,22] via the flexoelectric effect [START_REF] Megaw | Temperature changes in the crystal structure of barium titanium oxide[END_REF]24,25,26].

Although the foundations for the theoretical description of ferroelectrics have been established decades ago [27], understanding the complex physical properties of these material systems remains a challenge for fundamental research. Recent progress in this field, achieved to a large extent through advanced imaging techniques [28] and by employing modern numerical simulations [29], includes the discovery of highly complex polarization structures, such as flux closure [5,30] and bubble domains [31], meandering [32,[START_REF] Han | Mapping gradient-driven morphological phase transition at the conductive domain walls of strained multiferroic films[END_REF] and/or labyrinthine [11,[START_REF] Nahas | Inverse transition of labyrinthine domain patterns in ferroelectric thin films[END_REF] structures, non-Ising type chiral domain walls [START_REF] Cherifi-Hertel | Non-ising and chiral ferroelectric domain walls revealed by nonlinear optical microscopy[END_REF], polarization vortices in thin layers [START_REF] Rodriguez | Vortex polarization states in nanoscale ferroelectric arrays[END_REF][START_REF] Balke | Enhanced electric conductivity at ferroelectric vortex cores in BiFeO3[END_REF][START_REF] Yadav | Observation of polar vortices in oxide superlattices[END_REF], nanodots [START_REF] Liu | An isogeometric approach to flexoelectric effect in ferroelectric materials[END_REF] or nanopillars [START_REF] Li | Superelastic oxide micropillars enabled by surface tension-modulated 90° domain switching with excellent fatigue resistance[END_REF], or polar skyrmions [START_REF] Das | Observation of room-temperature polar skyrmions[END_REF][START_REF] Hlinka | Skyrmions in ferroelectric materials[END_REF].

While skyrmions and other chiral structures have dominated the past decade of research in magnetism [START_REF] Finocchio | Magnetic skyrmions: from fundamental to applications[END_REF], these topological states have received less attention by the ferroelectric community. Only recently a strong interest has emerged in chiral polarization structures, which can be attributed to the observation of skyrmion states in ferroelectrics [START_REF] Das | Observation of room-temperature polar skyrmions[END_REF][START_REF] Hlinka | Skyrmions in ferroelectric materials[END_REF]. However, the theoretical understanding of these structures is not as advanced as it is in the case of their magnetic counterparts, and the mechanism that underpins the formation of skyrmions in ferroelectrics is not fully understood. The fundamental interaction stabilizing the magnetic version of these structures in chiral ferromagnets [START_REF] Dzyaloshinsky | A thermodynamic theory of "weak" ferromagnetism of antiferromagnetics[END_REF][START_REF] Moriya | Anisotropic superexchange interaction and weak ferromagnetism[END_REF] is the Dzyaloshinky-Moriya Interaction (DMI). The DMI favors the formation of helical structures with a well-defined handedness as they occur, e.g., along the radial direction of skyrmions. As scientists working on ferroelectrics hope to replicate the success that chiral structures have witnessed in magnetism, the possibility of a "ferroelectric DMI" has recently been discussed [START_REF] Zhao | Dzyaloshinskii-Moriyalike interaction in ferroelectrics and antiferroelectrics[END_REF]. However, Erb and Hlinka [START_REF] Erb | Vector, bidirector, and Bloch skyrmion phases induced by structural crystallographic symmetry breaking[END_REF] showed that only very few exotic ferroelectrics could theoretically sustain an intrinsic DMI-type interaction since it requires particular symmetry properties of the crystal lattice. Here we discuss the flexoelectric coupling as an alternative mechanism that can generate chiral polarization states in ferroelectrics.

The thermodynamic description of the flexoelectric effect is given by the Lifshitz invariant in the free energy expansion [22]. It is known that, in magnetic materials, the occurrence of similar Lifshitz invariants converts directly into an antisymmetric coupling known as the DMI [START_REF] Bogdanov | Thermodynamically stable "vortices" in magnetically ordered crystals, The mixed state of magnets[END_REF][START_REF] Bogdanov | On the stability of localized states in nonlinear field models with Lifshitz invariants[END_REF], which favors the formation of helicoidal structures with a specific chirality.

The existence of a ferroelectric counterpart of the DMI was recently predicted by first-principles simulations [START_REF] Zhao | Dzyaloshinskii-Moriyalike interaction in ferroelectrics and antiferroelectrics[END_REF]. The ferroelectric analogue of the DMI was discussed in the context of Lifshitz invariants by Strukov and Levanyuk [START_REF] Strukov | Ferroelectric Phenomena in Crystals[END_REF], and more recently by Erb and Hlinka [START_REF] Erb | Vector, bidirector, and Bloch skyrmion phases induced by structural crystallographic symmetry breaking[END_REF], who argued that a ferroelectric DMI can exist. In addition to the remarkable similarity in the mathematical form of the flexoelectric Lifshitz invariant and DMI, the flexoelectric term appears to have a similar impact as the DMI in terms of the formation of chiral structures.

By means of the finite element modeling (FEM) based on the Landau-Ginzburg-Devonshire (LGD) theory, this paper shows that an anisotropic flexoelectric effect can give rise to a previously unexplored type of polarization state with distinct chiral properties. Remarkably, these homochiral properties are not induced by a DMI term. This finding suggests that the recently discussed DMI in ferroelectrics is not the only possible mechanism for the formation of homochiral polarization states, and that anisotropic flexoelectric effects offer an alternative pathway to stabilize such structures in ferroelectric nanostructures. We discuss common aspects of the DMI and the flexoelectric effect, which are both derived from Lifshitz invariants in the framework of the Landau theory of second-order phase transitions [22].

II. CONSIDERED PROBLEM AND MATERIAL PARAMETERS

Using a LGD phenomenological approach along with electrostatic equations and elasticity theory, we model the polarization, the internal electric field, and the elastic stresses and strains in a core-shell nanoparticle using FEM, where the ferroelectric core is made of BaTiO3 and has a cylindrical shape. The aspect ratio of the nanocylinder radius 𝑅 to its length ℎ is significantly higher than unity. The z-axis is parallel to the cylinder axis (Fig. 1). The shell is an elastically soft paraelectric or high-k semiconductor with a thickness ∆𝑅 ≪ 𝑅 and screening length Λ ≥ 1 nm. The coverage can be artificial (e.g., a soft organic semiconductor or vacancyenriched SrTiO3) or natural, and in the latter case it would originate from the polarization screening by surrounding media. The core-shell nanoparticle is placed in a very soft elastic medium. The LGD free energy functional G of the nanoparticle core includes a Landau energyan expansion on powers of 2-4-6 of the polarization (𝑃 𝑖 ), 𝐺 𝐿𝑎𝑛𝑑𝑎𝑢 ; a polarization gradient energy, 𝐺 𝑔𝑟𝑎𝑑 ; an electrostatic energy, 𝐺 𝑒𝑙 ; an elastic, electrostriction contribution 𝐺 𝑒𝑠 , a flexoelectric contribution, 𝐺 𝑓𝑙𝑒𝑥𝑜 ; and a surface energy, 𝐺 𝑆 . It has the form [START_REF] Eliseev | Flexo-Elastic Control Factors of Domain Morphology in Core-Shell Ferroelectric Nanoparticles: Soft and Rigid Shells[END_REF]:

𝐺 = 𝐺 𝐿𝑎𝑛𝑑𝑎𝑢 + 𝐺 𝑔𝑟𝑎𝑑 + 𝐺 𝑒𝑙 + 𝐺 𝑓𝑙𝑒𝑥𝑜 + 𝐺 𝑓𝑙𝑒𝑥𝑜 + 𝐺 𝑆 , (1a)
𝐺 𝐿𝑎𝑛𝑑𝑎𝑢 = ∫ 𝑑 3

𝐺 𝑒𝑠 = -∫ 𝑑 3 𝑟 𝑉𝑐 ( 𝑠 𝑖𝑗𝑘𝑙 2 𝜎 𝑖𝑗 𝜎 𝑘𝑙 + 𝑄 𝑖𝑗𝑘𝑙 𝜎 𝑖𝑗 𝑃 𝑘 𝑃 𝑙 ), (1d)
𝐺 𝑓𝑙𝑒𝑥𝑜 = -∫ 𝑑 3 𝑟 𝑉𝑐 𝐹 𝑖𝑗𝑘𝑙 2 (𝜎 𝑖𝑗 𝜕𝑃 𝑘 𝜕𝑥 𝑙 -𝑃 𝑘 𝜕𝜎 𝑖𝑗 𝜕𝑥 𝑙 ), (1e)
𝐺 𝑆 = 1 2 ∫ 𝑑 2 𝑟 𝑆 𝑎 𝑖𝑗 (𝑆) 𝑃 𝑖 𝑃 𝑗 . (1f)
Here 𝑉 𝐶 is the core volume. The coefficient 𝑎 𝑖 linearly depends on temperature T, 𝑎 𝑖 (𝑇) = 𝛼 𝑇 [𝑇 -𝑇 𝐶 ], where 𝛼 𝑇 is the inverse Curie-Weiss constant and 𝑇 𝐶 is the ferroelectric Curie temperature renormalized by surface tension/intrinsic surface stresses [START_REF] Ma | Surface tension and Curie temperature in ferroelectric nanowires and nanodots[END_REF][START_REF] Wang | Strain effect on phase transitions of BaTiO3 nanowires[END_REF][START_REF] Wang | Phase transitions and domain structures of ferroelectric nanoparticles: Phase field model incorporating strong elastic and dielectric inhomogeneity[END_REF] and surface

h shell soft medium z P3 P1,2 core 2R
bond contraction [START_REF] Huang | Surface bond contraction and its effect on the nanometric sized lead zirconate titanate[END_REF][START_REF] Huang | Grain-size effect on ferroelectric Pb(Zr1-xTix)O3 solid solutions induced by surface bond contraction[END_REF]. Tensor components 𝑎 𝑖𝑗 are regarded as temperature-independent. The tensor 𝑎 𝑖𝑗 is positively defined if the ferroelectric material undergoes a second order transition to the paraelectric phase and negative otherwise. The higher nonlinear tensor 𝑎 𝑖𝑗𝑘 and the gradient coefficients tensor 𝑔 𝑖𝑗𝑘𝑙 are positively defined and regarded as temperatureindependent. In Eq.(1e), 𝜎 𝑖𝑗 is the stress tensor, 𝑠 𝑖𝑗𝑘𝑙 is the elastic compliances tensor, and 𝑄 𝑖𝑗𝑘𝑙 is the electrostriction tensor. In the Lifshitz invariant, Eq.(1f), 𝐹 𝑖𝑗𝑘𝑙 is the flexoelectric tensor.

Landau-Khalatnikov equations [START_REF] Landau | On the anomalous absorption of sound near a second order phase transition point[END_REF][START_REF] Hlinka | Mobility of Ferroelectric Domain Walls in Barium Titanate[END_REF] obtained from a variation of the free energy (1), mathematical formulation of the electrostatic and elastic sub-problem (see e.g. [START_REF] Timoshenko | Theory of Elasticity[END_REF]), initial and boundary conditions (see e.g. [START_REF] Tagantsev | Short-range and long-range contributions to the size effect in metal-ferroelectric-metal heterostructures[END_REF][START_REF] Stengel | Origin of the dielectric dead layer in nanoscale capacitors[END_REF][START_REF] Stengel | First-principles modeling of ferroelectric capacitors via constrained displacement field calculations[END_REF]), sensitivity to the shape of the cylinder ends, polarization gradient coefficients, shell dielectric permittivity and semiconducting properties, and other details of FEM are given in Appendix A of Suppl. Mat. [63]. The ferroelectric, dielectric, and elastic properties of the BaTiO3 core are collected from Refs. [START_REF] Hlinka | Phenomenological model of 90-degree domain wall in BaTiO3 type ferroelectrics[END_REF][START_REF] Bell | Phenomenologically derived electric field-temperature phase diagrams and piezoelectric coefficients for single crystal barium titanate under fields along different axes[END_REF][START_REF] Berlincourt | Elastic and piezoelectric coefficients of single-crystal barium titanate[END_REF][START_REF] Ponomareva | Finite-temperature flexoelectricity in ferroelectric thin films from first principles[END_REF][START_REF] Shirane | Transition energy and volume change at three transitions in barium titanate[END_REF] and 69] and given in Table SI.

III. RESULTS OF FINITE ELEMENT MODELING

A. FEM Results at Room Temperature

Images in Figs. 2a and 3a are calculated without electrostriction (𝑄 𝑖𝑗 = 0) and flexoelectric (𝐹 𝑖𝑗 = 0) couplings between the electric polarization and elastic stresses. For the case a very prolate dipolar kernel oriented along z-axis appears inside the cylindrical core. The kernel has relatively thin 180-degree domain walls, which are mostly uncharged because they are parallel to the kernel axis and cylinder lateral surface. The bound charges appear at the walls only in a small spatial region near the kernel that is contact with the cylinder ends, where the 180-degree walls become counter head-to-head walls. The axial polarization 𝑃 3 inside the kernel is high, 𝑃 3 ~-(20 -25) C/cm 2 (this is very close to the bulk polarization of BaTiO3 ~26 C/cm 2 ), and the surrounding core has relatively small axial polarization of the opposite sign, 𝑃 3 ~(0 -5) C/cm 2 . The lateral components of polarization, 𝑃 1 and 𝑃 2 , form a twodimensional (2D) vortex without a central empty core, because a dipolar kernel evolves instead (Fig. S4 [63]). The two symmetrical Bloch points with 𝑷 = 0 are located at the junction of the dipolar kernel with the cylinder ends. The "up" or "down" orientation of polarization component SI. The extremal (maximal or minimal) value 𝑃 𝑒 in the center of the diffuse axial P3-domain is an odd function of 𝑓, which is zero at 𝐹 𝑖𝑗 = 0 (the green curve in Fig. 3f). Note that the 𝑃 𝑒 value frequently differs from 𝑃 𝑚𝑎𝑥 and 𝑃 𝑚𝑖𝑛 values due to the presence of the small sixteen 𝑃 3 -domains localized near the top and bottom junction of the sidewall with the cylinder ends (bottom row in Fig. S5 [63]).

For the remainder of the paper, we refer to the localized polarization structure near the wire ends as a "flexon" for the sake of brevity and to underline the flexoelectric nature of its for the integration over the cylinder cross-section {𝑥, 𝑦}. For the case of 𝑃 3 (𝜌 = 𝑅, 𝑧) → 0, z-dependence of the topological index is

𝑛(𝑧) = - 𝑃 3 (𝜌=0,𝑧) 2𝑃(𝜌=0,𝑧) ≅ - sign[𝑓]𝑧 2√1+(𝑧 2 𝐵 ⁄ ) . (3)
Here sign[𝑓] is the sign of the flexoelectric coefficients 𝐹 𝑖𝑗 , 𝐵 is a positive constant, which depends on the absolute value of |𝐹 𝑖𝑗 |. 𝑛(𝑧) is a normalized profile of 𝑃 3 (𝜌 = 0, 𝑧), and so 𝑛(𝑧) = 0 for 𝐹 𝑖𝑗 = 0, and its sign is defined by the sign of 𝐹 𝑖𝑗 . The dependence 𝑛(𝑧) is shown in Fig. 3g and Fig. S10a for zero, positive, and negative 𝐹 𝑖𝑗 . Since the value 𝑃(0, 𝑧) is very close to the 𝑃 3 (0, 𝑧) near the cylinder ends (Fig. S10a and S10b), and 𝑃 3 (0, 𝑧) vanishes in the central part of a nanoparticle, the topological index continuously changes from -½ to +½ with a zcoordinate change from one cylinder end to the other. The result clearly shows the localization of the chiral structuresthe flexonsat the ends of the wires. The topological index, which can be interpreted as the degree to which a structure is chiral, changes sign from one end to the other, and changes sign upon reversal of the sign of 𝐹 𝑖𝑗 . It also increases in magnitude with increasing absolute value of |𝐹 𝑖𝑗 |. These properties are evidence of an obvious correlation between the flexoelectric effect and the formation of chiral polarization structures.

The revealed type of isolated chiral polarization structures, i.e., flexons, display topological features of a three-dimensional meron. In this sense, the polarization vortex in the XY-plane can be interpreted as the Bloch-like transition region of a meron connecting polarization directions of opposite P3 sign in the core region and in the outer cylindrical shell (Fig. 2). The flexon polarization 𝑃 ⃗⃗ develops a characteristic drop-shape with a chiral structure localized near the surfaces of the cylinder that is reminiscent of the chiral-bobber state found in non-centrosymmetric magnetic films [START_REF] Rybakov | New type of stable particle-like states in chiral magnets[END_REF] and nanoparticles [START_REF] Pathak | Three-dimensional chiral magnetization structures in FeGe nanospheres[END_REF]. It is worth noting that similar, skyrmion-like configurations at the ends of cylindrical nanowires have also been predicted

analytically [START_REF] Braun | Nucleation in ferromagnetic nanowires-magnetostatics and topology[END_REF] and numerically [START_REF] Hertel | Magnetic drops in a soft-magnetic cylinder[END_REF] in the case of non-chiral ferromagnetic materials, but only in the form of transient configurations appearing during the dynamic magnetization reversal process. Here, the skyrmion-like polarization structures appear as stable states in the ferroelectrics, owing to a chiral-symmetry breaking effect of the flexoelectric coupling. In contrast to previous findings [START_REF] Yudin | Bichiral structure of ferroelectric domain walls driven by flexoelectricity[END_REF][START_REF] Eliseev | Structural phase transitions and electronic phenomena at 180-degree domain walls in rhombohedral BaTiO3[END_REF][START_REF] Gu | Non-Ising Character of a Ferroelectric Wall Arises from a Flexoelectric Effect[END_REF], the flexon structure is chiral [START_REF]different from its mirror image. Systems governed by chiral-symmetry breaking interactions, such as the DMI, favor the formation of chiral structures with one specific handedness. Those systems are sometimes called homochiral[END_REF] and almost uncharged because 𝑑𝑖𝑣𝑃 ⃗ ⃗ ≅ 0 (Fig. S9b [63]).

B. Temperature Behavior of the Flexon-Type Polarization Distribution

To define the temperature interval in which flexons exist as stable or meta-stable states, we performed FEM in the temperature range from 50 K to 400 K using different initial The axial counter domains inherent to flexons are the most pronounced feature over the narrower range 290 K < 𝑇 < 340 K.

IV. DISCUSSION

Any deviation from a four quadrants domain configuration in the flexon-type polarization distribution is found to be metastable. This is because the antiparallel bidomain-type polarization distribution (starting from a random noise) has a lower free energy in a BaTiO3 cylindrical nanoparticle. The derived energy values at room temperature are 𝐺 𝑓𝑙 = -3.6 • 10 -18 J and 𝐺 𝑏𝑑 = -4.0 • 10 -18 J in the flexon and the bidomain structure, respectively. The energy difference between these states, ∆𝐺 = 4 • 10 -19 J is much higher (about 100 𝑘 𝐵 𝑇) than the thermal energy barrier 𝑘 𝐵 𝑇 at room temperature. However, the ratio ∆𝐺 𝑘 𝐵 𝑇 strongly decreases with as the temperature increases. The linear relative dielectric permittivity in both states is about 110 at room temperature and it strongly increases with temperature. Furthermore, our numerical simulations show that a spontaneous off-field transition from the flexon to the bidomain polarization state does not occur, whereas the in-field transition is possible (corresponding hysteresis loops are shown in Appendix D [63]). Thus, the bidomain and flexon states of a prolate core-shell ferroelectric nanoparticle can be considered as the exited and ground state of a two-level system suitable for information recording. The two-level system can imitate qubits operating in the temperature range where 1 < ∆𝐺 𝑘 𝐵 𝑇 < 5. Furthermore, the bidomain polarization state corresponds to an antiferroelectric-like state of the nanoparticle polarization, which can be represented as two antiparallel nanoscopic dipoles. The flexon is a much more complex achiral vortex-like configuration containing two counter dipole nanodomains with diffuse relaxor-like polar properties. Thus, an ensemble of prolate core-shell ferroelectric nanoparticles, where a given nanoparticle is either in a flexon or a bidomain state, can be an alternative media for information processing. The media may exhibit unusual properties including antiferroelectric and/or relaxor-like polarization states, which can lead to additional functionalities. Note that the appearance of the antiferroelectric and relaxor-like glass states, as well as a newly discovered liquid glass state [START_REF] Roller | Observation of liquid glass in suspensions of ellipsoidal colloids[END_REF] with additional (anti)ferroelectric ordering and other cross-talk effects, are possible in a suspension of the prolate core-shell ferroelectric nanoparticles.

The relatively wide temperature range (about 100 K) corresponding to the stability or meta-stability of the flexon-type polarization distribution gives us the hope that the domain morphology can be observed experimentally. Specifically, the measurements of local vertical displacement by piezoresponse force microscopy (PFM) visualize the distribution of 𝑃 3 (𝑟 ⃗) at distances ~ 10 nm from the ends of a nanoparticle, but the resolution procedure for the local piezoresponse of diffuse domains under the surface is far not straightforward [START_REF] Kalinin | Nanoelectromechanics of piezoresponse force microscopy[END_REF][START_REF] Scrymgeour | Nanoscale piezoelectric response across a single antiparallel ferroelectric domain wall[END_REF]. This is because PFM is a near-field method. A complementary tool to probe chiral polar textures is farfield nonlinear optical microscopy [START_REF] Cherifi-Hertel | Shedding light on non-Ising polar domain walls: Insight from second harmonic generation microscopy and polarimetry analysis[END_REF], which has a comparatively much lower resolution than PFM, yet optimum focusing methods and the experimental geometry allow for overcoming the diffraction limit. For example, second-harmonic generation microscopy was successfully used by the community to precisely study semiconductor nanowires [START_REF] Bautista | Secondharmonic generation imaging of semiconductor nanowires with focused vector beams[END_REF]. This method should also be capable of providing complementary information on the 3D ferroelectric domain structure (being sensitive to 𝑃 1,2 and 𝑃 3 ) by using polarimetry analysis (see, e.g., [START_REF] Timofeeva | Polar second-harmonic imaging to resolve pure and mixed crystal phases along GaAs nanowires[END_REF]). Another promising method is resonant elastic soft X-ray scattering, a synchrotron-based method sensitive to chiral polar arrangements through dichroism effects [START_REF] Shafer | Emergent chirality in the electric polarization texture of titanate superlattices[END_REF][START_REF] Lovesey | Resonant x-ray diffraction from chiral electric-polarization structures[END_REF]. This method was successfully applied to detect different topological structures, including vortices [START_REF] Rodriguez | Vortex polarization states in nanoscale ferroelectric arrays[END_REF], skyrmions [START_REF] Das | Observation of room-temperature polar skyrmions[END_REF], and chiral domain walls [START_REF] Chauleau | Electric and antiferromagnetic chiral textures at multiferroic domain walls[END_REF].

Note that Liu et al. [START_REF] Liu | An isogeometric approach to flexoelectric effect in ferroelectric materials[END_REF] revealed that an axial polarization component of the vortex can appear in ferroelectric PbTiO3 nanodots due to the flexoelectric effect. Thus, the paper [START_REF] Liu | An isogeometric approach to flexoelectric effect in ferroelectric materials[END_REF] and this work predict different flexo-sensitive vortex-like states with an axial polarization in ferroelectrics nanoparticles of various geometry. At that this work reveals the critical influence of the flexoelectric coefficients sign and anisotropy on the appearance and properties of the axial polarization, and, most important, on the chirality of a ferroelectric vortex. Qualitatively, both works, [START_REF] Liu | An isogeometric approach to flexoelectric effect in ferroelectric materials[END_REF] and this one, illustrate that flexoelectricity can change the chiral state of a polarization texture, and this work studies the behavior of a topological index (in fact a skyrmion number) that quantifies the change (see Fig. 3g).

The main features characterizing polarization structures stabilized by DMI-type interactions are their breaking of chiral symmetry and their incommensurability, i.e., a longperiod modulation in space that is unrelated to the crystalline lattice parameter. The appearance of such chiral incommensurate phases can generally be attributed to specific energy terms in the Landau-Ginzburg potential, known as Lifshitz invariants [START_REF] Bogdanov | On the stability of localized states in nonlinear field models with Lifshitz invariants[END_REF][START_REF] Izyumov | Modulated, or long-periodic, magnetic structures of crystals[END_REF]. Only a few ferroelectrics have crystalline structures whose symmetry allows such Lifshitz invariants; therefore, an interaction directly analogous to the magnetic DMI is generally not found in ferroelectrics. However, although not identical in its mathematical form, the energy density due to the flexoelectric coupling [START_REF] Megaw | Temperature changes in the crystal structure of barium titanium oxide[END_REF] is similar to a DMI-type energy term related to a linear Lifshitz invariant.

Accordingly, we argue that the flexoelectric coupling can lead to polarization states with properties similar to those generated by a DMI-type interaction. Such a connection between flexoelectric coupling, Lifshitz invariants, and DMI has been discussed before in the case of liquid crystals [START_REF] Sparavigna | Electric-field effects on the spin-density wave in magnetic ferroelectrics[END_REF]. It was found that flexoelectricity in liquid crystals can play a central role in the development of modulated phases that are analogous to those known from chiral ferromagnets [START_REF] Sparavigna | Role of Lifshitz invariants in liquid crystals[END_REF]. Our study shows that -similar to liquid crystals, where elastic strains fields couple to electric fields-the flexoelectric effect through which strain fields couple to the electric polarization field can lead to comparable modulated phases with chiral properties in a ferroelectric nanoparticle.

V. SUMMARY AND CONCLUSIONS

Using FEM simulations based on the phenomenological LGD approach alongside electrostatic equations and elasticity theory, we identify a characteristic polarization structure developing between two oppositely oriented axial domains located near the cylinder ends. This polar structure, which we named "flexon", displays chiral features that are connected to the flexoelectric coupling. In the azimuthal plane, the flexon polarization forms a localized chiral structure resembling a meron, or a vortex with a central kernel. Analytical calculations and FEM prove that the flexon axial polarization, and thus its chirality, switches upon a change of the sign of the flexoelectric coefficients. We also observe that the anisotropy of the flexoelectric coupling critically influences the flexon formation and the related domain morphology. This observation corroborates the link between chirality and flexocoupling, and it identifies the flexoelectric effect as the driving force stabilizing these structures.

While in magnetic systems with strong DMI, similar localized chiral structures have been reported [START_REF] Rybakov | New type of stable particle-like states in chiral magnets[END_REF], the polarization state discussed here is formed without any ferroelectric counterpart of the DMI [START_REF] Zhao | Dzyaloshinskii-Moriyalike interaction in ferroelectrics and antiferroelectrics[END_REF]. We recall that, like the DMI, the flexoelectric coupling is derived as a Lifshitz invariant [Eq. (1f)] in the context of the Landau theory of phase transitions [22], and that such linear Lifshitz invariants generally play a key role in the formation of helical structures [START_REF] Bogdanov | On the stability of localized states in nonlinear field models with Lifshitz invariants[END_REF][START_REF] Izyumov | Modulated, or long-periodic, magnetic structures of crystals[END_REF]. The fact that both the DMI and the flexoelectric stabilize structures with a specific chirality demonstrates an analogy between these two interactions which appears to have been overlooked in the literature of ferroelectric solids. An important difference compared to the classical DMI is that the flexoelectric coupling is ubiquitous in ferroelectrics, whereas the ferroelectric DMI is forbidden by symmetry in most material types. Therefore, a chiral interaction mediated by the flexoelectric effect can potentially be found in all ferroelectrics. The coupling of the electric polarization and elastic strain gradients could thus be a much more commonly accessible alternative interaction for the formation of chiral and achiral structures [START_REF] Yudin | Bichiral structure of ferroelectric domain walls driven by flexoelectricity[END_REF][START_REF] Eliseev | Structural phase transitions and electronic phenomena at 180-degree domain walls in rhombohedral BaTiO3[END_REF]. This coupling could also open the possibility of generating and dissolving chiral polarization states through strain engineering [START_REF] Damodaran | New modalities of strain-control of ferroelectric thin films[END_REF].

We predict that the pronounced flexon-type polarization distribution with two axial counter domains exists in the temperature range 290 K < 𝑇 < 340 K. The relatively wide temperature range (about 50 K) corresponding to the stability or meta-stability of the flexon-type polarization distribution give us the hope that the flexons can be observed experimentally.

However, the analysis of the hysteresis loops leads to the conclusion that flexons and other domain configurations cannot be resolved from macroscopic measurements of the average polarization in a homogeneous electric field. We anticipate that flexons can be reliably observed, e.g., by the local methods using a strong gradient of electric field, such as PFM, which gives us the information about the distribution of polarization with a nanoscale resolution.

𝜎 𝑖𝑗 is the stress tensor, 𝑠 𝑖𝑗𝑘𝑙 is the elastic compliances tensor, 𝑄 𝑖𝑗𝑘𝑙 is the electrostriction tensor, and 𝐹 𝑖𝑗𝑘𝑙 is the flexoelectric tensor.

For cubic (m3m) point symmetry group of the parent phase the explicit form of the "half" Lifshitz invariant for the flexoeffect is 𝐹 11

∆𝐺 𝑓𝑙𝑒𝑥𝑜 = [𝜎 11
+ 𝜕𝑃 3 𝜕𝑥 2 )] (A.2)
Allowing for the Khalatnikov mechanism of polarization relaxation [12], minimization of the free energy (A.1) with respect to polarization leads to three coupled time-dependent Euler-Lagrange equations for polarization components inside the core,

𝛿𝐺 𝛿𝑃 𝑖 = -𝛤 𝜕𝑃 𝑖

𝜕𝑡

, where 𝑖 = 1, 2, 3. The explicit form of the equations for a ferroelectric crystal with m3m parent symmetry is: where n is the outer normal to the surface, i = 1, 2, 3. In our FEM studies, we use the so-called "natural" boundary conditions corresponding to 𝑎 𝑖𝑗 (𝑆) = 0.

𝛤 𝜕𝑃 1 𝜕𝑡 + 2𝑃 1 (
The electric field components 𝐸 𝑖 in Eq.(A.3) are derived from the electric potential 𝜑 in a conventional way, 𝐸 𝑖 = -𝜕𝜑 𝜕𝑥 𝑖 ⁄ . The potential 𝜑 𝑓 satisfies the Poisson equation in the ferroelectric cylinder (subscript "f"):

𝜀 0 𝜀 𝑏 ( 𝜕 2 𝜕𝑥 1 2 + 𝜕 2 𝜕𝑥 2 2 + 𝜕 2 𝜕𝑥 3 2 ) 𝜑 𝑓 = 𝜕𝑃 𝑖 𝜕𝑥 𝑖 , 𝑥 1 2 + 𝑥 2 2 ≤ 𝑅 ∩ 0 ≤ 𝑥 3 ≤ ℎ, (A.5a)
The electric potential 𝜑 𝑒 in the external region outside the core-shell nanoparticle satisfies the Laplace equation (subscript "e"):

𝜀 0 𝜀 𝑒 ( 𝜕 2 𝜕𝑥 1 2 + 𝜕 2 𝜕𝑥 2 2 + 𝜕 2 𝜕𝑥 3 2 ) 𝜑 𝑒 = 0, 𝑥 1 2 + 𝑥 2 2 > 𝑅 ∪ 𝑥 3 < 0 ∪ 𝑥 3 > ℎ, (
+ 𝜕𝑃 3 𝜕𝑥 2 ) (A.8f)
Equations (A.8) should be considered along with equations of mechanical equilibrium 𝜕𝜎 𝑖𝑗 (𝒙) 𝜕𝑥 𝑖 ⁄ = 0 and compatibility equations 𝑒 𝑖𝑘𝑙 𝑒 𝑗𝑚𝑛 𝜕 2 𝑢 𝑙𝑛 (𝒙) 𝜕𝑥 𝑘 𝜕𝑥 𝑚 ⁄ = 0, which are equivalent to the mechanical displacement vector Ui continuity [18]. The boundary conditions for elastic stresses are the virtual absence of their normal components at the nanoparticle surface, as the ambient media is regarded absolutely soft.

B. Finite Element Modelling Details

Electrical boundary conditions. We consider the case when the surface of the core is covered by an elastically soft ultrathin (its thickness ∆𝑅~(1 -4) nm) paraelectric or high-k semiconductor shell with a screening length Λ~1 nm. The coverage can be artificial (e.g., a soft organic semiconductor or vacancy-enriched SrTiO3) or natural, where in the latter case it would originate from the polarization screening by surrounding media. Note that a screening length larger than 0.1 nm weakly effects the core domain structure, in this case the screening acts as an electrical open-circuit condition. Thus, we

impose an electrical open-circuit condition at the cylinder ends to make the vortex-type polarization energetically favorable.

Initial conditions, shape, and gradient effects variability. We use a four 90-degree domain configuration in the XY-plane with superimposed small random noise as an initial distribution of polarization. These four domains determine the direction of the lateral polarization vorticity. When using a purely random noise as the initial distribution of polarization, the domain structure relaxes to structure resembling 180-degree domains in thin c + /c -films. This process occurs because the depolarization field favors a polarization orientation along the z-axis of the elongated cylinder; however, the vorticity of polarization components appears near the cylinder ends (an analog of the flux-closure a-domains in thin films).

FEM simulations are performed in COMSOL@MultiPhysics software, using electrostatics, solid mechanics, and general math (PDE toolbox) modules. The size of the computational region is not less than 4040160 nm 3 , and is commensurate with the cubic unit cell constant (about 0.4 nm) of

BaTiO3 at room temperature. The minimal size of a tetrahedral element in a mesh with fine discretization is equal to the unit cell size, 0.4 nm, the maximal size is (0.8 -1.2) nm in the core, 1.6 nm in the (1-4) nm thick elastically soft paraelectric or high-k semiconductor shell, and 4 nm in the dielectric medium (Fig. S1). The dependence on the mesh size is verified by increasing the minimal size to 0.8 nm. We verified that this only results in minor changes in the electric polarization, electric field, and elastic stress and strain, such that the spatial distribution of each of these quantities becomes less smooth (i.e., they contain numerical errors in the form of a small random noise). However, when using these larger cell sizes, all significant details remain visible, and more importantly, the system energy remains essentially the same with an accuracy of about 0.1%. Effect of the cylinder ends shape on the polarization distribution FEM, performed for the sharp and rounded ends of cylindrical core at room temperature with nominal values of the flexoelectric tensor components 𝐹 𝑖𝑗 (listed in Table SI), reveals that the distribution of the polarization component 𝑃 3 depends critically on the curvature radius of the cylinder ends (Fig. S2). The increase of the curvature leads to both a strong increase of the 𝑃 3 value and a simultaneous gradual transformation of the multiple edge domains into a bidomain c + /c -configuration, which contains flux-closure domains near the spherical ends (top row in Fig. S2). For the sharp and slightly rounded ends, the distribution 𝑃 3 becomes much more contrasting with increasing 𝐹 𝑖𝑗 ; but the dependence of 𝑃 3 on 𝐹 𝑖𝑗 virtually disappears as the curvature increases (compare the top and bottom rows in Fig. S2). FEM results presented below are performed for the cylindrical core with sharp ends, because this shape is the most sensitive to the flexoelectric coupling. Note, that the vortex stability in a ferroelectric nanoparticle depends on many physical factors including the particle geometry, electric screening and elastic boundary conditions. In particular, depending on the electric boundary conditions and aspect ratio of length to radius, the vortex-type polarization can be stable, metastable or unstable in a very prolate cylinder. Here we use not very prolate cylinders with an aspect ratio varying from 4:1 to 8:1 (see cylinders with an aspect ratio much more than 10 covered by the shell with a very small or zero effective screening length Λ (see e.g., Ref. [24]).

Effect of polarization gradient coefficients on the polarization distribution

We performed FEM for temperatures (200 -400) K using zero and nominal values of the -d.

The influence of small 𝑔 𝑖𝑗 values can be explained in the following way. According to the Ginzburg-Landau theory, the thickness of the domain wall is directly associated with the gradient coefficient. Thus, a decrease in the gradient coefficients leads to a decrease in the thickness of the domain walls, which affects the size of the domains in the transverse direction. This is exactly what is shown in Fig. S3, calculated for small values of the gradient coefficients. The small 𝑔 𝑖𝑗 values leads to the decrease of the domain wall energy, and, consequently, can lead to the formation of labyrinthine domain patterns [3,4] or other gradient-driven structural instabilities, such as domain wall meandering [25,26]. An increase in flexoelectric coefficients leads to an increase in effective gradient coefficients, as a consequence of this, the domain size increases. More importantly, the increase of |𝐹 𝑖𝑗 | to the value comparable with

√ |𝑔 𝑖𝑗 | 𝑐
, where 𝑐 is the geometry-dependent combination of elastic modulus, can lead to the appearance of spatially modulated domain patterns [27]. For the small 𝑔 𝑖𝑗 values used in Fig. S3b-d, the influence of 𝑔 𝑖𝑗 and 𝐹 𝑖𝑗 support the formation of spatially-modulated domain spots and Bloch points twisting (similar to meandering).

FEM performed for sharp and rounded ends of a cylindrical BaTiO3 core reveals that the distribution of the polarization components depends critically on the curvature radius of the cylinder ends and is highly sensitive to the values of polarization gradient coefficients.

Effect of the shell dielectric properties on the domain structure

The dominant vortex-type polarization 𝑃 ⃗ 1 + 𝑃 ⃗ 2 , which is directed tangentially to the lateral surface of the cylinder (Fig. S4, top row), appears much bigger than the axial component 𝑃 3 perpendicular to the cylinder ends. Actually, a comparison of the scales for 𝑃 3 and 𝑃 in Fig. S4 shows that 𝑃 ≫ |𝑃 3 | and 𝑃 ≈ |𝑃 ⃗ 1 + 𝑃 ⃗ 2 |, and so |𝑃 ⃗ 1 + 𝑃 ⃗ 2 | ≫ |𝑃 3 |. Consequently, the resulting depolarization field 𝐸 3 created by the component 𝑃 3 is very small inside the particle. The field 𝐸 3 is determined by the dielectric properties of the shell and, in particular, 𝐸 3 slightly decreases with the increase of shell dielectric permittivity. A shell with a very high relative dielectric permittivity (e.g., vacancyenriched paraelectric SrTiO3) can provide an effective screening of the axial polarization component in a prolate cylindrical multiaxial ferroelectric core, while the dielectric properties of the shell virtually do not effect on vortex-type polarization 𝑃 ⃗ 1 + 𝑃 ⃗ 2 . Being energetically preferable, the appearance of a shell-insensitive vortex-type polarization in a multiaxial ferroelectric core, such as BaTiO3, is very probable [28]. This qualitative statement is confirmed quantitatively by the very small 𝑑𝑖𝑣𝑃 ⃗ and depolarization electric field 𝐸 3 in the core (Fig. S9a andS9c). A strong influence of the shell dielectric properties on the core domain structure is possible for a uniaxial ferroelectric core [29], since a polarization rotation is impossible in this case. In particular, a weak screening can lead to the formation of labyrinthine domains in nanoparticles of uniaxial ferroelectrics CuInP2S6 [3] and Sn2P2S6 [4].

The elastic properties of the shell are important for spherical core-shell nanoparticles, where the vortices orientation is determined by an elastic anisotropy of the shell [30]. FEM calculations show that in the case of the cylindrical core-shell BaTiO3 nanoparticles, the vortex axis is oriented along the 1.9 • 10 -11 m 3 /C, corresponding to the isotropic symmetry of 𝐹 𝑖𝑗 (Fig. S6e). We revealed that an initial four-domain polarization distribution relaxes to a flexon-like domain structure in the temperature range 260 K < 𝑇 < 360 K (middle row in Figs. S7b-d). The structure of the P-vector is vortex-like, which weakly changes when approaching the surface in the same temperature range (as seen from the direction of arrows at the bottom image of Figs. S7b-d). A relatively small domain wall broadening exists near the cylinder ends, which can be seen from the polarization magnitude distribution at 260 K < 𝑇 < 360 K (top images of Figs. S7b-d).

The flexon gradually disappears at 𝑇 > 370 K (middle image in Figs. S7e). The ferroelectric polarization inside the core significantly decreases at 𝑇 > 370 K and completely disappears at 𝑇~400 K (top image in Figs. S7e). The lateral components 𝑃 1 and 𝑃 2 form a vortex-like structure and their distribution is nearly independent of the coordinate 𝑧 along the cylinder axis. Hence, the vortexlike structure of the P-vector is insensitive to the surface presence. The structure becomes faint with a temperature increase above 370 K (note that the arrows length decreases in the bottom image of The bottom row in Fig. S8 shows 𝑃 3 distributions for different voltages, from which it follows that different 𝑃 3 distributions averaged over the core volume correspond to the same average value, 𝑃 3 ̅̅̅ (𝑈). Thus, the analysis of the hysteresis loops leads to the conclusion that flexons cannot be distinguished from macroscopic (e.g., capacitance) measurements of the average polarization in a homogeneous electric field, they can only be registered by local probing methods using the strong gradient of electric field with a nanoscale resolution, such as piezoresponse force microscopy (PFM).

U = 0 U = -5 V U = +5 V -20 V +20 V 𝑃 ̅
300 𝐾, the distance between the electrodes is 60 nm. All other parameters are listed in .

In the virtual absence of a depolarization field, which is true with high accuracy for a vortex type polarization, the divergence 𝑑𝑖𝑣𝑃 ⃗ = 𝜕𝑃 𝑖 𝜕𝑥 𝑖 should be very small. Indeed, the approximation 𝑑𝑖𝑣𝑃 ⃗ ≈ 0 is valid inside the cylindrical core (Fig. S9b).

We also note the polarization magnitude 𝑃 = √𝑃 1 2 + 𝑃 2 2 + 𝑃 3 2 is very close to the constant value everywhere, except for the vortex core (Fig. S5, middle row), and its derivative

𝜕𝑃 𝜕х 3
is negligibly small everywhere, except in the immediate vicinity of the vortex core when in contact with the cylinder ends.

For this case we obtain that:

𝜕𝑇𝑟𝜎 ∂х ) is present in the Lifshitz invariant (1f). A similar invariant has been discussed in Ref. [31] in the context of incommensurate phases with a defined chirality. This is consistent with our finding of the formation of a flexon induced by the invariant.

F ij = 0 (a) ___6F ij ≫ 0 ___2 F ij > 0 ___ F ij > 0
The integration here is performed over the nanoparticle volume 𝑉 = 𝜋𝑅 2 ℎ. The toroidal moment (F.5) is independent on the coordinate origin of radius-vector 𝑟 , only if all three components of electric polarization have a zero average over the nanoparticle volume. This is true in the considered case for 𝑈 = 0 (no applied voltage).

Let us estimate the toroidal moment using the following approximation for polarization distribution:

FIGURE 1 .

 1 FIGURE 1. A cylindrical ferroelectric nanoparticle (core) of radius 𝑅, covered with an elastically soft semiconducting shell with a thickness ∆𝑅 ≪ 𝑅 and screening length Λ of 1 nm, placed in an isotropic elastically soft effective medium. The direction of axial polarization 𝑃 3 is shown by the straight orange arrow, and lateral components 𝑃 1,2 are shown by the curled red-blue arrow to highlight their vortex-type structure.

𝑃 3

 3 inside the kernel is determined by random noise in the initial conditions. Images in Figs. 2b-d and 3b-e are calculated for a nonzero electrostriction coupling (𝑄 𝑖𝑗 ≠ 0) and either negative, zero, or positive values of the flexoelectric coefficients 𝐹 𝑖𝑗 . In the presence of electrostriction coupling the dipolar kernel disappears completely (Figs. 2c and 3c). The flexoelectric effect induces an axial component of polarization consisting of two oppositely oriented diffuse 𝑃 3 -domains located near the cylinder ends and separated by a region with 𝑃 3 ≈ 0 (Figs. 2b, 2d and 3b, 3d).

FIGURE 2 .(

 2 FIGURE 2. Distribution of the polarization component 𝑃 3 (the top row) inside a cylindrical nanoparticle and a magnified view on the flexon structure (the bottom row). The arrows show the orientation of polarization vector 𝑷. The images are calculated without electrostriction (𝑄 𝑖𝑗 = 0) and flexoelectric (𝐹 𝑖𝑗 = 0) couplings (a); with electrostriction coupling (𝑄 𝑖𝑗 ≠ 0) and negative (b), or zero (c), or positive (d) values of flexoelectric coefficients 𝐹 𝑖𝑗 . The values of 𝐹 𝑖𝑗 and all other parameters are given in Table SI, 𝑇 = 300 K. Note the different scales for 𝑃 3 -distributions in plots (a)-(d) in order to maintain a contrast between the different regions.

  axial polarization. The main effect of a change of sign in the flexoelectric coefficients is the reorientation of the flexon axial polarization. The polarization structures at the wire ends shown in Fig. 2b-2d and Fig. 3b-2d display localized chiral structures with different chirality on opposite ends of the wire, and their chirality changes upon reversal of the sign of the flexoelectric coupling constant.

1 4𝜋

 1 𝜌, 𝜑, 𝑧), is the flexon.In order to quantify the chirality of the polarization structure and its variation along the cylinder axis, in Appendix F[63] we calculate the topological index 𝑛 =

FIGURE 3 . 4 .[

 34 FIGURE 3. Distribution of the polarization component 𝑃 3 (the top row) in the XZ cross-section of the nanoparticle. Images are calculated without electrostriction (𝑄 𝑖𝑗 = 0) and flexoelectric (𝐹 𝑖𝑗 = 0) coupling (a); with electrostriction coupling (𝑄 𝑖𝑗 ≠ 0) and negative (b), zero (c), positive (d), or high positive (e) values of flexoelectric coefficients 𝐹 𝑖𝑗 . The bottom part (f) is the dependence of the maximal (red curve, 𝑃 𝑚𝑎𝑥 ) and minimal (blue curve, 𝑃 𝑚𝑖𝑛 ) values of 𝑃 3 on the relative amplitude of the flexoelectric coupling strength 𝑓. The green curve is the extremal (maximal or minimal) value 𝑃 𝑒 in the center of the top axial 𝑃 3 -domain. Here 𝐹 𝑖𝑗 = 𝑓𝐹 𝑖𝑗 0 , the values of 𝐹 𝑖𝑗 0 and all other parameters are given in Table SI, 𝑇 = 300 K. The Z-profile of the polarization topological index 𝑛(𝑧) is shown in the inset (g) for zero (black line), positive (red curve), and negative (blue curve) 𝐹 𝑖𝑗 . Note the different scales for 𝑃 3 in the plots (a) and (b)-(e) in order to maintain a contrast between the different regions.

FIGURE 4 .

 4 FIGURE 4. Distribution of the polarization component 𝑃 3 in the XZ cross-section of cylindrical core (the top row). Images are calculated for the fixed values 𝐹 11 = 2.4 •10 -11 m 3 /C and 𝐹 12 = 0.5 • 10 -11 m 3 /C, while the value of 𝐹 44 varies from -6 to 6 (in 10 -11 m 3 /C) as indicated in the legends. The bottom part (f) is the dependence of the maximal (red curve, 𝑃 𝑚𝑎𝑥 ) and minimal (blue curve, 𝑃 𝑚𝑖𝑛 ) values of 𝑃 3 on the relative amplitude f of the flexoelectric coefficient 𝐹 44 in the core. The green curve is the extremal (maximal or minimal) value 𝑃 𝑒 in the center of the bottom axial 𝑃 3 -domain. Here 𝐹 44 = 𝑓𝐹 44 0 and 𝐹 44 0 = 0.06 • 10 -11 m 3 /C. The electrostriction coupling coefficients 𝑄 𝑖𝑗 and all other parameters are listed in Table SI, 𝑇 = 300 K. Note the different scales for 𝑃 3 in the plots (a)-(e) in order to maintain a contrast between the different regions.

6 P3

 6 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm 2 ) x-coordinate x-coordinate x-coordinate x-coordinate Flexoelectric coefficient F44 amplitude f (rel. units) distribution of polarization in a cylindrical core. Typical FEM results are shown in Fig. 5, where the columns (a)-(e) correspond to the temperature increase from 240 K to 370 K; the structure of the azimuthal components of the polarization vector, 𝑃 1 and 𝑃 2 , is vortex-like and shows weak variations when approaching the surface over the same temperature range (see the direction of arrows at the bottom image of Figs. S7 [63]).

FIGURE 5 .

 5 FIGURE 5. Distribution of polarization component 𝑃 3 in XZ cross-sections of the nanoparticle core. Different columns are calculated for the temperatures 𝑇 = 240, 270, 300, 340, and 370 K (a, b, c, d, e). All other parameters are listed in Table I. Note the different scales for 𝑃 3 in the plots (a) and (b)-(e) in order to maintain a contrast between the different regions.

  (d) T=340 K (e) T=370 K (b) T=270 K (c) T=300 K (a) T=240 K z-coordinate x-coordinatex-coordinate x-coordinate x-coordinate x-coordinate P 3 (μC/cm 2 ) P 3 (μC/cm 2 ) P 3 (μC/cm 2 ) P 3 (μC/cm 2 ) P 3 (μC/cm2

𝜕𝑥 2 ) + 𝐸 3 (

 23 A.3c)The temperature-dependent Khalatnikov coefficient 𝛤[13] determines the relaxation time of the polarization 𝜏 𝐾 = 𝛤 |𝛼| ⁄ . Consequently, 𝜏 𝐾 typically varies in the range (10 -9 -10 -6 ) seconds for temperatures far from TC. As argued by Hlinka et al.[14], we assumed that 𝑔 44 ′ = -𝑔 12 in Eqs.(A.3).The boundary condition for polarization at the core-shell interface 𝑟 = 𝑅 accounts for the flexoelectric effect:

  A.5b) Equations (A.5) are supplemented with the continuity conditions for electric potential and linear screening conditions for the normal components of the electric displacements at the cylinder surface S [4]: (𝜑 𝑒 -𝜑 𝑓 )| 𝑆 = 0, 𝒏(𝑫 𝑒 -𝑫 𝑓 )boundary condition (A.6a) corresponds to the surface of the core covered by an ultrathin semiconductor shell with an effective screening length Λ ≥1 nm [15, 16, 17]. Note that a screening length greater than 0.1 nm leads to the domain formation in the core. The case Λ → ∞ corresponds to an electrical open-circuit condition. We impose an electrical open-circuit condition at the cylinder ends to make the vortex-type polarization energetically favorable. Either charges are absent or the applied voltage is fixed at the boundaries of the computation region: state follow from the variation of the energy (A.1e) with respect to elastic stress, 𝛿𝐺 δ𝜎 𝑖𝑗 = -𝑢 𝑖𝑗 . In the oversimplified case

Figure S1 .

 S1 Figure S1. Examples of mesh for the sharp edges (a), rounded edges (b) and rounded caps (c) of the cylindrical nanoparticle. The color corresponds to the size of the mesh elements.

  (c) spherical ends 𝑅 𝑐ℎ = 5 nm (b) rounded ends with 𝑅 𝑐ℎ = 1 nm (a) "sharp" ends Element size (nm) Element size (nm) Element size (nm)

Figure S2 .

 S2 Figure S2. Distribution of the polarization component P3 in three different cross-sections of the cylindrical nanoparticle with sharp ends (a), rounded ends with curvature 𝑅 𝑐ℎ =1 nm (b) and 1.8 nm (c); and hemispherical ends (d). The values of the flexoelectric coefficients are 𝐹 11 = 2.4 •10 -11 m 3 /C, 𝐹 12 = 0.5 • 10 -11 m 3 /C, and 𝐹 44 = 0.06 • 10 -11 m 3 /C for the top row; 𝐹 11 = 144 •10 -11 m 3 /C, 𝐹 12 = 3 • 10 -11 m 3 /C, and 𝐹 44 = 0.36 • 10 -11 m 3 /C for the bottom row. Temperature 𝑇 = 300 𝐾. All other parameters are listed in Table SI.

  Typical distributions of the polarization component P3, and polarization magnitude 𝑃 calculated for small 𝑔 𝑖𝑗 and the different values of electrostriction and flexoelectric coupling coefficients are shown in Fig.S3. For 𝑄 𝑖𝑗 = 0 and 𝐹 𝑖𝑗 = 0, the reduction of 𝑔 𝑖𝑗 by a factor of ten leads to the decrease of the domain wall width only (Fig.S3a). For nonzero nominal values of 𝑄 𝑖𝑗 , the reduction of 𝑔 𝑖𝑗 by a factor of ten leads to the appearance of quasi-periodic spot-like 𝑃 3 -domains located near the lateral surface of the cylinder (shown by blue and red spots in the top row of Figs.S3b-d) instead of two diffuse axial domains shown in Fig.S2a. These spot-like domains are insensitive to the magnitude and sign of 𝐹 𝑖𝑗 , if its absolute value is less than 4 • 10 -11 m 3 /C (Figs.S3b-d). The sharp transformation of the quasi-stable spot-like 𝑃 3 -domains into a stable bidomain configuration with flux-closure domains at the cylinder ends appears at unrealistically high values of 𝐹 𝑖𝑗 (Fig.S3e). The distributions of polarization magnitude 𝑃 calculated for nonzero 𝑄 𝑖𝑗 and |𝐹 𝑖𝑗 | < 4 • 10 -11 m 3 /C reveal a twisted central line with quasi-periodically located multiple Bloch points, 𝑷 = 0, whose patterns are shown in the bottom row of Figs. S3b

Figure S3 .

 S3 Figure S3. Distribution of the polarization component P3 (top row) and polarization magnitude 𝑃 (bottom row) in three different cross-sections of the nanoparticle core. Different columns are calculated for the different values of electrostriction and flexoelectric coupling coefficients. The following values of gradient coefficients are used: 𝑔 11 = 1 •10 -10 m 3 /F, 𝑔 12 = -0.04 • 10 -10 m 3 /F, and 𝑔 44 = 0.04 10 -11 m 3 /F. Temperature 𝑇 = 300 𝐾. All other parameters are listed in Table SI.

Figure S4 .

 S4 Figure S4. Distribution of the polarization vector P inside a cylindrical nanoparticle (top row) and polarization magnitude 𝑃 (bottom row) in the XZ cross-section of the nanoparticle. The images are calculated without electrostriction (𝑄 𝑖𝑗 = 0) and flexoelectric (𝐹 𝑖𝑗 = 0) couplings (a). In the case where electrostriction coupling in non-zero (𝑄 𝑖𝑗 ≠ 0), the values of the flexoelectric coefficients 𝐹 𝑖𝑗 are varied: negative (b), zero (c), and positive (d). The arrow color corresponds to the value of 𝑃 3 (red-blue scale). The values of 𝐹 𝑖𝑗 and all other parameters are given in Table SI, 𝑇 = 300 K.

2 )Figure S5 .

 2S5 Figure S5. Distribution of the polarization component P3 (top row) and polarization magnitude (middle row) in the XZ cross-section of the nanoparticle. Bottom row: isosurfaces of zero polarization components P1 (gray), P2 (green), and nonzero 𝑃 3 =  0.1,  0.1,  0.2, and  0.4 µC/cm 2 (red and blue). Images are calculated without electrostriction (𝑄 𝑖𝑗 = 0) and flexoelectric (𝐹 𝑖𝑗 = 0) coupling (a); with electrostriction coupling (𝑄 𝑖𝑗 ≠ 0) and

Figure S6 .

 S6 Figure S6. Distribution of polarization component 𝑃 3 in three different cross-sections of the nanoparticle core (top row), and 𝑃 3 in the XZ cross-section of cylindrical core (middle row). Bottom row: isosurfaces of polarization components 𝑃 1 (gray), 𝑃 2 (green), and 𝑃 3 (red and blue). Zero values correspond 𝑃 1 and 𝑃 2 isosurfaces, while isosurfaces of 𝑃 3 correspond to the values 𝑃 3 =  0.1,  0.2, and  0.4 µC/cm 2 for the images (a)-(d), (e), and (f) respectively. The images are calculated for the fixed values 𝐹 11 = 2.4 •10 -11 m 3 /C and 𝐹 12 = 0.5 • 10 -11 m 3 /C, while the value of 𝐹 44 varies from -6 to 6 (in 10 -11 m 3 /C) as indicated in the legends. All other parameters are listed in Table SI. Temperature 𝑇 = 300 𝐾.

6 P3

 6 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm 2 ) Behavior of the Flexon-Type Polarization Distribution To define the temperature interval, where flexons exist being stable or meta-stable, we performed FEM in the temperature range from 50 K to 400 K using different initial distributions of polarization in a cylindrical core. Typical FEM results are shown in Fig.S7, where the columns (a)-(e) correspond to the temperature increase from 240 K to 370 K.A bidomain configuration of 𝑃 3 is stable at temperatures lower than 250 K (Figs.S7a). The bidomain structure has a relatively thin uncharged 180-degree domain wall inside the cylinder, which transforms into a flux-closure domain near the electrically-open cylinder ends (top and middle images in Figs. S7a). The origin of the flux-closure domain is the core tendency to minimize its electrostatic energy, because a flux-closure domain wall creates a much weaker depolarization field (in fact negligibly small) compared to the field that the charged 180-degree domain wall would create. The structure of the P-vector is vortex-like and changes moderately when approaching the surface at 𝑇 ≤ 240 K (as seen from the direction of arrows in the bottom image of Figs. S7a).

Figs. S7e).

  Figs. S7e). The flexon-type polarization distribution exists over a relatively wide temperature range 260 K < 𝑇 < 360 K. The axial counter domains inherent to flexons are the most pronounced feature over the narrower range 290 K < 𝑇 < 340 K. The relatively wide temperature range (about 100 K) corresponding to the stability or meta-stability of the flexon-type polarization distribution gives us the hope that the domain morphology can be observed experimentally.

Figure S7 .

 S7 Figure S7. Distribution of polarization magnitude 𝑃 (top row) and its component 𝑃 3 (middle row) in XZ crosssections of the nanoparticle core. Bottom row: distribution of polarization vector P shown by arrows colored according to the value of 𝑃 3 . Different columns are calculated for the temperatures 𝑇 = 240, 270, 300, 340, and 370 K (a, b, c, d, e). All other parameters are listed in Table SI.

Figure S8 .

 S8 Figure S8. Dependence of the average polarization component 𝑃 3 ̅̅̅ (black loops) and its minimal 𝑃 𝑚𝑖𝑛 (blue loops) and maximal 𝑃 𝑚𝑎𝑥 (red loops) values on the applied voltage 𝑈. Black dots connected with a vertical dashed line correspond to 𝑈 𝑖𝑛𝑡 . Plots (a-d) are calculated for different values of electrostriction and flexoelectric coupling coefficients, as written in the legends. The plots (a-c) are calculated for nonzero electrostriction coupling coefficients (𝑄 𝑖𝑗 ≠ 0) and negative (a), zero (b), or positive (c) flexoelectric coupling coefficients 𝐹 𝑖𝑗 . Plot (d) is calculated without electrostriction (𝑄 𝑖𝑗 = 0) and flexoelectric (𝐹 𝑖𝑗 = 0) coupling. Temperature 𝑇 =

Figure S9 .

 S9 Figure S9. Distribution of polarization component 𝑃 𝑧 (a), its divergency 𝑑𝑖𝑣𝑃 ⃗ (b), electric field component 𝐸 𝑧 (c), nonzero elastic stresses (d-i), their derivatives (j, k), strain components (l,m) and its gradients (n, o) in three different cross-sections of the cylindrical core. A cylindrical coordinate system {𝑟, 𝜑, 𝑧} is used for the mechanical stress, strain, electric field and polarization. Note that the stress component 𝜎 𝑧𝜑 = 0 (in spherical coordinates). Temperature 𝑇 = 300 𝐾. All other parameters are listed in Table SI.

2 ,

 2 Index and Toroidal MomentFirst, let us calculate the topological index 𝑛 of the unit polarization orientation[32] positive, and negative flexoelectric tensor coefficients. The integration is performed over the cylinder cross-section, and here 𝑝 = 𝑃 ⃗ 𝑃 . The flexon polarization has the following structure in Cartesian coordinates 𝑃 1 ≈ 𝑝(𝑟, 𝑧) cos[𝜑(𝛼, 𝑧)], 𝑃 2 ≈ 𝑝(𝑟, 𝑧) sin[𝜑(𝛼, 𝑧)], 𝑃 3 ≈ 𝜂(𝑟, 𝑧). (F.2a) Here 𝑝(𝑟, 𝑧) > 0, and 𝜑(𝛼, 𝑧) ≈ 𝛼 -𝜋 where 𝛼 is the polar angle, 𝑥 = 𝑟 cos 𝛼, and 𝑦 = 𝑟 sin 𝛼. The polarization magnitude is 𝑃(𝑟, 𝑧) = √𝑝 2 (𝑟, 𝑧) + 𝜂 2 (𝑟, 𝑧). Let us introduce the unit polarization as: 𝑝 ≈ {sin[𝜃(𝑟, 𝑧)] cos[𝜑(𝛼, 𝑧)] , sin[𝜃(𝑟, 𝑧)] sin[𝜑(𝛼, 𝑧)] , cos[𝜃(𝑟, 𝑧)]}, (F.2b) where the functions sin[𝜃(𝑟, 𝑧)] = 𝑝(𝑟,𝑧) √𝑝 2 (𝑟,𝑧)+𝜂 2 (𝑟,𝑧) and cos[𝜃(𝑟, 𝑧)] = 𝜂(𝑟,𝑧) √𝑝 2 (𝑟,𝑧)+𝜂 2 (𝑟,𝑧) . The gradients of the unit polarization can be written as we used 𝜕𝑟/𝜕𝑥 = 𝑥/𝑟 and 𝜕𝑟/𝜕𝑦 = 𝑦/𝑟. Below we use the following relations this relation remains unchanged even if one could take into account the radial dependence of polarization magnitude, 𝑃(𝑟, 𝑧) (F.2b) and reproducing the detailed calculations in Ref. [33], we obtain: 𝑛(𝑧) = ∫ sin[𝜃(𝑟, 𝑧)] 𝜃(𝑟, 𝑧)]] 𝜃(𝑟=0,𝑧) 𝜃(𝑟=𝑅,𝑧) [𝜑(𝛼, 𝑧)] 𝛼=0 𝛼=2𝜋 (F.3a) Here R is the cylinder radius. Substituting here 𝜑(𝛼, 𝑧) ≈ 𝛼 -𝜋 2 and cos[𝜃(𝑟, 𝑧)] = 𝜂(𝑟,𝑧) √𝑝 2 (𝑟,𝑧)+𝜂 2 (𝑟,𝑧) 𝜂(𝑅, 𝑧) = 0 and 𝑝(𝑅, 𝑧) > 0. The dependence 𝑛(𝑧) is shown in Fig. S10a for zero (green horizontal line), positive (3 solid curves), and negative (3 dashed curves) flexoelectric tensor coefficients 𝐹 𝑖𝑗 . The black curves 𝐹 𝑖𝑗 values are listed in TableSI; the red curves are calculated for twice the value of 𝐹 𝑖𝑗 (labeled as "2𝐹 𝑖𝑗 ") and the blue curves are calculated for six times the value of 𝐹 𝑖𝑗 (labeled as "6𝐹 𝑖𝑗 "). The Z-profile of the axial polarization 𝑃 3 (0, 𝑧) and polarization magnitude 𝑃(0, 𝑧) are shown in Fig.S10band Fig.S10c, respectively. Symbols are calculated by FEM for positive 𝐹 𝑖𝑗 (black diamonds), 2𝐹 𝑖𝑗 (red triangles) and 6𝐹 𝑖𝑗 (blue squires). An applied voltage is absent in Fig.S10.Solid and dashed curves are the interpolation functions.

  𝑓, 𝑔, 𝐴, 𝐵, z 0 , and z 𝑚 are the fitting parameters to FEM results, which are listed in TableSII.The length scale 𝐿 = 1 nm. The amplitude 𝑓 increases with the increase of flexoelectric coupling strength and saturates at high |𝐹 𝑖𝑗 |. Since the value 𝑃(0, 𝑧) is very close to the 𝑃 3 (0, 𝑧) near the cylinder end, but 𝑃 3 (0, 𝑧) vanishes in the nanoparticle center, the topological index continuously changes from -½ to +½ with a z-coordinate change from one cylinder end (𝑧 = -20 nm) to another (𝑧 = +20 nm).

Figure S10 .

 S10 Figure S10. (a) Z-dependence of the polarization orientation topological index 𝑛(𝑧) for zero (green horizontal line), positive (solid curves), and negative (dashed curves) flexoelectric tensor coefficients 𝐹 𝑖𝑗 . For the black curves 𝐹 𝑖𝑗 values are listed in Table S1, for the red curves we used 2𝐹 𝑖𝑗 , and blue curves, we used 6𝐹 𝑖𝑗 . (b-c) Z-profile of the axial polarization 𝑃 3 (b) and polarization magnitude 𝑃 (c) calculated at 𝑟 = 0. Symbols are calculated by FEM for positive 𝐹 𝑖𝑗 (black diamonds), 2𝐹 𝑖𝑗 (red triangles), and 6𝐹 𝑖𝑗 (blue squires). Solid and dashed curves are fitting functions listed in the text. Referenced values of 𝐹 𝑖𝑗 and all other parameters are given in Table SI. Temperature 𝑇 = 300 𝐾, 𝑈 = 0.

𝑃 1 ≈ 2 𝑅 2 ∫and

 122 𝑝(𝑟, 𝑧) sin[𝛼], 𝑃 2 ≈ -𝑝(𝑟, 𝑧) cos[𝛼], 𝑃 3 ≈ 𝜂(𝑟, 𝑧), (F.6a) where the coordinates are 𝑥 = 𝑟 cos[𝛼], 𝑦 = 𝑟 sin[𝛼], 𝑧 = 𝑧. (F.6a)The vectorial product is:[𝑃 ⃗ × 𝑟 ] = -𝑒 𝑥 ⃗⃗⃗ (𝑝(𝑟, 𝑧) cos[𝛼] 𝑧 + 𝜂(𝑟, 𝑧)𝑟 sin[𝛼]) + 𝑒 𝑦 ⃗⃗⃗⃗ (𝜂(𝑟, 𝑧)𝑟 cos[𝛼] -𝑝(𝑟, 𝑧) sin[𝛼] 𝑧) + 𝑒 𝑧 ⃗⃗⃗ 𝑝(𝑟, 𝑧)𝑟 (F.7a)After the integration we obtain that only the z-component of the TM is nonzero:𝑀 ⃗⃗ ≈ 𝑒 𝑧 ⃗⃗⃗ 𝑝(𝑟,𝑧)𝑟𝑑𝑟 𝑅 0 . (F.7b) Since the magnitude 𝑝(𝑟, 𝑧) is almost independent of the flexoelectric coupling (see e.g. Fig. S5, the middle row), the TM appears nearly the same for zero, positive, and negative flexoelectric tensor coefficients. To make an analytical estimate in Eq.(F.7b), one can use the following approximation for the magnitude 𝑝(𝑟, 𝑧) ≅ 𝑝 0 𝑡𝑎𝑛ℎ ( 𝑟 𝑟 0 ), where the 𝑝 0 and 𝑟 0 are temperature-dependent. This results in: 𝑀 ⃗⃗ ≈ 𝑒 𝑧 ⃗⃗⃗ 𝑝 0 {1 -𝑑 2 ( PolyLog[𝑥] is a polylogarithmic function. Figure S11 shows the dependence of the normalized TM 𝑀 𝑝 0 on the parameter d.

Figure S11 .

 S11 Figure S11. The dependence of the normalized toroidal moment 𝑀 𝑝 0 on the parameter d.

  + (𝜎 22 + 𝜎 33 )𝐹 12 ]

					𝜕𝑃 1 𝜕𝑥 1	+ [𝜎 22 𝐹 11 + (𝜎 11 + 𝜎 33 )𝐹 12 ]	𝜕𝑃 2 𝜕𝑥 2	+ [𝜎 33 𝐹 11 +
	(𝜎 11 + 𝜎 22 )𝐹 12 ]	𝜕𝑃 3 𝜕𝑥 3	+ 𝐹 44 [𝜎 12 ( 𝜕𝑥 2 𝜕𝑃 1	+	𝜕𝑃 2 𝜕𝑥 1	) + 𝜎 13 ( 𝜕𝑥 3 𝜕𝑃 1	+	𝜕𝑃 3 𝜕𝑥 1	) + 𝜎 23 (	𝜕𝑃 2 𝜕𝑥 3

  𝑎 1 -𝑄 12 (𝜎 22 + 𝜎 33 ) -𝑄 11 𝜎 11 ) -𝑄 44 (𝜎 12 𝑃 2 + 𝜎 13 𝑃 3 ) + 4𝑎 11 𝑃 1 𝑎 1 -𝑄 12 (𝜎 11 + 𝜎 22 ) -𝑄 11 𝜎 33 ) -𝑄 44 (𝜎 13 𝑃 1 + 𝜎 23 𝑃 2 ) + 4𝑎 11 𝑃 3

										3 + 2𝑎 12 𝑃 1 (𝑃 2 2 +
			𝑃 3 2 ) + 6𝑎 111 𝑃 1 5 + 2𝑎 112 𝑃 1 (𝑃 2 4 + 2𝑃 1 2 𝑃 2 2 + 𝑃 3 4 + 2𝑃 1 2 𝑃 3 2 ) + 2𝑎 123 𝑃 1 𝑃 2 2 𝑃 3 2 -𝑔 11	𝜕 2 𝑃 1 𝜕𝑥 1 2 -
			𝑔 44 ( 𝜕 2 𝑃 1 𝜕𝑥 2 2 +	𝜕 2 𝑃 1 𝜕𝑥 3 2 ) = -𝐹 11	𝜕𝜎 11 𝜕𝑥 1	-𝐹 12 ( 𝜕𝜎 22 𝜕𝑥 1	+	𝜕𝜎 33 𝜕𝑥 1	) -𝐹 44 ( 𝜕𝜎 12 𝜕𝑥 2	+	𝜕𝜎 13 𝜕𝑥 3	) + 𝐸 1
												(A.3a)
	𝛤	𝜕𝑃 2 𝜕𝑡	+ 2𝑃 2 (𝑎 1 -𝑄 12 (𝜎 11 + 𝜎 33 ) -𝑄 11 𝜎 22 ) -𝑄 44 (𝜎 12 𝑃 1 + 𝜎 23 𝑃 3 ) + 4𝑎 11 𝑃 2 3 + 2𝑎 12 𝑃 2 (𝑃 1 2 +
			𝑃 3 2 ) + 6𝑎 111 𝑃 2 5 + 2𝑎 112 𝑃 2 (𝑃 1 4 + 2𝑃 2 2 𝑃 1 2 + 𝑃 3 4 + 2𝑃 2 2 𝑃 3 2 ) + 2𝑎 123 𝑃 2 𝑃 1 2 𝑃 3 2 -𝑔 11	𝜕 2 𝑃 2 𝜕𝑥 2 2 -
			𝑔 44 ( 𝜕 2 𝑃 2 𝜕𝑥 1 2 +	𝜕 2 𝑃 2 𝜕𝑥 3 2 ) = -𝐹 11	𝜕𝜎 22 𝜕𝑥 2	-𝐹 12 ( 𝜕𝜎 11 𝜕𝑥 2	+	𝜕𝜎 33 𝜕𝑥 2	) -𝐹 44 (	𝜕𝜎 12 𝜕𝑥 1	+	𝜕𝜎 23 𝜕𝑥 3	) + 𝐸 2
												(A.3b)
	𝛤	𝜕𝑃 3 𝜕𝑡	+ 2𝑃 3 (3 + 2𝑎 12 𝑃 3 (𝑃 1 2 +
			𝑃 2 2 ) + 6𝑎 111 𝑃 3 5 + 2𝑎 112 𝑃 3 (𝑃 1 4 + 2𝑃 3 2 𝑃 1 2 + 𝑃 2 4 + 2𝑃 2 2 𝑃 3 2 ) + 2𝑎 123 𝑃 3 𝑃 1 2 𝑃 2 2 -𝑔 11	𝜕 2 𝑃 3 𝜕𝑥 3 2 -
			𝑔 44 ( 𝜕 2 𝑃 3 𝜕𝑥 1 2 +	𝜕 2 𝑃 3 𝜕𝑥 2 2 ) = -𝐹 11	𝜕𝜎 33 𝜕х 3	-𝐹 12 ( 𝜕𝜎 11 𝜕х 3	+	𝜕𝜎 22 𝜕х 3	) -𝐹 44 (	𝜕𝜎 13 𝜕𝑥 1	+	𝜕𝜎 23

  𝑠 𝑖𝑗𝑘𝑙 𝜎 𝑘𝑙 + 𝑄 𝑖𝑗𝑘𝑙 𝑃 𝑘 𝑃 𝑙 + 𝐹 𝑖𝑗𝑘𝑙 𝜎 𝑖𝑗 = 𝑐 𝑖𝑗𝑘𝑙 𝑢 𝑘𝑙 -𝑞 𝑖𝑗𝑘𝑙 𝑃 𝑘 𝑃 𝑙 -𝑓 𝑖𝑗𝑘𝑙 𝑢 𝑖𝑗 is the strain tensor components related to the displacement components 𝑈 𝑖 in the following way: 𝑢 𝑖𝑗 = (𝜕𝑈 𝑖 /𝜕𝑥 𝑗 + 𝜕𝑈 𝑗 /𝜕𝑥 𝑖 )/2. For cubic (m3m) point symmetry group of the parent phase the 𝑢 12 = 𝑠 44 𝜎 12 + 𝑄 44 𝑃 1 𝑃 2 + 𝐹 44 ( 𝑢 13 = 𝑠 44 𝜎 13 + 𝑄 44 𝑃 1 𝑃 3 + 𝐹 44 ( 𝑢 23 = 𝑠 44 𝜎 23 + 𝑄 44 𝑃 3 𝑃 2 + 𝐹 44 (

	𝜕𝑃 𝑙 𝜕𝑥 𝑘	, 0 < 𝑟 ≤ 𝑅, 0 < 𝑧 ≤ ℎ ,	(A.7b)
	where strain components are:
	𝑢 11 = 𝑠 11 𝜎 11 + 𝑠 12 (𝜎 22 + 𝜎 33 ) + 𝑄 11 𝑃 1 2 + 𝑄 12 (𝑃 2 2 + 𝑃 3 2 ) + 𝐹 11	𝜕𝑃 1 𝜕𝑥 1	+ 𝐹 12 ( 𝜕𝑥 2 𝜕𝑃 2	+	𝜕𝑃 3 𝜕𝑥 3	) (A.8a)
	𝑢 22 = 𝑠 11 𝜎 22 + 𝑠 12 (𝜎 11 + 𝜎 33 ) + 𝑄 11 𝑃 2 2 + 𝑄 12 (𝑃 1 2 + 𝑃 3 2 ) + 𝐹 11		𝜕𝑃 2 𝜕𝑥 2	+ 𝐹 12 ( 𝜕𝑥 1 𝜕𝑃 1	+	𝜕𝑃 3 𝜕𝑥 3	) (A.8b)
	𝑢 33 = 𝑠 11 𝜎 33 + 𝑠 12 (𝜎 11 + 𝜎 22 ) + 𝑄 11 𝑃 3 2 + 𝑄 12 (𝑃 2 2 + 𝑃 1 2 ) + 𝐹 11	𝜕𝑃 3 𝜕𝑥 3	+ 𝐹 12 ( 𝜕𝑥 1 𝜕𝑃 1	+	𝜕𝑃 2 𝜕𝑥 2	) (A.8c)
				𝜕𝑃 1 𝜕𝑥 2	+	𝜕𝑃 2 𝜕𝑥 1	)	(A.8d)
				𝜕𝑃 1 𝜕𝑥 3	+	𝜕𝑃 3 𝜕𝑥 1	)	(A.8e)
					𝜕𝑃 2
					𝜕𝑥 3
	𝜕𝑃 𝑙
	𝜕𝑥 𝑘

= 𝑢 𝑖𝑗 , 0 < 𝑟 ≤ 𝑅, 0 < 𝑧 ≤ ℎ , (A.7a)

Table SI .

 SI LGD coefficients and other material parameters of BaTiO3 nanocylinders

	Coefficient	Numerical value	References
	b, e		b = 7 (core background), e = 10 (surrounding)	14
	ai (in mJ/C 2 )	a1 = 3.34(T-381)10 5 , T = 3.3410 5 (a1 = -2.9410 7 at 298 K)	19
	aij (in m 5 J/C 4 )	a11 = 4.69(T-393)10 6 -2.0210 8 , a12 = 3.23010 8 ,	19
			(a11 = -6.7110 8 at 298 K)
	aijk (in m 9 J/C 6 )	(at 298 K a111 = 82.810 8 , a112 = 44.710 8 , a123 = 49.110 8 ) a111 = -5.52(T-393)10 7 +2.7610 9 , a112 = 4.4710 9 , a123 = 4.9110 9	19
	Qij (m 4 /C 2 )	Q11=0.11, Q12= -0.043, Q44=0.059	19
	sij (in 10 -12 Pa -1 )	s11=8.3, s12= -2.7, s44=9.24	20
	gij (in 10 -10 m 3 J/C 2 ) g11=5.0, g12= -0.2, g44= 0.2	14
	fij (in V) Fij (in 10 -11 m 3 /C)	The equality 𝐹 44 = 𝐹 11 -𝐹 12 is valid in the isotropic case. the values f11 = 5.1, f12 = 3.3, f44 = 0.065 V calculated in [21] F11 = 2.4, F12 = 0.5, F44 = 0.06 (these values are recalculated from	21
	vijklm		0 (since its characteristic values are unknown for BaTiO3 and other perovskites)
	ai (s)		0 (that corresponds to the so-called natural boundary conditions)
	𝛽 𝑇 (𝑐) (in 10 -6 K -1 )	9.8 (thermal expansion coefficient)	22, 23
	𝑎 𝑐𝑢𝑏𝑖𝑐 (𝑐)	(in Å)	4.035 Å lattice constant at 1000 C	Recalculated from [23]
	R and h (in nm)

R = 10 (vary from 2 to 20 nm), h = 80 (vary from 8 to 160 nm)

Table SI

 SI

); and the cylinder ends are electrically open, which corresponds to a weak electric screening. The electrical open-circuit condition at the cylinder ends makes the vortex-type polarization energetically favorable independently on the presence of the flexoelectric coupling. The single-domain state can be stable in a much more prolate P3 (μC/cm 2 ) (d) spherical ends with 𝑅 𝑐ℎ =5 nm (c) rounded ends with 𝑅 𝑐ℎ = 1.8 nm P3 (μC/cm 2 ) (b) rounded ends with 𝑅 𝑐ℎ = 1 nm (a) "sharp" ends P3 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm 2 ) P3 (μC/cm) P 3 (μC/cm) Nominal Fij High 6 Fij

  electrostriction tensor components 𝑄 𝑖𝑗 ; and zero, positive, and negative values of 𝐹 𝑖𝑗 . Three values of 𝑔 𝑖𝑗 were taken: a nominal value and a factor of ten larger and smaller than the nominal value. The nominal values of 𝑄 𝑖𝑗 , 𝐹 𝑖𝑗 , and 𝑔 𝑖𝑗 are listed in TableSI. Large values of 𝑔 𝑖𝑗 lead to an increase of the domain wall thickness and makes the domain structure insensitive to 𝐹 𝑖𝑗 values; very large values of 𝑔 𝑖𝑗 eventually prevent the domain formation. This result expected, since, according to the Ginzburg-Landau theory, the thickness of the uncharged 180-degree domain wall is proportional to √𝑔 44 for the considered geometry. The situation for small values of 𝑔 𝑖𝑗 is much more interesting.

Table SI .

 SI The bottom row shows 𝑃 3 distributions for different voltages U. Different distributions corresponding to the same U illustrate a rather small hysteresis effect. 𝜎 12 = 𝑐 44 𝑢 12 -𝑞 44 𝑃 1 𝑃 2 -𝑓 44 ( 𝜎 13 = 𝑐 44 𝑢 13 -𝑞 44 𝑃 1 𝑃 3 -𝑓 44 (

											𝜕𝑃 1 𝜕𝑥 2	+	𝜕𝑃 2 𝜕𝑥 1	),	(E.2d)
											𝜕𝑃 1 𝜕𝑥 3	+	𝜕𝑃 3 𝜕𝑥 1	),	(E.2e)
											APPENDIX E. Approximate Analytical Solution 𝜎 23 = 𝑐 44 𝑢 23 -𝑞 44 𝑃 3 𝑃 2 -𝑓 44 ( 𝜕𝑃 2 𝜕𝑥 3 + 𝜕𝑃 3 𝜕𝑥 2 ).	(E.2f)
	Since Here 𝑞 44 =	𝜕𝜎 𝑖1 𝜕𝑥 𝑖 𝑄 44 𝑠 44 , 𝑐 44 = = 𝜕𝜎 𝑖2 𝜕𝑥 𝑖 𝑠 44 = 1	𝜕𝜎 𝑖3 𝜕𝑥 𝑖 and 𝑓 44 = = 0 in accordance with mechanical equilibrium equations, the right-𝐹 44 𝑠 44
	hand sides in Eqs. (A.3) can be rewritten as:
		𝛤	𝜕𝑃 1 𝜕𝑡	+ 2𝑃 1 (𝑎 1 -𝑄 12 (𝜎 22 + 𝜎 33 ) -𝑄 11 𝜎 11 ) -𝑄 44 (𝜎 12 𝑃 2 + 𝜎 13 𝑃 3 ) + 4𝑎 11 𝑃 1 3 + 2𝑎 12 𝑃 1 (𝑃 2 2 +
							𝑃 3 2 ) + 6𝑎 111 𝑃 1 5 + 2𝑎 112 𝑃 1 (𝑃 2 4 + 2𝑃 1 2 𝑃 2 2 + 𝑃 3 4 + 2𝑃 1 2 𝑃 3 2 ) + 2𝑎 123 𝑃 1 𝑃 2 2 𝑃 3 2 -𝑔 11	𝜕 2 𝑃 1 𝜕𝑥 1 2 -
											𝑔 44 (	𝜕 2 𝑃 1 𝜕𝑥 2 2 +	𝜕 2 𝑃 1 𝜕𝑥 3 2 ) = -(𝐹 11 -𝐹 44 -𝐹 12 )	𝜕𝜎 11 𝜕х 1	-𝐹 12	𝜕 𝜕х 1	𝑇𝑟𝜎 ̂+ 𝐸 1 (E.1a)
		𝛤	𝜕𝑃 2 𝜕𝑡	+ 2𝑃 2 (𝑎 1 -𝑄 12 (𝜎 11 + 𝜎 33 ) -𝑄 11 𝜎 22 ) -𝑄 44 (𝜎 12 𝑃 1 + 𝜎 23 𝑃 3 ) + 4𝑎 11 𝑃 2 3 + 2𝑎 12 𝑃 2 (𝑃 1 2 +
							𝑃 3 2 ) + 6𝑎 111 𝑃 2 5 + 2𝑎 112 𝑃 2 (𝑃 1 4 + 2𝑃 2 2 𝑃 1 2 + 𝑃 3 4 + 2𝑃 2 2 𝑃 3 2 ) + 2𝑎 123 𝑃 2 𝑃 1 2 𝑃 3 2 -𝑔 11	𝜕 2 𝑃 2 𝜕𝑥 2 2 -
											𝑔 44 ( 𝜕 2 𝑃 2 𝜕𝑥 1 2 +	𝜕 2 𝑃 2 𝜕𝑥 3 2 ) = -(𝐹 11 -𝐹 44 -𝐹 12 )	𝜕𝜎 22 𝜕𝑥 2	-𝐹 12	𝜕 𝜕х 2	𝑇𝑟𝜎 ̂+ 𝐸 2 (E.1b)
		𝛤	𝜕𝑃 3 𝜕𝑡	+ 2𝑃 3 (𝑎 1 -𝑄 12 (𝜎 11 + 𝜎 22 ) -𝑄 11 𝜎 33 ) -𝑄 44 (𝜎 13 𝑃 1 + 𝜎 23 𝑃 2 ) + 4𝑎 11 𝑃 3 3 + 2𝑎 12 𝑃 3 (𝑃 1 2 +
							𝑃 2 2 ) + 6𝑎 111 𝑃 3 5 + 2𝑎 112 𝑃 3 (𝑃 1 4 + 2𝑃 3 2 𝑃 1 2 + 𝑃 2 4 + 2𝑃 2 2 𝑃 3 2 ) + 2𝑎 123 𝑃 3 𝑃 1 2 𝑃 2 2 -𝑔 11	𝜕 2 𝑃 3 𝜕𝑥 3 2 -
											𝑔 44 ( 𝜕 2 𝑃 3 𝜕𝑥 1 2 +	𝜕 2 𝑃 3 𝜕𝑥 2 2 ) = -(𝐹 11 -𝐹 44 -𝐹 12 )	𝜕𝜎 33 𝜕х 3	-𝐹 12	𝜕 𝜕х 3	𝑇𝑟𝜎 ̂+ 𝐸 3 (E.1c)
	The derivation is straightforward: -𝐹 11	𝜕𝜎 33 𝜕х 3	-𝐹 12 ( 𝜕𝜎 11 𝜕х 3	+	𝜕𝜎 22 𝜕х 3	) -𝐹 44 (	𝜕𝜎 13 𝜕𝑥 1	+	𝜕𝜎 23 𝜕𝑥 2	) = -(𝐹 11 -
	𝐹 44 )	𝜕𝜎 33 𝜕х 3	-𝐹 12 ( 𝜕𝜎 11 𝜕х 3	+	𝜕𝜎 22 𝜕х 3	) = -(𝐹 11 -𝐹 44 -𝐹 12 )	𝜕𝜎 33 𝜕х 3	-𝐹 12	𝜕 𝜕х 3	𝑇𝑟𝜎 ̂, where 𝑇𝑟𝜎 ̂= 𝜎 11 + 𝜎 22 +
	𝜎 33 and	𝜕𝜎 13 𝜕𝑥 1	+	𝜕𝜎 23 𝜕𝑥 2	≡ -	𝜕𝜎 33 𝜕х 3	.
		Elastic stresses existing in the system can be found from Eq.(A.6) as:
			𝜎 11 = 𝑐 11 𝑢 11 + 𝑐 12 (𝑢 22 + 𝑢 33 ) -𝑞 11 𝑃 1 2 -𝑞 12 (𝑃 2 2 + 𝑃 3 2 ) -𝑓 11	𝜕𝑃 1 𝜕𝑥 1	-𝑓 12 ( 𝜕𝑥 2 𝜕𝑃 2	+	𝜕𝑃 3 𝜕𝑥 3	), (E.2a)
			𝜎 22 = 𝑐 11 𝑢 22 + 𝑐 12 (𝑢 11 + 𝑢 33 ) -𝑞 11 𝑃 2 2 -𝑞 12 (𝑃 1 2 + 𝑃 3 2 ) -𝑓 11	𝜕𝑃 2 𝜕𝑥 2	-𝑓 12 ( 𝜕𝑥 1 𝜕𝑃 1	+	𝜕𝑃 3 𝜕𝑥 3	), (E.2b)
		𝜎 33 = 𝑐 11 𝑢 33 + 𝑐 12 (𝑢 11 + 𝑢 22 ) -𝑞 11 𝑃 3 2 -𝑞 12 (𝑃 2 2 + 𝑃 1 2 ) -𝑓 11	𝜕𝑃 3 𝜕𝑥 3	-𝑓 12 ( 𝜕𝑥 1 𝜕𝑃 1	+	𝜕𝑃 2 𝜕𝑥 2	) (E.2c)

  𝑄 44 (𝜎 13 𝑃 1 + 𝜎 23 𝑃 2 ) -(𝐹 11 -𝐹 44 -𝐹 12 )(𝑐 11 -𝑐 12 )It is easy to check that in the case 𝑑𝑖𝑣𝑃 ⃗ = 0 for arbitrary 𝑝(𝜌). For a 2D-vortex with an empty core we can assume that 𝑝(𝜌) ≈ 𝑡𝑎𝑛ℎ ( ) as a zero approximation we obtain from Eq.(E.2) that 𝜎 13 = 𝑐[START_REF] Dzyaloshinsky | A thermodynamic theory of "weak" ferromagnetism of antiferromagnetics[END_REF] 𝑢 13 , 𝜎 23 = 𝑐 44 𝑢 23 and 𝜎 33 = (𝑐 11 -𝑐 12 )𝑢 33 + 𝑐 12 Tr𝑢 ̂-𝑞 12 𝑝 2 𝑃 3 ≈ 𝑄 44 𝑐 44 𝑝(𝑢 13 𝑠𝑖𝑛𝜑-𝑢 23 𝑐𝑜𝑠𝜑)-(𝐹 11 -𝐹 44 -𝐹 12 )(𝑐 11 -𝑐 12 )Next, we consider several model cases. The first case is an absolutely rigid shell covering the ferroelectric core, when the maximal stresses evolved in the system can be roughly estimated from Eqs.(A.8) at 𝑢 𝑖𝑗 = 0. The second case corresponds to a very soft shell covering the rigid ferroelectric core, when the stresses are minimal and the strains (at least near the cylinder ends) can be roughly estimated from Eqs.(E.2) at 𝜎 𝑖𝑗 = 0. Furthermore, the compatibility conditions should be valid.Using the relations 𝑞 11 -𝑞 12 = (𝑄 11 -𝑄 12 )(𝑐 11 -𝑐 12 ) = The functions 𝑢 𝑖𝑗 (𝜌, 𝜑, 𝑧) are elastic strains, 𝑠 𝑖𝑗 are elastic compliances, 𝑄 𝑖𝑗 are electrostriction tensor components, 𝑔 𝑖𝑗 are polarization gradient coefficients written in Voight notations. From Eq.(E.7b) the axial part of the flexon polarization is proportional to 𝜎 𝑖𝑗 contains a contribution proportional to 𝑄 𝑖𝑗𝑘𝑙 𝑃 𝑘 𝑃 𝑘 , which originates from the electrostriction coupling. For the 𝑃 3 -component we determine that the term

						𝑃 𝑟 = 0, 𝑃 𝜑 = -𝑝(𝑟) and 𝑃 3 = 0. 𝐺 𝑓𝑙𝑒𝑥𝑜 ≈ 𝐹 11 -𝐹 44 -𝐹 12 2 (𝜎 𝑖𝑖 𝜕𝑃 𝑖 𝜕x 𝑖 -𝑃 𝑖 𝜕𝜎 𝑖𝑖 𝜕x 𝑖	).	(E.6b) (E.8b)
	Elementary, but cumbersome calculations lead to an "odd" flexo-field in the boundary conditions
	∫ 𝐺 𝑓𝑙𝑒𝑥𝑜 . 𝑉	≈ ∫ 𝐺 0 . 𝑉	+	𝐹 11 -𝐹 44 -𝐹 12 2	𝜌 𝐿 𝐶 𝑥 ), where 𝐿 𝐶 𝑥 is a transverse correlation length. 𝑄 11 -𝑄 12 3 ∫ 𝑃 𝑖 3 . 𝑆𝑖 .
	To verify the analogy with a Dyzaloshinskii-Moryia interaction (DMI), we attempt to convert the Lifshitz invariants into chiral interactions of a DMI type. Using 𝑑𝑖𝑣𝑃 ⃗ = 0 for the uncharged domain structures and making straightforward analytic manipulations for the cubic m3m point symmetry group Using Eqs.(E.6𝜕𝑢 33 𝜕𝑥 3 2[𝑎 1 -(𝑞 11 +2𝑞 12 )𝑝 2 -(𝑄 11 -𝑄 12 )(𝑐 11 -𝑐 12 )𝑢 33 +(𝑔 11 +(𝐹 11 -𝐹 44 -𝐹 12 )(𝑓 11 -𝑓 12 ))𝐿 𝐶 𝑧 +𝑔 44 𝐿 𝐶 , (E.7a) of the BaTiO3 parent phase, the explicit form of the Lifshitz invariant (1f) is 𝐺 𝑓𝑙𝑒𝑥𝑜 ≈ 𝑥 ] 𝐹 11 -𝐹 44 -𝐹 12 2 (𝜎 𝑖𝑖 𝜕𝑃 𝑖 𝜕𝑥 𝑖 -𝑃 𝑖 𝜕𝜎 𝑖𝑖 𝜕𝑥 𝑖 ), where a summation over "i" is performed [Eqs.(E.8)]. The elastic stress
						𝐹 11 -𝐹 44 -𝐹 12 2	𝑄 12 ((𝑃 1 2 + 𝑃 2 2 )	𝜕𝑃 3 𝜕𝑥 3	-
	𝑃 3	𝜕(𝑃 1 2 +𝑃 2 2 ) 𝜕𝑥 3
						𝑄 11 -𝑄 12 𝑠 11 -𝑠 12	, 𝑞 11 + 2𝑞 12 = (𝑄 11 +
	3 2𝑄 12 )(𝑐 11 + 2𝑐 12 ) = = 𝜕 𝜕х 3 [(𝑐 11 + 2𝑐 12 )𝑇𝑟𝑢 ̂-(𝑞 11 + 2𝑞 12 )𝑃 2 -(𝑓 11 + 2𝑓 12 )𝑑𝑖𝑣𝑃 ⃗ ] ≈ (𝑐 11 + 2𝑐 12 ) 𝑄 11 +2𝑄 12 𝑠 11 +2𝑠 12 , 𝑐 11 -𝑐 12 = 1 𝑠 11 -𝑠 12 , 𝑐 44 = 1 𝑠 44 and 𝑓 11 -𝑓 12 = (𝐹 11 -𝐹 12 )(𝑐 11 -𝜕𝑇𝑟𝑢 ∂х (E.3) 3 and 𝑐 12 ) = 𝐹 11 -𝐹 12 𝑠 11 -𝑠 12 , we obtain that
	𝜕𝜎 33 𝜕х 3 Assuming that = 𝜕 𝜕х 3 [𝑐 11 𝑢 33 + 𝑐 12 (𝑢 11 + 𝑢 22 ) + (𝑞 11 -𝑞 12 )𝑃 3 2 -𝑞 12 𝑃 2 -(𝑓 11 -𝑓 12 ) 𝜕 𝜕х 3 [(𝑐 11 -𝑐 12 )𝑢 33 + 𝑐 12 Tr𝑢 ̂+ (𝑞 11 -𝑞 12 )𝑃 3 2 ] -(𝑓 11 -𝑓 12 ) 𝜕𝑃 3 𝜕𝑥 3 -𝑓 12 𝑑𝑖𝑣𝑃 ⃗ ] ≈ 𝑃 3 ≈ 𝑄 44 𝑠 44 𝑝(𝑢 13 𝑠𝑖𝑛𝜑-𝑢 23 𝑐𝑜𝑠𝜑)-𝐹 11 -𝐹 44 -𝐹 12 𝑠 11 -𝑠 12 𝜕 𝜕z 𝑢 33 2[𝑎 1 -𝑄 11 +2𝑄 12 𝑠 11 +2𝑠 12 𝑝 2 -𝑄 11 -𝑄 12 𝑠 11 -𝑠 12 𝑢 33 +[𝑔 11 +(𝐹 11 -𝐹 44 -𝐹 12 ) 𝐹 11 -𝐹 12 ]𝐿 𝐶 𝑧 +𝑔 44 𝐿 𝐶 , (E.7b) 𝑥 ] 𝑠 11 -𝑠 12 𝜕 2 𝑃 3 𝜕𝑥 3 2 (E.4) 𝜕𝑇𝑟𝑢 where {𝜌, 𝜑, 𝑧} are cylindrical coordinates, the function 𝑝(𝜌) ≈ 𝑡𝑎𝑛ℎ ( 𝜌 𝑥 ), and 𝐿 𝐶 𝑥 and are lateral and 𝐿 𝐶 ∂х 3 ≈ 0 and using the smallness of 𝑃 3 , we can obtain the approximate linearized equation for 𝑃 3 : 2𝑎 1 * 𝑃 3 -𝑔 11 * 𝜕 2 𝑃 3 𝜕𝑥 3 2 -𝑔 44 ( 𝜕 2 𝑃 3 𝜕𝑥 1 2 + 𝜕 2 𝑃 3 𝜕𝑥 2 2 ) ≈ 𝜕𝑢 33 𝜕𝑥 3 + (E.5) axial correlation lengths. -𝐹 11 -𝐹 44 -𝐹 12 𝑠 11 -𝑠 12 𝜕 𝑢 33 (𝜌, 𝜑, 𝑧), and this proportionality along with Fig. S9o qualitatively describes the 𝜕z 𝐸 3 , curves' behavior in Fig. 3f and 4f in the main text.
	where	𝑎 1 * = 𝑎 1 -𝑄 12 𝑇𝑟𝜎 ̂-(𝑄 11 -𝑄 12 )𝜎 33 ≈ 𝑎 1 -𝑄 12 (𝑐 11 + 2𝑐 12 )𝑇𝑟𝑢 ̂-(𝑞 11 + 2𝑞 12 )𝑝 2 -Using Eqs.(E.1), we can "recover" an analog of Eq.(A.2) for a "full" Lifshitz invariant:
	(𝑄 11 -𝑄 12 )[(𝑐 11 -𝑐 12 )𝑢 33 + 𝑐 12 Tr𝑢 ̂] ≈ 𝑎 1 -(𝑞 11 + 2𝑞 12 )𝑝 2 -(𝑄 11 -𝑄 12 )(𝑐 11 -𝑐 12 )𝑢 33 𝑔 11 𝐺 𝑓𝑙𝑒𝑥𝑜 = 𝐹 11 -𝐹 44 -𝐹 12 2 (𝜎 𝑖𝑖 𝜕𝑃 𝑖 𝜕x 𝑖 -𝑃 𝑖 𝜕𝜎 𝑖𝑖 𝜕x 𝑖 ) + 𝐹 12 2 [𝑇𝑟(𝜎 ̂) 𝜕𝑃 𝑖 𝜕х 𝑖 -𝑃 𝑖 𝜕𝑇𝑟(𝜎 ̂) ], 𝜕х 𝑖 * = 𝑔 11 + (𝐹 11 -𝐹 44 -𝐹 12 )(𝑓 11 -𝑓 12 ). The simplest 2D polarization vortex can be modeled by the and (E.8a) functions 𝜕𝑃 𝑖 Assuming that 𝜕х 𝑖
						𝑃 1 = 𝑝(𝑟)𝑠𝑖𝑛𝜑, 𝑃 2 = -𝑝(𝑟)𝑐𝑜𝑠𝜑 and 𝑃 3 = 0.	(E.6a)
	These polarization components in cylindrical coordinates are

≈ 0 (as it should be for the mostly uncharged domain structures) and making an integration over parts in the second term in Eq.(E.8a), we obtain:

Table SII .

 SII Fitting parameters for Eqs.(F.4)

	Fitting				Flexoelectric coefficients 𝐹 𝑖𝑗
	parameter	-6𝐹 𝑖𝑗	-2𝐹 𝑖𝑗	-𝐹 𝑖𝑗	0	𝐹 𝑖𝑗	2𝐹 𝑖𝑗	6𝐹 𝑖𝑗
	𝑓 (μC/cm 2 )	-0.011	-0.0071	-0.0045	0	0.0045	0.0071	0.011
	𝑔 (μC/cm 2 )	0.0305	0.0224	0.0174	N/A	0.0174	0.0224	0.0305
	𝐴 (nm -2 )	0.0100	0.0125	0.0110	N/A	0.0110	0.0125	0.0100
	𝐵 (nm 2 )	5	10	15	N/A	15	10	5
	z 0 (nm)	1.8	2.2	2.3	N/A	2.3	2.2	1.8
	z 𝑚 (nm)	17.8	17.5	17.5	N/A	17.5	17.5	17.8

Acknowledgements

. A.N.M. acknowledges EOARD project 9IOE063 and related STCU partner project P751. R.H. and S.C.-H. acknowledge funding from the French National Research Agency through contract ANR-18-CE92-0052 "TOPELEC". V.Y.R. acknowledges the support of COST Action CA17139. A portion of FEM was conducted at the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility (CNMS Proposal ID: CNMS2021-B-00843).

Supplementary Materials to "Chiral Polarization Textures Induced by the Flexoelectric Effect in Ferroelectric

Nanocylinders" APPENDIX A. Mathematical Formulation of the Problem and FEM Details

A. Mathematical Formulation of the Problem

We use the Landau-Ginzburg-Devonshire (LGD) approach combined with electrostatic equations, because this method has proven to be successful in establishing the physical origin of anomalies in phase diagrams, determining polar and dielectric properties of ferroelectric nanoparticles [1,2], and calculating the changes of their domain structure morphology with size reduction [3,4]. The LGD approach allows for the consideration of various size and surface effects, such as correlation effects and depolarization fields arising in the case of incomplete polarization screening [5], surface bond contraction [6,7], and intrinsic surface stresses and strains [8,9,10].

The LGD free energy functional G additively includes a Landau expansion on powers of 2-4-6 of the polarization (𝑃 𝑖 ), 𝐺 𝐿𝑎𝑛𝑑𝑎𝑢 ; a polarization gradient energy contribution, 𝐺 𝑔𝑟𝑎𝑑 ; an electrostatic contribution, 𝐺 𝑒𝑙 ; the elastic, electrostriction, flexoelectric contributions, 𝐺 𝑒𝑠+𝑓𝑙𝑒𝑥𝑜 ; and a surface energy, 𝐺 𝑆 . It has the form [3,10,11]: )) (A.1e)

Here 𝑉 𝑐 and 𝑉 𝑠 are the core and shell volume, respectively. The coefficient 𝑎 𝑖 linearly depends on temperature T, 𝑎 𝑖 (𝑇) = 𝛼 𝑇 [𝑇 -𝑇 𝐶 ], where 𝛼 𝑇 is the inverse Curie-Weiss constant and 𝑇 𝐶 is the ferroelectric Curie temperature renormalized by electrostriction and surface tension. Tensor components 𝑎 𝑖𝑗 are regarded as temperature-independent. The tensor 𝑎 𝑖𝑗 is positively defined if the ferroelectric material undergoes a second order transition to the paraelectric phase and negative otherwise. The higher nonlinear tensor 𝑎 𝑖𝑗𝑘 and the gradient coefficients tensor 𝑔 𝑖𝑗𝑘𝑙 are positively defined and regarded as temperature-independent. The following designations are used in Eq.(A.1e):

cylinder axis and is nearly independent of the shell dielectric and elastic properties; therefore, role of the shell for a prolate multiaxial cylindrical core is much weaker than it is for a uniaxial spherical core.

APPENDIX B. Influence of the Flexoelectric Coupling Strength and Anisotropy on the

Domain Structure of a Cylindrical Nanoparticle

The distribution of the polarization component 𝑃 3 in three different cross-sections of the nanoparticle core, polarization magnitude 𝑃 in the XZ cross-section of the nanoparticle, and isosurfaces of polarization components 𝑃 1,2,3 are shown in Fig. S5. The images are calculated without electrostriction (𝑄 𝑖𝑗 = 0) and flexoelectric (𝐹 𝑖𝑗 = 0) coupling (Fig. S5a); with electrostriction coupling (𝑄 𝑖𝑗 ≠ 0) and negative (Fig. S5b), zero (Fig. S5c), positive (Fig. S5d), or high positive (Fig. S5e) values of the flexoelectric coefficients 𝐹 𝑖𝑗 .

The extremal (maximal or minimal) value 𝑃 𝑒 in the center of the diffuse axial P3-domain, shown by the green curve in Fig. 3f (main text), frequently differs from 𝑃 𝑚𝑎𝑥 and 𝑃 𝑚𝑖𝑛 values due to the presence of the small sixteen 𝑃 3 -domains localized near the top and bottom junction of cylindrical sidewall with the ends (bottom raw in Fig. S5). The polarization direction alternates in these 180degree domains, with four "up" and four "down" domains at each junction. These alternating 𝑃 3domains are not the part of flexon, because they are induced by the electrostriction coupling. Their structure depends weakly on the flexoelectric coupling strength 𝑓, if 𝑓 is small. However, the shape of each domain significantly changes and the size of each domain moderately increases as |𝐹 𝑖𝑗 | increases (compare Fig. S5d and Fig. S5e). The change of polarization direction occurs in each of the 16 alternating 𝑃 3 -domains at 𝐹 𝑖𝑗 → -𝐹 𝑖𝑗 (compare Fig. S5b and Fig. S5d).

In contrast to the topological index, the toroidal moment, 𝑀 ⃗⃗ = , where the magnitude 𝑝(𝜌, 𝑧) is nearly independent on the flexoelectric coupling (e.g., Fig. S5, the middle row).