Influence of domain walls and defects on the electrocaloric effect

The electrocaloric (EC) effect is the adiabatic temperature change of a material in a varying external electric field, which is promising for novel cooling devices. While the fundamental understanding of the caloric response of defect-free materials is well developed, there are important gaps in the knowledge about the reversibility and time-stability of the response. In particular, it is not settled how the time-dependent elements of microstructure that are always present in real materials act on the field-induced temperature changes. Ab initio based molecular dynamics simulations allow us to isolate and understand the effects arising from domain walls (DWs) and defect dipoles and to study their interplay. We show that DWs in cycling fields do not improve the response in either the ferroelectric (FE) phase or at the FE phase transition, but may result in irreversible heat losses. The presence of defect dipoles may be beneficial for the EC response for proper field protocols, and interestingly this benefit is not too sensitive to the defect configuration.


Introduction
Cooling technologies based on the electrocaloric (EC) effect-the temperature change in a material when an external electric field is applied under adiabatic condition-are heralded as potential energy efficient solid-state cooling devices [1]. This effect has been known since 1930s [2], however only since the 2006 report of a giant EC effect in thin ferroelectrics (FE) films [3] has the response been considered large enough for applications. This is reflected in the number of publications increasing from less than 100 in about 80 years to more than 1300 since 2006 3 . Despite this increased activity and development of EC prototypes [4], achieving real-world applications still requires scaling to industrial level and the fundamental understanding and optimization of EC materials themselves. Material properties that act as bottlenecks include too low efficiency (i.e., low response for small work input), narrow and often not suitable temperature windows with maximal responses, irreversibilities, aging and functional fatigue (i.e., the change and degradation of the response during subsequent cooling cycles).
The narrow temperature window with maximal EC response corresponding to the field induced paraelectric to FE phase transitions has been extensively studied in the literature [3,[5][6][7]. Furthermore, enhanced responses at lower temperatures have been reported at FE-FE transitions [8,9]. However, at these transitions a small operating range and irreversibilities [10,11] are potential obstacles. Furthermore, the phase boundary stabilizes multi-domain phases [12]. The response in a given FE phase is, typically, one order of magnitude smaller. It has been reported that fields applied in the opposite direction to the remnant polarization may enhance this response, however large heat losses are induced during field-induced switching [13].
It is established from the tuning of the piezoelectric response that walls may enhance functional responses via 'frustrated poling' [14]. Furthermore, based on a simple Landau picture, the caloric response of a material is typically linked to the field-induced macroscopic change of the polarization P as ∆S ∼ −∆P 2 , where ∆S denotes the field-induced isothermal entropy change. Within a given FE phase, the field induced change of P is typically larger in a multi-domain state. Based on this, one might expect that domain walls (DWs) may, in general, enhance the caloric response. However, as we show in this paper, this equation does not apply to inhomogeneous multidomain materials in which local and global orders differ. For a single domain material, it has been further shown that even the sign of the caloric response can be tuned by the relative direction of the field and polarization [15]. This raises the question of what other effects can be found in multidomain phase FE with multiple relative directions between field and polarization.
So far, additional extrinsic contributions to the temperature change by DWs have been reported at the transition between different multi-and single-domain states in strained films [16][17][18][19], in the presence of depolarization [20] or originating from complex polar vortex ordering within superlattices [21]. Furthermore, improved caloric responses coming from DWs themselves have been predicted, including positive and negative EC temperature changes at each side of the wall [22], enhancement of the adiabatic temperature change ∆T with the density and width of immobile elastic walls [23] and the enhancement of the EC on Bloch-like walls [24]. Many of these studies are based on thermodynamic models and do not answer the question of whether these local enhancements may be used in practical applications. In particular, wall motion has often been neglected. This, however, is particularly relevant, because it is known that field-induced switching by DW motion is related to irreversible heat losses [20,25,26].
The second element of microstructure that is often neglected in studies of the caloric response are internal bias fields, generated by the different dynamics of defect complexes. These bias fields are particularly relevant for deficient or hard doped materials where acceptor ions and oxygen vacancies form defect complexes such as defect dipoles [27]. These may act as restoring forces on the domain structure and thus reduce wall motion and polarization response [28] or mediate large reversible piezoelectric responses [29] and their arrangement and alignment within the surrounding domain structure play an important role in aging and fatigue [30]. Studies on the impact of defect dipoles on the EC effect reported reduced ∆P/∆E and reduced work losses [31] or improved reversibility and positive and negative responses, i.e., cooling under field removal or application [32][33][34]. So far, a fundamental understanding of the coupling between DWs, defects and caloric response is missing. This is because in experiments it is difficult to separate these elements of microstructure and, at the computational level, such a study requires large system sizes and multi-scale modeling.
In this paper, we study the adiabatic field response of multi-domain FEs in combination with symmetry conforming defect dipoles and the phase boundary between two FE phases. We use BaTiO 3 as a model system, as its phase and domain structure [12,35] and DW dynamics [36][37][38] are well established and as its experimental FE-FE phase transition from tetragonal (T) to orthorhombic (O) phase at 5 • C occurs at a suitable temperature range for applications. To separate the different effects, we first include only DWs and calculate the field response in the T and O phases. We then include defect dipoles along with the DWs and study the combined effect, and finally, we study the multi-domain system with defects near the phase transition.

Computational details
We perform molecular dynamics (MD) simulations within the feram code [39,40] utilizing the ab initio based effective Hamiltonian [41,42] which is specifically applicable to FE materials with a cubic perovskite structure at high temperatures. It uses coarse-graining to reduce the number of degrees of freedom as illustrated in figure 1(a). This method has been successfully used to determine phase transitions, electric field-response, DW properties and caloric response of ABO 3 pervoksites [9,15,[43][44][45][46][47]. We model BaTiO 3 using the parameters reported in [47]. We note that despite the high predictive power for qualitative trends, absolute transition temperatures are underestimated for the bulk material as tabulated in table 1, while the field strengths to activate DW motion or induce a phase transition are overestimated.
After thermalization at the temperature of interest, we directly monitor the change in temperature of the material by ramping on and off an electric field with the rate ∆E/∆t = 0.001 kV cm −1 fs −1 within the microcanonical ensemble [48]. Note that the system in the canonical ensemble shows temperature fluctuations with time and thus the initial temperature of the adiabatic step may deviate from the targeted temperature. The change in temperature may consist of two contributions, one from the EC effect and a second one from irreversible heat losses. In order to account for the reduced number of degrees of freedom in our model, we scale the temperature change by a factor of 1/5 [48].
We monitor the field-induced change of the polarization by the average dipole moments of N unit cells: where i is the Cartesian direction denoted by a, b or c. If the sum is taken over all of the unit cells, this corresponds to the global polarization P i of the system. For the multi-domain state, this quantity may be equal to zero and is not sufficient to characterize the state. The dipole moment per layer ⟨p i ⟩ x where the sum is restricted to layer x, resolves the change of the local polarization across the DW and its sign may The energy is given relative to the single domain (SD) material at the same temperature. Favorable wall positions at x0, x0 + 1 etc are higher in energy than the SD phase by the domain wall energy ∆EDW and are separated by the activation barrier ∆En. After each shift, this energy is dissipated as heat. If two walls annihilate the energy ∆En + ∆EDW is released. In the presence of defect dipoles, the energy minima shift by ∆E d for each plane with antiparallel defect dipoles crossed by the wall (purple line). In both cases, the mean polarization per layer along the y-direction ⟨p b ⟩ x changes sign across the walls. We initialize these walls by the pre-poling protocol from reference [38] i.e. thermalization of the system at a given temperature using the Nosé-Poincaré thermostat [49] under applied local fields, which are then successively removed. We use these converged configurations as input for subsequent cooling, heating or field-ramping simulations.
For the DW distance d considered here, the walls do not interact and are thermally stable [50]. Figure 1(e) sketches the energy landscape for DW motion: shifting of the walls is possible if the system overcomes the energy barrier ∆E n separating neighbouring favorable wall positions. Thus, field induced DW motion is possible above a temperature dependent activation field. We note that our aim is not to correctly pinpoint this activation field, and we estimate it by the change of ⟨p i ⟩ x with external field strength with a resolution of 20 kV cm −1 , only.
Defect dipoles with tetragonal symmetry can be microscopically realized by a pair of a doped acceptor ion and an oxygen vacancy, whereas the orthorhombic symmetry can be realized by the equal distribution of such tetragonal defects along both directions of the orthorhombic polarization or more complex defect structures. In order to study the impact of defect dipoles in an aged sample on a qualitative level, we follow the approach from [32] and include 1% of randomly distributed defect dipoles with a magnitude 0.1 eÅ fully aligned with the initial domain structure (SD, T180 or O90) which are frozen-in during the entire simulation cycle. We note that dipoles may change their orientation and symmetry with time, particularly in large fields or at phase transitions. However, at and below room temperature this re-alignment is slow, due to the large oxygen migration barrier [27] and can be ignored for a qualitative discussion.

Results and discussion
To answer the question how DWs modify the field induced temperature change within one FE phase, we first focus on the defect-free material, well apart from any transition in tetragonal and orthorhombic phases, and compare multi-and single-domain materials. We start with the single domain tetragonal phase at T i = 260 K and cycle a field from 0 to 80 kV cm −1 , to 0 and finally to −80 kV cm −1 , i.e., first parallel and then antiparallel to the polarization. The initial state has a remnant polarization of ⟨p b ⟩ = 28.8 µC cm −2 and in the interval ±40 kV cm −1 field ramping induces a reversible linear change of polarization of ±0.02 µC cm −1 kV −1 (see the red lines in figure 2(a)). For larger antiparallel fields (E < −40 kV cm −1 ), the system leaves the linear regime and the field-induced disorder increases. Here, the increase in the standard deviation of ⟨p b ⟩ by 6% during field ramping from 0 to −80 kV cm −1 serves as a measure of the decreasing order of the dipoles. With the change of ⟨p b ⟩ x = ⟨p b ⟩ the temperature changes reversibly with the field, reaching maximal heating and cooling of +0.5 K and −0.6 K at +80 kV cm −1 and −80 kV cm −1 , respectively, in good qualitative agreement with previous experimental and theoretical results in [13].
In the presence of T180 walls, the same field-protocol results in a different time evolution of polarization (compare solid lines in figure 2(a)). Increasing the field from zero induces an increase of ⟨p b ⟩ which is one order of magnitude larger compared to the SD phase, the polarization keeps increasing even during field removal and reaches a plateau around 130 ps, i.e., at 30 kV cm −1 . Inspection of the polarization profile ⟨p b ⟩ x (x) (not shown here) reveals two distinct field regimes: (I) Below a magnitude of the applied field of about 30 kV cm −1 in the first 30 ps and between 130 and 180 ps, the walls are pinned by the lattice. (II) If the system is exposed to fields with larger magnitude, the domains parallel to the field grow by DW motion.
In both regimes the changes of the macroscopic polarization ⟨p b ⟩ cannot be directly correlated to ∆T. For example, the presence of domains reduces ∆T during field application and the temperature starts to decrease in the plateau region (see panel (b)). Parts of these trends can be understood by the time evolution of the magnitude of the local polarization well in the domains |⟨p⟩ x | added as lines with symbols in panel (a); initially |⟨p⟩ x | in both domains is approximately equal to the SD case and the presence of domains reduces the slope of |⟨p⟩ x |(E), particularly in the −b domain, and in turn ∆T(∆E) is smaller in the presence of domains. The same holds between 130 ps and 180 ps, i.e., during ramping from 30 to −20 kV cm −1 ; the local polarization and the temperature decrease during field ramping but domains reduce the slopes compared to the SD material. Only within 1-2 f.u. directly at the walls ∆⟨|p|⟩ x is slightly enhanced. Zhang et al [24] correlated the same observation to an enhanced caloric response. However, the opposite signs of the response on either side of the wall in combination with the small volume of this region do not allow for an overall additional cooling. Now let us look at the temperature change in field regime (II): the maximal heating of about 0.6 K is reached after about 110 ps during field removal and after a full field-cycle at 160 ps the material is heated up by 0.4 K. This additional heating is related to irreversible heat losses due to DW motion rather than to the EC effect (see figure 1(e)). The field constantly does work to move the wall across the local energy barrier (∆E n ) which is dissipated as heat. This interpretation is also in agreement with the heating found for large negative fields; although ∆P and ∆|⟨p⟩ x | decrease, the dissipated energy does not depend on the direction of the wall shift but is always positive.
The heat dissipation depends not only on the magnitude of the field, but also on the duration that the system is exposed to a field larger than the activation field as the walls continuously shift through the system. For example, figures 2(c) and (d) show the field response if the field of 80 kV cm −1 is first ramped on and then kept constant. For a fixed strength of the field, the wall moves through the system with constant velocity and heat dissipation and the temperature increases linearly. The worst case scenario is visible after 140 ps: two neighboring walls meet and annihilate each other. Thereby, the energy ∆E n + ∆E DW is abruptly released as heat, resulting in a temperature jump of +0.3 K (cf. figure 1(e)). Both the constant heating under wall motion and the temperature increase under annihilation have been reported previously for elastically driven domains [26]. Note that this temperature change solely arises from irreversible heat losses and cannot be considered as a caloric response.
Can these findings be generalized to other types of walls? To answer this question, we picked the example of O90 walls with ab/a−b polarization (see figure 3). If the field is applied along b (blue lines) the response is analogous to T180 walls: in field interval (I) the walls are pinned and local changes of polarization as well as temperature are tiny, in interval (II) the field stabilizes the ab over the a−b domain and the macroscopic polarization increases due to DW motion, which is related to irreversible heat losses. Compared to the T180 case, the activation field for the wall motion is enhanced to about 120 kV cm −1 due to the finite angle between the polarization and the applied field and the reduced temperature, as the wall motion is a thermally activated process [50]. The reversible temperature change is smaller than the resolution of our data. Partly this is related to the already well saturated polarization at these low temperatures. Furthermore, it has been discussed that the conventional and inverse EC responses compete in the orthorhombic phase for fields applied along ⟨100⟩ [15].
Note that larger changes of |⟨p b ⟩ x | are again induced at the DW itself which broadens within the applied field. However, even at 200 kV cm −1 , this effect is restricted to 4 f.u. and does not improve the overall cooling. We conclude that, for the studied field direction, the walls in both phases reduce the EC response and not only the full switching of the polarization direction but already partial switching by wall motion induces irreversible work losses.
Can DWs improve the response for other field directions? To test this hypothesis, the black lines in figures 3(a) and (b) show the response of the O90 system for the same field protocol as above, but with the field applied along a. For this field direction, there is no driving force for wall motion as both domains incline the same angle with the field. Indeed, the walls do not move and the polarization component parallel to the field changes reversibly with the field strength. However, the walls are still not beneficial for caloric cooling, as the combination of small temperature and competing conventional and inverse EC effects hinder any caloric response up to about 100 kV cm −1 , while a reversible inverse EC of only 0.2 K is induced by 200 kV cm −1 . This response is related to a field-induced monoclinic distortion in both domains close to the FE phase transition. Indeed, by increasing the field strength further (see panels (c) and (d)), the system undergoes a transition to the tetragonal phase at a critical field strength of 240 kV cm −1 for the given temperature. However, the DWs disappear during the transition and further cycling of the field results in a sequence of phase transitions from SD T to SD O (445 ps, 145 KV cm −1 ) and finally from SD O to SD T (890 ps, 290 kV cm −1 ). A comparison of the response of the O90 system in the first field-cycle (between 200 and 300 ps) with the response of the O SD material (800-900 ps) shows that the DWs have a minor impact on the maximal cooling but reduce the critical field strength for the transition. Thus, for this field direction also the response in one phase is not improved by walls.
Could DWs improve the caloric response if they are pinned by defect dipoles? How do defect dipoles in general modify the caloric response and how do these effects depend on the change of the defect alignment with time? To answer these questions we compare the impact of defect dipoles aligned with the SD tetragonal phase and equidistant T180 domain structures. The red lines in figure 4 show the field response for SD defects: for 100 kV cm −1 applied parallel to the defects, the response is reversible and not considerably different than for the SD material. We note that further tuning of the field protocol to optimize the caloric response in the presence of such uniform defect dipoles was possible as discussed in [34] but is not our present focus.
Second, we monitor the field response of the system with defect dipoles aligned with the T180 domain structure. The blue and green lines in figure 4 compare the responses with and without defects. Without defects, we observe large heating even during field removal, as the walls move 24 f.u. between 30 and 150 ps and annihilate with the corresponding jump in temperature discussed before. Subsequent cooling comes from the normal EC response within the SD phase. In contrast, defects reduce the irreversible heating, and within the first 150 ps of the ramping cycle, the observed change in T is below the computational resolution. A closer inspection shows that in field regime (I) the tiny field response is not modified by defects nor is the activation field. Above the activation field in regime (II) the defects strongly modify the response. In the presence of defects, the walls move less than 15 f.u. (∆P < 10 µC cm −2 ), start to move back after 120 ps and the initial domain structure is fully restored after one field cycle and heating of about 0.1 K is observed during the ramping off from 50 kV cm −1 , only. The defect dipoles thus indeed reduce the irreversible heat losses, as the annihilation of walls is prevented and the energy barrier for wall motion is crossed less often. Note that the system can furthermore absorb heat in the first 120 ps as the energy increases by E d if a layer is polarized antiparallel to the defects (see figure 1(e)). This energy is, however, released during the back relaxation, resulting in overall heating. The reduction of heat losses by defects with T180 symmetry cannot be utilized for EC cooling as, at the same time, the reversible caloric response is strongly reduced. In the presence of defects, the field induced ∆|⟨p⟩ x | has opposite signs in both domains. In particular, a larger negative ∆|⟨p⟩ x | is induced in the layers where the polarization switches against the defect dipole. For the chosen field strength and equidistant domains the responses thus compensate each other.
We also simulated a range of defect concentrations and strengths (0%-2% and 0.05-0.2 eÅ), the field-induced shift of the wall and the irreversible heating scale with both the density and the strength. However, none of these configurations allows for a large reversible EC response. For all defect configurations, larger fields are able to overcome the restoring forces of the defects. Figures 4(c) and (d) show an example where, during application of 200 kV cm −1 to the system with 1% defects of 0.1 eÅ strength, the T180 walls annihilate and do not reappear during field removal, i.e., the whole system is finally poled parallel to the external field with one half of the volume with antiparallel defect dipoles. During the removal of 200 kV cm −1 , the system stays in the SD state and we find ∆|⟨p b ⟩ x | of −4 µC cm −2 and −5 µC cm −2 in the regions with parallel and antiparallel defects, resulting in a drop of temperature of about 1.3 K, i.e., comparable to the 1.2 K found above for the defect-free SD material (ramping from 80 to −80 kV cm −1 , cf figure 2). The defect dipoles aged within a multi-domain configuration are thus no obstacle for the EC response after full poling of the free dipoles.
After showing that defect dipoles can indeed act as a restoring force for the domain structure and that domain walls reduce the thermal hysteresis of the phase transition but may vanish at the transition [45], the next natural step is to study the effect of a combination of these two together near the phase transition. In figure 5, we plot the average mean polarization ⟨|p i |⟩ as a function of temperature for heating and cooling runs without an external field. Both types of defects promote the phase collinear to their alignment, i.e., in the presence of ±b defects the tetragonal phase is present at lower temperatures than in the SD case and vice versa for a±b defects (see transition temperatures in table 1). In the presence of b/−b defects, we furthermore find a reversible transition between T180 and O90 phases and a reduction of the thermal Figure 5. Change of the polarization components for zero-field cooling (blue) and heating (red) of the system with 1% defect dipoles fully aligned with the (a) T180 and (b) O90 domain structure. Note that the global polarization ⟨p b ⟩ is zero for all temperatures, and we instead plot the system average of the magnitudes of the components ⟨|p i |⟩. Figure 6. Response of the system with 1% defects aligned with O90 domains (a ± b): polarization (top panels) and adiabatic temperature change (bottom) as the field is varied along b (left) and a (right) directions are plotted as a function of simulation time at initial temperatures of 170 K (blue), 190 K (black) and 210 K (red). The field is increased from 0 to 100 kV cm −1 and then decreased to zero. Note that the system is close to its phase transition temperature of 195 K, i.e., at 170 K and 190 K the initial state is O90, while at 210 K it is in the monoclinic phase with |⟨px⟩| > |⟨py⟩|. hysteresis by 60%. Excitingly, defects aligned with ab/a−b domains completely close the thermal hysteresis of the phase transition within our temperature grid of 5 K, and we find a reversible transition with a fixed position of DWs from O90 to a monoclinic phase with ⟨|p a |⟩ ≫ ⟨|p b |⟩ ̸ = 0. Note that also a distribution of a and ±b defects with an overall orthorhombic symmetry which could be microscopically realized by the aging of acceptor-vacancy pairs in this phase, results in the same qualitative trends.
To answer whether the combination of O90 DWs and defects results in a large and reversible caloric response at the phase transition, we study the field response at temperatures close to the transition: for T < T c (170 K), O90 for T ⩽ T c (190 K) and for the monoclinic phase at T ⩾ T c (210 K) as shown in figure 6. In the latter case, ramping on the field along b induces a continuous polarization rotation from a towards b without any changes in the domain structure. For the configuration shown in figure 6 this rotation is fully reversible, reaches a monoclinic state ⟨|p a |⟩ ⩾ ⟨|p b |⟩ at 100 kV cm −1 and is related to an approximate reversible caloric response of 0.8 K, slightly larger than the SD response without defects. Note that increasing the field strength further (not shown) results in an increasing reversible ∆T until the transient state with ⟨|p a |⟩ = ⟨|p b |⟩ is reached. Even larger fields reduce ∆T again as the polarization rotates toward ⟨|p b |⟩ ≫ ⟨|p a |⟩. This is expected, due to the different entropies of orthorhombic and tetragonal phases [15]. Applying the field along the a-direction results in a smaller linear response. Again, the same qualitative trends are found for a combination of a and ±b defects. For 190 K, the system is initially in the state with |⟨p a ⟩ x | > |⟨p b ⟩ x | ̸ = 0, i.e. O90 phase with small monoclinic distortion. Applying the field along b disorders the dipoles in both domains with a larger decrease of ⟨|p b |⟩ x in the −b domain. Analogous to the continuous O to T transition discussed above (see figure 3(d)), the system shows an inverse response reaching about −0.4 K. But, no full transition to the tetragonal phase with polarization along a is possible, as the field is applied along b. Around 60 kV cm −1 the DW start to broaden and move slightly. Finally, the system reaches the SD monoclinic phase related to an increase in temperature by 0.8 K. This complex field response is approximately reversible under field removal. Changing the field direction to a results in an inverse EC effect with a maximal temperature change of about −0.4 K.
Reducing the temperature by 20 K to 170 K (blue lines in figure 6), the temperature change is similar to that observed for the defect free material at 120 K shown in figure 3(b) for both field directions. For the field along b, no significant change in temperature is observed during ramping until a small increase around 160 ps resulting from the DW motion. When the field is applied along the a-direction, an inverse temperature change associated with polarization rotation is observed. Again the same trends are found for a distribution of tetragonal defect dipoles with overall orthorhombic symmetry. In summary, the restoring forces of defect dipoles allow for a reversible EC response by polarization rotation at the phase boundary between the orthorhombic and tetragonal phases.

Conclusions
In this paper, we performed a systematic study on the caloric response of BaTiO 3 at and below room temperature phases and at the tetragonal-orthorhombic transition. Our MD simulations for an ab initio based Hamiltonian enable us to separately study the impact of DWs, defect dipoles and their interplay.
We show that the inclusion of DWs in the bulk material is overall not beneficial to the caloric response, because switching by wall motion results in irreversible heating. Even if the DWs are pinned by the lattice for small field strengths or additional defect dipoles, the presence of different domains reduces the overall EC response as the change of the local polarization is reduced. Also, the slightly larger change of the dipole ordering near the walls discussed in previous work does not improve the overall cooling, due to the small volume fraction and opposing signs of this response on either side of the wall.
On the other hand, defect dipoles may be beneficial to enhance the caloric response within one FE phase. In particular, defect dipoles aligned with the orthorhombic phase enhance the response at higher temperatures and defect dipoles aligned within the tetragonal phase are at least no obstacle for the response independent of whether the dipoles are poled in the SD material or have been aligned in a multi-domain alignment. The combination of DWs and defect dipoles allows us to reduce thermal hysteresis, particularly for defect dipoles aligned with an O90 domain structure. However, for moderate field strengths, it is not possible to convert the full latent heat of the transition into an adiabatic temperature change due to the complex interplay of external and internal electric fields.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.