Kinetics of 90° domain wall motions and high frequency mesoscopic dielectric response in strained ferroelectrics: a phase-field simulation.

The dielectric and ferroelectric behaviors of a ferroelectric are substantially determined by its domain structure and domain wall dynamics at mesoscopic level. A relationship between the domain walls and high frequency mesoscopic dielectric response is highly appreciated for high frequency applications of ferroelectrics. In this work we investigate the low electric field driven motion of 90°-domain walls and the frequency-domain spectrum of dielectric permittivity in normally strained ferroelectric lattice using the phase-field simulations. It is revealed that, the high-frequency dielectric permittivity is spatially inhomogeneous and reaches the highest value on the 90°-domain walls. A tensile strain favors the parallel domains but suppresses the kinetics of the 90° domain wall motion driven by electric field, while the compressive strain results in the opposite behaviors. The physics underlying the wall motions and thus the dielectric response is associated with the long-range elastic energy. The major contribution to the dielectric response is from the polarization fluctuations on the 90°-domain walls, which are more mobile than those inside the domains. The relevance of the simulated results wth recent experiments is discussed.

Substantial efforts have been devoted to ferroelectric (FE) materials and their applications in advanced electronic technologies12. The core ingredients in the physics of ferroelectrics include electric polarization P and its static/dynamic responses to electric field E or/and other stimuli. In the mesoscopic level, the dielectric and polarization responses of a ferroelectric are determined by the specific domain structure and realized by the domain wall motions. In particular, the high-frequency applications of ferroelectrics have received continuous attention. A number of FE materials in thin film form have been used for high-frequency devices in microwave and optical communications as well as computing applications because of their high dielectric permittivity34. A full understanding of the domain wall motion driven by either static or high-frequency dynamic force has been the core issue of the physics of ferroelectrics and the basis for their applications.

The FE domain wall motion is a complicated process and depends on a series of intrinsic particulars such as defects, domain structures, and strains567. The consequence of strain, a topic to be dealt with in this work, has been highly concerned. Besides conventional scope of ferroelasticity, an externally imposed strain in a paraelectric surprisingly enables remarkable ferroelectricity8. Externally imposed strain can modulate the FE phase transitions too91011. These strain effects raise particular attention to FE thin films deposited on rigid substrates, in which the induced strain can be controlled due to the lattice mismatch of the substrates with the films. It thus allows the performance improvement of the FE thin films by means of ‘strain-engineering’ the FE domains in the mesoscopic level.

Along this line, the most interested issue is the strain effect in perovskite FE oxide thin films with tetragonal lattice symmetry like BaTiO3 and PbTiO3. First, epitaxial BaTiO3 and PbTiO3 thin films deposited on substrates such as SrTiO3, LaAlO3, and MgO, are often explored not only for fundamental understanding but also for potential applications, while textured polycrystalline films are also investigated. Second, these thin films offer regular twin-like (stripe-like) 90°-domain structure in coexistence with 180°-domains due to the intrinsic ferroelastic effects. In most cases, the 90° ferroelastic domains are dominant, and the order parameter is strongly coupled to the strains. Furthermore, the substrate induced strain imposes a competitive or coherent coupling with the internal ferroelastic strain in the films, making the domain structure and domain wall motion complicated. These behaviors, particularly in the high frequency range, are much less studied. It is believed that externally imposed strains have strong impacts on the domain wall motion and thus the high frequency dielectric permittivity121314.

A prominent feature of the strain effects is the 90° domain walls vibration in response to an ac electric field of high frequency, while the static responses have been well addressed. The dynamics of the domain wall may be characterized by the variation in dielectric permittivity as a function of the ac electric field (frequency ω and amplitude E), as long as the domain structure is well defined. A recent experiment demonstrating this dielectric response was carried out on clamped Pb(Zr1-Ti)O3 (PZT) thin films deposited on Si substrates15. By creating cavities beneath the Pt/PZT/Si capacitors and cracking, these released PZT films show dramatic variation in the global dielectric nonlinearity and the frequency dependence as a function of mechanical clamping. The sequent piezoelectric examinations show that the increased local mobility of domain walls must be responsible for these behaviors, where “mobility” refers to the capability of domain walls in response to stimuli. This work among many others raises an important issue on how the 90°-domain structures in FE thin films couple dynamically with externally imposed strains16. This question motivates us to investigate the dynamics of domain walls in the 90° domain structure modulated by strains, and consequently the dielectric responses at the mesoccopic level.

