# Colossal reversible barocaloric effects in a plastic crystal mediated by
lattice vibrations and ion diffusion

Ming Zeng  Grup de Caracterizació de Materials, Departament de Física, EEBE
and Barcelona Research Center in Multiscale Science and Engineering
Universitat Politècnica de Catalunya, Av. Eduard Maristany 10-14, Barcelona
08019, Catalonia, Spain  Carlos Escorihuela-Sayalero  Grup de Caracterizació
de Materials, Departament de Física, EEBE and Barcelona Research Center in
Multiscale Science and Engineering Universitat Politècnica de Catalunya, Av.
Eduard Maristany 10-14, Barcelona 08019, Catalonia, Spain  Tamio Ikeshoji
Mathematics for Advanced Materials Open Innovation Laboratory (MathAM-OIL),
National Institute of Advanced Industrial Science and Technology (AIST), c/o
Advanced Institute for Materials Research (AIMR), Tohoku University, Sendai
980-8577, Japan  Shigeyuki Takagi  Institute for Materials Research (IMR),
Tohoku University, Sendai 980-8577, Japan  Sangryun Kim  Graduate School of
Energy Convergence, Gwangju Institute of Science and Technology (GIST), 123
Cheomdangwagi-ro, Buk-gu, Gwangju 61005, Republic of Korea  Shin-ichi Orimo
Advanced Institute for Materials Research (AIMR), Tohoku University, Sendai
980-8577, Japan  Institute for Materials Research (IMR), Tohoku University,
Sendai 980-8577, Japan  María Barrio  Grup de Caracterizació de Materials,
Departament de Física, EEBE and Barcelona Research Center in Multiscale
Science and Engineering Universitat Politècnica de Catalunya, Av. Eduard
Maristany 10-14, Barcelona 08019, Catalonia, Spain  Josep-Lluís Tamarit  Grup
de Caracterizació de Materials, Departament de Física, EEBE and Barcelona
Research Center in Multiscale Science and Engineering Universitat Politècnica
de Catalunya, Av. Eduard Maristany 10-14, Barcelona 08019, Catalonia, Spain
Pol Lloveras<EMAIL_ADDRESS>Grup de Caracterizació de Materials,
Departament de Física, EEBE and Barcelona Research Center in Multiscale
Science and Engineering Universitat Politècnica de Catalunya, Av. Eduard
Maristany 10-14, Barcelona 08019, Catalonia, Spain  Claudio Cazorla
<EMAIL_ADDRESS>Grup de Caracterizació de Materials, Departament de
Física, EEBE and Barcelona Research Center in Multiscale Science and
Engineering Universitat Politècnica de Catalunya, Av. Eduard Maristany 10-14,
Barcelona 08019, Catalonia, Spain  Kartik Sau<EMAIL_ADDRESS>Advanced
Institute for Materials Research (AIMR), Tohoku University, Sendai 980-8577,
Japan  Mathematics for Advanced Materials Open Innovation Laboratory (MathAM-
OIL), National Institute of Advanced Industrial Science and Technology (AIST),
c/o Advanced Institute for Materials Research (AIMR), Tohoku University,
Sendai 980-8577, Japan

###### Abstract

Solid-state methods for cooling and heating promise a more sustainable
alternative to current compression cycles of greenhouse gases and inefficient
fuel-burning heaters. Barocaloric effects (BCE) driven by hydrostatic pressure
($p$) are especially encouraging in terms of large adiabatic temperature
changes ($|\Delta T|\sim 10$ K) and colossal isothermal entropy changes
($|\Delta S|\sim 100$ J K-1 kg-1). However, BCE typically require large
pressure shifts due to irreversibility issues, and sizeable $|\Delta T|$ and
$|\Delta S|$ seldom are realized in a same material. Here, we demonstrate the
existence of colossal and reversible BCE in LiCB11H12, a well-known solid
electrolyte, near its order-disorder phase transition at $\approx 380$ K.
Specifically, for $\Delta p\approx 0.23$ (0.10) GPa we measured $|\Delta
S_{\rm rev}|=280$ (200) J K-1 kg-1 and $|\Delta T_{\rm rev}|=32$ (10) K, which
individually rival with state-of-the-art barocaloric shifts obtained under
similar pressure conditions. Furthermore, over a wide temperature range,
pressure shifts of the order of $0.1$ GPa yield huge reversible barocaloric
strengths of $\approx 2$ J K-1 kg-1 MPa-1. Molecular dynamics simulations were
carried out to quantify the role of lattice vibrations, molecular
reorientations and ion diffusion on the disclosed colossal BCE. Interestingly,
lattice vibrations were found to contribute the most to $|\Delta S|$ while the
diffusion of lithium ions, despite adding up only slightly to the accompanying
entropy change, was crucial in enabling the molecular order-disorder phase
transition. Our work expands the knowledge on plastic crystals and should
motivate the investigation of BCE in a variety of solid electrolytes
displaying ion diffusion and concomitant molecular orientational disorder.

solid-state refrigeration, barocaloric effects, orientational order-disorder
phase transition, lithium diffusion, molecular dynamics simulations

##  I Introduction

Solid-state methods for cooling and heating are energy efficient and
ecologically friendly techniques with potential for solving the environmental
problems posed by conventional refrigeration and heat pump technologies
relying on compression cycles of greenhouse gases and inefficient traditional
fuel-burning heaters Moya and Mathur (2020). Under moderate magnetic, electric
or mechanical field variations, auspicious caloric materials experience large
adiabatic temperature variations ($|\Delta T|\sim 1$–$10$ K) as a result of
phase transformations entailing large isothermal entropy changes ($|\Delta
S|\sim 10$–$100$ J K-1 kg-1) Mañosa  _et al._ (2013); Hou  _et al._ (2022).
Solid-state cooling and heat pumping capitalize on such caloric effects for
engineering refrigeration and heating cycles. From a practical point of view,
large and reversible $|\Delta T|$ and $|\Delta S|$ are both necessary for
achieving rapid and efficient devices under recursive application and removal
of the driving fields. In terms of largest $|\Delta T|$ and $|\Delta S|$,
mechanocaloric effects induced by uniaxial stress (elastocaloric effects) and
hydrostatic pressure (barocaloric effects –BCE–) are among the most promising
Mañosa and Planes (2017); Cazorla (2019); Lloveras and Tamarit (2021).

Recently, colossal and reversible BCE ($|\Delta S_{\rm rev}|\geq 100$ J K-1
kg-1) have been measured in several families of materials displaying order-
disorder phase transitions under pressure shifts of the order of $0.1$ GPa
Lloveras  _et al._ (2019); Li  _et al._ (2019); Aznar  _et al._ (2020, 2021);
Zhang  _et al._ (2022); Imamura  _et al._ (2020); Li  _et al._ (2021); Seo
_et al._ (2022); Salvatori  _et al._ (2022). On one hand, there are plastic
crystals like neopentane derivatives Lloveras  _et al._ (2019); Li  _et al._
(2019); Aznar  _et al._ (2020), adamantane derivatives Aznar  _et al._ (2021);
Salvatori  _et al._ (2022) and carboranes Zhang  _et al._ (2022) in which the
underlying phase transitions involve molecular orientational disorder
stabilized under increasing temperature. On the other hand, there are polymers
(e.g., acetoxy silicone rubber) Imamura  _et al._ (2020) and layered hybrid
organic-inorganic perovskites (e.g., [C10H21NH3]2MnCl4) Li  _et al._ (2021);
Seo  _et al._ (2022) in which the accompanying phase transformations entail
significant atomic rearrangements in the organic components. Another family of
disordered materials presenting also great barocaloric promise are solid
electrolytes (e.g., AgI, Li3N and Cu2Se) Aznar  _et al._ (2017); Sagotra  _et
al._ (2017, 2018); Min  _et al._ (2020), although in this latter case the
experimentally reported $|\Delta S_{\rm rev}|$ fall slightly below the
colossal threshold value of $100$ J K-1 kg-1 Aznar  _et al._ (2017).

Figure 1: Sketch of the order-disorder phase transition occurring in LCBH upon
increasing temperature. (a) Ball-stick representation of the low-$T$ ordered
(O) and high-$T$ disordered (D) phases. Lithium, carbon, boron and hydrogen
atoms are represented with red, brown, green and blue spheres, respectively.
In the high-$T$ phase, the Li+ cations diffuse throughout the crystalline
matrix while the [CB11H12]- anions reorient disorderly Tang  _et al._ (2015);
the volume increases significantly during the $T$-induced phase transition.
(b) Outline of the order-disorder phase transition in terms of Gibbs free
energies. The red dotted lines represent internal energies and the blue solid
lines Gibbs free energies; $T_{t}$ denotes the phase transition temperature.

In spite of these recent developments, finding barocaloric materials with
well-balanced and suitable features for developing thermal applications, e.g.,
$|\Delta T_{\rm rev}|\geq 20$ K and $|\Delta S_{\rm rev}|\geq 100$ J K-1 kg-1
driven by $\Delta p\lesssim 0.1$ GPa, is proving extremely difficult. From the
hundred of barocaloric materials known to date Lloveras and Tamarit (2021), to
the best of our knowledge only four fulfill the conditions specified above,
namely, the spin-crossover complex Fe3(bntrz)6(tcnset)6 ($|\Delta T_{\rm
rev}|=35$ K and $|\Delta S_{\rm rev}|=120$ J K-1 kg-1 for $\Delta p=0.26$ GPa)
Romanini  _et al._ (2021), the layered hybrid perovskite [C10H21NH3]2MnCl4
($|\Delta T_{\rm rev}|=27$ K and $|\Delta S_{\rm rev}|=250$ J K-1 kg-1 for
$\Delta p=0.19$ GPa) Li  _et al._ (2021); Seo  _et al._ (2022), the plastic
crystal 1-Br-adamantane ($|\Delta T_{\rm rev}|=20$ K and $|\Delta S_{\rm
rev}|=120$ J K-1 kg-1 for $\Delta p=0.10$ GPa) Aznar  _et al._ (2021), and the
elastomer acetoxy silicone ($|\Delta T_{\rm rev}|=22$ K and $|\Delta S_{\rm
rev}|=182$ J K-1 kg-1 for $\Delta p=0.17$ GPa) Imamura  _et al._ (2020).
Moreover, studies addressing a fundamental and quantitative understanding of
the atomistic mechanisms that bring on such colossal BCE are very scarce Li
_et al._ (2020a, 2022); de Oliveira (2023); Sau  _et al._ (2021a), thus
hindering the rational design of disordered materials with enhanced
barocaloric performances.

Figure 2: Experimental phase diagram of bulk LCBH and corresponding phase
transition entropy changes. (a) Volume per formula unit measured as a function
of temperature at normal pressure. (b) Isobaric heat flow data expressed as a
function of applied pressure and temperature; data collected during heating
(cooling) are represented in the positive (negative) y-axis. (c) Pressure and
temperature phase diagram; transition temperatures are determined from the
peaks in panel (b). (d) Phase transition entropy changes as a function of
pressure and transition path. $\Delta S_{\text{t}}$ remains practically
constant from atmospheric pressure all the way up to the triple point. At
$p\simeq 0.13$ GPa, $\Delta S_{\text{II}\rightarrow\text{I}}\approx\Delta
S_{\text{II}\rightarrow\text{III}}+\Delta S_{\text{III}\rightarrow\text{I}}$,
while above the triple point $\Delta
S_{\text{II}\rightarrow\text{III}}\approx\Delta
S_{\text{III}\rightarrow\text{I}}$. Straight lines at pressures above the
triple point are linear fits to $\Delta
S_{\text{II}\rightarrow\text{III}}+\Delta S_{\text{III}\rightarrow\text{I}}$.

In this work, we experimentally and theoretically demonstrate the existence of
colossal and reversible BCE in the monocarba-closo-dodecaborate LiCB11H12
(LCBH) near its order-disorder phase transition occurring at $T_{t}\approx
380$ K Tang  _et al._ (2015). LCBH is a well-known solid electrolyte in which
at temperatures above $T_{t}$ the lithium cations are highly mobile and the
molecular anions [CB11H12]- reorient disorderly Skripov  _et al._ (2015);
Mohtadi and Orimo (2016) (Fig. 1); thus, LCBH combines phase-transition
features of both plastic crystals and superionic compounds, two families of
materials for which colossal and giant BCE, respectively, have been previously
reported Lloveras  _et al._ (2019); Li  _et al._ (2019); Aznar  _et al._
(2020, 2017). In particular, we measured colossal values of $|\Delta T_{\rm
rev}|=32$ K and $|\Delta S_{\rm rev}|=280$ JK-1kg-1 for a pressure shift of
$0.23$ GPa, and large and reversible barocaloric strengths of $\approx 2$ J
K-1 kg-1 MPa-1 over a wide temperature interval of several tens of degrees.
Likewise, for a smaller pressure shift of $0.10$ GPa assuring values of
$|\Delta S_{\rm rev}|=200$ J K-1 kg-1 and $|\Delta T_{\rm rev}|=10$ K were
obtained. Atomistic molecular dynamics simulations were performed to reveal
key phase transition mechanisms and quantify the role played by the
vibrational, molecular orientational and ion diffusive degrees of freedom on
the disclosed BCE. Very interestingly, the contribution of the lattice
vibrations to $\Delta S$ was found to be the dominant at all pressures,
instead of the typically assumed one resulting from molecular reorientational
motion Li  _et al._ (2020a, 2022); de Oliveira (2023). Our results provide new
valuable insights into the physical behavior and functionality of plastic
crystals and suggest that colossal BCE similar to those reported here for LCBH
could also exist in other akin closo-borate materials like NaCB11B12 Tang  _et
al._ (2015); Skripov  _et al._ (2015), KCB11B12 Dimitrievska  _et al._ (2020),
and LiCB9H10 Kim  _et al._ (2019, 2020).

Figure 3: Experimentally measured colossal barocaloric effects in bulk LCBH.
(a)–(c) Isothermal entropy change, $\Delta S$, and (b)–(d) adiabatic
temperature change, $\Delta T$, obtained upon the application and removal of
pressure, $p$, considering (a)–(b) irreversible and (c)–(d) reversible
processes. Figure 4: Compendium of experimentally measured reversible BCE. The
size of the symbols represents the reversible barocaloric strength defined as
the ratio of $|\Delta S_{\rm rev}|$ by the corresponding pressure change
$\Delta p$. Material names are indicated near each symbol or in the right side
of the panel. NPG: neopentylglycol; PG: pentaglycerine; NPA: Neopentyl
alcohol; o-carb: orthocarborane; m-carb: metacarborane; p-carb: paracarborane;
1-Br-ada: 1-Bromoadamantane; 1-Cl-ada: 1-Chloroadamantane; 1ada-ol:
1-adamantanol; 2ada-ol: 2-adamantanol; 2m2ada-ol: 2-methyl-2-adamantanol; ASR:
Acetoxy Silicone Rubber. Numerical details and references can be found in the
Supplementary Table S1.

##  II Results and Discussion

###  II.1 LiCB11H12 general properties

