Colossal barocaloric effects in the complex hydride Li$_{2}$B$_{12}$H$_{12}$

Traditional refrigeration technologies based on compression cycles of greenhouse gases pose serious threats to the environment and cannot be downscaled to electronic device dimensions. Solid-state cooling exploits the thermal response of caloric materials to external fields and represents a promising alternative to current refrigeration methods. However, most of the caloric materials known to date present relatively small adiabatic temperature changes ($|\Delta T| \sim 1$ K) and/or limiting irreversibility issues resulting from significant phase-transition hysteresis. Here, we predict the existence of colossal barocaloric effects (isothermal entropy changes of $|\Delta S| \sim 100$ JK$^{-1}$kg$^{-1}$) in the energy material Li$_{2}$B$_{12}$H$_{12}$ by means of molecular dynamics simulations. Specifically, we estimate $|\Delta S| = 387$ JK$^{-1}$kg$^{-1}$ and $|\Delta T| = 26$ K for an applied pressure of $P = 0.4$ GPa at $T = 475$ K. The disclosed colossal barocaloric effects are originated by an order-disorder phase transformation that exhibits a fair degree of reversibility and involves coexisting Li$^{+}$ diffusion and (BH)$_{12}^{-2}$ reorientational motion at high temperatures.

Traditional refrigeration technologies based on compression cycles of greenhouse gases pose serious threats to the environment and cannot be downscaled to electronic device dimensions. Solid-state cooling exploits the thermal response of caloric materials to external fields and represents a promising alternative to current refrigeration methods. However, most of the caloric materials known to date present relatively small adiabatic temperature changes (|∆T | ∼ 1 K) and/or limiting irreversibility issues resulting from significant phase-transition hysteresis. Here, we predict the existence of colossal barocaloric effects (isothermal entropy changes of |∆S| ∼ 100 JK −1 kg −1 ) in the energy material Li 2 B 12 H 12 by means of molecular dynamics simulations.
Specifically, we estimate |∆S| = 387 JK −1 kg −1 and |∆T | = 26 K for an applied pressure of P = 0.4 GPa at T = 475 K. The disclosed colossal barocaloric effects are originated by an order-disorder phase transformation that exhibits a fair degree of reversibility and involves coexisting Li + diffusion and (BH) −2 12 reorientational motion at high temperatures.
Solid-state cooling is an environmentally friendly, energy efficient, and highly scalable technology that can solve most of the problems associated with conventional refrigeration methods based on compression cycles of greenhouse gases (i.e., environmental harm and lack of downsize scaling). Upon application of magnetic, electric or stress fields good caloric materials undergo noticeable temperature changes (|∆T | ∼ 1-10 K) as a result of induced phase transformations that involve large entropy variations (|∆S| ∼ 10-100 JK −1 kg −1 ) [1][2][3]. Solidstate cooling capitalizes on such caloric effects to engineer refrigeration cycles. From a performance point of view, that is, largest |∆T | and |∆S| (although these are not the only parameters defining cooling efficiency [4]), barocaloric effects driven by small hydrostatic pressure shifts appear to be the most promising [1][2][3]. * Corresponding Author FIG. 1. Low-T (ordered) and high-T (disordered) phases of bulk Li2B12H12. The low-T phase (α) presents cubic symmetry and space group P a3 [15]. In the high-T phase (β), cubic symmetry is preserved but the Li + ions are highly mobile and the (BH) −2 12 icosahedra present reorientational disorder [18]. The T -induced α → β phase transition is an order-disorder isosymmetric transformation [18]. Li, B, and H ions are represented with red, blue, and yellow colours, respectively.
Recently, colossal barocaloric effects (defined here as |∆S| ∼ 100 JK −1 kg −1 ) have been measured in two different families of materials that display intriguing orderdisorder phase transitions [5][6][7]. First, giant barocaloric effects have been theoretically predicted [8] and experimentally observed in the archetypal superionic compound AgI [5]. AgI exhibits a first-order normal (lowentropy) to superionic (high-entropy) phase transition that responds to both temperature and pressure [9] and which involves the presence of highly mobile silver ions in the high-T superionic state [10]. The entropy changes estimated for other normal to superionic phase transitions in general are large as well [11][12][13][14]. And second, colossal barocaloric effects have been reported for the molecular solid neopentylglycol [6,7], (CH 3 ) 2 C(CH 2 OH) 2 , and other plastic crystals [4]. In these solids molecules reorient almost freely around their centers of mass, which remain localized at well-defined lattice positions. Molecular rotations lead to orientational disorder, which renders high entropy. By using hydrostatic pressure, it is possible to block such molecular reorientational motion and thus induce a fully ordered state characterized by low entropy [14]. The barocaloric effects resulting from this class of order-disorder phase transition are huge and comparable in magnitude to those achieved in conventional refrigerators with environmentally harmful fluids [4,6,7].
Here, we report the prediction of colossal barocaloric effects (|∆S| ∼ 100 JK −1 kg −1 ) in the energy material Li 2 B 12 H 12 (LBH), a complex hydride that is already known from the fields of hydrogen storage [15][16][17] and solid-state batteries [18][19][20]. By using molecular dynamics simulations, we identify a pressure-induced isothermal entropy change of |∆S| = 387 JK −1 kg −1 and adiabatic temperature change of |∆T | = 26 K at T = 475 K. These colossal entropy and temperature changes are driven by moderate hydrostatic pressure shifts of P = 0.4 GPa, thus yielding huge barocaloric strengths of |∆S|/P = 968 JK −1 kg −1 GPa −1 and |∆T |/P = 65 K GPa −1 . The colossal barocaloric effects disclosed in bulk LBH are originated by simultaneous P -driven frustration and activation of Li + diffusion and (BH) −2 12 icosahedra reorientational motion. Thus, alkali-metal complex borohydrides (A 2 B 12 H 12 , A = Li, Na, K, Cs [21,22]) emerge as a promising new family of barocaloric materials in which the salient phase-transition features of fast-ion conductors and plastic crystals coexist.