While the strain coupling in polycrystalline PZT thin films may be too complicated, our purpose is to focus on tetragonal FE lattice (e.g. BaTiO3 or PbTiO3) and reveal the dynamics of the 90°-domain structure which couples the internal ferroelastic strains with externally imposed normal strains, in response to electric field E. Instead of addressing the dynamic motion equation171819, we perform the phase-field simulations by which a continuous distribution function of the domain orientations is used to characterize the domain wall motion2021. Such a phenomenological method enables us to directly examine the mesoscopic picture of these dynamic processes. We also investigate the underlying mechanism with which the strains modulate the domain wall motion in the energy landscape. In this phase-field model, the dipole-dipole and elastic interactions are considered22, since these long-range interactions play an important role in configuring the domain structure1723. The strain vector appears as an order parameter and can be modulated by external load.

We consider a tetragonal FE system as approached by an L × L two-dimensional (2D) lattice with periodic boundary conditions24, and start from the Ginzburg-Landau theory. On each site, an electric dipole P(r) = (P, P) normalized by a pre-factor P25 and an elastic displacement vector u(r) = (u, u) are imposed. We clarify that our simulation is not necessarily associated with a thin film, although the electric dipoles in our simulation can only relax on the in-plane configuration. An extension of the results to other geometry is plausible. For simplification, our simulation will not deal with the finite boundary problem since it will introduce complicated corrections to the mechanical balance25, while open-circuit boundary issue will only be discussed briefly. Given this assumption, the lattice can be mapped into a small region of real material. By these approximations, we can focus on the details of the local dielectric response in mesoscopic scale.

The electric field E = (E, E) is confined on the lattice plane. The total free energy for this FE lattice can be written as26: where F, F, F, F, F, and F are the Landau potential, gradient energy, dipole-dipole interaction, elastic energy, electrostrictive energy, and electrostatic energy, respectively. Term F extending to the sixth-order is written as: where A1, A11, A12, A111 and A112 are the Landau expansion coefficients and A1 = A10(T-T0). The lowest-order expression of term F is: where P and r = (x, y), and G11, G12, G44 and G′44 are the gradient energy coefficients20. Terms F and F can be written respectively as: where F includes two contributions. One is the anisotropic part F(P) and the other is the isotropic part F′(δP) if we treat P(r) = P+δP(r) where P is the spatially homogeneous average polarization and r is the spatial vector. It is noted that P can be defined by: where V is the volume of lattice.

Term F(P) describes the depolarization energy and shares the same form as F by difference of a 1/2 factor. In real system, the depolarization field is compensated by free charges for the cases with open boundaries. For the present case, no free charge is included. The effective electric field is independent of the depolarization field and simply equivalent to E. In experiments, the electric field is usually applied via the electrodes in close-circuit which introduces free charge onto the interface. This will cancel out the depolarization field. Term F counts an integration over the whole lattice, and a realistic calculation is done either by the Fourier transformation or by finite truncation treatment2728. For a 2D lattice, the finite truncation is a sufficiently accurate approximation as long as the truncating distance R is big (R = 8 in our simulation)242930.

The elastic energy F yields: where C is the elastic stiffness tensor which has only three independent elastic constants for a square lattice.