In a recent X-ray powder diffraction study Tang  _et al._ (2015), it has been
shown that at room temperature LiCB11H12 (LCBH) presents an ordered
orthorhombic structure (space group $Pca2_{1}$) in which the Li+ cations
reside near trigonal-planar sites surrounded by molecular [CB11H12]- anions
arranged in a cubic sublattice. An order-disorder phase transition occurs at
$T_{t}\approx 380$ K that stabilizes a disordered phase in which the Li+
cations are highly mobile and the molecular anions present fast
reorientational motion (Fig. 1a). At normal pressure, the lithium ion
conductivity measured just above $T_{t}$ exceeds values of $0.1$ S cm-1 Tang
_et al._ (2015) and the reorientational motion of the molecular anions can
reach frequencies of $10^{11}$ s-1 Tang  _et al._ (2015); Sau  _et al._
(2021b). Meanwhile, the $T$-induced order-disorder phase transition is
accompanied by a huge volume increase of the order of $\approx 10$% Sau  _et
al._ (2021b) that, based on the Clausius-Clapeyron (CC) equation $\Delta
S_{t}=\Delta V_{t}\frac{dp}{dT}$, suggests great barocaloric potential.

The described order-disorder phase transition can be qualitatively understood
in terms of the Gibbs free energy difference between the high-$T$ disordered
(D) and low-$T$ ordered (O) phases, $\Delta G\equiv G^{D}-G^{O}$ (Fig. 1b).
This free energy difference consists of an internal energy ($\Delta E$),
entropy ($-T\Delta S$), and volume ($p\Delta V$) terms. The internal energy
remains more or less constant during the phase transition while the volume
term disfavors the stabilization of the disordered phase since $\Delta V>0$.
Thus, the LCBH order-disorder phase transition appears to be governed by the
change in entropy, $\Delta S$, which in view of the ion conductivity and
molecular reorientational frequency measured above $T_{t}$ should be fairly
large.

Figure 5: Colossal BCE estimated for bulk LCBH with MD simulations. (a) Volume
change per formula unit across the phase transition expressed as a function of
temperature and pressure. (b) Total entropy curves expressed as a function of
pressure and temperature. Inset: theoretically calculated $p$–$T$ phase
diagram. (c) Isothermal entropy and (d) adiabatic temperature changes
expressed as a function of temperature and pressure. Results were obtained
from $NpT$-MD simulations. Figure 6: Atomistic insights into the order-
disorder phase transition in LCBH from MD simulations. (a) Lithium ion
diffusion coefficient, $D_{\rm Li}$. (b) Anionic reorientational frequency,
$\lambda_{\rm CBH}$. Solid lines correspond to Arrhenius law fits. (c)–(d)
Cumulative function of the vibrational entropy as a function of the phonon
energy and atomic species, calculated for the ordered ($T=400$ K) and
disordered ($T=412$ K) phases at zero pressure. Dashed lines indicate
analogous asymptotic values reached in the ordered phase. (e)–(f) Angular
probability density function estimated for the molecular (CB11H12)- anions
calculated in the ordered ($T=350$ K) and disordered ($T=550$ K) phases at
zero pressure, expressed as a function of the polar ($\theta$) and azimuthal
($\phi$) angles. Dark and bright areas represent low and high probability
regions, respectively.

###  II.2 Experimental barocaloric results

Conventional X-ray powder diffraction experiments performed at normal pressure
and under varying temperature confirmed the expected structures of the low-$T$
and high-$T$ phases (orthorhombic and cubic symmetry, respectively). Pattern
matching analysis of the obtained data yielded the temperature-dependent
volume of LCBH (see Fig. 2a), which shows a huge $\approx 13$% relative volume
increase at the endothermic transition corresponding to $\Delta V\approx
12\cdot 10^{-5}$ m3 kg-1.

High-pressure differential thermal analysis (HP-DTA) was carried out in the
pressure interval $0\leq p\leq 0.23$ GPa (Fig. 2b). At pressures below
$\approx 0.13$ GPa, a single peak in the heat flow was measured corresponding
to the aforementioned orthorhombic (ordered phase, II) $\leftrightarrow$ cubic
(disordered phase, I) first-order phase transition. At pressure above $\approx
0.13$ GPa, the HP-DTA signals exhibit two peaks thus indicating the appearance
of a new phase that we label here as III (high-pressure enantiotropy). To the
best of our knowledge, phase III has not been previously reported in the
literature and its specific crystalline structure remains unknown since we did
not resolve it. Interestingly, a broad peak was previously detected in
differential scanning calorimetry experiments Tang  _et al._ (2015) that hints
at the stabilization of phase III.

Transition temperatures were determined from the maximum of the HP-DTA peaks
(Fig. 2c) which allowed to estimate an upper threshold for the triple point at
$\approx$ (425 K, 0.13 GPa) given the width of the peaks obtained under the
chosen scanning rate. Considering only the data measured near atmospheric
pressure, the pressure dependence of the II$\rightarrow$I transition was
determined to be $\frac{dT}{dp}\approx 420$ K GPa-1, which slightly decreases
under increasing pressure due to the small convexity of the coexistence line.
For the II$\rightarrow$III and III$\rightarrow$I transitions, linear fits to
the obtained coexistence lines yielded $\frac{dT}{dp}\approx 135$ K GPa-1 and
$\frac{dT}{dp}\approx 310$ K GPa-1, respectively. Phase transition entropy
changes were calculated via integration of the $\frac{1}{T}\frac{dQ}{dT}$
function after baseline subtraction. As it was already expected, the $\Delta
S_{{\rm II}\rightarrow{\rm I}}$ values associated to the LCBH order-disorder
phase transition are noticeably large, namely, $\approx 208$ J K-1 kg-1 (Fig.
2d). By plugging the measured $\frac{dT}{dp}$ and $\Delta S_{{\rm
II}\rightarrow{\rm I}}$ values at atmospheric pressure in the CC equation we
obtain $\Delta V_{\text{CC}}\approx 9\cdot 10^{-5}$ m3 kg-1, which is in
reasonable agreement with the $\Delta V$ determined directly in the
experiments.

Above $p\approx 0.13$ GPa, due to the overlapping between the
II$\leftrightarrow$III and III$\leftrightarrow$I peaks, the contribution
associated to each phase transition was decided at the inflection point of the
cumulative entropy change function
$\int_{T_{1}}^{T}\frac{1}{T^{\prime}}\frac{dQ}{dT^{\prime}}dT^{\prime}$.
$\Delta S_{\text{t}}$ remains practically constant from atmospheric pressure
all the way up to the triple point. At $p\simeq 0.13$ GPa, we obtained $\Delta
S_{\text{II}\rightarrow\text{I}}\approx\Delta
S_{\text{II}\rightarrow\text{III}}+\Delta S_{\text{III}\rightarrow\text{I}}$,
as it is required by the condition of thermodynamic equilibrium, while above
the triple point $\Delta S_{\text{II}\rightarrow\text{III}}\approx\Delta
S_{\text{III}\rightarrow\text{I}}$. Splitting of the II$\rightarrow$I phase
transition into II$\rightarrow$III and III$\rightarrow$I might be associated
to the decoupling of the diffusive and orientational degrees of freedom right
at the stabilization of the high-$T$ phase, although further investigations
are necessary for a more conclusive assessment of phase III.

HP-DTA measurements along with experimental differential scanning calorimetry
(Supplementary Fig. S1), heat capacity (Supplementary Fig. S2) and theoretical
equations of state $V(T,p)$ (i.e., obtained from molecular dynamics
simulations, Sec. II.3) were used to determine the isobaric entropy curves
$S(T,p)$ (Supplementary Fig. S3), from which the BC effects can be directly
calculated (Methods). Figures 3a,b show representative isothermal entropy
changes, $|\Delta S|$, and adiabatic temperature changes, $|\Delta T|$,
obtained upon the first application and removal of the driving pressure shift.
It is worth noticing that a small $\Delta p\approx 0.03$ GPa already produced
colossal values of $|\Delta S|=100$ J K-1 kg-1 and $|\Delta T|=8$ K, and
similarly $\Delta p\approx 0.08$ GPa yielded $|\Delta S|=250$ J K-1 kg-1 and
$|\Delta T|=16$ K. For the largest pressure shift considered in this study,
namely, $\Delta p\approx 0.23$ GPa, the resulting $|\Delta S|$ and $|\Delta
T|$ amount to $300$ J K-1 kg-1 and $40$ K, respectively.

Operation of solid-state cooling and heating devices requires cyclic
application and removal of the driving fields, for which reversible caloric
effects, $|\Delta S_{\text{rev}}|$ and $|\Delta T_{\text{rev}}|$, must be
considered. By reversible caloric effects we mean acquitted of phase
transition hysteresis effects Aznar  _et al._ (2020). The obtained results are
shown in Figs. 3c,d. Colossal $|\Delta S_{\text{rev}}|$ were already obtained
for a minimum pressure shift of $\approx 0.08$ GPa. For instance, under a
moderate pressure change of $\approx 0.10$ GPa LCBH renders $|\Delta
S_{\text{rev}}|=200$ J K-1 kg-1 and $|\Delta T_{\text{rev}}|=10$ K. Meanwhile,
for the largest pressure shift considered in this study we measured
outstanding values of $|\Delta S_{\text{rev}}|=280$ J K-1 kg-1 and $|\Delta
T_{\text{rev}}|=32$ K.

Figure 4 compares most of the experimental $|\Delta S_{\text{rev}}|$ and
$|\Delta T_{\text{rev}}|$ reported thus far in the literature for barocaloric
materials. Additionally, the size of the symbols therein account for the
materials BC strength, which is defined as the ratio of $|\Delta
S_{\text{rev}}|$ by the corresponding pressure shift ${\Delta p}$. The best
performing barocaloric materials, therefore, should appear in the top right
side of the panel and with the largest possible symbol area. Each material has
been represented with one or two points that best illustrate their overall
barocaloric performance, while for LCBH we have selected a set of barocaloric
measurements.

Although LCBH is not the best performing material in terms of a single
quality, it displays an unprecedentedly well-balanced and accomplished
barocaloric portfolio consisting of colossal $|\Delta S_{\text{rev}}|$, large
$|\Delta T_{\text{rev}}|$ and large BC strength obtained under moderate
pressure shifts of the order of $0.10$ GPa. For instance, in terms of largest
$|\Delta S_{\text{rev}}|$ the plastic crystal neopentylglycol (NPG) emerges as
the clear winner since it holds a gigantic value of $\approx 400$ J K-1 kg-1
Aznar  _et al._ (2020); however, as regards $|\Delta T_{\text{rev}}|$ the same
material becomes a poor contestant in the presence of LCBH (that is, $\approx
8$ K versus $32$ K). Likewise, the $|\Delta T_{\text{rev}}|$ record holder,
namely, the spin-crossover complex Fe3(bntrz)6(tcnset)6 Romanini  _et al._
(2021), presents $|\Delta S_{\text{rev}}|$ and BC strength values that roughly
are halves of the LCBH maxima (for instance, $\approx 120$ J K-1 kg-1 versus
$280$ J K-1 kg-1). Therefore, LCBH can be deemed as one of the most thorough
and promising barocaloric materials reported to date owing to its unique
parity between sizable $|\Delta S_{\text{rev}}|$ and $|\Delta T_{\text{rev}}|$
obtained under moderate pressure shifts.

Figure 7: Partial contributions to the entropy change accompanying the order-
disorder phase transition in LCBH expressed as a function of pressure. Entropy
changes stem from the vibrational, $\Delta S_{\rm vib}$, molecular
orientational, $\Delta S_{\rm ori}$, and cation diffusive, $\Delta S_{\rm
diff}$, degrees of freedom. Results were obtained from comprehensive molecular
dynamics simulations and Gibbs free energy calculations (Methods).

###  II.3 Atomistic simulation of barocaloric effects

Figures 5a,b show the theoretical equation of state $V(T,p)$ and $p$–$T$ phase
diagram of bulk LCBH obtained from molecular dynamics (MD) simulations
(Methods). We determined the coexistence line of the high-$T$ (disordered) and
low-$T$ (ordered) phases by conducting numerous MD simulations at small
$p$–$T$ shifts of $0.025$ GPa and $12.5$ K. Each phase coexistence point in
Fig. 5b (inset) corresponds to sharp and simultaneous changes in the volume,
Li+ diffusion coefficient ($D_{\rm Li}$), and molecular [CB11H12]-
orientational frequency ($\lambda_{\rm CBH}$), as identified in the MD
simulations (Figs. 6a,b). At zero pressure, we estimated a huge volume
increase of about $11$% at the theoretical transition temperature
$T_{t}\approx 400$ K, along with the order parameter changes $\Delta D_{\rm
Li}=1.13\cdot 10^{-6}$ cm2 s-1 and $\Delta\lambda_{\rm CBH}=0.33\cdot 10^{11}$
s-1. It was found that the pressure dependence of the transition temperature
could be precisely reproduced by the second-order polynomial curve
$T_{t}(p)=412+438p-610p^{2}$ (red line in the inset of Fig. 5b), in which the
temperature and pressure are expressed in units of K and GPa, respectively.
The slight $\frac{dT}{dp}$ decrease under increasing compression is consistent
with the $p$-induced reduction of the transition volume change since $\Delta
S_{t}$ is roughly independent of pressure, in agreement with our experiments.
It is worth noting that phase-transition hysteresis effects cannot be
reproduced by the equilibrium MD approach employed in this study Sau  _et al._
(2021a).

The LCBH $p$–$T$ phase diagram obtained from MD simulations (Fig. 5b) is in
quantitative good agreement with the experiments performed below the triple
point found at $\approx 0.13$ GPa (Fig. 2b), although the transition
temperatures are slightly overestimated by theory. For example, at zero
pressure and $p=0.10$ GPa the MD simulations yielded $T_{t}=410\pm 15$ and
$440\pm 15$ K (Fig. 5), respectively, to be compared with the corresponding
experimental values $390\pm 10$ and $410\pm 10$ K (Fig. 2b). The agreement
between the predicted and measured volumes for the ordered and disordered
phases at zero pressure is also notable, finding only small relative
discrepancies of $\sim 1$% for the low-$T$ phase (Figs. 2a and 5a). Meanwhile,
the triple point observed in the experiments was not reproduced by the MD
simulations. It is worth noting, however, that under $p\neq 0$ conditions and
close to $T_{t}$ we observed pre-transitional effects in our simulations
consisting of few slowly diffusing Li ions in the ordered phase (Supplementary
Fig. S4).

Figures 5c,d show the theoretical barocaloric $|\Delta S|$ and $|\Delta T|$
deduced from the entropy curves $S(p,T)$ enclosed in Fig. 5b, which were
obtained from data generated in the MD simulations. The agreement between
these theoretical results and the corresponding experimental values is
remarkably good for pressures below the experimental triple point. For
example, for a pressure shift of $0.10$ GPa we estimated an isothermal entropy
change of $227$ J K-1 kg-1 and an adiabatic temperature change of $32$ K from
the MD simulations, to be compared with the corresponding experimental values
$250$ J K-1 kg-1 and $24$ K (Fig. 3a,b). In view of such a notable agreement,
we characterized with MD simulations the contributions to the phase transition
entropy change stemming from the vibrational, molecular orientational and
cation diffusive degrees of freedom, a highly valuable atomistic insight that
in principle cannot be obtained from the experiments.