RESULTS
At ambient conditions, lithium dodecahydrododecaborate (Li 2 B 12 H 12 ), LBH, presents an ordered cubic P a3 phase, referred to as α hereafter, which is characterized by Li + cations residing on near-trigonal-planar sites surrounded by three (BH) −2 12 icosahedron anions. In turn, each (BH) −2 12 anion resides in an octahedral cage surrounded by six Li + cations (Fig.1a) [15]. A symmetry preserving order-disorder phase transition occurs at high temperatures (∼ 500 K) that stabilises a disordered state, referred to as β hereafter, in which the Li + cations are mobile and the (BH) −2 12 anions present reorientational motion ( Fig.1b) [18]. The relative volume expansion that has been experimentally measured for such an orderdisorder phase transition is ≈ 8% [18]. This huge volume variation along with the accompanying, and pressumably also large, phase-transition entropy change could be propitious for barocaloric purposes if the involved phase transformation was responsive to moderate external pressures of ∼ 0.1 GPa. To the best of our knowledge, this possibility has not been hitherto explored. We performed classical molecular dynamics (MD) simulations based on a recently proposed LBH force field [23] to fill up such a knowledge gap (Methods and Supplementary Methods), which has clear implications for potential solid-state cooling applications. Figure 2a shows the P -T phase diagram that we estimated for bulk LBH using atomistic MD simulations. It was found that the temperature of the α → β phase transition is certainly sensitive to external pressure. Specifically, the dP/dT derivative of the corresponding phase boundary amounts to ≈ 0.008 GPa K −1 at zero pressure and to ≈ 0.02 GPa K −1 at P = 0.2 GPa. Likewise, the relative volume change ascribed to the α → β transformation is, according to our simulations, +4.6% at zero pressure and +3.4% at P = 0.2 GPa (Fig.2b). By using these thermodynamic data and the Clausius-Clapeyron relation [2], we roughly estimated an entropy change of ∆S ∼ 300 JK −1 kg −1 for the order-disorder transition occurring in LBH at P = 0.2 GPa. In view of these promising barocaloric descriptor values, we proceeded to accurately calculate the barocaloric isothermal entropy and adiabatic temperature changes, ∆S and ∆T , induced by pressures 0 ≤ P ≤ 0.4 GPa. To this end, we followed the numerical protocols described in the Methods section, which essentially involve the determination of the volume and heat capacity of bulk LBH (Fig.2c) as a function of pressure and temperature.
The results of our precise barocaloric calculations for temperatures and pressures in the intervals 450 ≤ T ≤ 525 K and 0 ≤ P ≤ 0.4 GPa are shown in Figs. 2d,e. The ∆S and ∆T values estimated for the α → β transformation in fact render colossal barocaloric effects. For example, at T = 490 K and P = 0.4 GPa (0.2 GPa) we calculated an isothermal entropy change of −365 JK −1 kg −1 (−135 JK −1 kg −1 ) and an adiabatic temperature change of +27 K (+10 K). The resulting barocaloric effects are direct, that is, ∆T > 0, because the low-entropy ordered state is stabilized under pressure (∆S < 0). A maximum |∆S| value of 387 JK −1 kg −1 was found at T = 475 K and P = 0.4 GPa (Fig.2d). For temperatures above ≈ 510 K, we estimated noticeably smaller |∆S| and ∆T values (e.g., 72 JK −1 kg −1 and 10 K for P = 0.4 GPa at T = 525 K), a trend that we link to some anomalous pressure-induced ionic diffusion (explained below). In the Discussion section, we will compare the barocaloric performance of LBH with those of other well-known barocaloric materials. In what follows, the atomistic mechanisms leading to the extraordinary ∆S and ∆T results just reported are unravelled.
There are two possible sources of large entropy variation in LBH, one stemming from the Li + ionic diffusion and the other from the (BH) −2 12 icosahedra reorientational motion. When hydrostatic pressure is applied on the disordered β phase at temperatures below ≈ 500 K, both the ionic diffusion and molecular orientational disorder are reduced and thus the crystal entropy diminishes significantly. This conclusion is straightforwardly deduced from the P -induced variation of the Li + diffusion coefficient, D Li , and reorientational (BH) −2 12 frequency, λ B12H12 , shown in Figs.3a,b (Methods). For instance, at T = 475 K and zero pressure D Li and λ B12H12 amount to 2.5·10 −6 cm 2 s −1 and 1.2·10 8 s −1 , respectively, whereas at P = 0.2 GPa both quantities are practically zero (Fig.3). The two resulting contributions to the system entropy variation are of the same sign and make |∆S| huge.
Which of these two P -induced order-restoring effects is most relevant for the barocaloric performance of bulk LBH? To answer this question, we performed constrained MD simulations in which we forced the Li + ions to remain localized around their equilibrium positions independently of temperature. This type of artificial condition in principle cannot be imposed in the experiments but can be easily enforced in the atomistic simulations.
The |∆S| values estimated in these constrained MD simulations were roughly half the value of the isothermal entropy changes obtained in the standard MD simulations. Therefore, we may conclude that at temperatures below ≈ 500 K the pressure-induced entropy changes stemming from the Li + ionic diffusion and (BH) −2 12 icosahedra reorientational motion variations play both an equally important role in the global barocaloric response of LBH. Figure 3a shows that at T 500 K the Li + diffusion coefficient increases under increasing pressure. For example, at T = 525 K and zero pressure we estimate D Li = 8.7 · 10 −6 cm 2 s −1 whereas at P = 0.4 GPa and the same temperature we obtain 17.2·10 −6 cm 2 s −1 . This ionic diffusion behaviour is highly anomalous because hydrostatic compression typically hinders ionic transport [9][10][11]. On the other hand, the reorientational motion of the (BH) −2 12 icosahedra behaves quite normally, that is, decreases under pressure [3,4]. For instance, at T = 525 K and zero pressure we estimate λ B12H12 = 1.4 · 10 8 s −1 whereas at P = 0.4 GPa and the same temperature we obtain 0.7 · 10 8 s −1 (Fig.3b). We hypothesize that the anomalous P -induced Li + diffusion behaviour observed in our MD simulations is due to the high anionic reorientational motion, which makes the (BH) −2 12 centers of mass to fluctuate and partially block the ionic current channels [24]. Consistently, when the frequency of the (BH) −2 12 rotations is reduced by effect of compression the ions can flow more easily throughout the crystal and Li + transport is enhanced. In this particular P -T region, the two contributions to the crystal entropy variation stemming from Li + ionic diffusion and (BH) −2 12 icosahedra reorientational motion have opposite signs hence |∆S| decreases significantly. The identified anomalous lithium diffusion behaviour, however, ceases at P ≈ 0.6 GPa since beyond that point D Li decreases systematically upon increasing pressure ( Supplementary Fig.1).
recently reported for the plastic crystal neopentylglycol by considering a similar pressure shift, 510 JK −1 kg −1 [6,7]. The rest of materials in Table I present isothermal entropy changes that are appreciably smaller, made the exception of the polar crystal (NH 4 ) 2 SO 4 which registers 130 JK −1 kg −1 . As regards |∆T |, the clear contestants of LBH are the fast-ion conductor AgI (36 K) and again the plastic crystal (CH 3 ) 2 C(CH 2 OH) 2 (45 K). The reason for the smaller |∆T | value estimated for LBH as compared to that of AgI is the significantly larger heat capacity of the former material, which results from a smaller molecular weight [12]. In terms of the barocaloric strengths defined as BSS ≡ |∆S|/P and BST ≡ |∆T |/P , LBH remains competitive with the best performers.
As it was mentioned in the Introduction, the magnitude of the |∆T | and |∆S| shifts are not the only parameters that define the barocaloric performance of a material. The degree of reversibility of the involved Pinduced phase transition, for instance, is another impor-tant barocaloric descriptor that provides information on the materials efficiency during successive pressure application/removal cycles. Specifically, the hysteresis of the transition makes the materials behaviour to depend on its cycling history and to increase the value of the external field that is required to bring the phase transition to completion [4]. As a consequence, the barocaloric performance of a hysteretic material can be significantly worse than that of its ideal non-hysteretic counterpart. In order to quantify the degree of reversibility associated with the α ↔ β phase transition in LBH, we performed a series of long MD simulations (∼ 2 ns) in which the pressure (temperature) was kept fixed while the temperature (pressure) was varied steadily first from 425 up to 625 K (from 0.0 up to 0.4 GPa) and subsequently from 625 back to 425 K (from 0.4 back to 0.0 GPa). The results of such field-changing simulations indicate that the degree of reversibility of the order-disorder α → β phase transition is quite acceptable (Supplementary Fig.2). For instance, by monitoring the variation of the system volume, we found that at zero pressure the difference between the transition temperatures observed during the heating and cooling stages was ∆T h ≡ T α→β − T β→α ≈ 50 K (Supplementary Fig.2a). The size of ∆T h , however, increases noticeably at higher pressures (≈ 100 K at 0.4 GPa). Meanwhile, at fixed temperature we found that the hysteresis of the phase transition as driven by pressure was practically null at T = 550 K (∆P h ≈ 0 GPa) and equal to 0.1 GPa at 475 K ( Supplementary Fig.2b).
Arguably the only weakness of bulk LBH in terms of barocaloric potential is that the critical temperature of the order-disorder α → β phase transition is significantly higher than room temperature. However, this practical problem can be efficiently solved by means of doping and alloying strategies. In fact, recently it has been experimentally shown that carbon-doped LBH, LiCB 11 H 12 , presents a much lower α → β transition temperature of ≈ 400 K [33], and that the disordered β phase is already stabilized at room temperature in Li(CB 9 H 10 )-Li(CB 11 H 12 ) solid solutions [34]. Moreover, the type of isosymmetric order-disorder phase transition underlying the exceptional barocaloric behaviour of LBH occurs also in analogous alkali-metal complex hydrides (A 2 B 12 H 12 , A = Na, K, Cs) [35] and other earth-abundant and nontoxic materials like KHPO 4 , NaAlSi 3 O 8 and KNO 3 [36]. Bulk KNO 3 , for example, displays a staggering volume collapse of ∼ 10% for a room-temperature phase transformation induced by a modest pressure of 0.3 GPa [37], which suggests great barocaloric potential as well.
In conclusion, we have predicted the existence of colossal barocaloric effects rendering isothermal entropy changes of ∼ 100 JK −1 kg −1 and adiabatic temperature shifts of ∼ 10 K in the complex hydride Li 2 B 12 H 12 , which are driven by moderate hydrostatic pressures of ∼ 0.1 GPa. The phase transition underlying such colossal barocaloric effects is remarkable as it combines key ingredients of fast-ion conductors (i.e., ionic diffusion) and molecular crystals (i.e., reorientational motion), ma-terials that individually have been proven to be excellent barocaloric materials. This same type of isosymmetric order-disorder phase transition is likely to occur also in other economically affordable and innocuous compounds (e.g., Cs 2 B 12 H 12 and KNO 3 ), thus broadening significantly the spectrum of caloric materials with commercial potential for solid-state cooling applications. We believe that our simulation study will stimulate experimental research on this new family of barocaloric materials, namely, alkali-metal complex hydrides, which are already known from other technological disciplines (e.g., hydrogen storage and electrochemical devices) and are routinely synthesized in the laboratory.