The electrostrictive energy F is: where q11 = C11Q11+2C12Q12, q12 = C11Q12+2C12(Q11+Q12), and q44 = 2C44Q44 are the effective electrostrictive coefficients, Q are the electrostrictive coefficients in external stress free state. Mathematically, we can also separate the total strain η into a homogenous component η, and a heterogeneous one δη which denotes the microscopic strain distribution at site (i, j)2225: here external strain η is imposed either by a substrate or an external load. Similar to generally accepted treatment202228, order parameter u(r) is reduced to a function of P(r)31. The step-by-step procedure of this derivation can be found in the Supplementary materials. To access the dielectric response over frequency domain, we compute dielectric permittivity ε(ω) over a broad range of frequency. The algorithms are presented in the Method section. The parameters for the simulations are listed in Table I, and for the details of the simulation procedure one can refer to earlier works3233.

Table 1

Physical parameters chosen for the simulation (τ −1 = |A1|D) [24, 25]. All these parameters appear in the dimensionless form

Parameter (unit)	Value	Parameter (unit)	Value	Parameter (unit)	Value
L	64~256	A*1 (|A1|)	−1.00	A*11 (A11P20/|A1|)	−0.24
A*12 (A12P20/|A1|)	2.50	A*111 (A111P40/|A1|)	0.49	A*112 (A112P40/|A1|)	1.20
G*11 (G11/A21|A1|)	1.60	G*12 (G12/A21|A1|)	0.00	G*44 (G44/A21|A1|)	0.80
G′*44 (G′44/A21|A1|)	0.80	C*11 (C11/|A1|P20)	2.75	C*12 (C12/|A1|P20)	1.79
C*44 (C44/|A1|P20)	0.543	q*11 (q11/|A1|)	0.143	q*12 (q12/|A1|)	−0.0074
q*44 (q44/|A1|)	0.0157	τ* (τ)	0.0004

Results

Domain wall kinetics under dc electric field

We first investigate the 90°-domain wall motion driven by a static (dc) E. To characterize the wall motion, we define the domain width l deviating from the equilibrium width l0 under zero electric field. Thus, Δl = l-l stands for the offset in response to E. In order to minimize the effect of domain wall irregularity, parameters l and l are calculated indirectly. A set of rectangle-like regions, whose two sides are on the domain walls to enable their areas as big as possible, are taken. We sum the areas of these rectangles. The local domain width l is obtained by dividing the area by its length and then performing statistics over sufficient number of rectangles in the same region.

In Figure 1 is shown the kinetics of field-driven domain wall motion, given different strains from slightly compressive ones to tensile ones. Both E and η are along the y-axis, thus the l indicates the width of a2-domain where the dipoles align along the y-axis. The evaluated Δl(t) data are plotted in Figure 1(a) where the applied strains are labeled numerically. For each strain, Δl(t) shows similar behavior and increases rapidly in the early stage and tends to be saturated at Δl in the late stage. Nevertheless, Δl depends on η, and the bigger η the smaller Δl. In Figure 1(b) we plot the evaluated l(η) for the a2-domain and corresponding Δl(η), showing the monotonous increasing of l(η) and decreasing of Δl(η). For definition of each symbol and motion pattern of domain wall, one can refer to the Supplement materials.

Figure 1

Evaluated parameter Δl as a function of time t given a series of η labeled numerically (a), l and Δl as a function of η respectively (b), parameters a and b as a function of η respectively (c), and initial domain wall motion speed v as a function of η (d).

E6|A1|P.

The kinetics of wall motion was investigated earlier34 and can be described by a simple model1723, yielding the kinetic equation Δl(t) = a(1-e−) where a = Δl and b are the fitting parameters. As shown in Figure 1(c), the monotonous decrease and increase of a and b as a function of η respectively imply that the a2-domain extension becomes tougher upon the transition from compressive strain to tensile one. It is also possible to evaluate the initial motion speed of the a2-domain walls by taking v = dΔl/dt|→0 = ab, as shown in Figure 1(d), consistent with the above argument. In other words, for compressive strain (η<0), the motion speed of domain walls is higher than that for tensile strain (η>0).

For the cases with strain η not parallel to E, the above analysis is qualitatively correct. The tensile strain makes the domains with P//η wider and the compressive one makes it narrower. This behavior allows a competition between the strain effect and electric field effect, and they can be coherent or cancelled depending on their orientation relationship. However, when the strain contains shear components, the situation becomes more complicated and extensive calculation is not discussed here.