Figures 6a,b reveal synchronized surges in $D_{\rm Li}$ and $\lambda_{\rm
CBH}$ at the order-disorder phase transition points. Thus, both ion diffusion
and molecular anion orientational disorder (Figs. 6e,f) contribute to the
transition entropy change and barocaloric effects disclosed in LCBH.
Nevertheless, there is a third possible source of entropy in the crystal which
is related to the lattice vibrations, $S_{\rm vib}$ (Supplementary Figs.
S5-S6). Figures 6c,d show examples of the cumulative $S_{\rm vib}$ function
expressed as a function of the vibrational phonon energy, calculated for LCBH
in the ordered and disordered phases at zero pressure and evaluated for each
atomic species. Therein, it is appreciated that the largest contribution to
the $S_{\rm vib}$ difference between the order and disordered phases comes
from the B atoms (followed by hydrogen). This outcome can be rationalized in
terms of the relative great abundance of this species in LCBH ($\approx 45$%)
and its larger mass as compared to that of H atoms ($10$ times heavier): B
ions have a predominant weight on the low-frequency vibrational modes (Fig.
6c-d) that most significantly contribute to $S_{\rm vib}$ near ambient
temperature.

Figure 7 shows the relative contributions of the vibrational, molecular
orientational and ion diffusion degrees of freedom to the phase transition
entropy change estimated at different pressures with MD simulations.
Interestingly, in all the analyzed cases the largest contribution stems from
changes in the lattice vibrations, $\Delta S_{\rm vib}$, followed by the
molecular reorientations, $\Delta S_{\rm ori}$, and finally ion diffusion,
$\Delta S_{\rm diff}$. For example, at zero pressure the vibrational,
molecular orientational and ion diffusive degrees of freedom respectively
contribute in $\approx 48$, $32$ and $20$% to $\Delta S_{t}$. The entropy
preeminence of the lattice vibrations can be rationalized in terms of (1) the
huge volume expansion accompanying the order-disorder phase transition ($\sim
10$%, Fig. 5a), which further curtails the frequency of the low-energy phonon
bands in the disordered phase (Supplementary Fig. S5), and (2) the
intensification and amplitude broadening of the molecular libration modes in
the disordered phase (inferred from the angular probability density variations
around the equilibrium positions in Figs. 6e-f). These outcomes are highly
valuable and insightful since thus far molecular reorientations were thought
to be the primary source of entropy variation in plastic crystals undergoing
order-disorder phase transitions Li  _et al._ (2020a, 2022); de Oliveira
(2023).

The vibrational and orientational entropy changes remain more or less constant
for pressures $\leq 0.1$ GPa, whereas $\Delta S_{\rm diff}$ significantly
decreases under compression. For instance, at $0.1$ GPa the diffusive degrees
of freedom contribute to $\Delta S_{t}$ in less than $4$%. These outcomes can
be understood in terms of the small fraction of diffusive ions in LCBH (i.e.,
one Li atom per formula unit) and the marked decline in $D_{\rm Li}$ induced
by pressure (Fig. 6a). The appearance of pre-transitional effects in our MD
simulations, specially under $p\neq 0$ conditions (Supplementary Fig. S4),
also contributes to the noticeable $\Delta S_{\rm diff}$ drop caused by
compression. Nonetheless, it is worth noting that despite the relative
minuteness of $\Delta S_{\rm diff}$, cation disorder was found to play a
critical role on triggering molecular orientational disorder, which by
contrast contributes very significantly to $\Delta S_{t}$. In particular, we
conducted constrained MD runs in which we fixed the positions of the lithium
ions so that they could not diffuse. It was found then that molecular
orientational disorder only emerged at temperatures well above $550$ K
(Supplementary Fig. S7). Therefore, it can be concluded that cation disorder
crucially assists on the realization of colossal BCE through the order-
disorder phase transition, a characteristic trait that differentiates LCBH
from other molecular plastic crystals bearing also great barocaloric promise.

##  III Conclusions

Colossal barocaloric effects (BCE) driven by pressure shifts of the order of
$0.10$ GPa were experimentally and theoretically disclosed in bulk LiCB11H12
(LCBH), a compound that at high temperatures presents disorder features
characteristic of both plastic crystals and superionic materials, namely,
molecular reorientational motion and ion diffusion. Reversible peaks of
$|\Delta S_{\rm rev}|=280$ J K-1 kg-1 and $|\Delta T_{\rm rev}|=32$ K were
experimentally measured around $400$ K for a pressure shift of $0.23$ GPa,
yielding huge and reversible barocaloric strengths of $\approx 2$ J K-1 kg-1
MPa-1 over tens of degrees intervals. Likewise, for a smaller pressure shift
of $0.10$ GPa we obtained very promising values of $|\Delta S_{\rm rev}|=200$
J K-1 kg-1 and $|\Delta T_{\rm rev}|=10$ K. These results place LCBH among the
best-known barocaloric materials in terms of huge and reversible isothermal
entropy and adiabatic temperature changes, two quantities that rarely are
found simultaneously in a same material.

Atomistic molecular dynamics simulations yielded theoretical $|\Delta S|$ and
$|\Delta T|$ in very good agreement with the experimental values, and allowed
to quantify the importance of vibrational, molecular orientational, and ion
diffusive degrees of freedom on the disclosed colossal BCE. It was found that
the contribution to the phase transition entropy change stemming from the
lattice vibrations was the largest, followed by that of molecular
reorientations and both being much superior than the entropy associated to
lithium diffusion alone. Nevertheless, cationic disorder was found to have a
critical influence on the stabilization of orientational disorder thus, in
spite of its small contribution to $\Delta S_{t}$, lithium diffusion appears
to be essential for the emergence of colossal BCE in bulk LCBH. These results
are of high significance since reveal the preeminence of the vibrational
degrees of freedom in the phase transition entropy change of a plastic
crystal, and demonstrate atomistic BCE mechanisms other than molecular
reorientational disorder (i.e., lattice vibrations and ion diffusion).

LCBH belongs to the family of closo-borate materials, a promising class of
solid electrolytes for all-solid-state batteries. Examples of akin compounds
that have been already synthesized in the laboratory and tested for
electrochemical energy storage applications are NaCB11H12 Tang  _et al._
(2015); Skripov  _et al._ (2015), KCB11H12 Dimitrievska  _et al._ (2020), and
LiCB9H10 Kim  _et al._ (2019, 2020). Colossal BCE could also exist in these
materials and in other similar compounds harboring both ion diffusion and
molecular orientational disorder at or near room temperature. Thus, the
present combined experimental-theoretical study opens new horizons in solid-
state cooling and heating and advances knowledge in the realization of
colossal BCE in plastic crystals.

## Methods

### Experimental techniques

Materials synthesis. LiCB11H12 was obtained by drying the hydrated compound
LiCB11H12$\cdot$xH2O (Katchem, Ltd.) under vacuum ($<5\times 10^{-4}$ Pa) at
160 ∘C for 12 h.
X-ray powder diffraction. High-resolution X-ray powder diffraction
measurements were performed using the Debye–Scherrer geometry and transmission
mode with a horizontally mounted cylindrical position-sensitive INEL detector
(CPS-120). Monochromatic Cu-K$\alpha_{1}$ radiation was selected by means of a
curved germanium monochromator. Temperature-dependent measurements were
performed using a liquid nitrogen 700 series Oxford Cryostream Cooler. Powder
samples were introduced into 0.5 mm diameter Lindemann capillaries. Volume was
obtained by pattern matching procedure.
Quasi-direct barocaloric measurements. A Q100 thermal analyzer (TA
Instruments) was used to perform differential scanning calorimetry experiments
at atmospheric pressure with $\sim 10$ mg of sample hermetically encapsulated
in Aluminum pans (Supplementary Fig. S1). The standard mode (at 3, 5 and 10 K
min-1) was used to determine the transition properties whereas the modulated
mode (isothermal conditions, modulation amplitude 1 ∘C, modulation period
$120$ s) was used to measure the heat capacity in each phase (Supplementary
Fig. S2).

Pressure-dependent calorimetry was performed with a custom-built high-pressure
differential thermal analyzer (from Irimo, Bellota Herramientas S.A.) that
uses Bridgman thermocouples as thermal sensors. The nominal operational
pressure range is from atmospheric to $0.3$ GPa and the temperature range is
from room temperature up to $473$ K. Heating ramps were performed at $3$ K
min-1 using a resistive heater whereas cooling were carried out at $\sim-2$ K
min-1 by an air stream. A few hundreds of mg of LiCB11H12 were mixed with an
inert perfluorinated fluid (Galden Bioblock Scientist) to remove air and
sealed within tin capsules. The pressure-transmitting fluid was Therm240
(Lauda).

Isobaric entropy functions $S(T,p)$ were determined with respect to a
reference temperature $T_{0}$ below the transition using the method explained
in Ref.Li  _et al._ (2020b) (Supplementary Fig. S3). The procedure is based on
the following thermodynamic equation:

$S(T,p)=S(T_{0},p)+\int_{T_{0}}^{T}\frac{1}{T}\left(C_{p}+\frac{dQ}{dT}\right)dT~{},$
(1)

where $\frac{dQ}{dT}$ is the heat flow in temperature due to the first-order
phase transition measured by pressure-dependent calorimetry.

In each phase, $C_{p}$ is the corresponding heat capacity and was considered
independent of pressure as indicated by the approximately linear behavior of
volume with temperature obtained in the two phases from MD simulations (Fig.
5) along with the thermodynamic equation:

$\left(\frac{\partial C_{p}}{\partial
p}\right)_{T}=-T\left(\frac{\partial^{2}V}{\partial T^{2}}\right)_{p}~{}.$ (2)

In the transition region $C_{p}$ was calculated as an average weighted
according to the fraction of each phase. To take into account the dependence
of the transition region with pressure, the overall $C_{p}$ function at
atmospheric pressure obtained in each phase and across the transition was
extrapolated to higher temperatures according to the experimental value of
$\frac{dT}{dp}\Delta p$, where $\Delta p$ is the pressure change applied in
each particular case. Experimental measurement of $C_{p}$ at atmospheric
pressure and the calculated curves at different pressures are shown in
Supplementary Fig. S2.

The pressure dependence of $S(T,p)$ was evaluated using the thermodynamic
equation:

$S(T,p)=S(T,p_{0})-\int_{p_{0}}^{p}\left(\frac{\partial V}{\partial
T}\right)_{T,p^{\prime}}dp^{\prime}~{},$ (3)

where $p_{0}$ was selected equal to $p_{\text{atm}}=1$ bar. Here, we make use
of the approximation $\left(\frac{\partial V}{\partial
T}\right)_{T,p}\simeq\left(\frac{\partial V}{\partial T}\right)_{T,p_{0}}$,
which is reasonable based on the $\left(\frac{\partial V}{\partial
T}\right)_{T,p}$ data obtained from the MD simulations (Fig. 5).

Once the entropy function $S(T,p)$ was determined for both heating and cooling
runs independently (Supplementary Fig. S3), BC effects obtained upon first
application or removal of the field were calculated as:

$\displaystyle\Delta S(T,p_{0}\rightarrow p_{1})=S(T,p_{1})-S(T,p_{0})~{}{\rm
and}$ (4) $\displaystyle\Delta T(T_{s},p_{0}\rightarrow
p_{1})=T(S,p_{1})-T_{s}(T,p_{0})~{},$ (5)

where $T_{s}$ is the starting temperature of the heating/cooling process.
Here, it must be considered that for materials with $\frac{dT}{dp}>0$ BC
effects on compression ($p_{0}=p_{\text{atm}}$, $p_{1}>p_{\text{atm}}$) and
decompression ($p_{0}>p_{\text{atm}}$, $p_{1}=p_{\text{atm}}$) are calculated
from $S(T,p)$ functions obtained on cooling and heating, respectively Aznar
_et al._ (2020). In turn, BC effects obtained reversibly on cyclic
compression-decompression processes were calculated from the $S(T,p)$ curves
obtained on heating at atmospheric pressure and cooling at high pressure.

### Simulation techniques

Molecular dynamics simulations. Force-field based molecular dynamics (MD)
simulations were performed using a previously reported interatomic potential
for LCBH Sau  _et al._ (2021b). This force field is a combination of Coulomb-
Buckingham (CB), harmonic bond, and angle-type potentials, namely:

$\displaystyle
U(r,\theta)=U_{\rm{CB}}(r)+U_{\rm{bond}}(r)+U_{\rm{angle}}(\theta)~{},$ (6a)
$\displaystyle
U_{\rm{CB}}(r)=\frac{q_{i}q_{j}}{4\pi\epsilon_{0}r}+A_{ij}\exp(-r/\rho)-\frac{C_{ij}}{r^{6}}~{},$
(6b) $\displaystyle U_{\rm{bond}}(r)=\frac{1}{2}k_{r}(r-r_{0})^{2}~{}{\rm
and}$ (6c) $\displaystyle
U_{\rm{angle}}(\theta)=\frac{1}{2}k_{\theta}(\theta-\theta_{0})^{2}~{},$ (6d)

where $q_{i}$ denotes the charge of the ion labeled $i$, $\epsilon_{0}$ the
vacuum permittivity, $A_{ij}$ and $\rho$ the short-range repulsive energy and
length scales for the pairs of atoms $ij$, and $C_{ij}$ the corresponding
dispersion interaction coefficient. $r_{0}$ and $\theta_{0}$ are an
equilibrium bond distance and angle, respectively, and $k_{r}$ and
$k_{\theta}$ the spring constants of the harmonic bond and angle potentials.
The numerical value of these potential parameters can be found in the
Supplementary Table S2.

We performed $NpT$-MD simulations in the temperature range $325\leq T\leq 525$
K at intervals of $12.5$ K, and pressure range $0\leq p\leq 0.15$ GPa at
intervals of $0.025$ GPa. The temperature and pressure in the system were
controlled with thermostating and barostating techniques, in which some
dynamic variables are coupled with the particle velocities and simulation box
dimensions. The simulation supercell comprised a total of $6400$ atoms. A time
step of $0.5$ fs was employed for integration of the atomic forces along with
the velocity Verlet algorithm. A typical $NpT$-MD run lasted for about $2$ ns
and the atomic trajectories were stored at intervals of $500$ fs. Detailed
analyses and statistical time averages were performed over the last $1$ ns of
such simulations. To guarantee proper convergence of the estimated
thermodynamic properties, in few instances longer simulation times of $10$ ns
were carried out. Periodic boundary conditions were applied along the three
Cartesian directions and the Ewald summation technique was used for evaluation
of the long-range Coulomb interactions with a short-range cut-off distance of
$13$ Å. All the $NpT$-MD simulations were carried out with the LAMMPS software
package Plimpton (1995).

