{'fraction_non_alphanumeric': 0.037795201652659746, 'fraction_numerical': 0.02089299967134607, 'mean_word_length': 5.023755656108597, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Labyrinthine domains in ferroelectric nanoparticles: a manifestation of gradient-driven morphological transition

In the framework of the Landau-Ginzburg-Devonshire (LGD) approach we studied finite size effects of the phase diagram and domain structure evolution in spherical nanoparticles of uniaxial ferroelectric. The particle surface is covered by a layer of screening charge characterized by finite screening length. The phase diagram, calculated in coordinates"particle radius vs. screening length"has a wide region of versatile poly-domain structures separating single-domain ferroelectric and nonpolar paraelectric phases. Unexpectedly, we revealed a region of stable irregular labyrinthine domains in the nanoparticles of uniaxial ferroelectric CuInP2S6 with the first order paraelectric-ferroelectric phase transition. We established that the origin of labyrinthine domains is the mutual balance of LGD, polarization gradient and electrostatic energies. The branching of the domain walls appears and increases rapidly when the polarization gradient energy decreases below the critical value. Allowing for the generality of LGD approach, we expect that the gradient-induced morphological transition can be the source of labyrinthine domains appearance in many spatially-confined ferroics with long-range order parameter, including relaxors, ferromagnetics, antiferrodistortive materials and materials with incommensurate ferroic phases.


I. INTRODUCTION
The ferroic materials described by Landau theory of symmetry-breaking have a substantial impact on fundamental science and various applications. Different types of topological defects in different ferroics (ferromagnets, ferroelecrics, ferroelastics) are even more numerous and enigmatic than different types of symmetry-breaking, and consequently, they become one of the key fundamental problems and hot topics in scientific community [1,2].
Complementary to the topological point defects [1], domain walls can be considered as extended 2D topological defects in ferroics (see e.g. chapter 8 in [2] and refs. therein). Vortices and vertices composed by the closure of four domain walls have been observed experimentally and described theoretically in a bulk and nanosized ferroelectrics [3,4,5,6,7]. Stable surface-induced labyrinthine domain structures were observed by Piezoresponse Force Microscopy (PFM) in ergodic ferroelectrics relaxors and explained by the presence of higher-order term in free-energy expansion that gives rise to the polarization modulations [8]. Fractal domain structures are sometimes observed in multiferroic thin films [9] and near the surface of relaxors close to relaxor-ferroelectric boundary [10], but the labyrinthine domains with a single characteristic length scale were observed by PFM in ergodic relaxors only [8]. These labyrinthine structures can coexist with classical ferroelectric domains closer to ferroelectric composition limit [11,12]. The labyrinthine domain structure was predicted theoretically in thin films of incommensurate ferroelectrics [13] and bi-layered ferroelectrics [14], being similar to those observed in ultrathin magnetic films [15].
However, we did not find any experimental observation or theoretical prediction of labyrinthine domains in the nanoparticles of ordered ferroelectrics, which intriguing polar and dielectric properties attract permanent attention of researchers. Classical examples are unexpected experimental results of Yadlovker and Berger [16,17,18], which reveal the enhancement of polar properties of cylindrical nanoparticles of Rochelle salt. Frey and Payne [19], Zhao et al [20], Drobnich et al [21], Erdem et al [22] and Golovina et al [23,24,25] demonstrate the possibility to control the phase transitions (including new polar phases appearance) for BaTiO3, S2P2S6, PbTiO3 and KTa1-хNbхO3 nanopowders and nanoceramics by finite size effects.
The theory of finite size effects in nanoparticles allows one to establish the physical origin of the polar and dielectric properties anomalies, and phase diagrams changes appeared under the nanoparticles sizes decrease. In particular, using the continual phenomenological approach Niepce [26], Huang et al [27,28], Ma [29], Eliseev et al [30] and Morozovska et al [31,32,33] have shown, that the changes of the transition temperatures, the enhancement or weakening of polar properties in a single-domain spherical and cylindrical nanoparticles are conditioned by the various physical mechanisms, such as correlation effect, depolarization field, flexoelectricity, electrostriction, surface tension and Vegard-type chemical pressure.
Notably depolarization field always decreases ferroelectric polarization and transition temperature, especially under the presence of imperfect screening [34,35,36]. For majority of models the particles were regarded covered with perfect electrodes and so their single-domain state would be stable. Only few models describing the imperfect screening effect in nanoparticles have been evolved [34][35][36].
To fill the gap in the knowledge, below we analyze the phase diagram and domain structure evolution in spherical nanoparticles of uniaxial ferroelectric CuInP2S6 (CIPS). We regarded that the particle surface is covered by a layer of screening charge characterized by finite screening length.
The imperfect screening and finite size effects are studied using the Landau-Ginzburg-Devonshire (LGD) approach combined with the electrostatic equations. We revealed that a regular stripe domain structure avalanche-like transforms into a labyrinth pattern with a gradient term decrease below the critical value and classified the event as a gradient-induced morphological transition.
The applicability of LGD approach for thin films and nanoparticles with radius less than (2-10) nm is corroborated by the fact, that the critical sizes of the long-range order appearance and properties calculated from atomistic [37 ,38, 39, 40, 41] and phenomenological [30-33, 42, 43] theories are in a good agreement with each other as well as with experimental results for nanosized ferromagnetics [44] and ferroelectrics [16-20, 22, 45]. Both atomistic simulations and LGDdescription are absent for CIPS nanoparticles.

II. THEORETICAL APPROACH
Let us consider a CIPS nanoparticle of radius R with a one-component ferroelectric polarization ( ) directed along the crystallographic axis 3 [ Fig.1(a)]. The particles are covered by a layer of screening charge with a surface charge density σ characterized by a nonzero screening length λ. The specific nature of the surface charge can be, e.g., Bardeen-type surface states [46]. For the case the screening charges can be localized at surface states caused by the strong band-bending via depolarization field [47,48,49,50,51], at that λ can be much smaller (≤0.1 nm) than a lattice constant (~0.5 nm) [34]. Concrete expression for λ can be derived in, e.g., Stephenson-Highland ionic adsorption model [52,53,54], by the linearization of σ , as  CIPS parameters are listed in Table I.
For a layered perovskite with layers plane (001) and ferroelectric dipoles directed in the outof-plane direction, we can assume that the dependence of the in-plane components of electric polarization on the inner field electric Ei is linear (see e.g. [30,35,56]): (1d) Here i E are electric field components related with electric potential ϕ as The coefficient β is temperature-independent and negative, since CIPS undergoes the first order transition to paraelectric phase. Coefficient γ and gradient coefficients g11 and g44 are positive and temperature independent. An isotropic approximation, 55 44 g g ≈ in (001) plane was taken for monoclinic CIPS structure. ij σ is the stress tensor in Eq.(1d). We omit the evident form of the flexo es G + for the sake of simplicity, it is listed in Refs. [57,58,59]. Since the values of the electrostriction and flexoelectric tensor components, ijkl Q and ijkl F , are unknown for CIPS, we performed numerical calculations with the coefficients varied in a physically reasonable range . Results proved the insignificant impact of electrostriction and flexoelectric coupling on domain morphology [60]. Other LGD parameters for a bulk ferroelectric CIPS were taken from Ref. [61] and are listed in Table I. Table I. LGD parameters for bulk ferroelectric CuInP2S6, used in calculations Corresponding Euler-Lagrange equation for ( ) The boundary condition for P at the spherical surface is natural, , n is the outer normal to the surface. The potential ϕ satisfies a Poisson equation inside the particle, and Laplace equation outside it, 3-D Laplace operator is denoted by the symbol ∆ . Equations (3) are supplemented by the condition of potential continuity at the particle surface, ( )

A. Main features of phase diagrams at different temperatures
Phase diagrams was studied at different temperatures in the range (300 -200) K.
Unfortunately we do not know the temperature dependence of λ, and so we perform all calculations regarding λ changing in the range (10 -3 -1) nm. Phase diagram of CIPS nanoparticles calculated at T=293 K and 200 K in coordinates "radius R -screening length λ" is shown in Fig. 2(a) and 2(b), respectively.
At room temperature the phase diagram has an unexpectedly wide region of stable poly-domain states (PDFE) separating single-domain ferroelectric (SDFE) and nonpolar paraelectric (PE) phases [see Fig. 2(a)]. The bottom row shows the typical changes of polarization distribution in the equatorial cross-section of the nanoparticle with R=5 nm, which happens with increase of λ. A single-domain state is stable at very small λ<0.01 nm, two-domain structure (electric quadrupole) is stable in the interval 0.01<λ<0.017 nm, three-domain structure (electric octupole) exist at 0.017<λ<0.019 nm, 2N-multipolar domain stripes are stable at 0.02<λ<0.035 nm. Coexistence of PDFE and PE phase when the nanoparticle consists of PE surface layer and ferroelectric domain stripes in the core appears at 0.035<λ<0.045 nm, and is followed by the size-induced phase transition into a stable PE at λ>0.045 nm. Unexpectedly, we revealed a region of stable "labyrinthine" domains (SLD) of irregular shape (yellow circles) inside the region of regular domain structures with quadruple two (purple circles) or multiple (magenta circles) domain stripes. SLD region is within a dashed parallelogram in Fig.2(a).
With the temperature decrease from 293 K (that is very close to the CIPS Curie temperature 302 K) to 200 K the region of a SDFE significantly increases towards smaller radii R (up to the very small R=1 nm for which LGD applicability becomes questionable) and higher λ (from e.g. λ=0.003 nm at 293 K to 0.02 nm at 200 K) [compare the size of SDFE regions in Fig.2(a) and 2(b)]. The wide region of PE phase (present at 293 K) almost disappears with the temperature decreasing to 200 K leading to the conclusion that PDFE state can be stable in ultra-small CIPS nanoparticles (with radius less than 2 nm) covered by a screening charge [compare the size of PE regions in Fig.2(a) and 2(b)].
The shift and increase of stable labyrinthine domains region(s) are evident with the temperature decrease from 300 to 200 K [compare the size and positions of SLD regions in Fig.2(a) and 2(b)].
The increase of SDFE, PDFE and SLD regions with the temperature T decrease stem from the well-   Table I. The effect of geometric catastrophe can be imagined from the images of SLD in the nanoparticles of radius 4, 3 and 2 nm for which the number of branches and sharp bendings of domain walls gradually decreases with the particle radius decrease from 4 nm to 2 nm [see right column in Fig. 2(b)]. Eventually SLD disappears for R=1 nm. Hence the effect of geometric catastrophe suppresses the SLD in small particles.
Note the validity of our prediction regarding SLD appearance and PDFE state conservation for nanoparticles of sizes more than 2R = 4 nm, because they corresponds to 10 lattice constants or more. It is general opinion that LGD approach can be valid only qualitatively for the sizes less than 10 l.c. [16-20, 22, 30-33], and must be approved by ab initio calculations.

B. Labyrinthine domains stability and evolution
Polarization distributions in the equatorial cross-section of the nanoparticle with radius   Table I.
An example of labyrinthine domains evolution with increase of the gradient coefficient g44 is illustrated in Fig.4(a)-(e). Stable labyrinths exist at g44 less than the critical value gcr≈2.5×10 -11 m 3 /F [ Fig. 4 Fig. 4(f)], which in turn disappears with g44 increase more than 2.8×10 -11 m 3 /F [ Fig. 4(g)]. Complementary we made sampling over 5 -20 different labyrinthine domain patterns for each g44, R and λ values, which emerged with computation time from different initial random distributions of polarization inside the particle (see Fig. S1 in Ref. [60]). From Fig.S1 we concluded the branching number Σ, defined as the total number of branched domain walls, dangling branches, separated stripes and loops, decreases sharply with g44 increase. Σ varies slightly for different samples far from the critical value gcr (top lines in Fig.S1), but the variation becomes bigger near the critical value gcr (bottom lines in Fig.S1).

Examples of Σ calculation are shown in Figs.5. A color image of complex labyrinthine pattern
with dangling branches, branch seedings, separated island and separated curved stripe is shown in Fig.5(a). Fig.5(b) shows the black domains with white walls corresponding to the structure (a).
Graphs (c, d) with numbered features have been drawn allowing for the connectivity between different domains and particle surface in plot 5(b). Proposed algorithm of Σ calculation counts all branching points, dangling branches and separated stripes ends, which do not cross the particle from one surface to another one. Meanwhile the straight or slightly curved stripes (even very small) that cross the particle from one surface to another one do not contribute to Σ. However, the algorithm is not ideal, because the criteria distinguishing "slightly curved" and "strongly curved" stripes are somehow voluntary. Actually, for some complex cases, like the one shown in Fig.5(a), visual recalculation of Σ lead to different results for peculiarities number corresponding to "red" (Σr=7) and "blue" (Σb=9) domains [compare Fig. 5

(b) and 5(c)].
To improve the situation we operate with the values averaged for "red" and "blue" domains, e.g. Σ =8 corresponds to Fig.5(a).  3 and 4), which do not cross the particle from one surface to another one, contribute to Σ number. Other four slightly curved stripes, which both ends at particle surfaces are marked with asterisks "*", do not contribute to Σ number.  Table I.   The sampling-averaged value Σ is not integer for fixed g44. From Fig.4(h), the dependence of Σ on the gradient term 44 g is described by the function , and so it continuously appears at cr g g = 44 . Hence we can associate Σ appearance with a rapid change of the domain walls connectivity.
We leave for further studies answer on the question how the threshold of labyrinthine domain appearance and the ranges of their stability at phase diagram can be derived analytically.
The theoretical prediction of labyrinthine domains requires urgent verification by PFM, that is an ideal tool for 3D visualization of the domain structure with nanoscale resolution (see e.g. [8,64,65,66,67] and refs therein). We are convinced by a numerical calculation that qualitatively similar SLD can be realized in other incompletely screened uniaxial ferroelectric nanoparticles, such as BiFeO3, however we leave a discussion of these results for further studies.