Domain wall vibrations and dielectric response: η0 = 0

Based on the above result, one understands that the 90°-domain structure clamped by external stress has different stability characteristics from the stress-free state. This difference can be discussed in the clamped 90°-domain structure driven by the ac-electric field, characterized by variation of dielectric permittivity as a function of η. In particular, the dielectric response in the frequency domain is related to the wall motion.

We first address the dielectric response in η = 0. Figure 2(a) and (b) show the real part ε(f) and imaginary part ε(f) at three θ angles. Here a small E = 0.6|A1|P is chosen, with the field direction defined by angle θ between the x-axis and E. In general, ε(f) decreases gradually with increasing f for all the three cases, while ε(f) shows two peaks at characteristic frequencies f ~ 0.1τ−1 and f ~ 7τ−1, which are respectively referred as the low-f and high-f dispersions. Here τ = 1/|A1|D is the characteristic time for electric dipole flip (see the Method section for details). The striking feature is the low-f dispersion anisotropy, i.e. the θ-dependence which is the most remarkable at θ = 90° (and θ = 0 too, four-fold symmetry). This θ-dependence is shown by ε(θ) at f = 0.01τ−1 as an example in Figure 2(c), characterized by the typical periodic variation.

Figure 2

Simulated dielectric permittivity spectra: (a) real part and (b) imaginary part over a wide range of frequency, at θ = 90°, 112°, and 135°, respectively. (c) Dielectric permittivity real part ε as a function of angle θ at a constant f = 0.01τ−1.

The above behaviors can be understood by investigating the instant evolution of domain structure. By turning E from θ = 0 to θ = 180°, one checks the domain evolution and dielectric dispersion. The low-f dispersion is related to the domain wall vibrations, while the high-f dispersion is attributed to the flip of individual dipoles. In fact, the dispersion peak around f is also observed in the mono-domain lattice.

Now we investigate the origin for the low-f dispersion anisotropy in the domain scale. In η = 0, the spatial distributions of ε at f = 0.05τ−1 in four θ angles are presented in Figure 3(a)~(d), respectively, where the color scales the intensity (ε = 0.0~1.0). It is immediately seen that the dielectric permittivity mainly comes from the contribution of wall vibrations, while those dipoles deeply inside the domains contribute little. This characteristic makes the dielectric response very specific.

Figure 3

Snapshoted patterns of dielectric permittivity real part ε at constant frequency f = 0.05τ−1 with angle θ = 45° (a), 90° (b), 112° (c), and 135° (d).

At θ = 45°, the domain structure and dielectric response are shown in Figure 3(a). Again, the dielectric response mainly comes from the contribution of electric dipoles on the walls and near regions. This feature can be understood from the electrostatic energy -P·E. The walls have the highest local mobility. No matter E is positive or negative, there is always one of the neighboring two domains, inside which all the dipoles take the 135° angle from E. This domain will shrink while the other will expand, making the wall move easily. This is also the reason why ε is the highest at θ = 45°, as seen in Figure 2(c).

At θ = 90° (or θ = 0), as shown in Figure 3(b), the wall mobility is slightly lower, and ε is lower too. In this case, E is parallel to the dipoles in one domain and perpendicular to those in the other. At θ = 112° and θ = 135°, as shown in Figure 3(c) and (d), respectively, the wall mobility falls down even further, resulting in even lower ε. The wall mobility becomes the lowest at θ = 135°. Figure 3(d) shows that almost the whole lattice has the identical ε except the very thin and dim lines on the walls. The reason is that term -P·E in two neighboring domains are almost equivalent and they compete with each other, hindering the wall motion. The dielectric response is weak with no frequency dispersion.

We finally check the dielectric response in the high-f range and one example is given in Figure 4(a) and (b) at f = 6τ−1 for two specific θ angles: θ = 90° and θ = 135°. The dielectric distribution over the whole lattice is roughly homogeneous. For details, one looks at the case in Figure 4(a) and finds only weak color contrast between the two neighboring domains a1 and a2. The domain a1 whose dipoles align along the x-axis but perpendicular to E shows slightly higher ε than that of domain a2. The reason is obvious that the dipoles in domain a1 are more fluctuating than those in domain a2.

