Ultrafast polarization switching in ferroelectrics

A method of ultrafast switching of ferroelectric polarization is suggested. The method is based on the interaction of a ferroelectric sample with the feedback field of a resonator in which the sample is inserted. The polarization reversal time can be of order of femtoseconds. The polarization switching produces a coherent electromagnetic pulse.

To regulate the processing of such devices, it is often necessary to be able to quickly vary the direction and magnitude of ferroelectric polarization. There exist two ways of polarization switching that can be called inhomogeneous (or incoherent) and homogeneous (or coherent).
First, the inhomogeneous way of polarization switching has been studied, being realized through the nucleation and growth of domains of opposite polarization, with moving domain walls under the influence of static electric fields [9][10][11][12][13]. This way, however, provides rather long switching times, on the order of nanoseconds, being limited by the domain recrystallization time that is typically hundreds of picoseconds [1,2,11,14,15]. Similar slow switching in the nanoscale volume of a ferroelectric can be realized by mechanical deformation of a ferroelectric sample [16]. Another slow mechanism of polarization switching is due to the chemical oxidization at the surface of a ferroelectric film [17]. Because of the principal restriction of the inhomogeneous switching caused by the limited domain recrystallization time, it has been necessary to find other ways that could provide much faster switching.
The other way that has been developed relatively recently is the homogeneous (or coherent) switching process realized by external alternating fields, applied perpendicular to the ferroelectric polarization, in the optical [18][19][20][21], terahertz [22][23][24][25][26][27][28], or infrared [29] regions. Under this process, the alternating field acts directly on all ions of a single-domain sample, and polarization switching occurs through a continuous homogeneous mechanism, without formation of new domains of opposite polarization. Homogeneous switching is facilitated in films of nanometer thickness, where inhomogeneous nucleation is strongly suppressed [30].
In order to realize a homogeneous switching, it is necessary that, first, all characteristic times of the process be much shorter than the domain nucleation-growth time and, second, that the sample be a single-domain ferroelectric. Under the ultrafast switching by means of alternating fields, the first condition is easy to accomplish, since the domain nucleation-growth time is sufficiently long, being of order of nanoseconds. And the preparation of single-domain ferroelectrics is a technical problem having several solutions [31][32][33][34]. For instance, single-domain states can be made stable by using strain [34] or doping with point defects [35]. Also, there are plenty of ferroelectric films of nanometer thickness, where domain nucleation is suppressed [30].
The homogeneous switching mechanism under the action of an alternating field, involving no nucleation and growth of oppositely polarized domains, can provide reversal times of order of picoseconds. However the present sources do not provide the strength of a pulse sufficient for completely switching the polarization. Experiments [36] have shown that the reversal can be only 40% of its equilibrium value. Although the reversal is quite fast, occurring in about 10 −13 s, but the reversed polarization rapidly, during the same 10 −13 s, returns to the initial state, similarly to the dynamics induced by terahertz pulses, when the reversal happens over a picosecond time scale, followed by its fast complete retrieval [37].
Finally, time-dependent density functional theory simulations show that by strongly exciting electrons via laser pulses it could be possible to change the underlying dynamical potential energy surface, which could result in the polarization switching within tens of picoseconds [38].
In the present paper, we consider the homogeneous way of polarization switching involving alternating fields, but with a rather different setup. The idea of the method is to put a ferroelectric, subject to an external constant electric field, into a resonator cavity. Then the polarization motion induces a resonator feedback field acting back on the ferroelectric. In such a way, there is no need in additionally imposed external electromagnetic pulses, but the ferroelectric produces the required pulse by itself through the feedback field. The polarization switching can be realized in femtoseconds. The suggested method of switching also uses alternating fields, similar to the techniques employing oscillating fields with fixed properties. However the principal difference is that here the alternating field is not imposed by external sources, but is self-organized, being created by moving polarization itself. Such a self-organized feedback field turns out to be essentially more effective than an externally imposed field. We use the system of units where the Planck constant is set to one.

