From Swiss-cheese to discrete ferroelectric composites: assessing the ferroelectric butterfly shape in polarization loops

We explore the polarization hysteretic behaviour and field-dependent permittivity of ferroelectric-dielectric 2D materials formed by random dispersions of low permittivity inclusions in a ferroelectric matrix, using finite element simulations. We show how the degree of impenetrability of dielectric inclusions plays a substantial role in controlling the coercive field, remnant and saturation polarizations of the homogenized materials. The results highlight the significance of the degree of impenetrability of inclusion in tuning the effective polarization properties of such ferroelectric composites: coercive field drops significantly as percolation threshold is attained and remnant polarization decreases faster than a linear decay.

The polarization in homogeneous ferroelectric materials under the action of a changing electric field possesses a hysteretic behaviour that has been usually modelled using the phenomenological Preisach [13,[29][30][31] or Jiles-Artherton (JA) models [32].These models were initially developed for describing ferromagnetic materials and later adapted for ferroelectrics.However, they have some inconvenient restraints.The former reproduces the switching behaviour of the system as a function of the local electric field and the maximum possible local applied field induced by a structure that has to be known in advance, while the latter uses Langevin function to model the polarization, resulting in its saturated behaviour at high fields that prohibits a simulation of hysteresis loops possessing a linear increase of polarization at high fields.Besides, there have been attempts to describe hysteretic behaviour of ferroelectrics by using the hyperbolic tangent function-based approach presented by Miller and coworkers [33,34], Monte Carlo method [35,36], density functional theory [37], and phase field model [38].Though all these models have proven to be essential and reliable in accurately describing the overall response of dense ferroelectrics, they cannot be directly applied for inhomogeneous ferroelectrics, due to the complex multipolar interactions inside such structures leading to highly inhomogeneous electric fields [17,18].
Therefore, a few recent works have analysed the hysteresis behaviour of ferroelectric-dielectric composites by combining the above-mentioned models with finite element analysis (FEM), which has proven to be a powerful technique to obtain the electrical properties of materials from first principles [39][40][41].These approaches present an efficient way to correctly estimate a field-dependent polarization of a range of two-phase random heterostructures [21][22][23][24][42][43][44].Previous studies have considered 2D microstructures generated for a given inclusion content of two-phase structures composed of randomly dispersed low-permittivity inclusions in a ferroelectric matrix [21][22][23]42].The current approach builds on these original papers based on hard inclusions, but we emphasize that there is no general form that was used to generate 2D structures where a degree of impenetrability is considered.Essentially, the degree of impenetrability controls the connectedness of the inclusion phase in the host matrix and can trigger the percolation transition at the critical fraction that dramatically influences the functional properties of composites [39][40][41]45].These issues with respect to the ferroelectric composites have been largely missing in the literature so far.
With that said, we propose a simple and powerful numerical approach to characterize the electric properties and polarization hysteresis behaviour of a random ferroelectric-dielectric composite.The approach integrates the phenomenological model proposed by Miller et al [34] for the description of the field-dependent polarization of pure ferroelectrics into a FEM analysis that solves continuum-electrostatics equations for local electric fields inside a composite.We apply our combined approach for a 2D random composite generated by using the Metropolis Monte Carlo algorithm and composed of monodisperse circular dielectric disks with specified surface fractions and degree of impenetrability randomly distributed in a continuous ferroelectric matrix.The effectiveness and robustness of our approach are proven by computing the electric field-dependent effective polarization hysteresis loop and permittivity of such ferroelectric composites with different compositions.The takeaway is that our results reveal a strong impact of the degree of disks impenetrability (and thus connectivity) on the butterfly shape of the hysteresis loop and the permittivity of the composite.Furthermore, simulations unveil a critical influence of the structural parameters and percolation threshold on the coercive field, remnant and saturation polarizations of the homogenized materials.