Figure 4

Snapshoted patterns of dielectric permittivity real part ε at constant frequency f = 6τ−1 with angle θ = 90° (a) and 135° (b). At constant frequency f = 0.05τ−1 with different tensile strains along the y-axis: (c) η = 0 and (d) η = 0.7%. The red lines below the snapshots indicate the alignment of ε along the . Hereafter, E0.6|A1|P.

Domain wall vibrations and dielectric response: η0>0

Now we discuss the cases with η≠0. We only present in details the results on normally strained lattice. When the lattice is normally strained, the domain structure is deformed. To clarify the similarity and difference between the strain-free and strained lattices, we consider the simplest situation: η = 0.7% under E with f = 0.05τ−1 and E = 0.6|A1|P, both aligned along the y-axis (θ = 90°). We present in Figure 4(c) and (d) the strain-free domain structure and strained structure respectively, as well as the spatial distribution of ε. For the two cases, the distributions are similar in amplitude but the high-ε spatial profiles across the domain walls are much wider for the strain-free lattice than those for the strained lattice. In the qualitative sense, the dielectric permittivity averaged over the whole strained lattice is lower than that over the strain-free lattice, at least in the low-f range. This η-dependence can be further illustrated in Figure 5(a), where ε(f) and ε(f) at three tensile strains (η = 0, 0.2%, 0.7%) are plotted. Both the real and imaginary parts are remarkably suppressed by tensile strain, and the f has a slight shift towards the high-f side. We also calculate the instant responses of the Δl and polarization component P against E at the three tensile strains, shown in Figure 5(b). The one-to-one correspondence between Δl, P, and E, respectively, is observed. The responses of Δl and P are synchronous with E in the low-f range, and delayed in the high-f range. Furthermore, the wall vibration amplitude is suppressed by the tensile strain.

Figure 5

(a) Simulated dielectric permittivity spectrum around f = 0.05τ−1 given η = 0, 0.2%, and 0.7% along the y-axis. (b) From top to bottom: parameter Δl as a function of time t given η = 0, 0.2%, and 0.7% at f = 0.05τ−1, the y component of total polarization as a function of time t, and the ac electric field E as a function of time t. The ac electric filed is along y-axis, E0.6|A1|P.

The above results refer to the simplest situation. When η and E are aligned along arbitrary directions independently, more complexity is seen in terms of the domain structure evolution and dielectric response. For some specific geometry, the strain and electric field may compete with each other. Some more discussion will be given below. However, in general, the calculated results are qualitatively similar to those for the simplest situation: the tensile strain suppresses the domain wall vibration, thus reducing the dielectric permittivity.

Domain wall vibrations and dielectric response: η0<0

Now we check the cases with η<0. Referring to the results under static (dc) E, as shown in Figure 1, one sees that the static compressive strain assists the 90° domain wall motion. It is thus expected that the dielectric permittivity in compressed lattice will increase. The η-dependences of ε(f) in three compressive strains (η = 0, −0.1%, −0.15%) along the y-axis are plotted in Figure 6(a), consistent with the expected results. Similar behaviors are observed for the strain and electric field both applied along the x-axis.

Figure 6

(a) Simulated dielectric permittivity spectrum around f = 0.05τ−1 given η = 0, −0.1%, and −0.15% along the y-axis with the ac electric along the y-axis too. The calculated real part ε at tensile strain (b) and compressive strain (c) perpendicular to the ac electric field along the y-axis. E0.6|A1|P.

We also check the situations where E is not parallel to η. One example is shown in Figure 6(b) for η>0 and Figure 6(c) for η<0, where E is along the y-axis and η is along the x-axis. Our extensive calculations establish the qualitatively similar behaviors: the tensile strain suppresses the 90° domain wall motion and thus the dielectric permittivity, while the compressive strain enhances the domain wall motion and the dielectric permittivity, no matter whether E is not parallel to η or not. We summarize spectrum ε(f, η) in Figure 7(a) and (b). At f0, while this tendency becomes negligible as f>f such as f = 0.5τ−1 since the single dipole response becomes dominant at this frequency. The overall evolution of ε(f, η) is plotted in Figure 7(b).

