Introduction

Hysteresis behavior of ferroelectric (FE) and antiferroelectric (AFE) crystals can be a footprint in evaluating the stability of the spontaneous dipoles and their possible configuration. A significant number of studies are published focusing on the phase stabilities of AFEs and AFE coupling at interfaces of multicomponent systems through the hysteresis shapes they exhibit as well as coexistence of the FE and AFE phases in a single composition [112]. In a practical sense, multilayers of FE and AFE components such as PbZr(1−x)TixO3–PbZrO3 (PZT–PZ) or PbTiO3–PbZrO3 (PT–PZ) have attracted interest as these structures were shown to exhibit high dielectric constants for critical compositional frequencies when in multilayer form [1317]. So, multilayers comprising FE and AFE components are gripping both from application to scientific point of views.

Theoretical studies that try to predict the hysteresis dynamics via adjustable parameters in the Hamiltonian also focus on the exchange coefficients and calculation of the dielectric anomalies for a set of chosen parameters and dipole arrangements [47, 9, 10]. Switching behavior and domain contributions to dynamics of FE films and superlattices have already been the focus of several research groups, highlighting the importance of interfaces [1822]. Prior to recent interest in FE films, studies have been extensively carried out for magnetic systems by applying the 2D and 3D Ising Model to these materials [23, 24], including attempts to select realistic interaction parameters for a Ising-type Hamiltonian extracted from experimental data [25]. Transverse Ising Model (TIM) [26] has been employed to some types of AFEs where an internal transverse field is present due to the distribution probability of a proton between the neighboring lattice sites. Moreover, a strong antiparallel exchange between sideways neighbors also produces well-defined AFE loops in the nearest neighbor 2D Ising limit [27], similar to FE domains or layers interacting with one another through a negative exchange coefficient at the domain interfaces [8]. With increasing temperature toward the Neel point, the nearly linear neck connecting the two FE parts of the loop gradually disappears. Before total destruction of the antiparallel configuration of the dipoles, loops reminiscent of a FE appear which again gradually disappears with increasing temperature-induced fluctuations. In theory, adjustable interaction parameters of a Hamiltonian could surely be shown to give rise to many types of hysteresis and phase transition behaviors to comment on real experimental observations. These studies are vital to assess and understand the internal energetics and competing mechanisms in systems undergoing order–disorder transitions, especially when departure from equilibrium is enforced due to the presence of an externally varying field. In spite of the interest in such systems, phase stabilities of AFEs sandwiched between layers of FE in the form of multilayers, have not been investigated in depth despite experimental studies reporting interesting results.

In this article, we carried out Monte–Carlo (MC) simulations on a 2D grid where we defined short-range and long-range interactions between lattice sites. The evolution of the grid was studied for cases of a pure FE, pure AFE, a bilayer, and a superlattice consisting of equal fractions of FE and AFE both with temperature and under applied field. Cooling runs under zero external field yielded information on the equilibrium order state of the systems as a function of temperature, emphasizing the great impact of the long-range dipolar interactions. By adjusting the strength of the short-range-to-long-range interactions and electrostatic-like coupling, various different configurations of the grid were obtained yielding various hysteresis loops under a triangular applied field with fixed frequency in all cases. During the hysteresis simulations of the bilayer and the superlattice, the effect of an electrostatic-like coupling between the layers that each spin feels was also studied. We include this coupling term with the reservation that it is not exactly corresponding to the behavior of real dipoles but is rather an energy term that has to be minimized in the simulations. The system’s strong tendency to minimize this term gives rise to very similar results obtained within the scope of the Landau–Ginzburg–Devonshire (LGD) theory that employs the electrostatic coupling in the presence of space dependent polarization variations. The depolarization effect either due to a dead layer near the film–electrode interfaces or due to imperfect screening of charges created by the spontaneous dipoles in the FE have been shown to be capable of suppressing ferroelectricity in very thin layers [28, 29]. In multilayers, the presence of such a term could then favor the FE layers to exist in a polydomain state where a similar effect was accounted for in our 2D hypothetical grids. With the inclusion of a depolarization-like effect, up-spins and down-spins with nearly equal fractions, reminiscent of 180° domains in real FEs, form in an alternating columnar fashion. Such a situation, of course, occurs when the strength of the depolarization term exceeds a critical value where it becomes energetically favorable to form strip-like domains at the expense of the domain-wall energy. The terms in the dipole–dipole (d–d) interaction that favor parallel alignment of spins with long-range coupling become a dominant parameter both in equilibrium configuration and switching under applied signal especially at low temperatures. For nearly equally stable FE and AFE layers, that is their equilibrium order is destroyed at approximately the same temperature, superlattices comprised of such layers with several interfaces have FE-like hysteresis loops while the bilayers display hysteresis loops that can be deconvoluted to a separate FE and an AFE component. We demonstrate that, in addition to different strain states of layers in superlattices, interfaces are just as important as regions exposing the components to one another where even a short-range penetration of one type of order of a component into the other can alter the equilibrium and dynamic phenomena.