II. EVOLUTION EQUATIONS
We consider a ferroelectric inserted into a cavity. Generally, if the sample is sufficiently large and especially when it is in contact with other media, say dielectric, then it can be separated into domains [39]. But we consider the case of a cavity containing no other materials inside it, except the ferroelectric itself, which is a singledomain sample.
Our main aim is to illustrate the idea of using a selforganized field of a resonator for accelerating polarization switching. We do not claim to treat a specific material, but we demonstrate the efficiency of the idea by a model. For this purpose, let us take the Hamiltonian in the pseudospin representation [40,41], having the form of an Ising-type model in a transverse field.
Here S α j is an α-component of the S = 1/2 spin operator characterizing an electric dipole at site j, Ω is tunneling frequency, J ij = J ji > 0 describes the strength of dipolar interactions, E tot is the total electric field acting on dipoles, and is a dipolar operator.
In what follows, the total electric field E tot will have two components that can be called longitudinal and transverse. It is important that the polarization would also have these two nonzero components.
This form of the Hamiltonian provides a good description of the so-called order-disorder ferroelectrics, although it can also be a reasonable approximation for other types of ferroelectrics [40,41]. Among orderdisorder ferroelectrics, it is possible to mention such as and their deuterated analogs, in which H 2 is replaced by D 2 . Similar Hamiltonians also are used for describing relaxor ferroelectrics [42].
Note that usually in the case of order-disorder ferroelectrics with spatially symmetric double wells, the vector of polarization possesses only the longitudinal zcomponent. However, we keep in mind the general case of nonsymmetric potentials (especially with respect to the inversion x → −x at the location of dipoles, because of which the polarization may have a transverse component. Hamiltonian (1) is the standard, widely used, Hamiltonian for describing macroscopic ferroelectric samples. For finite samples, in general, one should take into account the depolarizing field caused by the charges on the surfaces of the sample [43]. However, there are ways [44] of compensating surface charges, thus reducing or removing depolarizing fields.
Also, considering an external electric field E 0 , applied in the z-direction, leads to the appearance of the depolarizing field proportional to −P z , which is of the order of ρd 0 S, where ρ is the sample density. The energy, corresponding to the depolarizing field, is d o P z , which gives ρd 2 0 S. The latter expression equals the dephasing rate or transverse attenuation denoted by γ 2 = ρd 2 0 S. The magnitude of the energy, corresponding to the external field, is |d 0 E 0 |, which defines the dipole rotation frequency denoted as ω 0 = |d 0 E 0 |. For what follows, we will need a sufficiently strong external field, such that ω 0 be much larger than γ 2 . This is necessary for realizing coherent motion of dipoles, so that the reversal time be much shorter than the dephasing time. Under the condition ω 0 >> γ 2 , corrections, related to the depolarizing field, can be omitted in the Hamiltonian. At the same time, the attenuation rate γ 2 will be taken into account in the equations of motion.
The total electric field consists of two terms, where the first term is a field of the resonator cavity, in which the sample is inserted, and the second term is an external constant electric field. The resonator cavity is chosen such that it supports the TM 010 fundamental mode, whose electric field is directed along the cavity axis that is taken to be the x-axis. The resonator cavity electric field is a feedback field generated by the moving polarization with V being the sample volume. The equation for the feedback field can be derived in the standard way [45]. From the Maxwell equations inside the cavity with an inserted ferroelectric, it is straightforward to get the equation for the electric field generated by the ferroelectric polarization (4), where σ is conductivity and c is light velocity. It is possible to look for the solution to this equation in the form where e(r) is a cavity mode defined by the Helmholtz equation and normalized to the cavity volume V c , so that We are looking for the TM 010 fundamental mode, which, by definition, is the mode directed along the cavity axis that here is the axis x, which implies the conditions e y (r) = 0 , e z (r) = 0 .
The mode x-component is nonzero inside the cavity, while satisfies the boundary condition on the cavity cylindrical surface of radius R. The expression for the TM 010 fundamental mode is known [45] to be presented through the Bessel function of the first kind, with C 0 being the normalization constant. The boundary condition e x (R) = 0 corresponds to the first zero of the Bessel function J 0 (ωR/c) = 0, which defines the cavity natural frequency The normalization constant becomes Thus the TM 010 fundamental mode reads as where e x is a unit vector along the cavity axis x. Then we substitute expression (6), with mode (11), into Eq. (5), multiply the latter by form (11), take into account that P · e(r) = P x e x (r) , and integrate over the cavity volume. This leads to the equation in which γ = 2πσ is a cavity attenuation, ω is the cavity natural frequency (10), and is the filling factor corresponding to the TM 010 fundamental mode. The evolution equations, following from the Heisenberg equations of motion are The observable quantities are given by the statistical averages where N is the number of lattice sites and S = 1/2 is the spin value. The attenuation can be taken into account by employing the method of local fields [46,47], where the attenuation is caused by particle interactions acting in the local field formed by other particles, so that the dynamic variables are forced to relax to their local equilibrium values. In the present case, the latter are where S = 1/2 and the local equilibrium averages are expressed through variables (14) taken at the given moment of time. Since dipolar interactions are of long range, the mean-field approximation is applicable. In this way, the quantity S α j loc is defined as the average with the local equilibrium statistical operator for the ensemble of spins with the Hamiltonian Accomplishing explicit calculations for the local average (16), we keep in mind low temperatures, such that T ≪ J, where Thus we obtain the low-temperature local equilibrium values where we introduce the notation and define the frequency The positive value of ω 0 implies that the external electric field is directed downwards. Thus we come to the mean-field evolution equations for the pseudospin variables Here γ 1 is the longitudinal relaxation rate due to spinphonon interactions (see Blinc [41]), while γ 2 = ρd 2 0 S is the transverse attenuation caused by dipolar interactions.
By introducing the dimensionless feedback field where ρ = N/V , and taking account of the expression P x = ρd 0 Ss x , we get the feedback-field equation Here is the coupling rate characterizing the interaction between the ferroelectric sample and resonator.
Equations (21) and (23) define the dynamics of the variables s α and the feedback field h. The most interesting is the behavior of the dimensionless polarization s z = s z (t) as a function of time for different parameters, under the given initial polarization s z (0) = s 0 .

III. NUMERICAL SOLUTION
The polarization switching is realized in the following way. Suppose the ferroelectric sample is initially polarized along the axis z. The sample is placed inside a resonator cavity supporting the TM 010 fundamental mode and is subject to an external electric field directed opposite to the initial polarization. The polarization dynamics is governed by Eqs. (21) and (23). In order to precisely describe this dynamics, we have to fix realistic parameters typical for ferroelectrics [1,2,41].
First, let us notice that a ferroelectric characterized by Hamiltonian (1), without the last term containing electric fields, acquires spontaneous polarization below the critical temperature This temperature is positive for the tunneling frequency The interaction strength, due to dipolar forces, J ≈ ρd 2 0 . The electric dipole is d 0 = e 0 l 0 , with l 0 ∼ 10 −8 cm and the electric charge about a proton charge e 0 = 1.602×10 −19 C. Keeping in mind that one Coulomb 1C = 2.998 × 10 9 g 1/2 cm 3/2 /s, we find d 0 ∼ 10 −27 C cm ∼ 1D, where one Debye is 1D = 3.336 × 10 −28 C cm. For the density ρ ∼ 10 22 cm −3 , we obtain ρd 2 0 ∼ 10 −14 erg, that is ρd 2 0 ∼ 10 13 s −1 . Really, for typical ferroelectrics J ∼ 10 2 K, that is, J ∼ 10 13 s −1 . The dipolar forces induce the transverse attenuation γ 2 = ρd 2 0 S. Thus we have The longitudinal attenuation, caused by the interaction of pseudospins with phonons, is much smaller than the transverse attenuation, γ 1 ≪ γ 2 . The cavity is called resonant, since its natural frequency ω has to be tuned close to the dipole rotation frequency ω 0 , satisfying the quasi-resonance condition While the attenuations are to be smaller than ω 0 , so that Therefore, for parameter (19), we get Since JS = γ 2 ∼ 10 13 s −1 , to satisfy condition (29), we need that ω 0 be at least about 10 14 s −1 to 10 15 s −1 , which is in the near infrared or visible light range. This gives the wave vector k 0 ≡ ω 0 /c of the order 3 × 10 3 cm −1 to 3 × 10 4 cm −1 and the wavelength λ ∼ 10 −4 cm to 10 −3 cm. Resonator cavities in the range of visible light are widespread, and also there exist various cavities operating in the infrared region [48][49][50][51][52][53].
We solve numerically the system of Eqs. (21) and (23), concentrating our attention on the behavior of the polarization s z = s(t) as a function of time, for different parameters in the admissible range. In the figures, time is measured in units of 1/γ 2 and the frequency parameters are measured in units of γ 2 . As initial conditions, we need to fix the values s x (0) = 1 − s 2 0 , s y (0), s z (0) ≡ s(0) ≡ s 0 , h(0), and the time derivativeḣ(0). The initial polarization is positive, s 0 > 0. If some of other initial conditions, except s 0 , are not zero, the reversal begins immediately at t = 0. When all of them (except s 0 ) are zero, there is a time delay. In the figures, we show the results for the initial conditions s y (0) = 0, h(0) = 0, andḣ(0) = 0, while s 0 and, respectively, s x (0) = 1 − s 2 0 can be varied. The resonance is assumed, when ω = ω 0 . Figure 1 demonstrates the polarization reversal at different values of the resonator attenuation. For small γ, the polarization oscillates after the switching. Hence, in order to achieve a steady state after the polarization reversal, it is necessary to take larger γ. For γ = 10, the after-switching oscillations are suppressed. Thus, to avoid oscillations, it is preferable that the resonator ringing time τ ≡ 1/γ be shorter than the dephasing time T 2 ≡ 1/γ 2 . Figure 2 shows the dependence of the polarization switching on the tunneling frequency. The larger Ω, the shorter the delay time. As is clear from the evolution equations, this is because the tunneling triggers the polarization motion.
In Fig. 3, the role of the frequency ω is illustrated. The larger ω, the shorter the delay time and better the polarization inversion. This happens because the larger frequency makes stronger the coupling between the resonator cavity and the ferroelectric sample. Figure 4 shows the similar dependence of the polarization switching on the frequency ω, as in the previous figure, but for the initial polarization s 0 = 0.5, when s x (0) = 1 − s 2 0 is not zero. As is mentioned above, a nonzero s x (0) triggers the start of the polarization motion, so that there is no delay time, and the switching begins from t = 0.
In Fig. 5, the polarization switching for different initial polarizations is compared: s 0 = 0.5 (solid line) and s 0 = 1 (dashed-dotted line). For s 0 = 0.5, the initial transverse component s x (0) = 1 − s 2 0 is not zero, because of which the process of switching starts from the very beginning, practically at t = 0, without delay. Since the transverse component s x generates, by means of relation (23), the dimensionless cavity field h, hence the dimensional electric field inside the cavity E, the temporal behavior of h, as is seen from Fig. 7, is similar to that of the component s x .

IV. ANALYTICAL SOLUTION
Although the numerical solution of the previous section gives us an accurate description of the process of polarization switching, nevertheless it is desirable to have analytic solutions that would provide, at least approximately, explicit formulas allowing for the better understanding of the related physics and for straightforward estimates of characteristic quantities.
It is convenient to pass to the variables In terms of these variables, Eqs. (21) transform into the equation for the transverse component for the coherence intensity and for the polarization, Keeping in mind the case of resonance, defined by Eq. (27), and the existence of the small parameters described in Eqs. (26), (28), and (29), we notice that the variables u and h can be classified as fast, while the variables w and s, as slow. In that case, for solving the given system of equations, it is admissible to resort to the averaging techniques [54,55]. Below, we follow the variant of the method described in detail in Refs. [56,57]. First, we solve the equations for the fast variables, keeping there the slow variables as quasi-integrals of motion. Such a solution is straightforward, although cumbersome, since the equations for the fast variables become linear with respect to the latter, when the slow variables are kept fixed. Then the found solutions for the fast variables are substituted into the equations for the slow variables, with the averaging of the slow-variable equations over time. This yields the equations for the guiding centers, which can be analyzed. All this machinery has been thoroughly described in Refs. [56,57] and its use has been demonstrated for studying the dynamics of magnetic systems [58][59][60].
To present the solutions for the fast variables, under fixed slow variables, in a compact form, we take account of the small parameters and introduce several notations. We define the coupling parameter characterizing the strength of the coupling between the ferroelectric sample and resonator, the coupling function (35) describing the dynamics of the ferroelectric-resonator interaction, and the effective frequency Thus we obtain the transverse component and the feedback field Note that from Eq. (18), we have Substituting the fast variables into the equations for the slow variables, and averaging the resulting equations over time gives the guiding-center equations for the coherence intensity, and for the polarization where The parameter γ 3 is very small. However, it cannot be neglected, since it plays an important role in triggering the polarization motion at the initial stage. At the very beginning of the process, when t → 0, so that the coupling function (35) is close to zero. Then, keeping in mind that usually γ 1 ≪ γ 2 , the equations of motion become Their solutions are which, in view of inequalities (42), can be simplified to where w 0 ≡ w(0) = 1 − s 2 0 and s 0 ≡ s(0). These are the solutions at the initial stage, when the motion of individual polarizations is not mutually synchronized.
The coupling function (35) grows with time, implying the increase of the magnitude of the resonator feedback field, which collectivizes the individual polarizations, forcing them to move coherently. The influence of the feedback field becomes crucial after the coherence time t coh , when the coupling function grows so that This defines the coherence time with τ being the resonator ringing time.
If the ferroelectric-resonator coupling is strong, such that gs 0 ≫ 1, then the coherence time is At the coherence time, the solutions (44), that is can be written as where we assume that the coupling parameter g is sufficiently large, so that inequalities (42) are yet valid at t coh . After t coh , the coupling function (35) quickly grows reaching the value g(1 − As) ≃ g. At this stage, the parameters γ 1 and γ 3 , that are much smaller than γ 2 , and especially than gγ 2 , can be neglected. Then Eqs. (39) and (40) become The latter equations enjoy the exact solutions for the coherence intensity and polarization The quantities γ s and t 0 are the integration constants that are defined by sewing expressions (51) and (52) with Eqs. (49). This gives the switching time τ s ≡ 1/γ s , in which and the delay time The delay time shows the time, when the switching starts, while the switching time is the time during which the polarization reversal occurs. The delay time can also be written as which demonstrates its explicit dependence on w coh .

V. CHARACTERISTIC QUANTITIES
The derived analytic expressions provide a transparent illustration for the role of different system parameters and their combinations. To study more carefully these dependencies, let us keep in mind the case of strong ferroelectric-resonator coupling, when gs 0 ≫ 1. Then the coherence time (47) takes the form which shows that the larger the frequency ω 0 , the shorter this time. Quantities (53) read as γ s ≃ gγ 2 1 + 0.385 The switching time τ s ≡ 1/γ s becomes so that larger ω 0 makes the process of switching faster. The delay time (55) acquires the form which depends on the value w coh . The latter is connected with the initial polarization s 0 .
When the system at the initial moment of time is well polarized, with s 0 = 1 and w 0 = 0, then And the delay time turns into While when s 0 = 0.5 and w 0 = 0.75, then the delay time is To get concrete estimates, let us take the typical values of parameters as have been used when numerically solving the evolution equations: ω 0 ∼ 100γ 2 , Ω ∼ 0.1γ 2 , γ ∼ 10γ 2 , and s 0 ∼ 1. Then g ∼ 10, γ s ∼ γ g ∼ 10γ 2 , and γ 3 ∼ 10 −6 γ 2 . For the coherence time, we have t coh ∼ 10 −15 s and for the switching time we get τ s ∼ 10 −14 s. Diminishing s 0 decreases the delay time. Thus, if s 0 = 1 and w 0 = 0, then w coh ∼ 10 −8 and t 0 ∼ 10 −13 s. But if s 0 ≈ 0.5, so that w 0 ∼ 1, then t 0 ∼ 10 −14 s. These estimates are in good agreement with numerical calculations.
An important question is how the switching time is limited in realistic materials. In particular, what is the relation between the switching time τ s and the cavity ringing time (delay time τ ≡ 1/γ). From the above estimates, we find the ratio It looks that by varying the system parameters, it is possible to make this ratio rather small. However, there are limitations for the variation of the parameters. Thus, for a good quality cavity one has γ ≪ ω 0 . But γ cannot be arbitrarily small, since for γ < γ 2 , there appear oscillations in the polarization. In order to realize a stable switching without oscillations, it is necessary to take γ ≫ γ 2 . Therefore ratio (63) lies in the interval For the infrared region, where ω 0 /γ 2 ∼ 10, we have which actually means that the switching time τ s is of order of the ringing time τ . In the visible light region, when ω 0 /γ 2 ∼ 100, we find This tells us that again the switching time is correlated with the ringing time. For the visible light, it looks admissible to reach the shortest switching time of order τ s ∼ 10 −15 s.

VI. COHERENT RADIATION
The motion of electric dipoles has to produce electromagnetic radiation. If this motion is coherent, the produced radiation should also be coherent. Since the sample is inside a resonator cavity, the radiation can propagate only along the cavity axis, that is, along the axis x. The radiation intensity in the direction of n = e x consists of two terms describing incoherent and coherent radiation, I(n, t) = I inc (n, t) + I coh (n, t) .
Radiation, produced by moving dipoles can be described in the following way [45,57,61]. The incoherent radiation intensity reads as and the coherent radiation intensity is with the shape factor where L is the cavity length and γ 0 is the natural width For the frequency ω 0 ∼ 10 15 s −1 , we have the wavelength λ ∼ 10 −4 cm and the natural width γ 0 ∼ 10 4 s −1 . Then the radiation intensities at the maximum, are The number of dipoles that could radiate coherently can be estimated as N ∼ ρλ 3 . If we consider a small sample, with the length L smaller than the radiation wavelength λ, then the shape factor (68) is of order one. In that case, for the density ρ ∼ 10 22 cm −3 , we get N ∼ 10 10 . And for the radiation intensities, we find When the frequency is ω 0 ∼ 10 14 s −1 , then λ ∼ 10 −3 cm and γ 0 ∼ 10s −1 . The radiation intensities are Considering again coherently radiating dipoles, with the number N ∼ ρλ 3 , we have N ∼ 10 13 . This gives for the radiation intensities I inc (n, t) ∼ 10 −7 W , I coh (n, t) ∼ 10 5 W .
Such a level of radiation can be easily measured. The prevailing coherent component of radiation shows that this radiation is of the type of superradiance.

VII. CONCLUSION
A method of ultrafast polarization switching in ferroelectrics is suggested. The main idea is to place a ferroelectric sample into a resonator cavity. In the presence of a constant electric field, directed opposite to the ferroelectric polarization, the sample is in a nonequilibrium state. As soon as the polarization starts moving, it produces an electric field in the cavity. This field acts back on the sample forcing the polarization to move faster. Thus the ferroelectric itself generates a feedback field accelerating the polarization motion, so that there is no necessity of applying external alternating fields, as one usually does for realizing polarization switching. It turns out that the self-organized feedback field is essentially more effective for the polarization reversal than an externally imposed field.
The system of equations, describing the ferroelectric polarization and feedback field is solved numerically and also analytically by means of averaging techniques. This makes it possible to give a detailed description of the whole procedure, to study the role of the system parameters, and to estimate the characteristic quantities involved in the process. The polarization switching can be realized extremely fast: for the parameters of typical ferroelectrics, the switching time can reach femtoseconds. This ultrafast polarization reversal generates a coherent electromagnetic pulse.
We stress that the main goal of the paper is to attract attention to the possibility of accelerating the polarization switching in ferroelectrics by using the selfacceleration effect caused by the action of a resonatorcavity feedback field. This is the idea of principle that, to our knowledge, has not been considered for ferroelectrics before. It goes without saying that the method is not necessarily applicable to any particular material. Thus the method seems to be not applicable for the orderdisorder ferroelectrics with spatially symmetric double wells, where there is only the longitudinal polarization. However the considered model assumes the general case of asymmetric potentials for which the method is applicable. We also hope that, since there are various types of ferroelectric systems [41,[62][63][64], there can exist other materials for which the suggested idea could work.
In order that the suggested method could be realized, the existence of two spin components of polarization, longitudinal and transverse, is required. If one keeps in mind an order-disorder ferroelectric with lattice-site double wells that are ideally symmetric with respect to spatial inversion (especially with respect to the inversion x → −x), then there is only a longitudinal component, and the sample polarization is expressed through the zcomponent of the spin operator. But in the general case of an asymmetric double well, the sample polarization contains a term with the x-component of spin. The asymmetry can be induced by stress or by incorporating into the sample admixtures or vacancies. Thus, the inclusion of vacancies in order-disorder ferroelectrics is attributed to the breaking of spatial inversion symmetry along different directions [65]. The symmetry in order-disorder ferroelectrics can be distorted by the action of a transverse electric field [66].
To illustrate an idea, one has to consider some model. We considered an Ising-type model in a transverse field. This kind of models, employing the spin representation, is widely used for order-disorder ferroelectrics [41,[62][63][64] and relaxor ferroelectrics [42,67].
In addition to different ferroelectrics [41,[62][63][64]67] and multiferroics [68,69], ferroelectric-type spin mod-els are widely used for describing the systems of polar molecules, Rydberg atoms, Rydberg-dressed atoms, dipolar ions, vacancy centers in solids, and quantum dots [70][71][72][73][74][75]. These systems can form self-assembled lattice structures or can be loaded into external potentials imitating crystalline matter. The characteristics of these dipolar materials can be varied in a very wide range. Therefore the realization of the method considered in the paper looks feasible.