{'fraction_non_alphanumeric': 0.03420273897185328, 'fraction_numerical': 0.014038951207762714, 'mean_word_length': 4.643495145631068, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Direct molecular dynamics simulation of electrocaloric effect in BaTiO3

The electrocaloric effect (ECE) in BaTiO3 is simulated using two different first-principles based effective Hamiltonian molecular dynamics methods. The calculations are performed for a wide range of temperatures (30--900 K) and external electric fields (0--500 kV/cm). As expected, a large adiabatic temperature change, Delta-T, at the Curie temperature, T_C, is observed. It is found that for single crystals of pure BaTiO3, the temperature range where a large Delta-T is observed is narrow for small external electric fields (<50 kV/cm). Large fields (>100 kV/cm) may be required to broaden the effective temperature range. The effect of crystal anisotropy on the ECE Delta-T is also investigated. It is found that applying an external electric field along the [001] direction has a larger ECE than those along the [110] and [111] directions.


I. INTRODUCTION
The electrocaloric effect (ECE) is an adiabatic change in the temperature, ∆T , of a material upon applying an external electric field. In particular, if an electric field is applied to a ferroelectric material at just above its phase transition temperature, T C , and the field is then removed, a large reduction in temperature is expected. It is widely believed that this effect is applicable to solidstate refrigeration technologies.
In addition, recent developments in the techniques of vapor deposition enable the production of defect-free single-crystal ferroelectric thin films. Such high-quality films allow for the application of large external electric fields, which cannot be applied to bulk polycrystalline specimens. Consequently, this advance in processing allows for investigation of the ECE in these ferroelectric thin films. There has been much interest in this subject and experimental studies have shown the possibility of creating materials with a relatively large electrocaloric response [1][2][3] .
There have been several computational simulations of the ECE published in the literature. Ponomareva and Lisenkov have investigated the ECE of Ba 0.5 Sr 0.5 TiO 3 using Monte Carlo methods 4 . Rose and Cohen 5 have used molecular dynamics (MD) simulations and coreshell interatomic potentials to model the ECE in bulk LiNbO 3 . Also using this form of atomic potential, Chen and Fang 6 have simulated the ECE in BaTiO 3 nanoparticles. We have also used the so-called indirect MD method, discussed below, to calculate the ECE of bulk BaTiO 3 using a first-principles based effective Hamiltonian 7 . All of these simulations find a large ECE is observed just above T C due to the large change in entropy when transforming from the paraelectric to ferroelectric phase.
Following this introduction we present the methods used to calculate the ECE and in particular we introduce two MD methods that can be used to study the electrothermal coupling. The results are presented in Sec. III and in Sec. IV, the paper is summarized and conclusions are given. The methods and results presented here are a full detailed review that extends and explains the preliminary results given in the proceedings Ref. 8.

A. Effective Hamiltonian
The effective Hamiltonian, constructed from firstprinciples calculations, and used in the present MD simulations is essentially the same as that in Ref. 9 and 10 The true atomic structure has properties determined by the complex chemical bonding between the atoms, but in the model system the complexity is reduced; the collective atomic motion is coarse-grained by local soft mode vectors, u(R), and local acoustic displacement vectors, w(R), of each unit cell located at R in a simulation supercell as depicted in Fig. 1. Each term in the Hamiltonian bears a physical significance: η 4 = e 23 ), V elas, inho ({w}) is the elastic energy from inhomogeneous strains, V coup, homo ({u}, η 1 , . . . , η 6 ) is the coupling between the local soft modes and the homogeneous strain, V coup, inho ({u}, {w}) is the coupling between the soft modes and the inhomogeneous strains, and −Z * R E ·u(R) is the interaction between dipoles Z * u(R) and external electric field E, here Z * is the effective charge of the soft mode per unit cell. The coarse graining and calculation of V dpl ({u}) using fast-Fourier transform (FFT) enabled us to perform fast MD simulations in a large supercell, and was previously applied to bulk relaxor ferroelectrics 11,12 . A more detailed explanation of the effective Hamiltonian can be found in Refs. 13, 14, 10, and 15. The set of parameters for the effective Hamiltonian for BaTiO 3 are listed in Ref. 15. A temperature dependent effective negative pressure given by p = −0.005T GPa is applied to simulate the thermal expansion 15 .