Figure 7

(a) Simulated dielectric permittivity real part ε as a function of η with the ac electric field frequency f = 0.01τ−1, 0.05τ−1, and 0.5τ−1. (b) A 2D plot of ε as a function of η and f. The ac electric field and strains are all along the y-axis. E0.6|A1|P.

Discussion

It should be mentioned that the periodic boundary conditions allows the total strains in all directions to be exactly balanced out, and the strain effect is of long-range and coupled with mechanical boundary conditions. However, if other mechanical boundary conditions are considered, such as free boundary, the total strains can be partially relaxed through the free boundaries. The 90°-domain walls can be more moveable and thus show more significant response to electric field, as partially discussed in Ref. 15.

To understand the simulated results, one may give additional discussion by looking at the energy landscape and comparing the simulated results and experiments, although relevant experimental data are really rare.

Energy landscape

We calculate the elastic energy distribution associated with the 90° domain structures. For a single domain lattice, the total strain field is homogeneous, thus accommodating extremely large elastic energy. The single domain is decomposed into the 90° domain structure so that the total elastic energy can be relaxed. Owing to the lattice volume conservation, for one domain, the compressive strain along one direction (e.g. the x-axis) is always accompanied with the tensile strain along the other direction (the y-axis), and verse vice. One can refer to the Supplement materials for the strain distribution.

First, we consider the η = 0 case. Given a static electric field along the y-axis, the domain walls move into the a1-domain in compensation with the extension of the a2-domain width. The continuous shrinking of the a1-domain is accompanied with the increasing magnitude of elastic strain (e and e) inside the a1-domain. Since the total elastic energy is proportional to η2, the rapidly enhanced total elastic energy inside the shrunk a1-domain acts as the resistant force against the electric field, responsible for eventual termination of the wall motion when the a1-domain becomes sufficiently narrow, as shown in the left column of Figure 8(a) and Figure 8(b), where the total elastic energy F(x, y) is plotted. This explains why Δl tends to be saturated at Δl as time goes infinitive.

Figure 8

Evaluated spatial contours of elastic energy F in the 90°-domained lattice.

The left column (a) refers to η = 0 and the right one (b) refers to η = 0.4%. The top row refers to E = 0, the middle row refers to t = 720τ after E (dc) along the y-axis applies to the lattice, and the bottom row refers to t = 720τ given E (dc) along the x-axis. E6|A1|P.

This scenario applies to the lattice with externally imposed strain. Take η = 0.4% along the y-axis as an example. At t = 0 with E = 0, the whole elastic energy distribution on the right column of Figure 8 shifts upward with respect to the case of η = 0. However, one can observe a much higher energy distribution in the shrinking a1-domains. Therefore, the a1-domains are already under highly tensile state along the x-axis even at E = 0, when an externally imposed strain along the y-axis is applied. In this case, a static electric field along the y-axis may further force the wall moving into the a1-domain but will be highly resisted by the high strain energy, as shown in the right column of Figure 8. Therefore, additional shrinking of the a1-domain driven by the electric field would become even tougher, given that η = 0.4% already makes the a1-domain narrow.

If strain η is applied along other orientations, no matter whether it is parallel to E or not, the wall motion behaviors remain qualitatively similar. One example is given by applying E along the x-axis where E and η are normal to each other. As shown in Figure 8(c), a strain η along the y-axis enhances remarkably the total elastic energy in the a2-domain, making the further shrinking of the a2-domain more difficult than the case with η = 0. Therefore, the above-discussed results don't lose the generality.

Comparison with experiments