IV. CONCLUSION
In the framework of LGD approach combined with the equations of electrostatics, we studied the finite size effects of the phase diagrams and domain structure in spherical ferroelectric nanoparticles covered by a layer of a screening charge with finite screening length. The phase diagrams, calculated in coordinates "particle radius -screening length", has a wide region of versatile poly-domain states separating single-domain ferroelectric and nonpolar paraelectric phases. Quite unexpectedly we revealed that a regular stripe domain structure sharply transforms into a labyrinth pattern with a gradient term decrease below the critical value and named the event as a gradient- where ij σ is the stress tensor, electrostriction and flexoelectric tensor components are ijkl Q and ijkl F , respectively. Also equation of mechanical equilibrium should be solved with elastic boundary conditions. Electric boundary conditions should be modified by flexoeffect as explained in details in Refs. [2][3][4]. To the best of our knowledge the values of the electrostriction and flexoelectric tensor components are unknown for CIPS.
Our previous results [2][3][4] for PbTiO3 and BaTiO3 thin films with all known coefficients and preliminary numerical calculations made for nanoparticles in this work with the coefficients varied in a physically reasonable range (fleaxoelectric tensor ≤ ijkl F 10 11 m 3 /C and electrostriction tensor ≤ ijkl Q 0.1 m 4 /C 2 ) proved that that the impact of electrostriction leads to no more that 1% strain of the particle as well as to slightly visible decrease of the domain wall width. Hence the change of the nanoparticle spherical shape to ellipsoidal one (it indeed reveals the tiny elongation or contraction depending on the electrostriction coefficients sign in polar direction z) leads to rather small changes of electric field, which do not affect significantly on the morphology of domain structure, both labyrinthine and regular ones.
The flexoelectric impact is really small, but can be non-trivial for uniaxial and multiaxial ferroelectrics, since it can charge the nominally uncharged Ising-type walls [ 5 ]. However for the case of CIPS we did not see any significant changes of the domain morphology or walls charging with the flexoefffect included.