Electrical polarization hysteresis in ferroelectrics
We first briefly review polarization in ferroelectric materials under the action of a changing electric field in order to lay the groundwork for interpreting their hysteretic behaviour.In this work, the field-dependent polarization of such switchable dielectrics is modelled using the phenomenological approach presented by Miller and coworkers [33,34].For the sake of brevity, we will name it the Miller model through the manuscript.This model has a simple set of equations and employs a hyperbolic tangent function for description of the (dipole) polarization hysteresis loop.It is capable of describing the minor hysteresis loops (without the intrinsic contribution term), but accounts only the dipole switching contribution to the total polarization.In accordance with the Miller model, the slope of the major (saturated) polarization hysteresis curve, P m , calculated as the derivative of polarization with respect to the applied field, ∂P m /∂E, at any given field must be greater than the slope of minor loops, ∂P/∂E, at any specific electric field.If the polarization P at the specific electric field is inside the major loop, then it must asymptotically approach to that of the major loop with a change in the applied field.The basic equation is where the asymptotic nature of the minor loops is modelled using the positive function Γ as where P s is the saturated polarization and the direction indicator ξ = + 1 for increasing fields change and ξ = − 1 for decreasing fields.A value of the function Γ varies in the range from 0 to 1.When the polarization P in the material at a specific applied field is significantly different from the corresponding major loop value (P > P m ) or when P approaches P s (P < P s ), the slope coefficient Γ approaches zero.The polarization P stays almost constant against the change in the field until the difference with the major curve value P m is reduced.As the polarization P approaches the major curve value P m , the coefficient Γ approaches unity.Once the polarization meets the major curve, the coefficient value becomes unity, the slope of the minor loop matches the slope of the major loop and further minor loops follow the major loop (figure 1).The major polarization hysteresis loop in the material as a function of the applied field is modelled by using the hyperbolic tangent function as χ 00 is the electric susceptibility, E c is the coercive field at which the polarization is zero, and P r is the remnant polarization at zero field (E = 0).The minus and plus signs at E c correspond to the positive-going ( + P m ) and negative-going ( - P m ) branches of the polarization hysteresis loop, respectively.The electric susceptibility χ 00 is related to the high-field relative permittivity, ε 00 , of the material as χ 00 = ε 0 (ε 00 − 1) and characterizes the linear behaviour of the hysteresis loop at very high fields, after the saturation of the hyperbolic tangent function.Figure 1 summarizes the Millers model and shows the computed major (black curve) and family of minor polarization hysteresis loops (coloured curves) for a dense ferroelectric material under different electric field amplitudes.The hysteresis terms and parameters used in equations (1)-( 4) are also indicated.
This theoretical model has been proven to accurately predict the polarization hysteresis loop of real ferroelectrics [46][47][48].It will be used for the local description of a ferroelectric material inside a composite and incorporated in our FEM simulations described in the next section.The following assumptions are made regarding this study: the domain structures are not considered; the contribution of stress/strain to the total polarization is negligible, i.e. the electromechanical strain is not considered here; the system is free of charges and other sources of fields (e.g.interfacial charges due to stress); and the model does not account the effects associated with electrodes.