METHODS
Classical molecular dynamics simulations. Molecular dynamics (MD) (N, P, T ) simulations were performed with the LAMMPS code [38]. The pressure and temperature in the system were kept fluctuating around a set-point value by using thermostatting and barostatting techniques in which some dynamic variables are coupled to the particle velocities and simulation box dimensions. The interactions between atoms were modeled with the harmonic Coulomb-Buckingham interatomic potential reported in work [23], the details of which are provided in the Supplementary Methods. The employed interatomic potential reproduces satisfactorily the vibrational spectra, structure and lithium diffusion coefficients of bulk LBH [23] (Supplementary Discussion). We employed simulation boxes containing 6656 atoms and applied periodic boundary conditions along the three Cartesian directions. Newton's equations of motion were integrated using the customary Verlet's algorithm with a time-step length of 0.5 fs. The typical duration of a MD (N, P, T ) run was of 1 ns. A particle-particle particle-mesh k-space solver was used to compute long-range van der Waals and Coulomb interactions beyond a cut-off distance of 10Å at each time step.
Density functional theory calculations. Firstprinciples calculations based on density functional theory (DFT) [39] were performed to analyse the energy, structural, vibrational, and ionic transport properties of Li 2 B 12 H 12 . We performed these calculations with the VASP software [40] by following the generalized gradient approximation to the exchange-correlation energy due to Perdew et al. [41]. The projector augmented-wave method was used to represent the ionic cores [42], and the electronic states 1s-2s Li, 1s-2s-2p B and 1s H were considered as valence. Wave functions were represented in a plane-wave basis set truncated at 650 eV. By using these parameters and dense k-point grids for Brillouin zone integration, the resulting energies were converged to within 1 meV per formula unit. In the geometry relaxations, a tolerance of 0.01 eV·Å −1 was imposed on the atomic forces.
Ab initio molecular dynamics (AIMD) simulations based on DFT were carried out to assess the reliability of the interatomic potential model employed in the classical molecular dynamics simulations (Supplementary Fig.3 and Supplementary Discussion). The AIMD simulations were performed in the canonical (N, V, T ) ensemble considering constant number of particles, volume and temperature. The constrained volumes were equal to the equilibrium volumes determined at zero temperature, thus we neglected possible thermal expansion effects. Nevertheless, in view of previous first-principles work [43], it is reasonable to expect that thermal expansion effects do not affect significantly the estimation of lithium diffusion coefficients at the considered temperatures. The temperature in the AIMD simulations was kept fluctuating around a set-point value by using Nose-Hoover thermostats.
A large simulation box containing 832 atoms was employed in all the simulations, and periodic boundary conditions were applied along the three Cartesian directions. Newton's equations of motion were integrated by using the customary Verlet's algorithm and a time-step length of δt = 10 −3 ps. Γ-point sampling for integration within the first Brillouin zone was employed in all the AIMD simulations. The AIMD simulations comprised long simulation times of ∼ 100 ps.
Estimation of key quantities. The mean square displacement of lithium ions was estimated with the formula [43]: where r i (t j ) is the position of the migrating ion i at time t j (= j ·δt), τ represents a lag time, n τ = τ /δt, N ion is the total number of mobile ions, and N step the total number of time steps. The maximum n τ was chosen equal to N step /2, hence we could accumulate enough statistics to reduce significantly the fluctuations in MSD Li (τ ) at large τ 's. The diffusion coefficient of lithium ions then was obtained with the Einstein relation: by performing linear fits to the averaged MSD Li values calculated at long τ . The angular autocorrelation function of the closoborane (BH) 2− 12 icosahedra was estimated according to the expression [23]: φ B12H12 (τ ) = r(t) ·r(t + τ ) , wherer is a unitary vector connecting the center of mass of each closoborane unit with one of its edges and · · · denotes thermal average considering all the closoborane icosahedra. This autocorrelation function typically decays as ∝ exp [−λ B12H12 · τ ], where the parameter λ B12H12 represents a characteristic reorientational frequency. When the (BH) 2− 12 reorientational motion is significant, that is, λ B12H12 is large, the φ B12H12 function decreases rapidly to zero with time.
Isothermal entropy changes associated with the barocaloric effect were estimated with the formula [2,3]: where P represents the maximum applied hydrostatic pressure and V the volume of the system. Likewise, the accompanying adiabatic temperature shift was calculated as: where C P (T ) = dU dT P is the heat capacity of the crystal obtained at constant pressure and temperature conditions.
In order to accurately compute the ∆S(P, T ) and ∆T (P, T ) shifts induced by pressure, we calculated the corresponding volumes and heat capacities over dense grids of (P, T ) points spaced by δP = 0.1 GPa and δT = 25 K. Spline interpolations were subsequently applied to the calculated sets of points, which allowed for accurate determination of (∂V /∂T ) P and heat capacities. The ∆S and ∆T values appearing in Fig.2d-e were obtained by numerically integrating those spline functions with respect to pressure.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author (C.C.) upon reasonable request.