Density functional theory and ab initio molecular dynamics simulations. First-
principles calculations based on density functional theory (DFT) were
performed to analyze the energy, structural and vibrational properties of bulk
LCBH. The DFT calculations were carried out with the VASP code Kresse and
Hafner (1993) by following the generalized gradient approximation to the
exchange-correlation energy due to Perdew et al. (PBE) Perdew  _et al._
(1996). The projector augmented-wave method was used to represent the ionic
cores Blöchl (1994), and the electronic states $1s$–$2s$ Li, $2s$–$2p$ C,
$2s$–$2p$ B and $1s$ H were considered as valence. Wave functions were
represented in a plane-wave basis 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.005$ eV Å-1 was imposed in 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
MD simulations on the description of the vibrational degrees of freedom of
bulk LCBH (Supplementary Fig. S6). The AIMD simulations were performed in the
canonical ensemble $(N,V,T)$ considering constant number of particles, volume
and temperature. The constrained volumes were equal to the equilibrium volumes
determined at zero temperature, an approximation that has been shown to be
reasonable at moderate temperatures Sagotra  _et al._ (2019a). The temperature
in the AIMD simulations was kept fluctuating around a set-point value by using
Nose-Hoover thermostats. A large simulation box containing 800 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
$\delta t=10^{-3}$ ps. $\Gamma$-point sampling for integration within the
first Brillouin zone was employed in all the AIMD simulations. The AIMD
simulations comprised long simulation times of $\approx 200$ ps and
temperatures in the range $200\leq T\leq 500$ K.

Estimation of key quantities with MD simulations. The mean square displacement
of the lithium ions was estimated with the formula Sagotra  _et al._ (2019b):