S.2. Sampling over different labyrinthine domain patterns
We made sampling over 10 -20 different labyrinthine domain patterns, which emerged with computation time from different initial random distributions of polarization inside the particle (further named "initial seedings") for each g44, R, λ and T values. Six labyrinthine domain patterns corresponding to definite g44 values, which varies from 0.5×10 −11 m 3 /F (top row) to 2.5×10 −11 m 3 /F (bottom row), temperature T=300 K, particle radius R=10 nm, screening length λ=0.03 nm are shown in Fig. S1. The number of branches Σ of the labyrinthine domain structure is indicated below each pattern in the figure, and it is different for different "samples", while the width of the labyrinth pattern is independent on the sample (it depends on g44, R, λ and T and material parameters of CIPS). Σ varies slightly for different samples far from the critical value gcr (see top lines in Fig.S1), but the variation becomes bigger near the critical value gcr (see bottom lines in Fig.S1). Naturally, the average value of Σ calculated to the definite g44 is not integer. Other parameters are listed in Table I. The number of branches Σ of the labyrinthine domain structure is indicated below each figure. The algorithm of Σ calculation is illustrated in Fig. 5 and discussed in the main text.
Also we used different random seeding to generate regular domain stripes, and, as anticipated for the case, the results appeared independent of the initial seeding. We compared the energies for different samples, which appear close for different labyrinthine patterns (except for the immediate vicinity of gcr, where the fluctuations increase), and strongly differ from the energies of regular stripes. All points at the diagrams in Figs.2 correspond to minimal energy.

S.3. Clarifications concerning the charge of the domain walls
Earlier thin ferroelectric films with nominally uncharged 180-degree domain structure have been considered and it has been shown that in the region of the domain wall contact with the electricallyopen film surface the depolarization electric field appears and leads to the domain wall broadening near the surface [3,4,6]. The self-consistent calculations confirmed that broadening causes the decrease of depolarization field that in its turn leads to the decrease of the system electrostatic energy [3,4,6]. Notably, the broadened wall becomes curved and thus no more uncharged near the surface [3,4,6].
A similar effect was calculated by us for spherical particles and shown in the manuscript.
Namely the domain walls in both PDFE and SLD phases are uncharged in the central part of the particle (far from its poles) and become charged and broadened at the polar surfaces [see Fig. S3 and yellow-green regions near the poles in Fig.1(a)]. Far from the poles polarization vector (shown by white arrows in Fig.S5) is parallel to the wall planes and so the walls are uncharged. Inset shows the cross-section of the spatial region near the particle pole, where the polarization vector P in Inset shows the region near the particle pole, where polarization vector P in different domains (shown by red and blue arrows) are no more parallel to the curved surface of the broadened wall (shown by black curves).