B. MD methods and conditions
MD simulations for BaTiO 3 using the effective Hamiltonian in Eq. (1) are performed with the feram software.
feram is distributed freely under the GNU General Public License (GPL) and can be found at (http://loto.sourceforge.net/feram/). Details of the code can be found in Ref. 10. For the calculations presented here version feram-0.21.02 was used.
Examples of the input files are packaged within the source code under the feram-X.YY.ZZ/src/26example-BaTiO3-acoustic-MD/ directory.
As described in Ref. 7, for the so-called indirect approach, the pyroelectric coefficient is predicted as a function of applied electric field and temperature. Then, using the thermodynamic relationship the ECE ∆T is related to the temperature and external electric field. Here, P is polarization, and the external electric field, E, is switched from E 1 to E 2 . The specific heat capacity, C, is approximated as a constant and its value is taken from experimental observation 16 . However, it is well known that C is dependent on the external electric field and temperature, especially near the phase transition temperature T C . Moreover, under a weak external electric field, it is numerically difficult to determine ∂P /∂T at T C . Therefore, a better approach involves the so-called direct method to estimate the ECE ∆T . The procedure involves two steps: first constanttemperature MD is performed for a fixed external electric field, E, in the canonical ensemble using the velocityscaling thermostat. This allows the system to equilibriate. Next, the external electric field is switched off and the system is simultaneously changed to a constantenergy MD in the microcanonical ensemble that is allowed to evolve using the leapfrog method. The final state at the end of the constant-temperature MD is used as the initial state of the constant energy MD. A time step of ∆t = 2 fs is used in both ensembles.
There are two ways to determine the acoustic displacements, {w(R)}, at each time step of the MD simulation. One allows for the natural time evolution of the effective Hamiltonian Eq. (1) with MD. As discussed in Sec. II A, the number of degrees of freedom is reduced from 15 (5 atoms × 3 xyz-directions) to 6 (two 3-dimensional vectors) and consequently the specific heat capacity implicitly becomes 6/15 = 2/5 of the real system. Another approach is to optimize {w(R)} such that V elas, inho ({w})+V coup, inho ({u}, {w}) is minimized at each time step according to {u(R)}. In this case, {w(R)} is fully dependent on {u(R)}, and the degrees of freedom is further reduced from 6 to 3, meaning that the specific heat capacity of the calculated system is 1/5 of the real system. In both of these direct MD methods, {u(R)} evolves normally according to the effective Hamiltonian of Eq. (1).
The ECE response is calculated using the indirect method, direct MD method, and direct optimized method and are compared in Fig. 2. The values of ∆T , corrected for the underestimated heat capacity in the direct methods, are also given in the figure. The differing results for these three methods come primarily from two effects. One is that in the direct methods, the effective negative pressure, p = −0.005T GPa, is kept constant at the starting temperature, even though the constant energy evolution of the system causes a reduction in temperature. Another reason for the difference between the results is that the temperature and external electric field dependence of the specific heat capacity is automatically included in the direct methods whereas it must be approximated from experiment in the indirect methods.
The direct MD approach gives half of the raw temperature change ∆T during the constant-energy MD simulation as compared to the direct optimized approach due to the difference in the degrees of freedom. However, error is introduced because the T -dependent effective negative pressure is held constant during the constant energy temperature change. Therefore, in principle, the direct MD approach may give more accurate results than the direct optimized approach. In the study presented here, however, the direct optimized method is used because it is shown in Fig. 2 to produce nearly equivalent results and is almost twice as efficient computationally.
In Fig. 3, the ∆T is compared for two supercells of sizes L x × L y × L z = 48 × 48 × 48 and 96 × 96 × 96 unit cells. The greater fluctuations in the ∆T in the 48 × 48 × 48 system is likely due to the system-wide thermal fluctuations that evolve in the constant-energy MD. Therefore, a system size of 96 × 96 × 96 is employed for the simulations presented here. To accurately capture the temperature dependent phase transformations and account for the thermal fluctuations discussed above, a temperature step size of 0.5 K is used and numerous simulations are performed, e.g. 1,400 sets of MD simulations are performed for a sweep from 300.0 to 999.5 K.

III. RESULTS AND DISCUSSION
Using this method the impact of crystallographic anisotropy on the ECE is investigated by examining ∆T as a function of temperature for electric fields applied along the [001], [110], and [111] directions. As seen in Fig. 4 there exist singularities at the temperatures associated with the tetragonal-to-orthorhombic and orthorhombic-to-rhombohedral transformations; however, |∆T | is maximum at the cubic-to-tetragonal transition temperature T C for all applied fields. It is clearly seen that when the field is applied in the [001] direction ∆T is the largest and the available temperature range for the large ECE is the greatest. This is because of the strong coupling between the external electric field and the internal dipole moment. The first ferroelectric phase at T < T C is tetragonal and therefore the applied electric field acts to broaden the transformation, effectively increasing the transformation temperature. This is of practical importance for solid state cooling technologies because the engineering devices must be tuned to operate within a specific range of temperatures. For these applications ECE research should be focused on the cases where the applied electric field aligns with the intrinsic polarization at the Curie temperature, which in the case of BaTiO 3 is the [001] direction.
In Fig. 5, the temperature dependence for the ECE ∆T of BaTiO 3 , under various initial external electric fields is compared. It can be seen in Fig. 5(a), that even with a small initial external electric field (< 50 kV/cm), BaTiO 3 gives a large ∆T , but the temperature range where this large ∆T can be obtained is narrow. This result is consistent with previous theoretical studies of  other ferroelectric materials 4-6 in which |∆T | has a peak at T C . By increasing the applied fields (> 100 kV/cm) the range of applicable temperatures broadens as shown in Fig. 5(b). This is again due to the strong coupling between the tetragonal polarization and the external electric field, which both broadens the thermal range of the transformation and increases the pyroelectric response. For engineering applications it may be necessary to apply as large of field as possible to allow for an operating temperature range that is useful.

IV. SUMMARY AND CONCLUSIONS
In this work we have described two MD direct methods to evaluate the ECE and they are compared to the previous indirect approach. The reduction in the degrees of freedom associated with this effective Hamiltonian method are explained and its impact on the specific heat capacity, as applied to the direct MD and direct optimized methods, are observed. The three methods demonstrated here are found to produce roughly equivalent results; however, based on the computational errors discussed above it is concluded that the direct MD method is in principle more accurate although it is slower than the direct optimized method used in these calculations.
The effect of crystal anisotropy is found to be great. For BaTiO 3 , applying a field in the [001] direction results in the greatest electrocaloric response. This is good news for ferroelectric thin-films grown by vapor deposition methods that normally have the [001] polarization direction perpendicular to the film. It is also found that increasing the electric field that is switched broadens the range of applicable temperatures where a large ECE can be achieved.
These observations are physically intuitive. The applied electric field acts to elevate the transformation temperature so that for T > T C the crystal remains in the ordered ferroelectric phase. When the field is removed the crystal transforms to the disordered, paraelectric state and the change in the entropy results in an adiabatic change in temperature. The magnitude of the temperature change is directly proportional to the magnitude of the change in electric field. Although applying an electric field in arbitrary directions does act to stabilize the ordered phase above T C , the largest ∆T and widest range of applicable temperatures is achieved when the external field is applied parallel to the polarization direction of the first ferroelectric phase just below the Curie temperature. In the case of BaTiO 3 this is the [001] direction, but for other ferroelectric crystals the direction of applied field should be taken from knowledge of their phase diagrams. Because the electric field acts to raise the transformation temperature, having a larger switching field naturally broadens the range of applicable temperatures.  In frame (a) the external electric field switches from Ez = 5-50 to 0 kV/cm and in frame (b) it switches from Ez = 50-500 to 0 kV/cm. There is an ordinate for both the raw data ∆T and for the corrected data ∆T corrected = 1 5 ∆T that accounts for the underestimated specific heat capacity.