Theory and methodology

A 2D system consisting of sites that are strongly correlated to each other can be expressed within the 2D Ising limit as:

$$ H_{\text{SR}} = - J_{\text{HT}} \sum\limits_{1}^{{{N \mathord{\left/ {\vphantom {N 2}} \right. \kern-\nulldelimiterspace} 2}}} {S_{\text{i}} S_{{{\text{i}} \pm 1}} } - J_{\text{SW}} \sum\limits_{1}^{{{N \mathord{\left/ {\vphantom {N 2}} \right. \kern-\nulldelimiterspace} 2}}} {S_{\text{j}} S_{{{\text{j}} \pm 1}} } $$
(1)

where only nearest neighbor exchange is considered, subscript “SR” in HSR stands for “short-range” and JHT is the exchange coefficient for FE order between the sites along y-axis (head-to-tail) with Si = ±1/2 being the local spin of the site i and JSW is the sideways exchange coefficient between neighbor spins, imposing an antiparallel state when negative; N represents the total number of sites in the grid. HSR is sufficient to induce a long-range order in a 2D system with nearest neighbor exchange below the characteristic Curie point that is determined by the randomization effect of temperature on spin states. The Curie point of the system is basically the kT value, with k being the Boltzmann constant, and T temperature, above which probability of having aligned spin couples is around 50%, meaning disorder in the system within the current algorithm. In an ensemble of local spins (and dipoles), we should also incorporate the dipole–dipole (d–d) interaction that has a long-range nature:

$$ H_{\text{DD}} = A\sum\limits_{{{\text{i}} \ne {\text{j}}}} {\frac{{S_{\text{i}} \cdot S_{\text{j}} - 3(S_{\text{i}} \cdot n)(S_{\text{j}} \cdot n)}}{{\left| r \right|^{3} }}} $$
(2)

where A is the interaction constant and \( r = \sqrt {(x_{\text{i}} - x_{\text{j}} )^{2} + (y_{\text{i}} - y_{\text{j}} )^{2} } \) is the distance between sites i and j (\( {\text{i}} \ne {\text{j}},\,r_{\text{i}} - r_{{{\text{i}} \pm 1}} \,{\text{and}}\,r_{\text{j}} - r_{{{\text{j}} \pm 1}} \) is taken as unity), n is a unit vector directed along the line joining the two sites considered, Si · Sj and Si(j) · n are the dot products. Note that a parallel configuration is favored for head-to-tail spins while a full antiparallel state has lower energy only for sideways dipoles. The Si · Sj term in (2) can only be minimized by antiparallel alignment of dipoles for two degrees of freedom of spins in this work. We plot the magnitude of the cos2θ arising from the 3(Si · n)(Sj · n) term considering a central dipole interacting with others in a 2D grid had all dipoles been pointing up as given in Fig. 1a. The regions of zero (0) value for cos2θ are the regions where antiparallel alignment, dictated by the Si · Sj of (2) term, is dominant while one (1) favors parallel alignment. For a clear representation, Fig. 1b provides a map of favored interaction type as a function of position with respect to a dipole at the center of the map. As we will discuss in the forthcoming section, the information contained in Fig. 1a and b will prove very useful in clarifying the trends in the systems considered. For values of A comparable to JHT and JSW in magnitude, (2) has a great impact on the properties and equilibrium states of the considered systems owing to its long-range influence. The electrostatic energy of an applied field is added to the system in the form:

$$ E_{\text{EL}} = 2\mu_{\text{i}} (E_{\text{App}} - E_{0} )\sum\limits_{1}^{N} {S_{\text{i}} } $$
(3)