The implication of the above simulation data can be checked by comparing the simulated results with experimental data on the contribution of domain wall vibrations to the dielectric response in BaTiO3 (BTO) and Pb(Zr1-xTix)O3 (PZT) ceramics353637, although these data have not yet received confirm from the direct detecting of the mesoscopic scale dielectric distribution across the 90° domain walls. Therefore, such comparisons may not be so direct in quantitative sense. For PZT, experiments35 showed that above 1011 Hz the real part of dielectric constant begins to drop, due to the frozen dipole oscillators themselves. In our simulation, the corresponding frequency is assigned as f = 7τ−1 in Figure 2(a). Thus we have the characteristic inverse time scale τ−1 ~1010 s−1, and the calculated f ~109 Hz.

Regarding the PZT ceramics, experimental data36 also indicate the dielectric dispersion around 1.0 GHz, which was argued to arise from the vibrations of the frozen 90° domain walls. For the dielectric permittivity magnitude, we take the dimensionless factor (|A1|ε*)−1~103, where A1 is 3.8(T-479) × 105C−2m2N taken from PZT and ε* is the vacuum permittivity. The calculated real part of the dielectric permittivity at θ = 45° is ~200, as shown in Figure 2(c), and the maximum of the imaginary part is ~100. Indeed, experimental measurements37 on ceramics PZT (x = 0.48–0.52) gave the difference in the real part between the extremely high frequency and the very low frequency, which is 250–500. The peak of the imaginary part is 150–265. These values agree roughly with our simulation results here. The above comparisons allow us to claim that the present calculations are reliable even in quantitative sense.

We also find some experimental results about the effects of tensile and compressive strains on the overall dielectric constant. In microwave frequency range, it was shown that dielectric constant and tunability of BTO films grown on MgO gradually decrease as the in-plane strain goes from the compressive type to tensile type38. Similar behaviors were observed in polycrystalline Ba0.6Sr0.4TiO3 thin films upon increasing tensile strain39. Those results are consistent with our calculation in a qualitative sense. Unfortunately, none of those experiments establishes a clear logic between the dielectric response and the 90° domain wall motion, although the frequency range is properly related to the 90° domain wall vibration. Establishing this logic is a technical challenge for experimentalists. To reveal the domain wall response to external electric field, one may map the local polarization response in real-space. Nevertheless, this task becomes difficult for such a high frequency. In this sense, the present simulation seems to be unique for bridging the relationship between the dielectric response in microwave frequency and the strain via the microscopic domain scale.

Methods

Calculation of ac dielectric permittivity

The temporal evolution of the dipole lattice is tracked by solving the Ginzburg-Landau (TDGL) equation which takes the following form: where t is time scaled in unit τ with τ −1 = |A1|D, and D is the kinetic coefficient. We also calculate the dielectric permittivity as a function of E. For a dielectric system that cannot polarize instantaneously in response to an electric field, the total electric polarization as a function of t can be described as: where P(t) is a convolution of electric field E(t) at previous times with time-dependent permittivity ε(t) = ε(t)+iε(t). Therefore, its Fourier transformation can be directly written as where P(ω) and E(ω) are the Fourier transformations of P(t) and E(t) respectively32. For the case of ac sine electric field such as E(t) = E·sin(ω0t) with ω0 = 2πf and ω = 2πf, E(ω) can be directly written as: where E′0 is the coefficient of Fourier transformation and ω0 is the frequency. Polarization P(ω) can be calculated by Fourier-transforming the temporal evolution spectrum P(r, t) in Eq.(10). In details, the real-spatial spectrum of P(r, t) at site r and its spatial average P(t) over sufficient number of time periods is calculated by solving Eq.(10). The Fourier-transformed frequency spectrum P(r, ω) can be expressed as: which is used to transform P(r, t) and P(t) into frequency domain to obtain P(r, ω) and P(ω). Subsequently, one can compute the corresponding dielectric permittivity ε(r, ω0) at site r and its spatial average ε(ω0) at frequency ω0 using Eq.(12).

Author Contributions

J.M.L. conceived the research project and P.C. performed the computations. D.P.C., Y.L.W., Y.L.X., Z.B.Y., J.G.W. and J.Y.L. commented the modeling and discussed the results. P.C. and J.M.L. wrote the paper.