Modelling concepts and numerical implementation
3.1.Generation of random two-phase microstructures Here, using a general procedure to generate two-phase ferroelectric materials we consider 2D random microstructures of equilibrium distributions of monodisperse circular dielectric disks at a specified disk surface fraction f d and degree of impenetrability λ in a ferroelectric matrix at surface fraction f f = 1 − f d , as illustrated in figure 2 [17,21,[39][40][41].In this continuum scale analysis, each disk of radius R is composed of an impenetrable core of radius λR, which is encompassed by a penetrable concentric shell of thickness (1 − λ)R.The extreme limits λ = 0 and λ = 1 correspond, respectively, to fully penetrable disks (Swiss-cheese model) and totally impenetrable disks.Thus, by continuously varying λ between 0 and 1 (i.e.decreasing the allowed overlap between disks), one can vary the exclusion-area effects and, hence, the level of connectedness of the inclusions.We employed the traditional Metropolis Monte Carlo (MC) algorithm to produce statistically homogeneous and isotropic equilibrated realizations of such composite [21].The explanation and details of the key steps of the algorithm can be found in our previous papers [39][40][41].Figure 2 demonstrates the effect of the impenetrability parameter λ and disk surface fraction f d on the generated realizations of heterostructure.Thus, the impenetrability parameter λ provides us with an additional degree of freedom for simulating a broad range of realistic composites.In particular, these configurations can represent a variety of porous composite materials as well as a large class of ferroelectric thin films fabricated experimentally [1].
Interestingly, figure 2 also reveals the formation of connected clusters of particles as disk surface fraction increases.Indeed, for a small number of disks in the unit cell (i.e.small f d ) most of the disks are well isolated from each other.Placing more disks one by one (i.e.increasing f d ) results in more disks being into close proximity with each other and, thus, the probability of overlaps between the disks increases, leading to the formation of small clusters of connected particles initially.These clusters become larger as the number of disks is increased and, eventually results at some content (f dc ) in a large cluster with size extending from one side to the opposite side of the unit cell, achieving percolation [18,45].As will be shown later, the percolation threshold significantly impacts the effective properties of our composite.Note that due to the finite size of our system, each generated random realization of the unit cell at a particular λ evidences a slight variation in percolation threshold.Therefore, the effective dielectric properties studied here show an abrupt change when the average value of percolation threshold (f dc ) is considered.Typical values of the percolation threshold at which a continuous cluster of connected disks is formed in these systems range from f dc = 0.67 (λ = 0) to f dc = 0.84 (λ = 1) [17,[39][40][41].

Model of the hysteresis behaviour of ferroelectric-dielectric media
The Miller model presented in section 2 can be directly used to describe the polarization hysteresis loop of pure ferroelectrics.However, for highly disordered inhomogeneous composites that generate complex local field fluctuations under an applied field, the use of numerical methods is essential [21].Therefore, we combine the Miller model used for the local description of the polarization of a ferroelectric phase with a FEM analysis that accurately evaluates the local electric field in a composite.In what follows, we describe a robust and effective method for calculation of the polarization characteristics, as the macroscopic polarization hysteresis loop P(E app ) and field-dependent differential permittivity ε(E app ), of the ferroelectric-dielectric heterostructures with varying composition.
The 2D composite model is analysed in a parallel plate capacitor arrangement with capacitor plates on the horizontal edges of the composite sample, as shown in figure 3 [21].The capacitor plates separated by a distance h are applied with a potential difference (across the y-axis) to simulate electrodes, j bottom = 0V at the bottom plate and j top = j app at the top plate, where ( ) , j 0 is the maximal applied voltage, ω is the angular frequency, and t is the time.The sine wave voltage is applied in order to create a time-varying electric field across the composite structure.The fringing effect at the edges is eliminated by periodic boundary condition on the left and right edges, i.e., j left = j right .The FEM analysis is performed at equal time steps to calculate the local electric field inside the composite at different values of the applied voltage and, subsequently, the corresponding local polarization employing the Miller model for each generated realization of the composite.Specifically, the local electric fields were calculated by solving continuum-electrostatics equations using COMSOL Multiphysics ® software [49] with the aid of a Java environment.For that, a differential form of Gauss law for the electric field displacement D(r) with no free electric charges was solved using the equation The constitutive relation for the low-permittivity disks is defined as where ε 0 is the permittivity of vacuum, ε d is the relative permittivity of disks, and E(r) is the local electric field calculated from the spatial distribution of the electrostatic potential j(r) inside the capacitor as E(r) = − ∇j(r).
Otherwise, the ferroelectric phase (matrix) is modelled using the constituent equation where P(r) is the local polarization vector that depends on the local electric field acting at the corresponding position r and is locally calculated using equation (1) of the Miller model.This means that the local polarization will be scattered in values in the ferroelectric phase depending on the degree of local field inhomogeneity.Due to the ferroelectric nonlinearity, the field-dependent macroscopic (effective) polarization P(E app ), where E app = − j app /h, was obtained using the iterative procedure by successively solving equation (5) under the applied sinusoidal alternating field across the y-axis of the capacitor, starting with zero applied field (E app = 0) and assuming a zero spontaneous polarization as an initial condition for the ferroelectric phase [21].The time range t was chosen to be 1.25 cycles of the sine wave applied field which produces the complete (closed) polarization hysteresis loop (one cycle of the sine curve), as shown in figure 4 (red curve).The first quarter cycle of the sine wave starting from zero applied field increases the composite polarization from zero to negative maximum corresponding to = -E E app app max , where E app max is the maximum applied field.This part of the curve is not shown in figure 4. The next complete cycle of the sinusoidal applied field from -E app max to +E app max (the lower part of the hysteresis loop) and back to -E app max (the upper part of the hysteresis loop) corresponds to the redcoloured hysteresis loop in figure 4. At each field step, the local polarization P(r) in the ferroelectric phase only is recomputed by equation (1) using the local electric field E(r) calculated at the previous field step [21].Once equation ( 5) is solved, the calculated projection of the local polarization on the applied field direction is averaged over the capacitor area to obtain the macroscopic polarization P(E app ) of the composite realization.In fact, the macroscopic P − E hysteresis loop of the ferroelectric composite is a combination of local minor and major P − E hysteresis loops formed due to the local electric field inhomogeneity.Finally, the field-dependent differential permittivity, ε(E app ) was calculated from the macroscopic polarization hysteresis loop as [

Further computational details, materials, and model validation
The composite studied in this work is a mixture consisting of ferroelectric and dielectric phases.Any simulation of the homogenized response of such composite requires knowledge of material properties of these phases.We aim to simulate a realistic composite by using experimental material values for the phases.The randomly distributed circular disks are assigned to be made of a linear dielectric material described by a relative permittivity which is arbitrarily set to ε d = 1, but any other nominal value can be chosen.For the pure ferroelectric material (matrix phase), PbZr 0.53 Ti 0.47 O 3 (PZT 53/47) is considered as an illustrating example.Its hysteresis characteristics are extracted from the experimentally measured polarization hysteresis loop provided in [47] (their figure 12) and reproduced in our figure 4 (open circles): saturated polarization P s = 21 μCcm −2 , remnant polarization P r = 13 μC cm −2 , and coercive field E c = 65 kV/cm.These values are required as inputs of the Miller model employed here.
Our approach for the calculation of polarization hysteresis behaviour of the random ferroelectric-dielectric composite described above is summarised as follows [21].First, random configurations of the composite for a specified disk surface fraction f d and degree of impenetrability λ are generated.In all calculations, we use a square unit cell of size L = 1cm and the disks radius R = 0.04 cm.The surface fraction occupied by the dielectric disks f d varies in this study between 0 and 0.9.The typical generated realizations have from ∼20 to ∼550 disks in the unit cell depending on values of the disk surface fraction, degree of impenetrability, and statistical distribution of disks in the unit cell.Second, the field-dependent macroscopic polarization P(E app ) and differential permittivity ε(E app ) for each generated realization of a composite are calculated by combining the Miller model with FEM analysis.These characteristics are evaluated at different values of the applied external field in the range from -E app max to +E app max with = E 255 kV app max /cm, j 0 = 255 kV, and ω = 6.3 rad/s.Finally, the field-dependent macroscopic polarization and differential permittivity are averaged over the results of 200 statistically independent realizations of the composite at specified f d and λ.
To verify the correctness of the implementation, the polarization hysteresis loop of a pure PZT 53/47 ferroelectric material (f d = 0) was numerically calculated by employing the Miller-FEM model presented above and hysteresis characteristics extracted from the experimentally measured polarization loop.The saturated permittivity was chosen to be ε f = ε 00 = 600 that corresponds to the electric susceptibility χ 00 = 5.3 × 10 −9 F/ m. Figure 4 demonstrates that the numerical (solid curve) and experimental (open circles) hysteresis loops are in excellent agreement with each other, proving the validity of the developed model.It should be noted here that the resulted coercive field of the numerical hysteresis loop (E c = 53 kV/cm) is slightly lower than that extracted from the experimental hysteresis loop (E c = 65 kV/cm).This is mainly owing to the asymmetry of the experimental hysteresis curve with respect to the polarization axis (y-axis), as can be clearly observed in figure 4 (open circles).As another comment, we would like to note that the Miller-FEM model is required here since the shape of hysteresis curve at each position of the two-phase random composite depends on the local electric field that is inhomogeneous inside the composite.

Results and discussion
Our results are discussed from two perspectives.First, we discuss the impact of the electric field on the shape of the hysteresis loop as a function of the disk surface fraction f d and degree of impenetrability λ.Secondly, we are interested in the coercive field, remnant and saturation polarizations of the ferroelectric-dielectric composite.Figure 5 presents the effective polarization hysteresis loops, P(E app ), and the field-dependent differential permittivity, ε(E app ), at selected values of the disk surface fraction and degree of impenetrability.For all λ, the magnitude of polarization and hysteresis area gradually decrease as well as the hysteresis loop becomes more tilted toward the E app axis as f d increases.This is due to the reduction of the ferroelectric material in the composite and the presence of the additional depolarization effect caused by a random dispersion of the dielectric phase, leading to the reduction of the electric field in the ferroelectric phase [12,27,50].Both P(E app ) and ε(E app ) are very similar for all values of f d when the degree of impenetrability λ < 0.6.The degree of impenetrability of dielectric inclusions λ has a significant effect on the shape and area of the polarization hysteresis loop but this effect is more pronounced for λ 0.6 and f d 0.4.For an arbitrary value of f d , both the magnitude of polarization and hysteresis loop area of the composite gradually increase with increasing λ.Furthermore, the shape of hysteresis loop changes by getting a more visible upright stance when λ is increased.Thus, the overlapping between inclusions significantly influences the local electric field, and consequently the polarization inside the composite.It turns out that the averaged polarization in the ferroelectric phase is larger for the structures with larger values of λ possessing the smaller allowed overlap between disks.This is because the probability of formation of connected clusters of disks is higher in structures with smaller λ and, as will be shown later, the occurrence of such clusters has for effect to reduce the local electric field inside the ferroelectric matrix.For the same reason, the field-dependent differential permittivity, ε(E app ) increases when increasing λ for an arbitrary value of f d (figure 5(b)).For all cases, the butterfly shape of the differential permittivity remains symmetric [3].Importantly, a notably different behaviour of P(E app ) and ε(E app ) is observed below and above a specific value f dc that depends on λ (f dc > 0.6).For all structures with f d < f dc , the effective polarization and permittivity monotonically reduce with increasing f d , while for f d > f dc both P(E app ) and ε(E app ) fall sharply with a different degree depending on λ and the former possesses a tendency towards linearization.This abrupt change in behaviour is an evident manifestation of a critical transition in the morphology of composite related to the appearance of a physically infinite cluster of connected disks spanning the unit cell, i.e. the percolation threshold f dc .For these 2D composites, it was reported to be f dc = 0.67, 0.71, 0.76, 0.84 for the degree of impenetrability λ = 0, 0.7, 0.9, 1, correspondingly, that is consistent with our observations [17,40,45].Actually, the dramatic impact of percolation on the properties of ferroelectric composites is often disregarded in many studies.
To explain observations above and gain further insight into the effect of disk penetrability and percolation, we carried out a comparative study of the spatial electric field for different configurations of the composite taken at E app = 320 V/cm (figure 6).For a small surface fraction (f d = 0.1), the majority of circular-shaped disk are well isolated from each other and exhibit much higher electric field values than the ferroelectric matrix.The field strength inside the particular disk is quite homogeneous and mainly governed by the permittivity of the ferroelectric and dielectric phases if there are no other disks located in close proximity.The intensity of the electric field is also very high on the ferroelectric side close to disk-matrix interfaces especially at the parts of interfaces oriented in the direction perpendicular to the applied field.Far enough from disks, the electric field intensity of the ferroelectric matrix is quite uniform and comparable with the value of the applied external field [42].With increasing the surface fraction, more disks come into close contact with each other, resulting in strong interactions between them manifested in a high field concentration within the host matrix in between the disks.When two or more particles form a disconnected or connected chainlike structure oriented in the direction roughly perpendicular to the applied field, the electric field inside disks and outside disks along the chain dramatically increases.This effect is much stronger for the structures with smaller values of λ where the probability of formation of the fully connected chains is higher.This behaviour is accompanied by a decrease of the electric field in the overall ferroelectric phase as the surface fraction of disks increases.Besides, the inhomogeneity degree of the local electric field in disks and ferroelectric matrix increases.Eventually, the first physically connected chain of disks spanning the composite from the left to the right side is formed at the surface fraction near the percolation threshold f dc (depending on λ), breaking in such a way a continuity of the ferroelectric matrix (see white arrows in figure 6 for λ = 0.5 and f d = 0.7).This leads to a very large enhancement of the electric field inside such chain and a significant redaction of the field intensity in the ferroelectric phase.In this case, a high variation of the electric potential difference tends to be located in the lowpermittivity chainlike cluster (causing the high electric field) and avoids the high-permittivity regions of the ferroelectric matrix (high concentration of isopotential lines in disks, not shown here).This explains the abrupt linearization and decrease of the polarization above the percolation threshold demonstrated in figure 5(a).
So far, we have discussed the impact of the surface fraction and connectivity of disks on the polarization, permittivity, and associated electric fields.Now, we take a closer look at the dependence of the coercive field, remnant polarization, and saturation polarization with respect to the impenetrability parameter, as displayed in figure 7.For the entire range of λ values explored, we observe that the coercivity gradually increases from 53 kV/ cm to the maximal value of 65.6 kV/cm for the surface fraction well below the percolation threshold.Though most of the curves are quite close to each other, we can assume that the higher degree of disk penetrability (i.e.small λ), the higher value of the coercive field for a specific f d .In other words, the ferroelectric phase in structures with smaller λ is overall subjected to the lower field that, in turn, requires a higher applied electric field for domain switching and saturation.Over the last two decades, the impact of a dielectric phase on the effective coercive field of composites has been debated.In agreement with our results, increasing (or slightly increasing) E c with increasing porosity in ferroelectric composites has been observed in many experimental [27,50,51] and numerical [13,27] investigations.For example, Wersing et al [52] experimentally showed a nonlinear increase of E c as a function of porosity volume fraction in a porous PZT ceramic.However, other works demonstrated either no change [53,54] or a small decrease [14,55,56] of the coercive field as a function of porosity; whereas some studies reported either a decrease or increase of E c depending on size [57] and shape [9] of pores.Interestingly, near the percolation threshold, the coercivity becomes more nonlinear and, eventually, followed by a sharp decrease for f d > f dc , an important effect that must be always taken into account in interpreting results but is rather largely ignored in the description of the ferroelectric properties.d [50].The grey area corresponds to the percolation thresholds: λ = 0 (f dc = 0.67) and λ = 1 (f dc = 0.84) [17,[39][40][41].
The decrease of the remnant polarization with increasing disk surface fraction is consistent with many experimental observations [14,53,55,58].It has been also argued that the remnant polarization of a porous material reduces disproportionately with the reduction of the ferroelectric material density expressed by a simple linear equation ( ) is the remnant polarization of the dense material [24,50].This is in accordance with the results shown in figure 7(b).Besides, our results unveil a highly nonlinear behaviour of the remnant polarization at concentrations of disks near the percolation threshold and its sharp decrease above f dc .We believe our results complement nicely those described by Zhang and coworkers [24].A clue to the interpretation of their results is based on a fit of experimental data of the remnant polarization by considering the following form ( ) ( ) , where A denotes the depolarization factor describing the shape of the inclusion.However, this linear form does not explain the nonlinearity of P r as well as the distinctive λ-dependent polarization with increasing values of the dielectric disk surface fraction.Taking a nonlinear ansatz for the normalized remnant polarization ( ) ( ) 2 leads to better fits for the entire range of f d values as presented in figure 8, but it remains to explain more precisely the physical significance of the fit parameters.In any event, it would appear that the details thereof have not been disentangled.One has no rigorous calculation of the characteristic length scale on which polarization ought to proceed.The above linear form of P r (f d ) has another questionable aspect since it relies on simple proportionality to the amount of ferroelectric material in the system, whereas we know that polarization is an intensive variable [59][60][61].Additionally, it remains unclear how a strong variance of the depolarization factor of the inclusions can be consistent with what one would expect from an electrostatics perspective [62].The computation of A relies heavily on the geometrical features of the dielectric phase: this requires an understanding of the local electric field distribution in the complex morphologies shown in figure 2   Figure 7(c) illustrates one more noteworthy property of the polarization: the saturation polarization follows the same decay variation versus f d as the remnant polarization [13,24].The rectangularity of the hysteresis loop, P r /P s , continuously reduces with an increase of the surface fraction below f dc (figure 7(d)).For a specific f d , the rectangularity is smaller for the smaller values of λ, namely, for disks with a higher degree of penetrability.This is consistent with figure 5(a) demonstrating the increase in the tilting degree of the hysteresis loop when increasing f d and decreasing λ.For the values of f d > f dc , corresponding to the grey area, the rectangularity strongly decreases for penetrable disks (λ < 1) in sharp contrast to the case of totally impenetrable disks (λ = 1).Finally, we note a distinct behaviour of the non-percolating structure (impenetrable disks) with parameters λ = 1 and f d = 0.9 that corresponds to a periodic hexagonal array of disks in a ferroelectric matrix (see also the corresponding electric field in figure 6).
Overall, figures 5 and 7 reveal non-trivial composition-structure-property relationships in the ferroelectricdielectric composite.The coercive field, remnant and saturation polarizations demonstrate a nonlinear behaviour that is more pronounced at higher contents of the dielectric phase and, eventually, a critical transition associated with the percolation threshold.Furthermore, these ferroelectric properties depend on the degree of disk connectivity and associated clustering characteristics of microstructures.All these effects can alter the form of an effective polarization significantly.Finally, our results provide evidence on the prospects of utilizing the connectivity and percolation properties for tailoring different combinations of the reduced permittivity and other ferroelectric characteristics for tunability applications.

Final remarks
At the end, it is worth noting several comments concerning this study.Firstly, we have considered a composite material having intrinsic permittivity of disks ε d = 1, although our approach is robust to any real or complex value of disk permittivity.Secondly, though the generated microstructures are composed of randomly distributed circular disks, by varying the phase concentration and degree of disk penetrability, one can vary the connectedness of the disks and, hence, the disk clustering characteristics, and degree of phases connectivity.The interaction between such diverse particles within the host matrix involves complicated interfacial phenomena.Thirdly, we note that our results presented for 2D composites, are also directly valid for their 3D counterparts consisting of parallel infinitely long, partially penetrable cylindrical disks randomly distributed throughout a ferroelectric matrix.Moreover, any 2D material can be also considered as a 3D material with a lateral dimension much smaller than the other dimensions.However, we would like to comment on a distinct behaviour between 2D and 3D composites.Ferroelectric properties of our 2D random composite, such as the polarization, permittivity, and coercivity, show an abrupt decrease above some critical concentration of inclusions.However, these properties would decrease more smoothly up to the high volume fractions of inclusions (spheres) for 3D ferroelectric composites, despite having percolation threshold at considerably lower volume fractions than that in 2D composites.This is mainly because the ferroelectric matrix becomes discontinuous when the dielectric phase is sample-spanning (percolating) in 2D composites, while 3D composites above the percolation threshold can be bicontinuous, allowing high electric field concentration in both the inclusion and ferroelectric phases, in contrast to 2D structures.This can explain the differences in behaviour of 2D and 3D random composites.

Conclusions
This study adds an important dimension to our understanding of irregularly shaped low permittivity disks described by a degree of impenetrability by revealing the role of connectedness on the effective polarization characteristics of ferroelectric composites.We found that the effect of disk concentration can be subtle, with large deviations in the coercive field, and remnant and saturation polarizations compared to homogenized dense materials: coercive field drops significantly as percolation threshold is attained and remnant polarization decreases faster than a linear law.We demonstrate these deviations in careful computation of the effective permittivity for 2D two-phase dielectric-ferroelectric composites.These results illustrate that any effective modelling of the polarization properties of such composites behaves very differently than one might naively assume, particularly for large surface fraction, when a degree of impenetrability of inclusion forming the twophase composite is specified.Nonetheless it remains an interesting problem to solve the case of 3D dielectricferroelectric composite materials and allow for polarization space charges [62,63] as well, getting a more complete picture of polarization loops.A recurrent theme in ferroelectricity is related to local symmetry breaking.One way our analysis could apply in that case is considering strain gradients in ferroelectric materials containing anisotropic inhomogeneities.This is a critical question that we would like to see checked in future work.We anticipate that such issue is also challenging for an abundance of applications.The inclusion of porosity in ferroelectric thin films would enhance the surface area and reactivity, leading to a potential improvement of the photoelectrochemical performance [66].Due to their nonlinear response to an electric field, dielectric-ferroelectric materials are also of great interest for microwave devices such as tunable capacitors, controlled filters, and phase shifters for a new generation of antennas.

Figure 1 .
Figure 1.Major (saturated) polarization hysteresis loop (black curve) and minor (unsaturated) hysteresis loops (coloured curves) evolving under different applied fields generated by the Millers model for a dense ferroelectric.Arrows indicate the direction taken around the loops.

Figure 2 .
Figure 2. Typical equilibrium sample realizations illustrating the local microstructure of the two-phase composite consisting of monodisperse circular partially penetrable dielectric disks, at a specified disk surface fraction f d (0 f d 1) and degree of impenetrability λ (0 λ 1), randomly distributed within a square matrix (ferroelectric material).

Figure 3 .
Figure 3. Schematic illustration of the numerical setup for the calculation of a polarization hysteresis loop for a typical realization of a random ferroelectric-dielectric composite.A parallel plate capacitor is filled with a composite consisting of monodisperse circular disks, of linear dielectric permittivity ε d , surface fraction f d and with a degree of impenetrability λ (0 λ 1), randomly dispersed within the nonlinear ferroelectric host material of field-dependent permittivity ε f and surface fraction f f = 1 − f d .Each disk of diameter D has a mutually impenetrable core of diameter λD (region inside the inner black circle) encompassed by a perfectly penetrable shell of thickness (1 − λ)D/2.A sinusoidal input voltage with a frequency ω is applied on capacitor plates.Periodic boundary conditions are imposed on the left and right edges.

Figure 4 .
Figure 4. Polarization hysteresis loop of the PZT 53/47 ferroelectric ceramic simulated by using the Miller-FEM model (solid red curve) and hysteresis characteristics extracted from the experimentally measured polarization loop taken from [47] (black circles).

Figure 5 .
Figure 5. (a) Effective polarization-field hysteresis loops P(E app ) simulated by Miller-FEM model as a function of disk surface fraction f d and impenetrability parameter λ of the ferroelectric composite.(b) Effective differential permittivity ε(E app ) as a function of electric field reconstructed from the curves of hysteresis loops in (a).

Figure 6 .
Figure 6.Illustration of the spatial electric field distribution |E| at E app = 3.2 × 10 4 V/m for the different configurations shown in figure 2. White arrows in the map for λ = 0.5 and f d = 0.7 indicate the percolating cluster of disks spanning the composite from the left to the right side.The colour bar denotes the intensity of the electric field in V/m.

Figure 7 .
Figure 7. (a) Coercive field, (b) remnant polarization, (c) saturation polarization, and (d) rectangularity as functions of surface fraction f d and impenetrability parameter λ extracted from the curves of hysteresis loops P(E app ) shown in figure 5(a).The inset in (a) shows details of the data.The dashed line in (b) represents a linear dependence of the form ( ) ( ) f f = -P P 1 r d r 0 [63][64][65].

Figure 8 .
Figure 8.A fit (red curve) of numerical data of the normalized remnant polarization (blue dots) by considering a nonlinear ansatz of the form ( ) ( ) f f f ~--P a b 1 r d d d 2 as functions of surface fraction f d and impenetrability parameter λ.The values of parameters a and b for each λ are shown in the corresponding inset. 47,48]