with EApp being the externally applied field, and

$$ E_{0} = \beta \bar{S} $$
(4)

β is a coefficient that establishes the strength of the electrostatic-like coupling, \( \bar{S} \) is the average spin of the entire system and μi being the dipole magnitude of site i. The term (4) stands for the overall electrostatic-like coupling similar to the terms used in the LGD theory of FEs. In the latter, the term scaling with \( 4\pi \bar{P} \) (in Gaussian units, \( \bar{P}/\varepsilon_{0} \) in SI units, \( \bar{P} \) is total polarization of the system) stabilizes 180° electrical domains [2830] to compensate for the internal depolarization field induced by the variation of the order parameter near interfaces and due to imperfectly screening electrodes. One must keep in mind that the depolarization term in a FE is a function of sample shape and depends on the thickness for the case of a film that is infinite in the plane. In our study, we arbitrarily adjust this term to demonstrate the effect of electrostatic coupling between the layers due to the different intrinsic order the FE and AFE layers tend to attain. The total energy of the 2D grid becomes:

$$ \bar{H} = H_{\text{SR}} + H_{\text{DD}} + E_{\text{EL}} $$
(5)

Defining a periodic structure such as a bilayer or a superlattice will clearly be through assigning alternating values of \( J_{\text{HT}}^{1,2} \) and \( J_{\text{SW}}^{1,2} \) as demonstrated schematically in Fig. 2, with superscripts 1, 2 denoting the FE and the AFE layer, respectively. The layer fractions are taken equal with interface layers assumed to have JHT, JSW = 0. Minimization of the \( \bar{H} \) was done with a MC approach where the system was allowed to evolve toward its equilibrium configuration at a given temperature, T. A grid size of 70 × 70 was constructed with free boundary conditions. Simulations were run for a variety of cases for comparison such as a pure FE grid, AFE, bilayers, and superlattices. Throughout the hysteresis simulations, the “order state” of the systems considered were tracked via the average spin value of the system given as

$$ {\text{Avr}}(S) = \frac{1}{N}\sum\limits_{1}^{N} {S_{\text{i}} } $$
(6)

A Markov chain was constructed with random spin-site selection for flipping, and the kinetic Glauber formula [31] was used to decide if sign change of a spin at a chosen site would be accepted in the form of a probability, P:

$$ P = \left\{ {\begin{array}{*{20}c} 1 & {{\text{If }}\Updelta \bar{H} < 0} \\ {\frac{1}{{2\tau_{\text{MC}} }}\left[ {1 - \tanh \left( {\frac{{\Updelta \bar{H}}}{2kT}} \right)} \right]} & {{\text{If }}\Updelta \bar{H} > 0} \\ \end{array} } \right. $$
(7)

where τMC is a time step taken as unity. \( \Updelta \bar{H} \) is the energy difference between two consecutive configurations of the system that differ by only one single flip. Average spin versus applied field hysteresis loops of the structures were obtained by applying a triangular electric field (amplitude varying from zero to EMAX where \( E_{\text{MAX}} \gg 4J_{\text{HT}}^{1,2} \)) with a total of 100 incremental steps. At each field-step the grid was allowed to relax for 20 MC steps (MCS) and the resulting configuration constituted the initial state for the next incremental applied field. One MCS represents 702 flip attempts on randomly chosen sites throughout the grid. For a fully parallel oriented system of spins, the field at which switching will occur can be found from \( \Updelta \bar{H} \) approaching to zero, meaning the sum of the first two terms of (5) will become equal to the last term for a given external field. This will result in a very high acceptance rate of spin flip attempts to reduce energy in just a few numbers of grid-sweeps until the full-parallel configuration is attained. Note that \( \Updelta \bar{H} \) attains near zero values in the AFE component when a critical field is reached followed by stabilization of a field-induced FE alignment.

Fig. 1
figure 1

a The value of \( \cos^{2} \theta \) in the dot product \( S_{{{\text{i}}\left( {\text{j}} \right)}} \cdot n \) as a function of position for a pair of interacting dipoles one of which is fixed as the central dipole in reference (bold gray) and b Map of parallel and antiparallel alignment of dipoles interacting with a central reference dipole for spin-up/spin-down degree of freedom to reduce dipole-dipole interaction energy. Both plots are for 34 × 34 sites corresponding to ±17 distance units along x-axis and y-axis around the central reference dipole (x = 0, y = 0)

Fig. 2
figure 2

The schematic of the superlattice and the bilayer grids used in this study (black: AFE, white: FE, gray: interfaces)

Results and discussion

Single component, FE, and AFE grids

Before going on to the simulations of bilayers and superlattices, we reproduced well-defined hysteresis loops for single component FE and AFE grids, resembling room temperature experiments in real systems, which we will keep as reference systems in the rest of the study (see Fig. 3). In all our simulations, we limited the maximum distance between spin couples undergoing long range d–d interactions to 8 units where each unit represents one lattice parameter. As the strength of the dipolar interaction scales with 1/83 for sites separated by r = 8, the term in (2) becomes negligible at larger r and it helps us save computational time without any sacrifice from the actual trends of the lattice grid. Although a very well-known energy contribution to both electrically and magnetically ordered systems, it is rather hard to separately judge the impact of the d–d term in FEs experimentally. Dipolar interaction is at its minimum value for a head-to-tail column of spins along y-axis, stabilizing a FE arrangement. Taking into account the degree of freedom for spins in this study and Eq. 2, we find that around 65% of the sites in the grid will have the tendency to exist in an all-parallel state at equilibrium as given in Fig. 1b. The black area corresponds to the case of d–d energy with negative values for parallel spins whereas the white area has negative energy for antiparallel spins. We extracted this ratio by summing the term Si · Sj − 3(Si · n)(Sj · n) with normalization of the distance between considered sites in (2) at each point and deciding which spin configuration minimizes the local d–d energy simply by comparing FE and AFE alignment possibilities. There will thus be a competition between the d–d term in (2) and short-range interaction that imposes an AFE phase depending on the \( J_{\text{HT}}^{\text{i}} /J_{\text{SW}}^{\text{i}} \) ratio for a layer i when approaching equilibrium. Note that the AFE hysteresis in Fig. 3 has A = JHT.

Fig. 3
figure 3

The reference loops used in the study of the bilayers and the superlattices, \( A = J_{\text{HT}} = 8kT \) in both the FE and the AFE, \( J_{\text{SW}}^{2} = - 2J_{\text{HT}}^{2} = - 16kT \)

Stability of the anti-parallel dipole configuration of an AFE within the presence of long-range, FE-favoring interactions remains somewhat an intriguing subject. For example, the very well-known double-loops of PZ [3234] corroborate the fact that the energy-minimizing mechanism that promotes the antiparallel distortions occuring in the crystal must be quite dominant over the long-range d–d term [35]. This must especially be the circumstance when the dipoles are constrained by the lattice to a certain orientation. The long, linear neck in between the two field-induced loops published for PZ in several studies support this argument [3234]. In a most basic view, the long-range d–d energy mostly favors FE-order except for immediate sideway dipoles in the system regardless of the degree of freedom for the dipoles at low fields. That the AFE-to-FE transition requires quite high fields also signals the strength of the AFE favoring mechanism. Looking at Fig. 1, the full antiparallel alignment of dipoles is only favorable for sideways neighbors within the limit of (2) when 3(Si · n)(Sj · n) = 0 while Si · Sj has its maximum.

Within the scope of a Hamiltonian constructed around a short-range exchange and a long-range d–d term as often encountered in literature, we particularly conclude that the sideways exchange must be the dominant contribution in stabilization of an AFE phase. Field-induced transition of the loops from the AFE to the FE also occurs at fields that are comparable or larger than the coercive fields in most FE systems. Hence, the values chosen for JSW were also adjusted accordingly where a clear double-loop AFE hysteresis was obtained (Fig. 3). Following this short discussion, one should also consider the impact of strain on the stability of such systems in addition to the intrinsic terms. Nearly all hysteresis data published for PZ were acquired in relatively thick films, at the order of a few hundred nanometers. It is clear that such structures will undergo a misfit strain relaxation on misfitting substrates. Regardless of whether epitaxial or polycrystalline, very distinct AFE loops were observed for PZ, indicating that the AFE phase can be stabilized in thin films with varying crystal orientations, a situation that one would not strongly anticipate. In a recent report published by our group, such clear AFE loops tend to disappear when epitaxially grown PZ is relatively thin and is in the form of layers sandwiched between PZT 80/20 layers [17]. For the extreme case of a many-layer superlattice consisting of nearly equal fractions of PZT80/20 and PZ, no trace of an AFE behavior was observed, the possible reasons for what are discussed in Sects. “Phase transition behavior of the bilayer and superlattice grids” and “Hysteresis loops of the superlattice and the bilayer”. In our simulations we qualitatively observed the same trend in the comparison of the hysteresis of bilayers and superlattices without altering any of the exchange and d–d coefficients.

Phase transition behavior of the bilayer and superlattice grids

Using the systems whose hysteresis curves are given in Fig. 3 as components comprising the layers, we created the bilayer and the superlattice grids whose schematic were already given in Fig. 2. To shed light on the phase transition behavior and the Curie points of the bilayer and the superlattice with \( J_{\text{SW}}^{1} = J_{\text{HT}}^{1} = J_{\text{HT}}^{2} \) and \( J_{\text{SW}}^{2} = - 2J_{\text{HT}}^{2} \) with \( J_{HT}^{1} ,J_{HT}^{2} = 8kT \) and \( A = J_{\text{HT}}^{1,2} \), we carried out cooling runs. We cooled both the superlattice and the bilayer slowly starting from 2.5 kT/JHT to 0 kT/JHT where each system was kept at chosen temperature intervals for 200 MCS and the same procedure was repeated four times for statistical integrity. We essentially noted that such a relaxation allowance at each temperature is sufficient as the systems reach a level-off and do not evolve into further configurations especially at low kT/JHT. This picture, of course, changes with increasing kT/JHT where thermal fluctuations cause a large variation of average spin values as expected. The results of the cooling runs are plotted in Fig. 4. A very interesting behavior is displayed where the superlattice has a strong FE order at very low temperatures that drops with a quite steep slope toward net zero spin value. The bilayer, on the other hand, apparently has a much higher Curie point.

Fig. 4
figure 4

Average spin as a function of temperature (cooling curves). Points 1 and 2 denoted on the dashed vertical line at 0.85 kT/JHT are the configurations of the bilayer and the superlattice after 200 MCS, given under the plot. The dashed curves are guides for the eye. \( A = J_{\text{HT}} = 8kT \) in this figure

To check whether the superlattice indeed undergoes a phase transition or if the apparent disappearance of the net spin value is a consequence of the spin clustering, we give the configurations of the grid at the kT/JHT where net spin in the systems are denoted by 1 and 2 in Fig. 4. Following the relaxation of both systems at 1500 K, we note that the net average spin of the superlattice system approaches zero not due to a phase transition (total destruction of the order with thermal fluctuations) but due to the domain-like formations in the FE layers. It is quite straightforward to explain this behavior as \( J_{\text{SW}}^{2} = - 2J_{\text{HT}}^{2} \) and that \( J_{\text{HT}}^{1} ,\,J_{\text{HT}}^{2} ,\,A \cong kT \), meaning that thermal fluctuations are dominating and the all-parallel alignment in the FE layers is somewhat destroyed under the influence of the more stable anti-parallel configuration in the AFE layers. However, just the opposite trend occurs at low temperatures where the long-range FE ordering prevails in the entire system despite \( J_{\text{SW}}^{2} = - 2J_{\text{HT}}^{2} \). The antiparallel spin arrangement in the AFE layer of the bilayer sustains stability at low temperatures for \( J_{\text{SW}}^{2} = - 2J_{\text{HT}}^{2} \) and \( A = J_{\text{HT}}^{1,2} \). Values of \( A > J_{\text{HT}}^{1,2} \) can still permit an antiparallel arrangement in the AFE of the bilayer due to the fact that about 35% of spins will exist in an antiparallel state owing to the form of (2) and cosθ in the second dot product of (2). The latter are true for the bilayer as there is just one interface with the FE component leading to less influence of the FE–AFE interaction through the d–d energy.

One also must note that these are generalized discussions and that in real samples the layers are often under different strains with probably quite different phase transition behavior. In this article, we show that a phase transition in a multilayer system may not be strictly or only related to possible different strain states of the layers. For the sake of demonstration, we give the cooling curves of a superlattice and a bilayer excluding the long-range interactions in Fig. 5. The absence of the long-range dipolar term significantly impacts the Curie point of the two systems, and the interfaces in the superlattice give rise to a decrease in the net spin as these regions were considered “transition regions \( \left( {J_{\text{HT}}^{1,2} = 0,\,J_{\text{SW}}^{1,2} = 0} \right) \)” in between the layers. Still, we should add here that the interface susceptibilities remain somewhat insignificant compared to the d–d interaction unless the number of interfaces approach that of the individual layers.

Fig. 5
figure 5

Average spin as a function of temperature (cooling curves) in the same temperature range when \( A = 0 \). The dashed curves are guides for the eye. Note that the bilayer and the superlattice have nearly the same Curie point

Hysteresis loops of the superlattice and the bilayer

In principle, the cooling curves already provide the evidence that at low temperatures away from the Curie point, the FE arrangement is dominant. To see how the switching of grids occurs under applied field, we give the hysteresis of the bilayer and the superlattice at \( {{kT} \mathord{\left/ {\vphantom {{kT} {J_{\text{HT}} }}} \right. \kern-\nulldelimiterspace} {J_{\text{HT}} }} = 0.125 \) in Fig. 6. The FE layer has JHT = JSW and the AFE is characterized by JSW = −2JHT in Fig. 6a and the same constants and applied signal frequency and amplitude were used to get the plots in Fig. 3 were employed. The two structures behave very differently where the bilayer exhibits both FE and AFE switching while both tend to merge into one loop with decreasing \( J_{\text{SW}}^{2} \) (in AFE layer). Clearly, the magnitude of |JSW| determines the stability of the AFE layer. The single component AFE hysteresis loop for JSW = −JHT is given in Fig. 6d. There is still a very clear double loop with a smaller AFE-to-FE transition field compared to the one given in Fig. 3. The loop for a bilayer and a superlattice using JHT = JSW for the FE and JSW = −JHT for the AFE is in Fig. 6b. These loops evident that near all-parallel alignment is taking place, with the exception of some AFE clusters still remaining in the AFE part of the bilayer that switch to all-parallel alignment at a normalized field of about 0.3. In the AFE model of Kittel and Cross [3638], the variation of JSW identically corresponds to the adjustments of the term R in RPa Pb, where Pa Pb is the product of the sublattice polarizations, which is determining the strength of the “antiparallel remanence” of spins until an electrostatic energy overcomes this barrier to induce FE alignment via an externally applied field.

Fig. 6
figure 6

ac Various hysteresis loops of a bilayer (shaded) and a superlattice (solid black line) for the given Hamiltonian coefficients obtained using a triangular field signal at fixed frequency. Note that the loops tend to merge into a single one in the superlattice with decreasing \( \left| {J_{\text{SW}} } \right| \) . A = JHT in all plots. d The single component AFE loop when \( J_{\text{SW}}^{1} = J_{\text{HT}}^{1} ,\,J_{\text{SW}}^{2} = - J_{\text{HT}}^{2} \) (solid black line). The AFE loop previously provided in Fig. 3 is given for comparison (shaded)

Overall, it is clear that the field-induced AFE-to-FE transition in the sandwiched AFE layers occurs at lower applied fields than in a pure AFE structure, a consequence of the 3(Si · n)(Sj · n) in (2) term acting to impose all-parallel alignment for angles higher than π/6 (with respect to the horizontal axis) between interacting spin couples near the interfaces. Also note that β is taken as zero until now and cases where β > 0 will be discussed in the forthcoming section. In literature, very wide range of hysteresis shapes have often been reproduced in theoretical studies but mostly by varying the Hamiltonian parameters for bilayers and multilayers. We here demonstrate that structural or compositional periodicity is just as important, particularly when long-range interactions are taken into account. Therefore, through more exposure of the layers to each other at the interfaces, different phase transition behavior can be exhibited, excluding any strain-related arguments or varying the Hamiltonian parameters. This is in good qualitative agreement with experiments where FE-like behavior has often been encountered in real multilayer structures with high periodicity. For bilayers, the regions that are far from the interface tend to behave rather independently, reflected in the extensions of the hysteresis loops but this happens at low fields due to the assistance provided by the already switched spins, especially those in the FE layer. Not exactly knowing if the very same mechanism is the reason, we had observed loops of AFE character in epitaxial FE–AFE bilayers reported in one of our recent articles [17].

For a complementary view and in order to display the behavior of the systems in weak d–d interaction energy, we provide the hysteresis runs of the bilayer and the superlattice with \( A = 0.5J_{\text{HT}}^{1,2} \) in Fig. 7. With decreasing A, the layers start switching independently and it turns out that the superlattice and the bilayer have identical loops when A = 0.5 (Fig. 7). It is very fortunate that this picture is in total contradiction with the one provided in Fig. 6a especially where the superlattice switches with significantly merged double loops at both applied field polarities. The components in the bilayer of Fig. 6b, however, switch independently in a relative sense but not as in a manner profound as in Fig. 7.

Fig. 7
figure 7

Hysteresis loop when \( J_{\text{HT}}^{1,2} = J_{\text{SW}}^{1} = 8kT,\,J_{\text{SW}}^{2} = - 2J_{\text{HT}}^{2} \,A = 0.5J_{\text{HT}}^{1,2} \)

Effect of a depolarization-like field (β > 0)

Throughout our simulations, the superlattice has a significantly larger remanence at zero applied fields than the bilayer but one additional remark we would like to emphasize is the effect of an electrostatic-like coupling on the hysteresis loops of FE–AFE layers. To account for a depolarization-field effect, we assumed an arbitrary \( \beta = 5J_{\text{HT}}^{1,2} \). Inclusion of this term into the hysteresis simulations leads to slimmer and tilted hysteresis loops as it reduces the effective field that the spins are exposed to in addition to the domain formation to minimize the related term \( \beta \bar{S} \) when EApp < E0 in (3) (See Fig. 8). Thus, switching occurs gradually in a range of applied field values, i.e., in a rather diffuse fashion. Another prominent effect is the loss of the remanence near zero-field as 180° domains start to form in the FE layer in the form of spin clusters with flat interfaces (at low enough temperatures). For very large β, there is no hysteresis but just a linear response for all systems without any apparent “remnant order parameter” where the applied field only changes the up-spin domain/down-spin domain fraction. A similar result was recently reported using the time-dependent Landau–Ginzburg equation for BaTiO3 thin films with thick dead layers (corresponding to strong depolarization field) by Ahluwalia and Srolovitz [39].

Fig. 8
figure 8

Hysteresis loop when \( A = J_{\text{HT}}^{1,2} = J_{\text{SW}}^{1} = 8kT \) and \( J_{\text{SW}}^{2} = - J_{\text{HT}}^{2} \) with \( \beta = 5J_{\text{HT}}^{1,2} \) (electrostatic-like coupling)

Interestingly, the superlattice grid hysteresis in Fig. 8 has a qualitatively very similar shape to a recently reported result for an epitaxial bilayer [17]. As mentioned earlier, the superlattice grids in this study have a high FE–AFE interface-to-volume ratio where the components feel each other’s presence. Such a scenario could also be true when one thinks of coexistence of the FE and AFE phases in the same layer and the hysteresis will exhibit both characteristics [40]. Considering that the bilayers are relatively more relaxed than the superlattices and the presence of just one interface [1417], it is possible to expect that the components of the bilayer will display a relatively independent behavior. The effect of varying strain in the layers in this study can be incorporated by choosing appropriate coefficients in the Hamiltonian and we would like to state that no double loops are of consequence when a strong FE (high Curie point) and a weak AFE (low Neel point) are thought to comprise the grid in the presence of d–d interactions. This is not so when d–d interactions are excluded in the simulations, thereby indicating the large impact of the d–d term in (2) in addition to any possible strain arguments.

Zero-field near-equilibirum configuration for β > 0

In order to examine the conditions stabilizing a FE or AFE type ordering under electrostatic-like coupling, we also carried out zero-applied field runs where we observed the evolution of the grid under the influence of β and A while kT/JHT = 0.125 is fixed. Each run had 1000 MCS. In Fig. 9, to emphasize the contribution of the long-range dipolar energy, we show the case when \( A = J_{\text{HT}}^{1,2} \) in Fig. 9a and b, \( J_{\text{HT}}^{1} ,J_{\text{HT}}^{2} = 8kT \) and \( J_{\text{SW}}^{1} = 8kT \), \( J_{\text{SW}}^{2} = - 16kT \) in all. For β = 0, the spins remain in an all-parallel state due to the short range exchange and dipole-dipole interaction, which turns out to be dominant when the 3(Si · n)(Sj · n) term in Eqn. (2) attains values for angles equal or larger than π/6 in the cross product (see Fig. 1a).

Fig. 9
figure 9

The bilayer configurations after 1000 MCS for: a \( \beta = 0,\,A = J_{\text{HT}}^{1,2} = 8kT \), b \( \beta = 100,\,A = J_{\text{HT}}^{1,2} = 8kT \), and c \( \beta = 100,\,A = 0 \) (No d–d interaction)

Figure 9b and c reveal the configurations for \( \beta = 10J_{\text{FE}} \bar{S} \) that is sufficient to force the FE component of the grids to evolve into spin-up and spin-down domains. In the presence of \( A = J_{\text{HT}}^{1,2} \) in (2), the FE layer transforms into a periodic one having fine laths with flat interfaces. The flat domain walls are a consequence of the competition where the sign of the sideways exchange and the part of the dipole-dipole interaction are effective with respect to FE ordering, coming mainly from the 3(Si · n)(Sj · n) term that is zero when \( n \bot \left\langle {S_{\text{i,j}} } \right\rangle \), leaving (2) effective with the \( S_{\text{i}} \cdot S_{\text{j}} \) product. The domain size depends on the strength of the dipolar interaction term where the sideways anti-parallel alignment and head-to-tail parallel alignment compete. In the absence of any long-range contribution in Fig. 9c (A = 0), the cost of the 180° domain formation is only a slightly perturbed interface, hence it is not surprising to observe that the FE layer splits into two equal halves of up-spins and down-spins to minimize the term in (4) after 1000 MCS. We would like to remind here that a needle-like FE crystal does not suffer from a depolarization field due to the small area where the polarization vector pointing along the needle terminates. Such a geometric effect is not accounted for in this study, i.e., we externally introduce this term that the systems try to minimize.

The superlattice splits into domains for the same values of the coefficients used in the bilayer runs but at much smaller β as provided in Fig. 10. As much as there is the imposition for parallel alignment of the FE layers on the AFE ones, there is the AFE influence on the FE layers and vice versa near the interfaces. Therefore, the interfaces act as nucleation centers for the domains in the FE even at relatively small values of β. This could also be a qualitative explanation for the slimmer hysteresis loops often observed for multilayers in experiments where the structure switches much easier than bulk at lower coercive fields. The influence of layer periodicity has also been investigated for FE–PE multilayer structures by Stephanovic et al. [41] using an analytical method to demonstrate the influence of FE layers on each other through the electrostatic effects without any strain-related parameters. They conclude that above a critical layer frequency, the FE layers start interacting with one another and the entire structure has a minimally varying polarization profile within each domain.

Fig. 10
figure 10

Superlattice configuration for \( J_{\text{HT}}^{1,2} = J_{\text{SW}}^{1} = 8kT,\,J_{\text{SW}}^{2} = - 2J_{\text{HT}}^{2} ,\,A = J_{\text{HT}}^{1,2} \,{\text{and }}\,\beta = 5J_{\text{HT}}^{1} \) after 1000 MCS

Conclusions

We carried out Monte–Carlo simulations where the phase transition behavior, hysteresis characteristics, and equilibrium configurations of FE, AFE, and multilayers consisting of both were investigated. The long-range interactions have a substantial influence on the phase transition behavior and configurational order of the system employed in this study. For the case of superlattices with several interfaces, quite different hysteresis behaviors were displayed when a stable FE and a stable AFE were thought to form the structure. The AFE characteristics tend to totally disappear in hysteresis loops of superlattices while the bilayer can still exhibit independently switching AFE clusters characterized by extensions of the hysteresis loops. Nevertheless, the multilayer and the bilayer loops are much slimmer, implying that lower applied fields suffice to switch the systems compared to the loops obtained from the single component FE and AFE grids, in good qualitative agreement with real experimental observations. The electrostatic-like coupling does not induce an AFE-to-FE transition in the AFE layer as the FE layer splits into clusters of spins with flat interfaces similar to 180° domains extensively treated in numerous studies particularly using the LGD formalism.

In short, using hypothetical order–disorder systems, we analyzed the influence of long-range dipolar interactions for various behavior of multilayers that often yield FE-like loops for FE–AFE structures with high component periodicity often observed in experiments. There can certainly be other influences such as the coexistence of FE and AFE phases in a strained AFE layer and these formations have to be examined under the knowledge of individual strain states and relaxation behavior of the layers.