$\displaystyle{\rm MSD_{\rm Li}}(\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{N_{\rm ion}\left(N_{\rm step}-n_{\tau}\right)}\times$
$\displaystyle\sum_{i=1}^{N_{\rm ion}}\sum_{j=1}^{N_{\rm step}-n_{\tau}}|{\bf
r}_{i}(t_{j}+\tau)-{\bf r}_{i}(t_{j})|^{2}~{},$

where ${\bf r}_{i}(t_{j})$ is the position of the migrating ion $i$ at time
$t_{j}$ ($=j\cdot\delta t$), $\tau$ represents a lag time,
$n_{\tau}=\tau/\delta t$, $N_{\rm ion}$ is the total number of mobile ions,
and $N_{\rm step}$ the total number of time steps. The maximum $n_{\tau}$ was
chosen equal to $N_{\rm step}/2$, hence we could accumulate enough statistics
to reduce significantly the fluctuations in ${\rm MSD_{\rm Li}}(\tau)$ at
large $\tau$’s. The diffusion coefficient of lithium ions was calculated with
the Einstein’s relation:

$D_{\rm Li}=\lim_{\tau\to\infty}\frac{{\rm MSD_{Li}}(\tau)}{6\tau}~{},$ (8)

by performing linear fits to the averaged ${\rm MSD_{Li}}$ values calculated
at long $\tau$.

The angular autocorrelation function of the molecular [CB11H12]- anions was
estimated using the expression Sau  _et al._ (2021a):

$\phi_{\rm CBH}(\tau)=\langle\hat{{\bf r}}(t)\cdot\hat{{\bf
r}}(t+\tau)\rangle~{},$ (9)

where $\hat{{\bf r}}$ is a unitary vector connecting the center of mass of
each closoborane unit with one of its edges and $\langle\cdots\rangle$ denotes
statistical average in the $(N,p,T)$ ensemble considering all the molecular
anions. This autocorrelation function typically decays as
$\propto\exp{[-\lambda_{\rm CBH}\cdot\tau]}$, where the parameter
$\lambda_{\rm CBH}$ represents a characteristic reorientational frequency. For
significant anion reorientational motion, that is, large $\lambda_{\rm CBH}$,
the $\phi_{\rm CBH}$ function decreases rapidly to zero with time.

The temperature dependence of the lithium diffusion coefficient was assumed to
follow an Arrhenius law at any pressure of the form:

$D_{\rm Li}(T)=D_{0}\cdot e^{-(\frac{E_{a}}{k_{B}T})}~{},$ (10)

where $D_{0}$ and $E_{a}$ are parameters that depend on $p$ and $k_{B}$
represents the Boltzmann constant. The reorientational frequency of
closoborane units, $\lambda_{\rm CBH}$, was assumed to follow a similar
dependence on temperature.

The entropy of each phase was calculated as a function of temperature and
pressure, $S(p,T)$, by fully considering the vibrational, molecular
orientational and ion diffusive degrees of freedom:

$S(p,T)=S_{\rm vib}(p,T)+S_{\rm ori}(p,T)+S_{\rm diff}(p,T)~{}.$ (11)

In the low-$T$ phase, $S_{\rm ori}$ and $S_{\rm diff}$ are null while in the
high-$T$ phase are finite and positive.

The vibrational density of states (VDOS), $g(\omega)$, was calculated via the
Fourier transform of the velocity-velocity autocorrelation function obtained
directly from the $NpT$-MD simulations, namely:

$g(\omega)=\frac{1}{N_{ion}}\sum_{i}^{N_{ion}}\int_{0}^{\infty}\langle{\bf
v}_{i}(\tau)\cdot{\bf v}_{i}(0)\rangle e^{i\omega\tau}d\tau~{},$ (12)

where ${\bf v}_{i}(t)$ represents the velocity of the atom labeled $i$ at time
$t$, and $\langle\cdots\rangle$ denotes statistical average in the $(N,p,T)$
ensemble. The vibrational entropy was subsequently estimated with the formula
Togo  _et al._ (2010):

$\displaystyle S_{\rm vib}(p,T)$ $\displaystyle=$
$\displaystyle-\int_{0}^{\infty}k_{B}\ln{\left[2\sinh{\left(\frac{\hbar\omega}{2k_{B}T}\right)}\right]}\hat{g}(\omega)d\omega+$
(13)
$\displaystyle\int_{0}^{\infty}\frac{\hbar\omega}{2T}\tanh^{-1}{\left(\frac{\hbar\omega}{2k_{B}T}\right)}\hat{g}(\omega)d\omega~{},$

where $\hat{g}(\omega)$ is the normalized vibrational density of states
($\int_{0}^{\infty}\hat{g}(\omega)d\omega=3N_{ion}$) and the dependence on
pressure (and also temperature) is implicitly contained in $\hat{g}(\omega)$.

The orientational entropy of the molecular anions, $S_{\rm ori}$, was directly
calculated from the angular probability density, $\rho(\theta,\phi)$, like
Takagi  _et al._ (2020):

$S_{\rm
ori}(p,T)=-k_{B}\int_{0}^{\pi}\int_{0}^{2\pi}\rho(\theta,\phi)\ln{\rho(\theta,\phi)}~{}d\theta
d\phi~{},$ (14)

where $\rho(\theta,\phi)$ was obtained from the $NpT$-MD simulation runs in
the form of average histograms (Fig. 6).

The ion diffusive entropy difference was estimated at the phase transition
points via equalization of the Gibbs free energies of the low-$T$ (O) and
high-$T$ (D) phases, namely, $G^{D}(p,T_{t})=G^{O}(p,T_{t})$, thus leading to
the expression:

$\Delta S_{\rm diff}(p,T_{t})=\frac{\langle\Delta
E\rangle}{T_{t}}+p\frac{\langle\Delta V\rangle}{T_{t}}-\Delta S_{\rm
vib}-\Delta S_{\rm ori}~{},$ (15)

where $\Delta X\equiv X^{D}-X^{O}$ and $E$ represents the internal energy of
the system. For any pressure, $\Delta S_{\rm diff}$ was assumed to be constant
at temperatures $T_{t}\leq T$.

## References

  * Moya and Mathur (2020) X. Moya and N. D. Mathur, Science 370, 797 (2020).
  * Mañosa  _et al._ (2013) L. Mañosa, A. Planes, and M. Acet, Journal of Materials Chemistry A 1, 4925 (2013).
  * Hou  _et al._ (2022) H. Hou, S. Qian, and I. Takeuchi, Nature Reviews Materials 7, 633 (2022).
  * Mañosa and Planes (2017) L. Mañosa and A. Planes, Advanced Materials 29, 1603607 (2017).
  * Cazorla (2019) C. Cazorla, Applied Physics Reviews 6, 041316 (2019).
  * Lloveras and Tamarit (2021) P. Lloveras and J.-L. Tamarit, MRS Energy & Sustainability 8, 3 (2021).
  * Lloveras  _et al._ (2019) P. Lloveras, A. Aznar, M. Barrio, P. Negrier, C. Popescu, A. Planes, L. Mañosa, E. Stern-Taulats, A. Avramenko, N. D. Mathur, X. Moya, and J.-L. Tamarit, Nature Communications 10, 1803 (2019).
  * Li  _et al._ (2019) B. Li, Y. Kawakita, S. Ohira-Kawamura, T. Sugahara, H. Wang, J. Wang, Y. Chen, S. I. Kawaguchi, S. Kawaguchi, K. Ohara, K. Li, D. Yu, R. Mole, T. Hattori, T. Kikuchi, S.-i. Yano, Z. Zhang, Z. Zhang, W. Ren, S. Lin, O. Sakata, K. Nakajima, and Z. Zhang, Nature 567, 506 (2019).
  * Aznar  _et al._ (2020) A. Aznar, P. Lloveras, M. Barrio, P. Negrier, A. Planes, L. Mañosa, N. D. Mathur, X. Moya, and J.-L. Tamarit, Journal of Materials Chemistry A 8, 639 (2020).
  * Aznar  _et al._ (2021) A. Aznar, P. Negrier, A. Planes, L. Mañosa, E. Stern-Taulats, X. Moya, M. Barrio, J.-L. Tamarit, and P. Lloveras, Applied Materials Today 23, 101023 (2021).
  * Zhang  _et al._ (2022) K. Zhang, R. Song, J. Qi, Z. Zhang, Z. Zhang, C. Yu, K. Li, Z. Zhang, and B. Li, Advanced Functional Materials 32, 2112622 (2022).
  * Imamura  _et al._ (2020) W. Imamura, É. O. Usuda, L. S. Paixão, N. M. Bom, A. M. Gomes, and A. M. G. Carvalho, Chinese Journal of Polymer Science 38, 999 (2020).
  * Li  _et al._ (2021) J. Li, M. Barrio, D. J. Dunstan, R. Dixey, X. Lou, J.-L. Tamarit, A. E. Phillips, and P. Lloveras, Advanced Functional Materials 31, 2105154 (2021).
  * Seo  _et al._ (2022) J. Seo, R. D. McGillicuddy, A. H. Slavney, S. Zhang, R. Ukani, A. A. Yakovenko, S.-L. Zheng, and J. A. Mason, Nature Communications 13, 2536 (2022).
  * Salvatori  _et al._ (2022) A. Salvatori, P. Negrier, A. Aznar, M. Barrio, J. L. Tamarit, and P. Lloveras, APL Materials 10, 111117 (2022).
  * Aznar  _et al._ (2017) A. Aznar, P. Lloveras, M. Romanini, M. Barrio, J.-L. Tamarit, C. Cazorla, D. Errandonea, N. D. Mathur, A. Planes, X. Moya, and L. Mañosa, Nature Communications 8, 1851 (2017).
  * Sagotra  _et al._ (2017) A. K. Sagotra, D. Errandonea, and C. Cazorla, Nature Communications 8, 963 (2017).
  * Sagotra  _et al._ (2018) A. K. Sagotra, D. Chu, and C. Cazorla, Nature Communications 9, 3337 (2018).
  * Min  _et al._ (2020) J. Min, A. K. Sagotra, and C. Cazorla, Phys. Rev. Materials 4, 015403 (2020).
  * Tang  _et al._ (2015) W. Tang, A. Unemoto, W. Zhou, V. Stavila, M. Matsuo, H. Wu, S.-I. Orimo, and T. Udovic, Energy and Environmental Science 8, 3637 (2015).
  * Romanini  _et al._ (2021) M. Romanini, Y. Wang, K. Gürpinar, G. Ornelas, P. Lloveras, Y. Zhang, W. Zheng, M. Barrio, A. Aznar, A. Gràcia-Condal, B. Emre, O. Atakol, C. Popescu, H. Zhang, Y. Long, L. Balicas, J. Lluís Tamarit, A. Planes, M. Shatruk, and L. Mañosa, Advanced Materials 33, 2008076 (2021).
  * Li  _et al._ (2020a) F. B. Li, M. Li, X. Xu, Z. C. Yang, H. Xu, C. K. Jia, K. Li, J. He, B. Li, and H. Wang, Nature Communications 11, 4190 (2020a).
  * Li  _et al._ (2022) F. Li, M. Li, C. Niu, and H. Wang, Applied Physics Letters 120, 073902 (2022).
  * de Oliveira (2023) N. de Oliveira, Acta Materialia 246, 118657 (2023).
  * Sau  _et al._ (2021a) K. Sau, T. Ikeshoji, S. Takagi, S.-i. Orimo, D. Errandonea, D. Chu, and C. Cazorla, Scientific Reports 11, 11915 (2021a).
  * Skripov  _et al._ (2015) A. V. Skripov, R. V. Skoryunov, A. V. Soloninin, O. A. Babanova, W. S. Tang, V. Stavila, and T. J. Udovic, The Journal of Physical Chemistry C 119, 26912 (2015).
  * Mohtadi and Orimo (2016) R. Mohtadi and S.-i. Orimo, Nature Reviews Materials 2, 16091 (2016).
  * Dimitrievska  _et al._ (2020) M. Dimitrievska, H. Wu, V. Stavila, O. A. Babanova, R. V. Skoryunov, A. V. Soloninin, W. Zhou, B. A. Trump, M. S. Andersson, A. V. Skripov, and T. J. Udovic, The Journal of Physical Chemistry C 124, 17992 (2020).
  * Kim  _et al._ (2019) S. Kim, H. Oguchi, N. Toyama, T. Sato, S. Takagi, T. Otomo, D. Arunkumar, N. Kuwata, J. Kawamura, and S.-i. Orimo, Nature Communications 10, 1081 (2019).
  * Kim  _et al._ (2020) S. Kim, K. Kisu, S. Takagi, H. Oguchi, and S.-I. Orimo, ACS Applied Energy Materials 3, 4831 (2020).
  * Sau  _et al._ (2021b) K. Sau, T. Ikeshoji, S. Kim, S. Takagi, and S. I. Orimo, Chemistry of Materials 33, 2357 (2021b).
  * Li  _et al._ (2020b) J. Li, D. Dunstan, X. Lou, A. Planes, L. Mañosa, M. Barrio, J.-L. Tamarit, and P. Lloveras, Journal of Materials Chemistry A 8, 20354 (2020b).
  * Plimpton (1995) S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Perdew  _et al._ (1996) J. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Blöchl (1994) P. Blöchl, Physical Review B 50, 17953 (1994).
  * Sagotra  _et al._ (2019a) A. K. Sagotra, D. Chu, and C. Cazorla, Phys. Rev. Mater. 3, 035405 (2019a).
  * Sagotra  _et al._ (2019b) A. K. Sagotra, D. Chu, and C. Cazorla, Physical Review Materials 3, 035405 (2019b).
  * Togo  _et al._ (2010) A. Togo, L. Chaput, I. Tanaka, and G. Hug, Physical Review B 81, 174301 (2010).
  * Takagi  _et al._ (2020) S. Takagi, T. Ikeshoji, T. Sato, and S.-i. Orimo, Applied Physics Letters 116, 173901 (2020).

## Acknowledgements

K.S. thanks the JSPS International Fellowship. This work was supported by JSPS
KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas
“Hydrogenomics”, No. JP18H05513 and JSPS Fellowship grant (21F21345). We
gratefully acknowledge the Center for Computational Materials Science of
Institute for Materials Research, Tohoku University for permitting the use of
MASAMUNE- IMR (MAterial science Supercomputing system for Advanced MUltiscale
simulations toward Next-generation Institute of Material Research) (project
no. 202203-SCKXX-0408). This work was also supported by MINECO Project No.
PID2020-112975GB-I00 (Spain) and DGU Project No. 2021SGR-00343 (Catalonia).
C.C. acknowledges financial support from the Spanish Ministry of Science,
Innovation and Universities under the “Ramón y Cajal” fellowship
RYC2018-024947-I and the TED2021-130265B-C22 project, and computational
support from the Red Española de Supercomputación (RES) under the grants
FI-2022-1-0006, FI-2022-2-0003 and FI-2022-3-0014. M.Z. (CSC No. 202008310208)
thanks the China Scholarship Council.

## Conflict of Interest

The authors declare no conflict of interest.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.


#  2D discrete Hodge-Dirac operator on the torus

Volodymyr Sushch
Koszalin University of Technology
Sniadeckich 2, 75-453 Koszalin, Poland
<EMAIL_ADDRESS>
###### Abstract

We discuss a discretisation of the de Rham-Hodge theory in the two-dimensional
case based on a discrete exterior calculus framework. We present discrete
analogues of the Hodge-Dirac and Laplace operators in which key geometric
aspects of the continuum counterpart are captured. We provide and prove a
discrete version of the Hodge decomposition theorem. Special attention has
been paid to discrete models on a combinatorial torus. In this particular
case, we also define and calculate the cohomology groups.

Key words and phrases: discrete exterior calculus, discrete operators, Hodge-
Dirac operator, discrete Laplacian, Hodge decomposition, combinatorial torus,
cohomology groups

Math. Subj. Clas.: 39A12, 39A70, 58A14

##  1 Introduction

The choice of technique to approximate the solution of partial differential
equations depends on a discretisation scheme. The interest for discrete models
which preserve a geometric structure of continuum counterparts has grown in
the computing community [1, 2, 3, 10, 12]. It is known that a geometric
discretisation scheme would be ideal if this had the same properties as the
continuum. However, there are difficulties usually in definitions of discrete
counterparts of the Hodge star and the wedge product between forms. Various
approaches with geometric discretisation have been proposed in the literature.
See, for example, [4, 5, 8, 9, 11, 19, 20, 22]. Recently, in [6], it was
described a quite general framework of discrete calculus based on a new type
of discrete geometry called script geometry. The proposed approach in this
article was introduced by Dezin [7] and later developed in the author’s
previous papers [13, 14, 15, 16, 17, 18].

Let
$\Lambda(\mathbb{R}^{2})=\Lambda^{0}(\mathbb{R}^{2})\oplus\Lambda^{1}(\mathbb{R}^{2})\oplus\Lambda^{2}(\mathbb{R}^{2})$
denote the graded vector space of smooth differential forms on
$\mathbb{R}^{2}$, where $\Lambda^{r}(\mathbb{R}^{2})$ denotes the subspace of
$r$-forms, $r=0,1,2$. Let
$d:\Lambda^{r}(\mathbb{R}^{2})\rightarrow\Lambda^{r+1}(\mathbb{R}^{2})$ be the
exterior derivative. The codifferential
$\delta:\Lambda^{r}(\mathbb{R}^{2})\rightarrow\Lambda^{r-1}(\mathbb{R}^{2})$
is defined by $\delta=-\ast d\ast$, where $\ast$ is the Hodge star operator
such that
$\ast:\Lambda^{r}(\mathbb{R}^{2})\rightarrow\Lambda^{2-r}(\mathbb{R}^{2})$ and
$\ast^{2}=(-1)^{r}$. The operator

$d+\delta:\Lambda(\mathbb{R}^{2})\rightarrow\Lambda(\mathbb{R}^{2})$ (1.1)

is called the Hodge-Dirac operator on $\mathbb{R}^{2}$. The Laplacian
$\Delta:\Lambda^{r}(\mathbb{R}^{2})\rightarrow\Lambda^{r}(\mathbb{R}^{2})$ is
defined by

$\Delta=(d+\delta)^{2}=d\delta+\delta d.$ (1.2)

Our purpose is to develop a satisfactory discrete model of de Rham-Hodge
theory on manifolds which are homeomorphic to the torus. We consider a chain
complex as a combinatorial model of $\mathbb{R}^{2}$. When we add to this
discrete analogues of the exterior derivative, the Hodge star operator, and
the exterior product acting on cochains, we have all of the basic ingredients
for the calculus of discrete counterparts of differential forms. We show that
discrete analogues of the operators (1.1) and (1.2) have properties like those
in the continual case. We formulate and prove a discrete version of the Hodge
decomposition theorem. We give an example illustrating how cohomology groups
are calculated for our discrete model. Note that our construction of discrete
versions of the Hodge-Dirac and Laplace operators is very close to the
construction in Section 5 of [6]. Matrix forms of these discrete operators on
the torus are the same in both of case. The difference between our approach
and one in [6] is in the definitions of discrete Hodge-Dirac and Laplace
operators. In [6], the definitions are given in terms of the exterior
derivative and boundary operators. Meanwhile, as in the continual theory, we
define these operators in terms of the exterior derivative and its adjoint.

##  2 Discrete model

In this section, we briefly review the construction of a discrete exterior
calculus framework, which was initiated in [7] and developed in e.g., [14,
15].

The starting point for consideration is a two-dimensional chain complex (a
combinatorial model of $\mathbb{R}^{2}$). Let the sets $\\{x_{k}\\}$ and
$\\{e_{k}\\}$, $k\in{\mathbb{Z}}$, be the generators of free abelian groups of
zero-dimensional and one-dimensional chains of the one-dimensional complex
$C=C_{0}\oplus C_{1}$. The free abelian group is understood as the direct sum
of infinity cyclic groups generated by $\\{x_{k}\\}$, $\\{e_{k}\\}$. The
boundary operator $\partial:C_{1}\rightarrow C_{0}$ is the homomorphism
defined by $\partial e_{k}=x_{k+1}-x_{k}$ and the boundary of every zero chain
is defined to be zero. Geometrically we can interpret the zero-dimensional
basis elements $x_{k}$ as points of the real line and the one-dimensional
basis elements $e_{k}$ as open intervals between points. We call the complex
$C$ a combinatorial real line. Let the tensor product

$C(2)=C\otimes C$

be a combinatorial model of the two-dimensional Euclidean space
$\mathbb{R}^{2}$. The two-dimensional complex $C(2)=C_{0}(2)\oplus
C_{1}(2)\oplus C_{2}(2)$ contains the free abelian groups of $r$-chains,
$r=0,1,2$, generated by the basic elements

$\displaystyle x_{k,s}=x_{k}\otimes x_{s},\qquad e_{k,s}^{1}=e_{k}\otimes
x_{s},\qquad e_{k,s}^{2}=x_{k}\otimes e_{s},\qquad V_{k,s}=e_{k}\otimes
e_{s},$

where $k,s\in{\mathbb{Z}}$. It is convenient to introduce the shift operators
$\tau,\sigma$ in the set of indices by

$\tau k=k+1\qquad\sigma k=k-1.$ (2.1)

The boundary operator $\partial:C_{r}(2)\rightarrow C_{r-1}(2)$ is given by

$\displaystyle\partial x_{k,s}=0,\qquad\partial e_{k,s}^{1}=x_{\tau
k,s}-x_{k,s}\qquad\partial e_{k,s}^{2}=x_{k,\tau s}-x_{k,s},$
$\displaystyle\partial V_{k,s}=e_{k,s}^{1}+e_{\tau k,s}^{2}-e_{k,\tau
s}^{1}-e_{k,s}^{2}.$ (2.2)

The definition (2) is extended to arbitrary chains by linearity.

Let $K(2)$ be a complex of cochains with real coefficients. The cochain
complex $K(2)$ with a coboundary operator defined in it is the dual object to
the chain complex $C(2)$. It has a similar structure to $C(2)$ and consists of
cochains of dimension 0, 1 and 2. Then $K(2)$ can be expressed by

$K(2)=K^{0}(2)\oplus K^{1}(2)\oplus K^{2}(2),$

where $K^{r}(2)$ is the set of all $r$-cochains. We will call cochains forms
(or discrete forms) emphasizing their relationship with differential forms.
Then the complex $K(2)$ is a discrete analogue of the grade algebra of
differential forms $\Lambda(\mathbb{R}^{2})$. Denote by $\\{x^{k,s}\\}$,
$\\{e^{k,s}_{1},\ e^{k,s}_{2}\\}$ and $\\{V^{k,s}\\}$ the basis elements of
$K^{0}(2)$, $K^{1}(2)$ and $K^{2}(2)$ respectively. The pairing is defined
with the basis elements of $C(2)$ by the rule

$\langle x_{k,s},\ x^{i,j}\rangle=\langle e_{k,s}^{1},\
e^{i,j}_{1}\rangle=\langle e_{k,s}^{2},\ e^{i,j}_{2}\rangle=\langle V_{k,s},\
V^{i,j}\rangle=\delta_{k}^{i}\delta_{s}^{j},$ (2.3)

where $\delta_{k}^{i}$ is the Kronecker delta. The operation (2.3) is linearly
extended to arbitrary chains and cochains. Let $\overset{r}{\omega}\in
K^{r}(2)$, then we have

$\overset{0}{\omega}=\sum_{k,s}\overset{0}{\omega}_{k,s}x^{k,s},\quad\overset{1}{\omega}=\sum_{k,s}(\omega^{1}_{k,s}e_{1}^{k,s}+\omega^{2}_{k,s}e_{2}^{k,s}),\quad\overset{2}{\omega}=\sum_{k,s}\overset{2}{\omega}_{k,s}V^{k,s},$
(2.4)

where $\overset{0}{\omega}_{k,s}$, $\omega^{1}_{k,s}$, $\omega^{2}_{k,s}$ and
$\overset{2}{\omega}_{k,s}$ are real numbers for any $k,s\in{\mathbb{Z}}$.

The coboundary operator $d^{c}:K^{r}(2)\rightarrow K^{r+1}(2)$ is defined by

$\langle\partial a_{r+1},\ \overset{r}{\omega}\rangle=\langle a_{r+1},\
d^{c}\overset{r}{\omega}\rangle,$ (2.5)

where $a_{r+1}\in C_{r+1}(2)$. The operator $d^{c}$ is an analog of the
exterior differential. From the above it follows that

$d^{c}\overset{2}{\omega}=0\quad\mbox{and}\quad
d^{c}d^{c}\overset{r}{\omega}=0\quad\mbox{for any}\quad r=0,1.$

By (2) and (2.3) we can calculate

$d^{c}\overset{0}{\omega}=\sum_{k,s}(\Delta_{k}\overset{0}{\omega}_{k,s})e_{1}^{k,s}+(\Delta_{s}\overset{0}{\omega}_{k,s})e_{2}^{k,s},$
(2.6)
$d^{c}\overset{1}{\omega}=\sum_{k,s}(\Delta_{k}\omega_{k,s}^{2}-\Delta_{s}\omega_{k,s}^{1})V^{k,s},$
(2.7)

where $\Delta_{k}$ and $\Delta_{s}$ are the difference operators defined by

$\Delta_{k}\overset{r}{\omega}_{k,s}=\overset{r}{\omega}_{\tau
k,s}-\overset{r}{\omega}_{k,s},\qquad\Delta_{s}\overset{r}{\omega}_{k,s}=\overset{r}{\omega}_{k,\tau
s}-\overset{r}{\omega}_{k,s}.$ (2.8)

Here $\overset{r}{\omega}_{k,s}$ is a component of $\overset{r}{\omega}\in
K^{r}(2)$ and $\tau$ is given by (2.1). Note that
$\overset{1}{\omega}_{k,s}=\\{\omega^{1}_{k,s},\ \omega^{2}_{k,s}\\}$.

We now consider a multiplication of discrete forms which is an analogue of the
exterior multiplication for differential forms. Denote by $\cup$ this
multiplication. For the basis elements of $K(2)$ the $\cup$-multiplication is
defined as follows

$x^{k,s}\cup x^{k,s}=x^{k,s},\qquad x^{k,s}\cup e^{k,s}_{1}=e^{k,s}_{1},\qquad
x^{k,s}\cup e^{k,s}_{2}=e^{k,s}_{2},$ $x^{k,s}\cup V^{k,s}=V^{k,s},\qquad
e^{k,s}_{1}\cup x^{\tau k,s}=e^{k,s}_{1},\qquad e^{k,s}_{2}\cup x^{k,\tau
s}=e^{k,s}_{2},$ $V^{k,s}\cup x^{\tau k,\tau s}=V^{k,s},\qquad e^{k,s}_{1}\cup
e^{\tau k,s}_{2}=V^{k,s},\qquad e^{k,s}_{2}\cup e^{k,\tau s}_{1}=-V^{k,s},$

supposing the product to be zero in all other cases. The operation is extended
to arbitrary forms by linearity. It is important to note that this definition
leads to the following discrete counterpart of the Leibniz rule for
differential forms.

######  Proposition 2.1.

Let $\overset{r}{\omega}\in K^{r}(2)$ and $\overset{p}{\varphi}\in K^{p}(2)$.
Then

$d^{c}(\overset{r}{\omega}\cup\overset{p}{\varphi})=d^{c}\overset{r}{\omega}\cup\overset{p}{\varphi}+(-1)^{r}\overset{r}{\omega}\cup
d^{c}\overset{p}{\varphi}.$ (2.9)

This was proved by Dezin [7].

Define the operation $\ast:K^{r}(2)\rightarrow K^{2-r}(2)$ by the rule

$\ast x^{k,s}=V^{k,s},\quad\ast e^{k,s}_{1}=e^{\tau k,s}_{2},\quad\ast
e^{k,s}_{2}=-e^{k,\tau s}_{1},\quad\ast V^{k,s}=x^{\tau k,\tau s}.$ (2.10)

Again, the operation is extended to arbitrary forms by linearity. This
operation is a discrete analogue of the Hodge star operator. It is true that
for any $\overset{r}{\omega}\in K^{r}(2)$ we have

$\overset{r}{\omega}\cup\ast\overset{r}{\omega}=\sum_{k,s}(\overset{r}{\omega}_{k,s})^{2}V^{k,s}.$
(2.11)

Consider the two-dimensional finite chain $V\in C_{2}(2)$ with unit
coefficients of the form

$V=\sum_{k,s}V_{k,s},\qquad k=1,2,...,N,\quad s=1,2,..,M.$ (2.12)

This chain imitates a rectangle. Using (2) we have

$\partial
V=\sum_{k=1}^{N}e^{1}_{k,1}+\sum_{s=1}^{M}e^{2}_{N,s}-\sum_{k=1}^{N}e^{1}_{k,M}-\sum_{s=1}^{M}e^{2}_{1,s}.$
(2.13)

Then for forms $\overset{r}{\varphi},\ \overset{r}{\omega}\in K^{r}(2)$ of the
same degree $r$ the inner product over the set $V$ is defined by the rule

$(\overset{r}{\varphi},\ \overset{r}{\omega})_{V}=\langle V,\
\overset{r}{\varphi}\cup\ast\overset{r}{\omega}\rangle.$ (2.14)

For forms of different degrees the product (2.14) is set equal to zero. From
(2.11) and (2.3) we have

$(\overset{r}{\varphi},\
\overset{r}{\omega})_{V}=\sum_{k=1}^{N}\sum_{s=1}^{M}\overset{r}{\varphi}_{k,s}\overset{r}{\omega}_{k,s}.$

######  Proposition 2.2.

Let $\overset{r}{\varphi}\in K^{r}(2)$ and $\overset{r+1}{\omega}\in
K^{r+1}(2)$, $r=0,1$. Then we have

$(d^{c}\overset{r}{\varphi},\ \overset{r+1}{\omega})_{V}=\langle\partial V,\
\overset{r}{\varphi}\cup\ast\overset{r+1}{\omega}\rangle+(\overset{r}{\varphi},\
\delta^{c}\overset{r+1}{\omega})_{V},$ (2.15)

where

$\delta^{c}\overset{r+1}{\omega}=(-1)^{r+1}\ast^{-1}d^{c}\ast\overset{r+1}{\omega}$
(2.16)

is the operator formally adjoint of $d^{c}$.

Here $\ast^{-1}$ is the inverse of $\ast$, i.e., $\ast\ast^{-1}=1$. By (2.10)
for the basic elements we have

$\ast^{-1}x^{k,s}=V^{\sigma k,\sigma s},\quad\ast^{-1}e^{k,s}_{1}=-e^{k,\sigma
s}_{2},\quad\ast^{-1}e^{k,s}_{2}=e^{\sigma
k,s}_{1},\quad\ast^{-1}V^{k,s}=x^{k,s},$

where $\sigma$ is given by (2.1).

###### Proof.

By Definitions (2.5), (2.14) and Formula (2.9) we obtain

$\displaystyle(d^{c}\overset{r}{\varphi},\ \overset{r+1}{\omega})_{V}=\langle
V,\ d^{c}\overset{r}{\varphi}\cup\ast\overset{r+1}{\omega}\rangle=\langle V,\
d^{c}(\overset{r}{\varphi}\cup\ast\overset{r+1}{\omega})-(-1)^{r}\overset{r}{\varphi}\cup
d^{c}\ast\overset{r+1}{\omega}\rangle$ $\displaystyle=\langle\partial V,\
\overset{r}{\varphi}\cup\ast\overset{r+1}{\omega}\rangle+(-1)^{r+1}\langle V,\
\overset{r}{\varphi}\cup\ast\ast^{-1}d^{c}\ast\overset{r+1}{\omega}\rangle$
$\displaystyle=\langle\partial V,\
\overset{r}{\varphi}\cup\ast\overset{r+1}{\omega}\rangle+\langle V,\
\overset{r}{\varphi}\cup\ast(\delta^{c}\overset{r+1}{\omega})\rangle.$

∎

The operator $\delta^{c}:K^{r+1}(2)\rightarrow K^{r}(2)$ given by (2.16) is a
discrete analogue of the codifferential $\delta$. For the 0-form
$\overset{0}{\omega}\in K^{0}(2)$ we have $\delta^{c}\overset{0}{\omega}=0$.
It is obvious from (2.16) that $\delta^{c}\delta^{c}\overset{r}{\omega}=0$ for
any $r=1,2.$ Using (2.6)–(2.8), (2.10) and (2.16) we can calculate

$\delta^{c}\overset{1}{\omega}=\sum_{k,s}(-\Delta_{k}\omega_{\sigma
k,s}^{1}-\Delta_{s}\omega_{k,\sigma s}^{2})x^{k,s},$ (2.17)
$\delta^{c}\overset{2}{\omega}=\sum_{k,s}(\Delta_{s}\overset{2}{\omega}_{k,\sigma
s})e_{1}^{k,s}-(\Delta_{k}\overset{2}{\omega}_{\sigma k,s})e_{2}^{k,s}.$
(2.18)

In the particular case $r=0$, the equality (2.15) can be expressed as

$\displaystyle\sum_{k=1}^{N}\sum_{s=1}^{M}((\Delta_{k}\overset{0}{\varphi}_{k,s}){\omega}_{k,s}^{1}+(\Delta_{s}\overset{0}{\varphi}_{k,s}){\omega}_{k,s}^{2})$
$\displaystyle=\sum_{k=1}^{N}(\overset{0}{\varphi}_{k,\tau
M}\omega^{2}_{k,M}-\overset{0}{\varphi}_{k,1}\omega^{2}_{k,0})+\sum_{s=1}^{M}(\overset{0}{\varphi}_{\tau
N,s}\omega^{1}_{N,s}-\overset{0}{\varphi}_{1,s}\omega^{1}_{0,s})$
$\displaystyle+\sum_{k=1}^{N}\sum_{s=1}^{M}\overset{0}{\varphi}_{k,s}(-\Delta_{k}\omega_{\sigma
k,s}^{1}-\Delta_{s}\omega_{k,\sigma s}^{2}),$

where $\overset{1}{\omega}_{k,s}=\\{\omega^{1}_{k,s},\ \omega^{2}_{k,s}\\}$.
The similar equality holds in the case $r=1$.

It should be noted that the relation (2.15) includes not only the forms with
the components $\overset{r}{\varphi}_{k,s}$ and $\overset{r+1}{\omega}_{k,s}$,
where the subscripts $k,s$ would run only over the values from (2.12), but
also the components $\overset{r}{\varphi}_{0,s}$, $\overset{r}{\varphi}_{\tau
N,s}$, $\overset{r+1}{\omega}_{0,s}$, $\overset{r+1}{\omega}_{\tau N,s}$,
$\overset{r}{\varphi}_{k,0}$, $\overset{r}{\varphi}_{k,\tau M}$,
$\overset{r+1}{\omega}_{k,0}$ and $\overset{r+1}{\omega}_{k,\tau M}$.

Let us set

$\displaystyle\overset{r}{\omega}_{0,s}=\overset{r}{\omega}_{N,s},\qquad\overset{r}{\omega}_{\tau
N,s}=\overset{r}{\omega}_{1,s},\qquad s=1,2,..,M,$
$\displaystyle\overset{r}{\omega}_{k,0}=\overset{r}{\omega}_{k,M},\qquad\overset{r}{\omega}_{k,\tau
M}=\overset{r}{\omega}_{k,1},\qquad k=1,2,...,N.$ (2.19)

For $r$-forms satisfying conditions (2) the inner product (2.14) generates the
finite-dimensional Hilbert spaces $H^{r}(V)$. Now we consider the operators

$d^{c}:H^{r}(V)\rightarrow H^{r+1}(V),\qquad\delta^{c}:H^{r+1}(V)\rightarrow
H^{r}(V).$

######  Proposition 2.3.

Let $\overset{r}{\varphi}\in H^{r}(V)$ and $\overset{r+1}{\omega}\in
H^{r+1}(V)$, $r=0,1$. Then we have

$(d^{c}\overset{r}{\varphi},\
\overset{r+1}{\omega})_{V}=(\overset{r}{\varphi},\
\delta^{c}\overset{r+1}{\omega})_{V}.$ (2.20)

###### Proof.

In fact, by use of conditions (2) one has

$\langle\partial V,\
\overset{0}{\varphi}\cup\ast\overset{1}{\omega}\rangle=\sum_{k=1}^{N}(\overset{0}{\varphi}_{k,\tau
M}\omega^{2}_{k,M}-\overset{0}{\varphi}_{k,1}\omega^{2}_{k,0})+\sum_{s=1}^{M}(\overset{0}{\varphi}_{\tau
N,s}\omega^{1}_{N,s}-\overset{0}{\varphi}_{1,s}\omega^{1}_{0,s})=0,$
$\langle\partial V,\
\overset{1}{\varphi}\cup\ast\overset{2}{\omega}\rangle=\sum_{k=1}^{N}(\varphi^{1}_{k,1}\overset{2}{\omega}_{k,0}-\varphi^{1}_{k,\tau
M}\overset{2}{\omega}_{k,M})+\sum_{s=1}^{M}(\varphi^{2}_{\tau
N,s}\overset{2}{\omega}_{N,s}-\varphi^{2}_{1,s}\overset{2}{\omega}_{0,s})=0,$

where $\overset{1}{\omega}_{k,s}=\\{\omega^{1}_{k,s},\ \omega^{2}_{k,s}\\}$
and $\overset{1}{\varphi}_{k,s}=\\{\varphi^{1}_{k,s},\ \varphi^{2}_{k,s}\\}$.
Hence by Proposition 2.2 it follows (2.20). ∎

##  3 Discrete Hodge decomposition

In this section, we discuss the properties of discrete analogues of the
Laplacian and Hodge-Dirac operators using the concepts of the previous
section. We also present a discrete version of the Hodge decomposition
theorem, emphasizing that it provides an exact counterpart to the continuum
theory.

Let us consider the operator

$\Delta^{c}=d^{c}\delta^{c}+\delta^{c}d^{c}:H^{r}(V)\rightarrow H^{r}(V).$

This is a discrete analogue of the Laplacian (1.2).

######  Proposition 3.1.

For any $r$-form $\varphi\in H^{r}(V)$ we have

$(\Delta^{c}\varphi,\ \varphi)_{V}\geq 0$

and $(\Delta^{c}\varphi,\ \varphi)_{V}=0$ if and only if
$\Delta^{c}\varphi=0$.

###### Proof.

By Proposition 2.3 one has

$\displaystyle(\Delta^{c}\varphi,\ \varphi)_{V}=(d^{c}\delta^{c}\varphi,\
\varphi)_{V}+(\delta^{c}d^{c}\varphi,\ \varphi)_{V}$
$\displaystyle=(\delta^{c}\varphi,\ \delta^{c}\varphi)_{V}+(d^{c}\varphi,\
d^{c}\varphi)_{V}=\|\delta^{c}\varphi\|^{2}+\|d^{c}\varphi\|^{2},$

where $\|\cdot\|$ denotes the norm and

$\|\varphi\|^{2}=(\varphi,\
\varphi)_{V}=\sum_{k=1}^{N}\sum_{s=1}^{M}({\varphi}_{k,s})^{2}.$

From this if $(\Delta^{c}\varphi,\ \varphi)_{V}=0$ then
$\|\delta^{c}\varphi\|^{2}=0$ and $\|d^{c}\varphi\|^{2}=0$. It gives
$\delta^{c}\varphi=0$ and $d^{c}\varphi=0$. Hence

$\Delta^{c}\varphi=d^{c}\delta^{c}\varphi+\delta^{c}d^{c}\varphi=0.$

∎

######  Corollary 3.2.

$\Delta^{c}\varphi=0$ if and only if $\delta^{c}\varphi=0$ and
$d^{c}\varphi=0$.

######  Proposition 3.3.

The operator $\Delta^{c}:H^{r}(V)\rightarrow H^{r}(V)$ is self-adjoint, i.e.

$(\Delta^{c}\varphi,\ \omega)_{V}=(\varphi,\ \Delta^{c}\omega)_{V}.$

###### Proof.

By (2.20) it is obvious. ∎

Consider the spaces

$R_{d^{c}}^{r}=\\{d^{c}\varphi\in H^{r}(V):\ \varphi\in H^{r-1}(V)\\},$
$R_{\delta^{c}}^{r}=\\{\delta^{c}\omega\in H^{r}(V):\ \omega\in
H^{r+1}(V)\\},$

and

$N_{\Delta^{c}}^{r}=\\{\psi\in H^{r}(V):\ \Delta^{c}\psi=0\\}.$

By analogy with the continuum case the discrete $r$-form $\omega$ is called
closed if $d^{c}\omega=0$ and exact if $\omega\in R_{d^{c}}^{r}$.

######  Proposition 3.4.

For each $r=0,1,2$ we have the direct sum decomposition

$H^{r}(V)=R_{d^{c}}^{r}\oplus R_{\delta^{c}}^{r}\oplus N_{\Delta^{c}}^{r}.$

###### Proof.

The space $H^{1}(V)$ decomposes into

$H^{1}(V)=R_{d^{c}}^{1}\oplus N_{\delta^{c}}^{1},\qquad
H^{1}(V)=R_{\delta^{c}}^{1}\oplus N_{d^{c}}^{1},$ (3.1)

where $N_{\delta^{c}}^{1}$ and $N_{d^{c}}^{1}$ are the orthogonal complements
of the corresponding spaces. For any $\omega\in R_{\delta^{c}}^{1}$ we have
$\omega=\delta^{c}\psi$ and

$(d^{c}\varphi,\ \omega)_{V}=(d^{c}\varphi,\
\delta^{c}\psi)_{V}=(d^{c}d^{c}\varphi,\ \psi)_{V}=0$

for any $\varphi\in H^{0}(V)$. Therefore $\omega$ is orthogonal to
$R_{d^{c}}^{1}$. It follows that

$R_{\delta^{c}}^{1}\subset N_{\delta^{c}}^{1}.$

Similarly we find $R_{d^{c}}^{1}\subset N_{\delta^{c}}^{1}$. Hence (3.1)
becomes

$H^{1}(V)=R_{d^{c}}^{1}\oplus R_{\delta^{c}}^{1}\oplus N^{1},$

where

$N^{1}=N_{\delta^{c}}^{1}\cap N_{d^{c}}^{1}=\\{\varphi\in H^{1}(V):\
d^{c}\varphi=0,\ \delta^{c}\varphi=0\\}.$

By Corollary 3.2 we have $N^{1}=N_{\Delta^{c}}^{1}$.

Similar reasonings apply to the spaces $H^{0}(V)$ and $H^{2}(V)$. Thus we have

$H^{0}(V)=R_{\delta^{c}}^{0}\oplus N_{\Delta^{c}}^{0},\qquad
H^{2}(V)=R_{d^{c}}^{2}\oplus N_{\Delta^{c}}^{2},$

since $R_{d^{c}}^{0}=\\{0\\}$ and $R_{\delta^{c}}^{2}=\\{0\\}$. ∎

The Proposition 3.4 is a discrete version of the well-known Hodge
decomposition theorem (see, e.g., [21]).

Let $\Omega$ be an inhomogeneous discrete form, i.e.
$\Omega=\overset{0}{\omega}+\overset{1}{\omega}+\overset{2}{\omega}$, where
$\overset{r}{\omega}\in H^{r}(V)$. The inner product (2.14) can be extend to
an inner product of inhomogeneous discrete forms by the rule

$(\Omega,\ \Phi)_{V}=\sum_{r=0}^{2}(\overset{r}{\omega},\
\overset{r}{\varphi})_{V},$ (3.2)

where $\Phi=\overset{0}{\varphi}+\overset{1}{\varphi}+\overset{2}{\varphi}$.
The inner product (3.2) generates the finite-dimensional Hilbert space $H(V)$.
It is true that

$H(V)=H^{0}(V)\oplus H^{1}(V)\oplus H^{2}(V).$

By Proposition 3.4 the following holds

$H(V)=R_{d^{c}}\oplus R_{\delta^{c}}\oplus N_{\Delta^{c}},$ (3.3)

where

$R_{d^{c}}=R_{d^{c}}^{1}\oplus R_{d^{c}}^{2},\qquad
R_{\delta^{c}}=R_{\delta^{c}}^{0}\oplus R_{\delta^{c}}^{1},$

and

$N_{\Delta^{c}}=N_{\Delta^{c}}^{0}\oplus N_{\Delta^{c}}^{1}\oplus
N_{\Delta^{c}}^{2}.$

The discrete Hodge-Dirac operator is defined as

$d^{c}+\delta^{c}:H(V)\rightarrow H(V).$ (3.4)

######  Proposition 3.5.

Operator (3.4) is self-adjoint with respect to the inner product (3.2), i.e.

$((d^{c}+\delta^{c})\Phi,\ \Omega)_{V}=(\Phi,\ (d^{c}+\delta^{c})\Omega)_{V}.$

###### Proof.

By (2.20) it is obvious. ∎

######  Proposition 3.6.

$(d^{c}+\delta^{c})\Omega=0\quad\mbox{if and only if}\quad\Omega\in
N_{\Delta^{c}}.$

###### Proof.

The equation $(d^{c}+\delta^{c})\Omega=0$ can be written as

$\displaystyle\delta^{c}\overset{1}{\omega}=0,\quad
d^{c}\overset{1}{\omega}=0,$ $\displaystyle
d^{c}\overset{0}{\omega}+\delta^{c}\overset{2}{\omega}=0.$

Let $\Omega\in N_{\Delta^{c}}$. This means that $\overset{0}{\omega}\in
N_{\Delta^{c}}^{0}$, $\overset{1}{\omega}\in N_{\Delta^{c}}^{1}$ and
$\overset{2}{\omega}\in N_{\Delta^{c}}^{2}$ By Corollary 3.2,
$\overset{1}{\omega}\in N_{\Delta^{c}}^{1}$ if and only if
$\delta^{c}\overset{1}{\omega}=0$ and $d^{c}\overset{1}{\omega}=0$. It is easy
to show that

$\|d^{c}\overset{0}{\omega}+\delta^{c}\overset{2}{\omega}\|^{2}=(d^{c}\overset{0}{\omega}+\delta^{c}\overset{2}{\omega},\
d^{c}\overset{0}{\omega}+\delta^{c}\overset{2}{\omega})_{V}=\|d^{c}\overset{0}{\omega}\|^{2}+\|\delta^{c}\overset{2}{\omega}\|^{2}$
(3.5)

for any $\overset{0}{\omega}\in H^{0}(V)$ and $\overset{2}{\omega}\in
H^{2}(V)$. For $\overset{0}{\omega}\in N_{\Delta^{c}}^{0}$ and
$\overset{2}{\omega}\in N_{\Delta^{c}}^{2}$, by (2.20) and (3.5), we have

$0=(\delta^{c}d^{c}\overset{0}{\omega},\overset{0}{\omega})_{V}+(d^{c}\delta^{c}\overset{2}{\omega},\overset{2}{\omega})_{V}=(d^{c}\overset{0}{\omega},\
d^{c}\overset{0}{\omega})_{V}+(\delta^{c}\overset{2}{\omega},\
\delta^{c}\overset{2}{\omega})_{V}=\|d^{c}\overset{0}{\omega}+\delta^{c}\overset{2}{\omega}\|^{2}$

and thus $d^{c}\overset{0}{\omega}+\delta^{c}\overset{2}{\omega}=0$. The
converse is trivially true. ∎

######  Proposition 3.7.

For any inhomogeneous form $F\in R_{d^{c}}\oplus R_{\delta^{c}}$, there exists
a unique form $\Omega\in R_{d^{c}}\oplus R_{\delta^{c}}$ which is a solution
to the equation

$(d^{c}+\delta^{c})\Omega=F.$ (3.6)

###### Proof.

Since the operator (3.4) is self-adjoint and $H(V)$ is a finite-dimensional
Hilbert space, the existence of the solution is a consequence of the
uniqueness of the solution and vice versa. By Proposition 3.6,
$(d^{c}+\delta^{c})\Omega=0$ implies $\Omega\in N_{\Delta^{c}}$. From this and
by (3.3), the uniqueness of the solution $\Omega\in R_{d^{c}}\oplus
R_{\delta^{c}}$ of Eq. (3.6) follows immediately. ∎

######  Corollary 3.8.

If $\Omega\in H(V)$ is a solution of Eq. (3.6) then the following holds

$\|\Omega\|^{2}\leq
c(\|d^{c}\Omega\|^{2}+\|\delta^{c}\Omega\|^{2})+\|\Omega_{\Delta^{c}}\|^{2},$

where $c$ is a constant,
$d^{c}\Omega=d^{c}\overset{0}{\omega}+d^{c}\overset{1}{\omega}$,
$\delta^{c}\Omega=\delta^{c}\overset{1}{\omega}+\delta^{c}\overset{2}{\omega}$
and $\Omega_{\Delta^{c}}$ is the projection of $\Omega$ onto $N_{\Delta^{c}}$.

##  4 Combinatorial torus

In this section, we consider an example of a discrete model of the torus in
detail. We recall that the torus can be regarded as the topological space
obtained by taking a rectangle and identifying each pair of opposite sides
with the same orientation. Let consider the partitioning of the plane
$\mathbb{R}^{2}$ by the straight lines $x=k$ and $y=s$, where
$k,s\in\mathbb{Z}$. Denote by $V_{k,s}$ an open square bounded by the lines
$x=k,\ x=\tau k,\ y=s$ and $y=\tau s$, where $\tau$ is given by (2.1). Denote
the vertices of $V_{k,s}$ by $x_{k,s},\ x_{\tau k,s},\ x_{k,\tau s}$, $x_{\tau
k,\tau s}$. Let $e_{k,s}^{1}$ and $e_{k,s}^{2}$ be the open intervals
$(x_{k,s},\ x_{\tau k,s})$ and $(x_{k,s},\ x_{k,\tau s})$, respectively. Such
introduced geometric objects can be identified with the combinatorial objects
we have considered in the previous sections. We identify the collection
$V_{k,s}$ with $V$ given by (2.12) and let $N=M=2$. In this case the
conditions (2) take the form

$\overset{r}{\omega}_{0,s}=\overset{r}{\omega}_{2,s},\qquad\overset{r}{\omega}_{3,s}=\overset{r}{\omega}_{1,s},\qquad\overset{r}{\omega}_{k,0}=\overset{r}{\omega}_{k,2},\qquad\overset{r}{\omega}_{k,3}=\overset{r}{\omega}_{k,1},$
(4.1)

where $k=1,2$ and $s=1,2$. If we identify the points and the intervals on the
boundary of $V$ in the following way

$\displaystyle x_{1,1}=x_{3,1}=x_{1,3}=x_{3,3},\qquad x_{1,2}=x_{3,2},\qquad
x_{2,1}=x_{2,3},$ $\displaystyle e_{1,1}^{1}=e_{1,3}^{1},\qquad
e_{2,1}^{1}=e_{2,3}^{1},\qquad e_{1,1}^{2}=e_{3,1}^{2},\qquad
e_{1,2}^{2}=e_{3,2}^{2},$ (4.2)

we obtain the geometric object which is homomorphic to the torus (see Figure
1). Denote by $C(T)$ the complex $C(2)$ which corresponds to the introduced
geometric object. We call $C(T)$ a combinatorial torus.

$x_{1,1}$$x_{2,1}$$x_{1,1}$$x_{2,2}$$e^{1}_{1,1}$$e^{1}_{1,1}$$e^{1}_{2,1}$$e^{1}_{2,1}$$e^{2}_{1,1}$$e^{2}_{1,2}$$e^{2}_{2,2}$$e^{2}_{1,2}$$e^{2}_{1,1}$$e^{2}_{2,1}$$e^{1}_{1,2}$$e^{1}_{2,2}$$x_{1,2}$$x_{1,1}$$x_{1,2}$$x_{1,1}$$x_{2,1}$$V_{1,1}$$V_{2,1}$$V_{2,2}$$V_{1,2}$
Figure 1: Combinatorial torus

It is clear that by (4) the conditions (4.1) hold for any $r$-form on the
combinatorial torus. The Hilbert space $H(V)$ considered in the previous
sections can be regarded as a space of cochains of the complex $K(T)$ dual to
$C(T)$.

Let now consider the forms $\varphi\in K^{0}(T)$, $\omega\in K^{1}(T)$ and
$\psi\in K^{2}(T)$, that is

$\varphi=\sum_{k=1}^{2}\sum_{s=1}^{2}\varphi_{k,s}x^{k,s},\quad\omega=\sum_{k=1}^{2}\sum_{s=1}^{2}(u_{k,s}e_{1}^{k,s}+v_{k,s}e_{2}^{k,s}),\quad\psi=\sum_{k=1}^{2}\sum_{s=1}^{2}\psi_{k,s}V^{k,s}.$

By (4) for these forms the formulas (2.6), (2.7), (2.17) and (2.18) become

$\displaystyle
d^{c}\varphi=(\varphi_{2,1}-\varphi_{1,1})e_{1}^{1,1}+(\varphi_{1,1}-\varphi_{2,1})e_{1}^{2,1}+(\varphi_{2,2}-\varphi_{1,2})e_{1}^{1,2}$
$\displaystyle+(\varphi_{1,2}-\varphi_{2,2})e_{1}^{2,2}+(\varphi_{1,2}-\varphi_{1,1})e_{2}^{1,1}+(\varphi_{1,1}-\varphi_{1,2})e_{2}^{1,2}$
$\displaystyle+(\varphi_{2,2}-\varphi_{2,1})e_{2}^{2,1}+(\varphi_{2,1}-\varphi_{2,2})e_{2}^{2,2},$
(4.3) $\displaystyle
d^{c}\omega=(u_{1,1}-u_{1,2}+v_{2,1}-v_{1,1})V^{1,1}+(u_{2,1}-u_{2,2}-v_{2,1}+v_{1,1})V^{2,1}$
$\displaystyle+(u_{1,2}-u_{1,1}+v_{2,2}-v_{1,2})V^{1,2}+(u_{2,2}-u_{2,1}+v_{1,2}-v_{2,2})V^{2,2},$
(4.4)
$\displaystyle\delta^{c}\omega=(u_{2,1}-u_{1,1}+v_{1,2}-v_{1,1})x^{1,1}+(u_{1,1}-u_{2,1}-v_{2,1}+v_{2,2})x^{2,1}$
$\displaystyle+(u_{2,2}-u_{1,2}+v_{1,1}-v_{1,2})x^{1,2}+(u_{1,2}-u_{2,2}+v_{2,1}-v_{2,2})x^{2,2},$
(4.5)
$\displaystyle\delta^{c}\psi=(\psi_{1,1}-\psi_{1,2})e_{1}^{1,1}+(\psi_{2,1}-\psi_{2,2})e_{1}^{2,1}+(\psi_{2,2}-\psi_{1,2})e_{2}^{1,2}$
$\displaystyle+(\psi_{2,1}-\psi_{1,1})e_{2}^{1,1}+(\psi_{1,2}-\psi_{1,1})e_{1}^{1,2}+(\psi_{2,2}-\psi_{2,1})e_{1}^{2,2}$
$\displaystyle+(\psi_{1,2}-\psi_{2,2})e_{2}^{2,2}+(\psi_{1,1}-\psi_{2,1})e_{2}^{2,1}.$
(4.6)

It should be noted that the formulas above can be expressed in matric form.
Let introduce the following row vectors

$[\varphi]=[\varphi_{1,1}\ \varphi_{2,1}\ \varphi_{1,2}\
\varphi_{2,2}],\quad[\omega]=[u_{1,1}\ u_{2,1}\ v_{1,2}\ v_{1,1}\ u_{1,2}\
u_{2,2}\ v_{2,2}\ v_{2,1}],$ $[x]=[x^{1,1}\ x^{2,1}\ x^{1,2}\
x^{2,2}],\quad[e]=[e_{1}^{1,1}\ e_{1}^{2,1}\ e_{2}^{1,2}\ e_{2}^{1,1}\
e_{1}^{1,2}\ e_{1}^{2,2}\ e_{2}^{2,2}\ e_{2}^{2,1}],$ $[\psi]=[\psi_{1,2}\
\psi_{2,2}\ \psi_{1,1}\ \psi_{2,1}],\qquad[V]=[V^{1,2}\ V^{2,2}\ V^{1,1}\
V^{2,1}].$

Denote by $[\cdot]^{T}$ a corresponding column vector. Then we have

$d^{c}\varphi=A[\varphi]^{T}[e],\quad
d^{c}\omega=B[\omega]^{T}[V],\quad\delta^{c}\omega=A^{T}[\omega]^{T}[x],\quad\delta^{c}\psi=B^{T}[\psi]^{T}[e],$

where

$A=\begin{bmatrix}-1&1&0&0\\\ 1&-1&0&0\\\ 1&0&-1&0\\\ -1&0&1&0\\\ 0&0&-1&1\\\
0&0&1&-1\\\ 0&1&0&-1\\\ 0&-1&0&1\end{bmatrix},\quad
B=\begin{bmatrix}-1&0&-1&0&1&0&1&0\\\ 0&-1&1&0&0&1&-1&0\\\
1&0&0&-1&-1&0&0&1\\\ 0&1&0&1&0&-1&0&-1\end{bmatrix}$

and $A^{T}$, $B^{T}$ are the transpose of $A$, $B$.

The discrete Hodge-Dirac operator (3.4) on the combinatorial torus can be
represented by the following block matrix

$\begin{bmatrix}\mbox{\large 0}&A^{T}&\mbox{\large 0}\\\ A&\mbox{\large
0}&B^{T}\\\ \mbox{\large 0}&B&\mbox{\large 0}\end{bmatrix},$

where 0 is a zero square matrix of the corresponding size.

In the same way the discrete Laplacian on the combinatorial torus can be
written as

$\Delta^{c}\varphi=D[\varphi]^{T}[x],\quad\Delta^{c}\omega=D_{1}[\omega]^{T}[e],\quad\Delta^{c}\psi=D[\psi]^{T}[V],$

where

$D=\begin{bmatrix}4&-2&-2&0\\\ -2&4&0&-2\\\ -2&0&4&-2\\\
0&-2&-2&4\end{bmatrix},\ D_{1}=\begin{bmatrix}4&-2&0&0&-2&0&0&0\\\
-2&4&0&0&0&-2&0&0\\\ 0&0&4&-2&0&0&-2&0\\\ 0&0&-2&4&0&0&0&-2\\\
-2&0&0&0&4&-2&0&0\\\ 0&-2&0&0&-2&4&0&0\\\ 0&0&-2&0&0&0&4&-2\\\
0&0&0&-2&0&0&-2&4\end{bmatrix}.$

Let us define analogues of the cohomology groups $\mathcal{H}^{r}(T)$ of the
combinatorial torus $C(T)$. The quotient space of the linear space of closed
$r$-forms

$N_{d^{c}}^{r}(T)=\\{\omega\in K^{r}(T):\ d^{c}\omega=0\\}$

modulo the subspace of exact $r$-forms

$R_{d^{c}}^{r}(T)=\\{\omega\in K^{r}(T):\ \exists\varphi\in
K^{r-1}(T)\quad\omega=d^{c}\varphi\\}$

is called the $r$-th cohomology group of $C(T)$, that is,

$\mathcal{H}^{r}(T)=N_{d^{c}}^{r}(T)/R_{d^{c}}^{r}(T).$

Two closed $r$-forms $\omega^{1}$ and $\omega^{2}$ are cohomologous,
$\omega^{1}\sim\omega^{2}$, if and only if they differ by an exact form, i.e.,

$\omega^{1}\sim\omega^{2}\Leftrightarrow\omega^{1}-\omega^{2}\in
R_{d^{c}}^{r}(T).$

An element of $\mathcal{H}^{r}(T)$ is thus an equivalence class $[\omega]$ of
closed $r$-forms $\omega+d^{c}\varphi$, defined by the equivalence relation
$\sim$. These equivalence classes endow $\mathcal{H}^{r}(T)$ with a group
structure.

Calculation of $\mathcal{H}^{0}(T)$. Since there are no $-1$-forms, a $0$-form
$\varphi$ can never be exact, i.e. $R_{d^{c}}^{0}(T)=\\{0\\}$. If $\varphi\in
N_{d^{c}}^{0}(T)$ then $d^{c}\varphi=0$. By (4) it follows immediately that
$\varphi_{1,1}=\varphi_{2,1}=\varphi_{1,2}=\varphi_{2,2}$. Hence

$\varphi=c(x^{1,1}+x^{2,1}+x^{1,2}+x^{2,2}),$

where $c\in\mathbb{R}$ is a constant. Thus $\mathcal{H}^{0}(T)$ is isomorphic
to the group generated by one independent generator which is isomorphic to
$\mathbb{R}$ and we write $\mathcal{H}^{0}(T)\cong\mathbb{R}$.

Calculation of $\mathcal{H}^{1}(T)$. Let $\omega^{1}=\\{u^{1},v^{1}\\}\in
R_{d^{c}}^{1}(T)$. Then $\omega^{1}=d^{c}\varphi$ for some $\varphi\in
K^{0}(T)$. From (4) we have

$\displaystyle u^{1}_{1,1}=\varphi_{2,1}-\varphi_{1,1},\qquad
u^{1}_{1,2}=\varphi_{2,2}-\varphi_{1,2},$ $\displaystyle
u^{1}_{2,1}=\varphi_{1,1}-\varphi_{2,1},\qquad
u^{1}_{2,2}=\varphi_{1,2}-\varphi_{2,2},$ $\displaystyle
v^{1}_{1,1}=\varphi_{1,2}-\varphi_{1,1},\qquad
v^{1}_{2,1}=\varphi_{2,2}-\varphi_{2,1},$ $\displaystyle
v^{1}_{1,2}=\varphi_{1,1}-\varphi_{1,2},\qquad
v^{1}_{2,2}=\varphi_{2,1}-\varphi_{2,2}.$

It follows that

$u^{1}_{2,1}=-u^{1}_{1,1},\qquad u^{1}_{1,2}=-u^{1}_{2,2},\qquad
v^{1}_{1,2}=-v^{1}_{1,1},\qquad v^{1}_{2,1}=-v^{1}_{2,2}.$

Hence any form $\omega^{1}\in R_{d^{c}}^{1}(T)$ can be written as

$\displaystyle\omega^{1}=u^{1}_{1,1}e_{1}^{1,1}-u^{1}_{1,1}e_{1}^{2,1}-u^{1}_{2,2}e_{1}^{1,2}+u^{1}_{2,2}e_{1}^{2,2}$
$\displaystyle+v^{1}_{1,1}e_{2}^{1,1}-v^{1}_{1,1}e_{2}^{1,2}-v^{1}_{2,2}e_{2}^{2,1}+v^{1}_{2,2}e_{2}^{2,2}$
$\displaystyle=u^{1}_{1,1}(e_{1}^{1,1}-e_{1}^{2,1})+u^{1}_{2,2}(e_{1}^{2,2}-e_{1}^{1,2})+v^{1}_{1,1}(e_{2}^{1,1}-e_{2}^{1,2})+v^{1}_{2,2}(e_{2}^{2,2}-e_{2}^{2,1}).$

Let now $\omega=\\{u,v\\}$ be a closed 1-form, i. e., $\omega\in
N_{d^{c}}^{1}(T)$ and $d^{c}\omega=0$. By (4) the requirement $d^{c}\omega=0$
means that

$\displaystyle u_{1,1}-u_{1,2}+v_{2,1}-v_{1,1}=0,$ $\displaystyle
u_{2,1}-u_{2,2}-v_{2,1}+v_{1,1}=0,$ $\displaystyle
u_{1,2}-u_{1,1}+v_{2,2}-v_{1,2}=0,$ $\displaystyle
u_{2,2}-u_{2,1}+v_{1,2}-v_{2,2}=0.$

From this we obtain

$\displaystyle u_{1,1}-u_{1,2}+u_{2,1}-u_{2,2}=0,$ $\displaystyle
v_{1,1}-v_{2,1}+v_{1,2}-v_{2,2}=0.$

Hence any form $\omega\in N_{d^{c}}^{1}(T)$ can be written as

$\displaystyle\omega=u_{1,1}e_{1}^{1,1}+(u_{1,2}-u_{1,1}+u_{2,2})e_{1}^{2,1}+u_{1,2}e_{1}^{1,2}+u_{2,2}e_{1}^{2,2}$
$\displaystyle+v_{1,1}e_{2}^{1,1}+(v_{1,1}+v_{1,2}-v_{2,2})e_{2}^{2,1}+v_{1,2}e_{2}^{1,2}+v_{2,2}e_{2}^{2,2}.$

This yields

$\displaystyle\omega=u_{1,1}(e_{1}^{1,1}-e_{1}^{2,1})+u_{1,2}(e_{1}^{2,1}+e_{1}^{1,2})+u_{2,2}(e_{1}^{2,2}+e_{1}^{2,1})$
$\displaystyle+v_{1,1}(e_{2}^{1,1}+e_{2}^{2,1})+v_{1,2}(e_{2}^{2,1}+e_{2}^{1,2})+v_{2,2}(e_{2}^{2,2}-e_{2}^{2,1})$
$\displaystyle=(u_{1,2}+u_{2,2})(e_{1}^{2,1}+e_{1}^{1,2})+(v_{1,2}+v_{1,1})(e_{2}^{2,1}+e_{2}^{1,2})+\omega^{0},$

where

$\omega^{0}=u_{1,1}(e_{1}^{1,1}-e_{1}^{2,1})+u_{2,2}(e_{1}^{2,2}-e_{1}^{1,2})+v_{1,1}(e_{2}^{1,1}-e_{2}^{1,2})+v_{2,2}(e_{2}^{2,2}-e_{2}^{2,1})$

and note that $\omega^{0}\in R_{d^{c}}^{1}(T)$. Hence

$\omega\sim(u_{1,2}+u_{2,2})(e_{1}^{2,1}+e_{1}^{1,2})+(v_{1,2}+v_{1,1})(e_{2}^{2,1}+e_{2}^{1,2})\in\mathcal{H}^{1}(T).$

This means that $\mathcal{H}^{1}(T)$ has two independent generators. Thus
$\mathcal{H}^{1}(T)\cong\mathbb{R}^{2}$.

Calculation of $\mathcal{H}^{2}(T)$. A 2-form

$\psi^{1}=\psi^{1}_{1,1}V^{1,1}+\psi^{1}_{2,1}V^{2,1}+\psi^{1}_{1,2}V^{1,2}+\psi^{1}_{2,2}V^{2,2}$

is an element of $R_{d^{c}}^{2}(T)$ if $\psi^{1}=d^{c}\omega$ for some 1-form
$\omega=\\{u,v\\}$. Using (4) $\psi^{1}=d^{c}\omega$ gives rise to

$\displaystyle\psi^{1}_{1,1}=u_{1,1}-u_{1,2}+v_{2,1}-v_{1,1},\quad\psi^{1}_{2,1}=u_{2,1}-u_{2,2}-v_{2,1}+v_{1,1},$
$\displaystyle\psi^{1}_{1,2}=u_{1,2}-u_{1,1}+v_{2,2}-v_{1,2},\quad\psi^{1}_{2,2}=u_{2,2}-u_{2,1}+v_{1,2}-v_{2,2}.$

Adding these equations we obtain

$\psi^{1}_{1,1}+\psi^{1}_{2,1}+\psi^{1}_{1,2}+\psi^{1}_{2,2}=0.$

Hence any element $\psi^{1}\in R_{d^{c}}^{2}(T)$ can be written as

$\psi^{1}=\psi^{1}_{1,1}V^{1,1}+\psi^{1}_{2,1}V^{2,1}+\psi^{1}_{1,2}V^{1,2}+(-\psi^{1}_{1,1}-\psi^{1}_{2,1}-\psi^{1}_{1,2})V^{2,2}.$

Since $d^{c}\psi=0$ for any 2-form $\psi$, $N_{d^{c}}^{2}(T)=K^{2}(T)$. Any
element of $N_{d^{c}}^{2}(T)$ can be expressed as

$\psi=(\psi_{1,1}+\psi_{2,1}+\psi_{1,2}+\psi_{2,2})V^{2,2}+\psi^{0},$

where

$\psi^{0}=\psi_{1,1}V^{1,1}+\psi_{2,1}V^{2,1}+\psi_{1,2}V^{1,2}+(-\psi_{1,1}-\psi_{2,1}-\psi_{1,2})V^{2,2}.$

Thus

$\psi\sim(\psi_{1,1}+\psi_{2,1}+\psi_{1,2}+\psi_{2,2})V^{2,2},$

since $\psi^{0}\in R_{d^{c}}^{2}(T)$. This means that $\mathcal{H}^{2}(T)$ has
only one independent generator. So that $\mathcal{H}^{2}(T)\cong\mathbb{R}$.

Thus the cohomology groups are exactly the same as in the continuum case.

## References

  * [1] D.N. Arnold, R.S. Falk, and R. Winther. Finite element exterior calculus, homological techniques, and applications, Acta Numer. 15 (2006), 1-155.
  * [2] D.N. Arnold, R.S. Falk, and R. Winther. Finite element exterior calculus: from Hodge theory to numerical stability, Bull. Amer. Math. Soc. 47(2) (2010), 281-354.
  * [3] R. Ayoub, A. Hamdouni, and D. Razafindralandy. A new Hodge operator in discrete exterior calculus. Application to fluid mechanics. Communications on Pure & Applied Analysis 20(6) (2021), 2155–2185.
  * [4] V. de Beaucé, S. Sen and J.C. Sexton. Chiral dirac fermions on the lattice using geometric discretisation, Nucl. Phys. B (Proc. Suppl.) 129–130 (2004), 468–470.
  * [5] P. Becher and H. Joos. The Dirac-Kähler equation and fermions on the lattice, Z. Phys. C 15 (1982), 343–365.
  * [6] P. Cerejeiras, U. Kähler, F. Sommen , and A. Vajiac. Script Geometry. In: Bernstein S., Kähler U., Sabadini I., Sommen F. (eds) Modern Trends in Hypercomplex Analysis. Trends in Mathematics. Birkhäuser, Cham. (2016), 79–110.
  * [7] A.A. Dezin. Multidimensional analysis and discrete models, CRC Press, Boca Raton, 1995.
  * [8] J. Dodziuk. Finite-difference approach to Hodge theory of harmonic forms, Amer. J. Math. 98 (1976), 79–104.
  * [9] J.M. Rabin. Homology theory of lattice fermion doubling, Nucl. Phys. B 201 (2) (1982), 315–332.
  * [10] I. Nitschke, S.Reuther, and A. Voigt. Discrete Exterior Calculus (DEC) for the Surface Navier-Stokes Equation. In: Bothe D., Reusken A. (eds) Transport Processes at Fluidic Interfaces. Advances in Mathematical Fluid Mechanics. Birkhäuser, Cham. (2017), 177–197.
  * [11] S. Sen, S. Sen, J.C. Sexton, and D. Adams. A geometric discretisation scheme applied to the Abelian Chern-Simons theory, Phys. Rev. E 61 (2000), 3174–3185.
  * [12] A. Stern, P. Leopardi. The abstract Hodge-Dirac operator and its stable discretization, SIAM Journal on Numerical Analysis 54(6) (2016), 3258-–3279.
  * [13] V. Sushch. Difference Poisson equation on a curvilinear mesh. Differ. Equations 32(5) (1996), 693–697; translation from Differ. Uravn. 32(5) (1996), 684–688.
  * [14] V. Sushch. On some finite-difference analogs of invariant first-order hyperbolic systems, Diff. Equ. 35(3) (1999), 414–420; translation from Differ. Uravn. 35(3) (1999), 411–417.
  * [15] V. Sushch. Green function for a two-dimensional discrete Laplace-Beltrami operator. Cubo 10(2) (2008), 47-59.
  * [16] V. Sushch. A discrete model of the Dirac-Kähler equation. Rep. Math. Phys. 73(1) (2014), 109–125.
  * [17] V. Sushch. A discrete Dirac-Kähler equation using a geometric discretisation scheme. Adv. Appl. Clifford Alg. 28:72, 1–17 (2018)
  * [18] V. Sushch. A discrete version of plane wave solutions of the Dirac equation in the Joyce form. Adv. Appl. Clifford Alg. 30, No. 3, Paper No. 46, 20 p. (2020).
  * [19] F.L. Teixeira, W.C. Chew. Lattice electromagnetic theory from a topological viewpoint, J. Math. Phys. 40(1) (1999), 168–187.
  * [20] S. Watterson. The chiral and flavour projection of Dirac-Kähler fermions in the geometric discretization, Int. J. Geom. Methods Mod. Phys. 5(3) (2008), 345-362.
  * [21] F. W. Warner. Foundations of Differentiable Manifolds and Lie Groups, Springer-Verlag, New York, 1983.
  * [22] S.O. Wilson. Differential forms, fluids, and finite models, Proc. of Amer. Math. Soc. 139 (2011), 2597–2604.