# System Support for Environmentally Sustainable Computing in Data Centers

Fan Chen
Luddy School of Informatics, Computing, and Engineering, Indiana University,
Bloomington, IN
E-mail<EMAIL_ADDRESS>

###### Abstract

Modern data centers suffer from a growing carbon footprint due to insufficient
support for environmental sustainability. While hardware accelerators and
renewable energy have been utilized to enhance sustainability, addressing
Quality of Service (QoS) degradation caused by renewable energy supply and
hardware recycling remains challenging: (1) prior accelerators exhibit
significant carbon footprints due to limited reconfigurability and inability
to adapt to renewable energy fluctuations; (2) integrating recycled NAND flash
chips in data centers poses challenges due to their short lifetime, increasing
energy consumption; (3) the absence of a sustainability estimator impedes data
centers and users in evaluating and improving their environmental impact. This
study aims to improve system support for environmentally sustainable data
centers by proposing a reconfigurable hardware accelerator for intensive
computing primitives and developing a fractional NAND flash cell to extend the
lifetime of recycled flash chips while supporting graceful capacity
degradation. We also introduce a sustainability estimator to evaluate user
task energy consumption and promote sustainable practices. We present our
preliminary results and recognize this as an ongoing initiative with
significant potential to advance environmentally sustainable computing in data
centers and stimulate further exploration in this critical research domain.

###### Index Terms:

Environmental sustainability, embodied carbon, operational carbon, hardware
accelerators, storage systems.

##  I Introduction

Modern data centers contribute significantly to an unsustainable carbon
footprint, from both operational and embodied sources. Operational emissions
arise from routine hardware usage, while embodied emissions result from
manufacturing and construction processes. Global data center operational
energy is projected to reach 4.5% of global demand by 2025 [1]. Additionally,
Over the past decade, global data center infrastructure capacity, reflecting
embodied carbon footprint, has surged $6\times$ [2]. Renewable energy and
hardware recycling are two critical methods to enhance the environmental
sustainability of future data centers. Several cloud companies [3] have
collectively invested in 22 GW of renewable energy generation to reduce the
operational carbon footprint. Meanwhile, recent work considers reusing
obsolete hardware, e.g., mobile CPUs [4], mobile DRAMs [4], and FPGAs [5], to
reduce the embodied carbon footprint in data centers. However, mitigating
Quality of Service (QoS) degradation in future sustainable data centers caused
by renewable energy supply and hardware recycling poses a complex challenge.

On one side, renewable energy supply significantly degrades QoS due to its
intermittent nature. Solar and wind generate 47% and 34% of the renewable
energy [6], leading to hourly and seasonal fluctuations in energy supply.
While large energy storage [3, 7] can mitigate these fluctuations, their
addition significantly increases the embodied carbon footprint of data
centers. Consequently, servers may need to be slowed down or shut down during
low energy supply periods. Software-based techniques such as carbon-aware
scheduling [3] and virtual transient servers [7] are proposed to schedule non-
critical tasks among CPUs without significant QoS degradation in the face of
energy supply fluctuations. However, emerging workloads such as big data, deep
learning, and cryptography heavily use accelerators [8, 9] and large data
storage [10]. There is a near total lack of hardware (particularly accelerator
and storage) support for maintaining the QoS in data centers powered by
renewable energy supply.

On the other side, while hardware recycling reduces the embodied carbon
footprint, it compromises data center QoS due to the limited performance and
reliability of obsolete hardware. Recent attempts to integrate outdated
computing components like mobile processors [4] and DRAMs [4] fall short in
handling demanding tasks in emerging workloads. Thus, future data centers will
require reconfigurable application-specific hardware acceleration with low
embodied carbon footprints. Although few works [11] explore reusing obsolete
NAND flash chips in data center storage systems, utilizing aging flash chips
with about-to-worn-out blocks poses challenges in supporting long-term
workloads in future data centers.

Below, we provide a detailed overview of these challenges across the computing
and storage layers of system stacks:

  * •

Hardware Acceleration: Previous GPU-, FPGA- [5, 12], ASIC- [8] and PIM-based
[9] accelerators have been utilized to mitigate performance loss caused by
hardware recycling. However, none effectively addresses energy supply
fluctuations or maintains a low embodied carbon footprint. Specifically, GPU-
and application-specific-integrated-circuit (ASIC)-based accelerators lack the
ability for rapid wake-up or shutdown in response to energy supply
fluctuations. While FPGA- [13, 14] and PIM-based [9, 15] accelerators adopt
nonvolatile memory arrays for data retention, critical components such as
SRAM-based switches in FPGAs and analog-to-digital converters (ADCs) in PIMs
still lose information without power, hindering prompt workload resumption
after power loss. ASIC-based accelerators [8] focusing on accelerating only
one type of computing kernels suffer from high embodied carbon footprint, due
to their lack of reconfigurability.

  * •

NAND Flash: To reduce embodied carbon footprint, old NAND flash chips will be
adopted by future sustainable data centers. However, flash suffers from short
cell endurance [16], where about-to-worn-out blocks limit chip lifetime.
Despite error correcting code (ECC) [17, 18] and wear leveling [19]
assistance, no prior technique significantly extends the lifetime of about-to-
worn-out blocks, hindering widespread adoption of recycled NAND flash chips in
future data centers. Moreover, due to fluctuations in renewable energy supply,
DRAM snapshots are frequently stored in and loaded from storage systems.
Unfortunately, most recycled flash chips use multi-level cells (MLCs) [20],
suffering from long program and read operations.

  * •

Sustainability Estimator: Without an environmental sustainability estimator,
data centers struggle to assess the operational and embodied energy of user
workloads, leading to a lack of visibility regarding the sustainability costs
of user actions. Additionally, the absence of such metrics inhibits the
development of flexible billing policies that could reward users who opt for
recycled hardware and encourage broader adoption of environmentally friendly
practices.

##  II Approaches

This project aims to develop system support for environmentally sustainable
data centers, with three main objectives: (1) Facilitate energy-efficient
hardware acceleration for memory-intensive tasks such as big data, deep
learning, and cryptography, ensuring steady progress even in the face of
renewable energy fluctuations; (2) Extend the lifetime of about-to-worn-out
blocks in recycled NAND flash chips, thereby mitigating the embodied energy
consumption; and (3) Provide an accurate evaluation of both operational and
embodied energy consumption for data center user tasks and incentivize users
to make their computing more environmentally sustainable. The primary approach
entails constructing an environmentally sustainable computing system. This
system encompasses a nonvolatile reconfigurable PIM accelerator, a storage
system based on recycled NAND flash, and an accurate environmental
sustainability estimator.

###  II-A Amoeba - A Reconfigurable Nonvolatile Accelerator

Existing PIMs in data centers encounter substantial manufacturing costs and
operational inefficiencies, particularly when faced with intermittent
renewable energy sources, resulting in significant embodied and operational
carbon emissions. We propose Amoeba, depicted in Figure 1, as a reconfigurable
FeFET-based PIM architecture aimed at addressing these challenges. Our
approach involves profiling compute- and resource-intensive computational
kernels in data centers, including NTT for lattice-based cryptography, SHA3
for blockchains, and convolution for machine learning, leveraging insights
from prior studies. Categorizing basic operations and supporting computation
primitives with three types of processing engines (PE) based on the FeFET
crossbar structure, Amoeba offers fine-grained PE-level reconfiguration to
markedly minimize embodied carbon emissions. Furthermore, serving as a fully
nonvolatile accelerator, it ensures consistent forward progress even amidst
renewable energy fluctuations, thereby contributing to a notable reduction in
the operational carbon footprint.

Figure 1: The Amoeba architecture.

PE-Level Reconfigurability for Reducing Embodied Carbon Footprint. Crossbars
in Amoeba can be logically configured into different modes: (1) Associative PE
(APE) for LUT and ADD. Each crossbar acts as an embedded CAM block for
associative computing. As in previous work [21], each unit consists of
complementary cells (H and L) for 1’ and 0’, respectively. The search word is
encoded by complementary voltages (V/0 for 1’ and 0/V for 0’) and broadcasted
onto vertical search lines. After a parallel search, the horizontal match line
identifies matches and outputs the vector for post-processing. LUT is directly
implemented. For ADD, bitwise search-based additions are cascaded for multi-
bit addition, allowing for parallel execution. (2) Multiplicaiton PE (MPE) for
MVM and SHIFT. The crossbar-based MVM unit, as in [22, 23], maps dataflow in a
weight-stationary manner. Weights are pre-coded onto columns, and input
feature maps are converted to read voltages and applied to horizontal
wordlines. The ADC aggregates current at vertical bitlines, representing dot-
products. SHIFT operations comprise $>$40% of NTT workloads, but they aren’t
inherently MVM operations and can’t be directly mapped to crossbars. We
propose transferring SHIFT to MVM by pre-coding the transformation matrix onto
the crossbar and applying the input as a multiplicand. This scheme, based on
linear algebra, can be generalized to support any cyclic permutation. (3)
Computing PE (CPE) for LOGIC operations. We perform Logic operations like AND
and XOR within the crossbar array without using CMOS ALUs. In contrast to
prior methods [24, 25], which necessitate a larger crossbar array, our
approach only requires a 2$\times$N crossbar, where $N$ represents the input
width. Take ADD as an example, we adjust ADC sensing levels to output a
logical ’1’ only when both cells are in a low-resistance state. (4) Combining
APE and MPE for MUL. Prior algorithms [26, 27] broke an N-bit MUL to SHIFT and
ADD of partial products. [28] implemented such MUL with bitwise operations
based ADD and implicit data selection based SHIFT. Instead, we combine the APE
and MPE to implement MUL. The benefits are twofold: first, the search-based
ADD reduces expensive read/write operations in [28]; second, the MVM based
explicit SHIFT reduces the complex control for memory block selection in [28].

Figure 2: A Fraction NAND flash cell (FRAC)

Fully Nonvolatile Accelerator for Reducing Operational Carbon Footprint. To
reduce the operational carbon footprint, we propose a fully nonvolatile
accelerator with FeFET-based True Random Number Generator (TRG) and ultra-low
power FeFET-based ADC. (1) FeFET-based TRG. We leverage the stochastic
switching behavior in scaled FeFET devices [29] as an on-chip entropy source
for TRG designs, same as prior work [30]. Our baseline TRG employs a
0.2$\mu{m}$$\times$0.2$\mu{m}$ FeFET device with a read voltage set at 0.2V.
Initial findings indicate a bias towards ‘0’s in the output bits. To mitigate
this bias, we propose a tracking scheme involving an 8-bit counter to record
output probabilities in consecutive 256-bit segments. The counter’s output
then acts as the control signal for adjusting write voltages in generating
subsequent segments. (2) Low-Power ADC with Tunable Precision. We develop a
precision-scalable, low-power ADC using FeFET technology for faster analog-to-
digital conversion. Figure 3(d) illustrates a four-bit ADC. We utilize four
partially polarized 1$\mu{m}$$\times$1$\mu{m}$ FeFET devices (①$\sim$④). The
write line toggles between 0$V$ and $V_{DD}$ for enabling/disabling. Inputs
are applied to the data line. Due to variations in $V_{th}$, read currents
differ. For instance, with a 0.9V input in Figure 3(d), higher-than-threshold
current is detected on ① and ②, while lower current is observed on ③ and ④.
This results in a 4-bit output ”1100” with an appropriately set current
threshold. The ADC’s precision can be easily adjusted by disabling ① and ③.
Threshold currents for each device can be dynamically programmed.

Figure 3: A fully nonvolatile accelerator. Figure 4: An environmental
sustainability estimator (ESE).

###  II-B FRAC - Recycling Used NAND Flash Chips

To reduce embodied energy, repurposing outdated NAND flash chips in
sustainable data centers is crucial. These chips, manufactured using older
process technologies, may have undergone varying levels of write operations
before recycling, leading to differing remaining lifetimes. While existing
work [31, 32, 33] salvages already-worn-out blocks, the remaining lifetime of
recycled chips is predominantly influenced by about-to-worn-out blocks. We
propose FRAC to improve the lifetime of about-to-worn-out blocks by exploring
the trade-off between chip capacity and cell lifetime. FRAC supports graceful
storage capacity degradation by gradually decreasing the number of $V_{th}$
states in a cell.

FRAC: A NAND Flash Fraction Cell. A conventional flash cell (e.g., SLC, MLC,
TLC, or QLC) uses $2^{n}$ $V_{th}$ states to store $n$ bits respectively, with
$n\in[1,4]$.Figure 2(a) illustrates a TLC cell. In contrast, FRAC uses $m$
$V_{th}$ states to store data, where $m\in[2,2^{n}]$, e.g., $m=3$. However, a
fraction cell having 3 $V_{th}$ states can represent only 1 bit, wasting 1
$V_{th}$ state. To use more $V_{th}$ states in a cell, we propose to use
$\alpha$ FRAC cells, each of which has $m$ $V_{th}$ states, to store
$\lfloor\log_{2}(m^{\alpha})\rfloor$ bits. For instance, two 3-state cells can
store 3 bits. The truth table of two 3-state cells is shown in Figure 2(b). By
looking up this table, the value of two 3-state cells is translated to 3 bits.
We define cell utilization as the number of $V_{th}$ states representing data
divided by the total number of $V_{th}$ states of FRAC cells. We explore the
cell utilization of 3-, 5-, and 7-state cells in Figure 2(c). To achieve the
highest cell utilization, we can store 11 bits in seven 3-state cells, 16 bits
in ten 5-state cells, and 16 bits in five 7-state cells.

Tradeoff Between Capacity and Lifetime. The endurance of a flash cell ($L$)
has a power-dependence on the number of program/erase (P/E) cyclings
($N_{PE}$) [34, 35], i.e., $L\propto N_{PE}^{\beta}$, where $\beta\geq 0.3$.
Producing more $V_{th}$ states in a flash cell requires more program
operations [36], thereby greatly decreasing the cell endurance. Compared to a
3-bit TLC, a 2-state cell shown in Figure 2(a) prolongs the cell endurance by
$10\times$ [37]. Although prior work [38] improves the flash lifetime by
converting 2-bit MLC to 1-bit SLC, the 50% capacity loss does not maintain a
gradual and graceful capacity degradation. In contrast, FRAC can gracefully
explore the tradeoff between the chip endurance and capacity, as shown in
Figure 2(d). By gradually reducing the number of $V_{th}$ states from 8 to 2,
FRAC gracefully degrades the capacity of a page from 4KB to 1.3KB and prolongs
the cell endurance from $1\times$ to $10\times$.

Read and Write. To sense data from a TLC, as Figure 2(e) exhibits, totally
$\log_{2}(8)=3$ iterations are required. The reading reference of this
iteration is decided by the result of the previous iteration. For example,
after comparing against $r_{3}$, if the $V_{th}$ state of the cell is larger
than $r_{3}$, $r_{5}$ will be selected in the next iteration. Otherwise,
$r_{1}$ will be chosen. The same as conventional MLC, TLC, and QLC, a read
operation on a $m$-state faction cell also requires $\lceil\log_{2}(m)\rceil$
iterations. After the data are sensed out of the cells in a page, the data in
multiple cells is translated to one value. After an erase operation, to
program a TLC, the incremental step pulse scheme [39] issues multiple
programming pulses, each of which is stronger than the previous one, as shown
in Figure 2(f). Between two programming pulses, there is a verify that checks
whether the cell reaches the target $V_{th}$ or not. On the contrary, to
program an $m$-state cell, instead of a small pulse, the incremental step
pulse scheme can directly start with a larger pulse. In this way, FRAC
requires fewer pulses and thus prolongs the cell endurance, when $m<8$.

###  II-C ESE - An Environmental Sustainability Estimator

We introduce the Environmental Sustainability Estimator (ESE), an innovative
tool facilitating accurate assessment of CO2 emissions linked with user
activities. ESE integrates environmental sustainability into data center
billing policies, incentivizing users to adopt more sustainable computing
practices. Illustrated in Figure 4(a), ESE comprises a hardware estimator, a
data center energy model, and an energy source predictor. First, the hardware
estimator analyzes user source code and expected Quality of Service (latency),
partitions tasks, and distributes them across CPUs, GPUs, and PIMs to ensure
latency requirements are met. It also assesses storage capacity and network
bandwidth based on task data access patterns. Second, the data center energy
model will use the hardware estimator output to produce both operational
energy consumption ($E_{ope}$) and embodied energy consumption ($E_{emb}$).
Third, the user inputs the task’s start time to access weather forecasts. The
energy source predictor utilizes this information along with calendar data to
determine the net energy demand, representing the disparity between total data
center operational energy and current renewable energy generation. Finally,
based on the values of $E_{ope}$, $E_{emb}$, and net energy demand, the data
center uses different billing policies to decide the user charge.

Hardware estimator. The hardware estimator minimizes $E_{ope}$ and $E_{emb}$
for each user task within a latency constraint on a heterogeneous hardware
platform. It takes source code and expected latency as inputs. The source code
is compiled into binaries for an initial hardware configuration, and static
features are extracted directly. Runtime features are collected by executing
the binaries on real hardware. A latency model processes both types of
features to estimate CPU/GPU/PIM/net/SSD latency. The task’s estimated latency
is the sum of these values. If the estimated latency exceeds the expected
latency, kernels are moved to GPUs or PIMs, and the estimator iterates until
the expected latency is achieved. Otherwise, the minimum estimated latency is
returned.

Traditional task partition methods [40, 41] rely on static features such as
branch count and arithmetic operations, extracted from source code using
intermediate representations (IRs). However, these handcrafted features often
overlook inter-kernel data flows, treating each kernel independently. Our
hardware estimator employs a transformer architecture [42] to directly extract
more informative static features from the source code, capturing interactions
among all kernels. Leveraging a pre-trained transformer [43] tailored for code
generation eliminates the need for training from scratch. Runtime features
encompass CPU cache misses, GPU core utilization, PIM computing unit
utilization, and host-device data transfers, among others. Measured using
performance profiling tools like Intel VTune and NVIDIA visual profiler, these
features are assessed after deploying compiled binaries on real hardware.

To expedite latency evaluation, we train a CNN as a latency model utilizing
static and runtime features alongside hardware configuration to predict
latency. We create a latency dataset by generating varied partition results
and measuring latency values on CPUs, GPUs, and PIMs for training. Post-
training, partition results won’t require execution on real hardware; instead,
latency models will directly provide values.

Data Center Energy Model. The energy model comprises an operational and an
embodied energy model. The energy model inputs are the static and runtime
features, hardware configuration, and hardware latency values, while its
outputs are operational energy and embodied energy. We will train a CNN as the
operational energy model. Similar to latency models in the hardware estimator,
we will randomly generate different partition results and measure the energy
consumption on CPUs, GPUs, PIMs, networks, and storage systems to build an
energy dataset for training the energy model. Besides the dynamic and leakage
energy consumption of various hardware units, we will also consider the
cooling energy cost, idle equipment energy, power delivery loss, and other
operation-related energy in the operational energy model. The embodied energy
model is a linear equation. We will compute the total embodied energy of a
hardware unit by considering the energy consumption during production and
manufacture, transport, use and maintenance, recycling, and other stages in
the whole hardware lifetime. The embodied energy ($E_{emb}$) of a user task
can be computed as $\sum_{i\in
X}TBE_{i}\times\frac{latency_{i}}{lifetime_{i}}$, where $X$ indicates the
hardware units used by the task; $TBE_{i}$ is the total embodied energy of the
hardware unit $i$; $latency_{i}$ means the latency of the hardware unit $i$;
and $lifetime_{i}$ represents the lifetime of the hardware unit $i$.

Figure 5: Carbon footprints and progresses under energy fluctuations compared
to prior designs.

Figure 6: The Raw Bit Error Rate (RBER) of FRAC-based pages in a recycled
flash chip.

Figure 7: The prediction on the wind energy generation in the California grid.

Energy Source Predictor. The energy source predictor aims to characterize the
distribution of net energy demand forecast errors under various weather
conditions and data center workload capacities, as shown in Figure 4(d). The
energy source predictor will be built as a long short-term memory (LSTM)
network to output simultaneous quantile forecasts of net energy demand and
generated renewable energy. The network inputs include an array of calendar
data and weather information. The network will be trained on historical
records and predicts the forecasts and forecast errors of $T_{0}+5$-minute,
$T_{0}+10$-minute, and $T_{0}+15$-minute at $T_{0}$. After training, the
network will ingest near-past 5-minute forecasts as well as near-past
10-minute and near-past 15-minute historical data available at $T_{0}$. More
specifically, the network predicts the conditional quantiles of forecast
errors for net energy demand and renewable energy production. We will consider
seven target quantiles, i.e., P2.5, P5, P25, P50, P75, P95, and P97.5.

##  III Results

We offer preliminary results of our ongoing research. These findings highlight
the potential of our work to advance environmentally sustainable computing in
data centers.

Amoeba Results. We evaluated Amoebaby implementing a simple prototype. We
compare Amoeba against prior accelerators designs, i.e., FPGA [44], CMOS [45],
RRAM [46], and FeFET [47]. We use three benchmarks including a 32k NTT using
Montgomery reduction with a fixed $q$ =12289 [28], SHA3 with a 1088 blocksize
and 1600-bit state sizes, and AlexNet [45]. For embodied carbon emissions, we
mainly focused on manufacturing costs and used the total cost for three
different domain accelerators as the overall footprint. We consider California
grid [48] as the renewable power supply and adopt its historical data, taking
into account dynamic intermittency and fluctuations. We adopt the recently
proposed Carbon Explorer framework [48] to study the Pareto Frontier analysis
and report the design solution for each scheme in Figure 7(left). In general,
the Amoeba achieves state-of-the-art carbon minimization. Figure 7(right)
compares the breakdown of simulation run times and computational progress for
AlexNet inference at different times of a week under California grid [48]
supply. For all architectures, Amoeba achieves the highest forward progress
since the fully nonvolatile accelerator can work below threshold power (i.e.,
below Thld) with reduced operational energy consumption. In addition, we can
also see that the fluctuation of the power supply imposes large rollover
penalties for nonvolatile processors and CMOS circuitries in existing RRAM and
FeFET accelerators.

FRAC Results. We evaluated the viability of FRAC using a prototype built upon
a Zynq FPGA [49] with a NAND flash daughter broad. Recycled flash chips can be
plugged into the flash daughter board, while we implement the control logic
for FRAC read, program, and erase operations on the FPGA that can access those
chips. Our prototype is a minimal storage system built upon recycled flash
chips and can be used to study the lifetime improvement achieved by FRAC. Our
prototype primarily translates high-level read/write I/O requests to device-
level NAND flash commands. We tested the RBER of FRAC cells with different
numbers of $V_{th}$ states using our prototype. Figure 7 highlights the RBER
of FRAC-based pages subjected to 6k programming and erasure cycles in an aged
flash chip. A higher RBER indicates a page nearing the end of its lifetime
[10]. When the cells in the recycled flash chip have only two $V_{th}$ states,
the RBER is only $0.6\%$ on average. In contrast, three $V_{th}$ states in a
cell increase the RBER to $0.9\%$ on average. Four $V_{th}$-state-cells
further enlarge the RBER of the same page to 1.4% on average. These data
suggest that a full-fledged FRAC design can fully explore the trade-off
between cell endurance and chip capacity.

ESE Results. We evaluated ESE using a simplified prototype of the energy
source predictor. Our prototype predicts the average value of renewable energy
generated by wind in the California grid every 30 minutes, utilizing a 2-layer
LSTM network with forget, input, and output gates. Historical wind energy
generation data from the California grid served as ground truth, sourced from
California ISO [50], while wind data from Apple Weather [51] provided input.
Training utilized 70% of the data, with 10% for validation and the remaining
20% for testing. Figure 7 illustrates the average wind energy generation
prediction. While our prototype’s predictions capture the general trend of the
ground truth data, the 30-minute prediction interval results in larger
fluctuations compared to the ground truth, suggesting the need for shorter
intervals (e.g., $5\sim 15$ minutes). These findings imply that a
comprehensive ESE could accurately evaluate operational and embodied energy
consumption in data centers.

##  IV Conclusion

Modern data centers face significant environmental sustainability challenges
throughout their lifecycle. The project aims to tackle these challenges by
developing novel technologies including a reconfigurable hardware accelerator
architecture, a fraction NAND flash cell and its system integration, and an
environmental sustainability estimator. These advancements lay a solid
foundation for improving environmental sustainability in data centers,
enabling both users and data centers to conduct data-intensive applications in
a more sustainable and scalable manner. This not only enhances computing
capabilities but also maximizes the societal benefits derived from computing.

## References

  * [1] Y. Liu et al., “Energy consumption and emission mitigation prediction based on data center traffic and pue for global data centers,” Global Energy Interconnection, vol. 3, no. 3, pp. 272–282, 2020.
  * [2] E. Masanet et al., “Recalibrating global data center energy-use estimates,” Science, vol. 367, no. 6481, pp. 984–986, 2020.
  * [3] B. Acun et al., “Carbon explorer: A holistic framework for designing carbon aware datacenters,” in ASPLOS, pp. 118–132, 2023.
  * [4] J. Switzer et al., “Junkyard computing: Repurposing discarded smartphones to minimize carbon,” in ASPLOS, pp. 400––412, 2023.
  * [5] H. Dogan et al., “Aging analysis for recycled fpga detection,” in IEEE International Symposium on Defect and Fault Tolerance in VLSI and Nanotechnology Systems, pp. 171–176, 2014.
  * [6] S. Nalley et al., “Annual energy outlook 2022,” Energy Information Agency, p. 23, 2022.
  * [7] R. Singh et al., “Yank: Enabling green data centers to pull the plug,” in 10th USENIX Conference on Networked Systems Design and Implementation, (USA), pp. 143––156, USENIX Association, 2013.
  * [8] Y. Chen et al., “Dadiannao: A machine-learning supercomputer,” in International Symposium on Microarchitecture, 2014.
  * [9] A. Shafiee et al., “Isaac: A convolutional neural network accelerator with in-situ analog arithmetic in crossbars,” in ACM/IEEE International Symposium on Computer Architecture, pp. 14––26, 2016.
  * [10] L. M. Grupp et al., “The bleak future of nand flash memory,” in USENIX Conference on File and Storage Technologies, 2012.
  * [11] S. Boyd et al., “Life-cycle assessment of nand flash memory,” IEEE Transactions on Semiconductor Manufacturing, 2011.
  * [12] Y. Isaka et al., “Systematic unsupervised recycled field-programmable gate array detection,” IEEE Transactions on Device and Materials Reliability, vol. 22, no. 2, pp. 154–163, 2022.
  * [13] Y. Y. Liauw et al., “Nonvolatile 3d-fpga with monolithically stacked rram-based configuration memory,” in IEEE International Solid-State Circuits Conference, pp. 406–408, 2012.
  * [14] Y. Chen et al., “3d-nonfar: Three-dimensional non-volatile fpga architecture using phase change memory,” in ACM/IEEE International Symposium on Low Power Electronics and Design, pp. 55––60, 2010.
  * [15] K. Qiu et al., “Resirca: A resilient energy harvesting reram crossbar-based accelerator for intelligent embedded processors,” in IEEE International Symposium on High Performance Computer Architecture, pp. 315–327, 2020.
  * [16] Y. Cai et al., “Error patterns in mlc nand flash memory: Measurement, characterization, and analysis,” in Design, Automation & Test in Europe Conference & Exhibition, pp. 521–526, 2012.
  * [17] T. Nakamura et al., “Adaptive artificial neural network-coupled ldpc ecc as universal solution for 3-d and 2-d, charge-trap and floating-gate nand flash memories,” IEEE Journal of Solid-State Circuits, vol. 54, no. 3, pp. 745–754, 2019.
  * [18] T. Tanzawa et al., “A compact on-chip ecc for low cost flash memories,” IEEE Journal of Solid-State Circuits, vol. 32, no. 5, pp. 662–669, 1997.
  * [19] X. Jimenez et al., “Wear unleveling: Improving NAND flash lifetime by balancing page endurance,” in USENIX Conference on File and Storage Technologies, pp. 47–59, 2014.
  * [20] Y. Cai et al., “Program interference in mlc nand flash memory: Characterization, modeling, and mitigation,” in IEEE International Conference on Computer Design, pp. 123–130, 2013.
  * [21] F. Chen et al., “Parc: A processing-in-cam architecture for genomic long read pairwise alignment using reram,” in IEEE Asia and South Pacific Design Automation Conference, pp. 175–180, 2020.
  * [22] F. Chen et al., “Regan: A pipelined reram-based accelerator for generative adversarial networks,” in IEEE Asia and South Pacific Design Automation Conference, pp. 178–183, 2018.
  * [23] F. Chen et al., “Zara: A novel zero-free dataflow accelerator for generative adversarial networks in 3d reram,” in Design Automation Conference, p. 133, 2019.
  * [24] Y. Zha et al., “Reconfigurable in-memory computing with resistive memory crossbar,” in International Conference on Computer-Aided Design, 2016.
  * [25] Y. Zha et al., “Liquid silicon: A data-centric reconfigurable architecture enabled by rram technology,” in ACM/SIGDA International Symposium on Field-Programmable Gate Arrays, 2018.
  * [26] A. Haj-Ali et al., “Efficient algorithms for in-memory fixed point multiplication using magic,” in IEEE International Symposium on Circuits and Systems, 2018.
  * [27] A. Haj-Ali et al., “Imaging: In-memory algorithms for image processing,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 65, no. 12, pp. 4258–4271, 2018.
  * [28] H. Nejatollahi et al., “Cryptopim: In-memory acceleration for lattice-based cryptographic hardware,” in Design Automation Conference, 2020.
  * [29] S. Deng et al., “A comprehensive model for ferroelectric fet capturing the key behaviors: Scalability, variation, stochasticity, and accumulation,” in IEEE symposium on VLSI technology, pp. 1–2, 2020.
  * [30] F. Chen, “Puffin: An efficient dnn training accelerator for direct feedback alignment in fefet,” in International Symposium on Low Power Electronics and Design, pp. 1–6, 2021.
  * [31] C. Wang et al., “Extending the lifetime of nand flash memory by salvaging bad blocks,” in Design, Automation & Test in Europe Conference & Exhibition, pp. 260–263, 2012.
  * [32] J.-W. Hsieh et al., “Implementation strategy for downgraded flash-memory storage devices,” ACM Transaction Embedded Computing System, vol. 12, mar 2013.
  * [33] H.-Y. Lin et al., “Revive bad flash-memory pages by hlc scheme,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 38, no. 5, pp. 860–873, 2019.
  * [34] H. Yang et al., “Reliability issues and models of sub-90nm nand flash memory cells,” in International Conference on Solid-State and Integrated Circuit Technology, pp. 760–762, 2006.
  * [35] V. Mohan et al., “How i learned to stop worrying and love flash endurance,” in USENIX Conference on Hot Topics in Storage and File Systems, p. 3, 2010.
  * [36] D.-H. Kim et al., “A 1tb 4b/cell nand flash memory with tprog=2ms, tr=110µs and 1.2gb/s high-speed io rate,” in IEEE International Solid- State Circuits Conference, pp. 218–220, 2020.
  * [37] B. Gu et al., “Dynamic file cache optimization for hybrid ssds with high-density and low-cost flash memory,” in IEEE 39th International Conference on Computer Design, pp. 170–173, 2021.
  * [38] X. Jimenez et al., “Phoenix: Reviving mlc blocks as slc to extend nand flash devices lifetime,” in Design, Automation & Test in Europe Conference & Exhibition, pp. 226–229, 2013.
  * [39] K.-D. Suh et al., “A 3.3 v 32 mb nand flash memory with incremental step pulse programming scheme,” IEEE Journal of Solid-State Circuits, vol. 30, no. 11, pp. 1149–1156, 1995.
  * [40] Z. Wang et al., “Partitioning streaming parallelism for multi-cores: a machine learning based approach,” in IEEE International Conference on Parallel Architectures and Compilation Techniques, pp. 307–318, 2010.
  * [41] Y. Wen et al., “Smart multi-task scheduling for opencl programs on cpu/gpu heterogeneous platforms,” in IEEE International conference on high performance computing, pp. 1–10, 2014.
  * [42] A. Svyatkovskiy et al., “Intellicode compose: Code generation using transformer,” in ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering, 2020.
  * [43] Z. Feng et al., “CodeBERT: A pre-trained model for programming and natural languages,” in Findings of the Association for Computational Linguistics: EMNLP 2020, 2020.
  * [44] C. Zhang et al., “Optimizing fpga-based accelerator design for deep convolutional neural networks,” in ACM/SIGDA international symposium on field-programmable gate arrays, 2015.
  * [45] Y.-H. Chen et al., “Eyeriss: A spatial architecture for energy-efficient dataflow for convolutional neural networks,” in ACM/IEEE International Symposium on Computer Architecture, 2016.
  * [46] F. Wang et al., “Reram-based processing-in-memory architecture for blockchain platforms,” in IEEE Asia and South Pacific Design Automation Conference, 2019.
  * [47] R. Agrawal et al., “Security and privacy of blockchain-based single-bit cache memory architecture for iot systems,” IEEE Access, 2022.
  * [48] B. Acun et al., “Carbon explorer: A holistic framework for designing carbon aware datacenters,” in ACM International Conference on Architectural Support for Programming Languages and Operating Systems, 2023.
  * [49] eeNews Europe, “Fpga-based nvme data storage solution.” https://www.eenewseurope.com/en/fpga-based-nvme-data-storage-solution/, 2019.
  * [50] C. I. S. Operator, “California iso open access same-time information system.” http://oasis.caiso.com/mrioasis/logon.do, 2021.
  * [51] D. Sky, “Dark sky and apple weather.” https://darksky.net/apple-weather, 2021.


Present address: ]Walmart Global Tech, Sunnyvale, California 94086, USA
Present address: ]Institute of Physics, Chinese Academy of Sciences, Beijing
100190, China

# Gapless spinons and a field-induced soliton gap in the hyper-honeycomb Cu
oxalate framework compound [(C2H5)3NH]2Cu2(C2O4)3

C. Dissanayake  Department of Physics, University of Central Florida, Orlando,
Florida 32816, USA  A. C. Jacko  School of Mathematics and Physics, The
University of Queensland, Brisbane, Queensland 4072, Australia  K.
Kumarasinghe  R. Munir  H. Siddiquee  Department of Physics, University of
Central Florida, Orlando, Florida 32816, USA  W. J. Newsome  Department of
Chemistry, University of Central Florida, Orlando, Florida 32816, USA  F. J.
Uribe-Romo  Department of Chemistry, University of Central Florida, Orlando,
Florida 32816, USA  REACT: Renewable Energy and Chemical Transformations
Cluster, University of Central Florida, Orlando, Florida 32816, USA  E. S.
Choi  National High Magnetic Field Laboratory, Florida State University,
Tallahassee, Florida 32310, USA  S. Yadav  [  Department of Physics,
University of Florida, Gainesville, Florida 32611, USA  X.-Z. Hu  [
Department of Physics, University of Florida, Gainesville, Florida 32611, USA
Y. Takano  Department of Physics, University of Florida, Gainesville, Florida
32611, USA  S. Pakhira  Ames National Laboratory, Iowa State University, Ames,
Iowa 50011, USA  Institute for Quantum Materials and Technologies, Karlsruhe
Institute of Technology, D-76021 Karlsruhe, Germany  D. C. Johnston  Q.-P.
Ding  Y. Furukawa  Ames National Laboratory, Iowa State University, Ames, Iowa
50011, USA  Department of Physics and Astronomy, Iowa State University, Ames,
Iowa 50011, USA  B. J. Powell  School of Mathematics and Physics, The
University of Queensland, Brisbane, Queensland 4072, Australia  Y. Nakajima
<EMAIL_ADDRESS>Department of Physics, University of Central
Florida, Orlando, Florida 32816, USA

###### Abstract

We report a detailed study of the specific heat and magnetic susceptibility of
single crystals of a spin liquid candidate: the hyper-honeycomb Cu oxalate
framework compound [(C2H5)3NH]2Cu2(C2O4)3. The specific heat shows no anomaly
associated with a magnetic transition at low temperatures down to $T\sim$ 180
mK in zero magnetic field. We observe a large linear-in-$T$ contribution to
the specific heat $\gamma T$, $\gamma=98(1)$ mK/mol K2, at low temperatures,
indicative of the presence of fermionic excitations despite the Mott
insulating state. The low-$T$ specific heat is strongly suppressed by applied
magnetic fields $H$, which induce an energy gap, $\Delta(H)$, in the spin-
excitation spectrum. We use the four-component relativistic density-functional
theory (DFT) to calculate the magnetic interactions, including the
Dzyaloshinskii-Moriya antisymmetric exchange, which causes an effective
staggered field acting on one copper sublattice. The magnitude and field
dependence of the field-induced gap, $\Delta(H)\propto H^{2/3}$, are
accurately predicted by the soliton mass calculated from the sine-Gordon model
of weakly coupled antiferromagnetic Heisenberg chains with all parameters
determined by our DFT calculations. Thus our experiment and calculations are
entirely consistent with a model of [(C2H5)3NH]2Cu2(C2O4)3 in which
anisotropic magnetic exchange interactions due to Jahn-Teller distortion cause
one copper sublattice to dimerize, leaving a second sublattice of weakly
coupled antiferromagnetic chains. We also show that this model quantitatively
accounts for the measured temperature-dependent magnetic susceptibility. Thus
[(C2H5)3NH]2Cu2(C2O4)3 is a canonical example of a one-dimensional spin-1/2
Heisenberg antiferromagnet and not a resonating-valence-bond quantum spin
liquid, as previously proposed.

###### pacs:

##  I Introduction

Frustrated quantum magnets have attracted great interest because of their rich
ground states [1, 2]. For instance, one of the most interesting ground states
is the Kitaev spin liquid, which harbors itinerant Majorana fermions and
localized visons [3]. This intriguing spin-liquid state can be stabilized in
tricoordinate lattices with bond-dependent Ising interactions, including
honeycomb and hyper-honeycomb lattices with strong spin-orbit coupling.
Extensive theoretical and experimental studies have been conducted to discover
materials exhibiting quantum spin-liquid states [4, 5].

Recently, metal-organic frameworks (MOFs) have been pointed out to be an
excellent platform for such exotic spin-liquid states. MOFs are crystalline
materials consisting of inorganic building units containing metal ions and
organic linkers [6]. The versatile choice of the building units and linkers
potentially provides MOFs with a significant advantage over inorganic
materials in stabilizing quantum spin liquids, including the Kitaev spin
liquid. Indeed, recent theoretical work predicts that the Ru oxalate
Ru2(C2O4)3 [7, 8] with the 82.10a-structure [9], or the lig net in O’Keeffe’s
three-letter reticular chemistry structure resource (RCSR) codes [10], becomes
a Kitaev spin liquid.

[(C2H5)3NH]2Cu2(C2O4)3 is a promising material for hosting an exotic spin-
liquid ground state. This Cu oxalate framework compound crystallizes in a
monoclinic structure with space group $P2_{1}/c$ [11]. The Cu ions form a
hyper-honeycomb structure, also known as the (10.3)b or ths net, connected by
oxalate linkers, as shown in Fig. 1 [12]. The Kitaev model has been exactly
solved on the (10.3)b net, and has a quantum spin-liquid ground state with
gapless Majorana modes and gapped vison excitations arising from the emergent
$\mathbb{Z}_{2}$ gauge field [13]. However, in [(C2H5)3NH]2Cu2(C2O4)3 there
are two distinct Cu2+ $S=1/2$ ions; the Cu1 ions form chains along the $c$
axis [Fig. 1(a)] while the Cu2 ions form chains along the $a$ axis [Fig.
1(b)]. Due to Jahn-Teller distortion, the magnetic interactions between the Cu
ions are anisotropic. Based on the orbital analysis and the distance between
the Cu sites [14], the magnetic interaction between the Cu2 ions with the
shorter interatomic distance, $J_{4}$, is expected to be a strong
antiferromagnetic (AF) interaction ($J_{4}/\mathrm{k_{B}}$> 200 K). The
interaction between Cu1 ions, $J_{1}$, is expected to be smaller and also AF
($J_{1}/\mathrm{k_{B}}\sim$ 60 K). The interaction between Cu1 and Cu2 ions,
$J_{2}$, is expected to be weakly AF, with magnitude
$J_{2}/\mathrm{k_{B}}\sim$ 25 K. The interaction $J_{3}$ between Cu2 ions with
the larger interatomic distance is much weaker than the other interactions.

The ground state of [(C2H5)3NH]2Cu2(C2O4)3 is quite exotic. The material is
experimentally confirmed to be a magnetic insulator, whereas density
functional theory calculations find a metal [15]. This discrepancy suggests
that the electronic correlations between electrons on the Cu2+ ions lead to a
Mott insulating phase and play an important role in the low-energy physics of
this material [15]. Muon spin rotation ($\mu$SR) [14] and nuclear magnetic
resonance (NMR) [16] detect no magnetic order down to 50 mK. Despite the Mott
insulating state, low-temperature ($T$) specific heat $C_{\mathrm{p}}$ down to
2 K reveals a large linear-in-$T$ contribution, suggestive of charge-neutral
fermionic excitations in this compound [14]. To explain this exotic ground
state, Zhang et al. proposed that a quantum spin liquid with a resonating-
valence-bond (RVB) state may be realized in this system [14]. However, it was
later argued, on the basis of a combined density-functional theory (DFT) and
quantum many-body theory [15], that the quantum disorder observed in
[(C2H5)3NH]2Cu2(C2O4)3 arises because the Cu2 ions form dimers along the
$J_{4}$ bonds [Fig. 1(c)], leaving weakly coupled chains of Cu1 ions along the
$J_{1}$ bonds. Indeed it was proposed [15] that the effective interchain
coupling should be unfrustrated but extremely weak, $\sim J_{2}^{2}/2J_{4}$,
leading to an extremely low Néel temperature even though the system is
bipartite.

In this paper, we probe the low-energy physics of single crystals of the MOF
oxalate [(C2H5)3NH]2Cu2(C2O4)3 through a systematic study of specific heat
$C_{\mathrm{p}}$ in applied magnetic fields $H$. We confirm the absence of
magnetic order down to 180 mK in $H$ = 0 and observe a sizable zero-field
linear-in-$T$ contribution to $C_{\mathrm{p}}$ at low $T$. Applied magnetic
fields strongly suppress the low-temperature specific heat, with
characteristic gapped behavior. We calculate the four-component relativistic
band structure and thence estimate the Dzyaloshinskii–Moriya interaction
(DMI). We find that the component of the DMI perpendicular to the $b$ axis
alternates along the $J_{1}$ bonds. The low-energy physics of a quasi-one-
dimensional antiferromagnetic Heisenberg chain (AFHC) in a magnetic field with
an alternating DMI perpendicular to the field is described by the sine-Gordon
model [17]. We show that the measured field dependence and magnitude of the
gap are in good agreement with those expected for solitons, calculated from
the sine-Gordon model using the parameters estimated from our DFT
calculations. We also find that the $T\rightarrow 0$ K Wilson ratio estimated
from our measurements is $R_{\mathrm{W}}=2.1(1)$, which agrees with the value
$R_{\mathrm{W}}=2$ expected for spinons in spin-1/2 AFHCs [18]. This strongly
suggests that [(C2H5)3NH]2Cu2(C2O4)3 is a quasi-one-dimensional Heisenberg
antiferromagnet and not an RVB quantum spin liquid.

Figure 1:  Crystal structure of [(C2H5)3NH]2Cu2(C2O4)3 [12]. Views are shown
along (a) the $a$ axis and (b) the $c$ axis. Cu1 ions form chains along the
$c$ axis, while Cu2 ions form chains along the $a$ axis. (c) Schematic diagram
of Cu network in [(C2H5)3NH]2Cu2(C2O4)3 forming a hyper-honeycomb structure,
also known as the ths net in O’Keeffe’s three-letter RCSR codes [10]. Due to
Jahn-Teller distortion, magnetic interactions between Cu2+ ions
($J_{1},J_{2},J_{3}$, and $J_{4}$) are anisotropic (see the main text). The
inset shows a photo image of a single crystal of [(C2H5)3NH]2Cu2(C2O4)3.
Figure 2:  Specific heat $C_{\mathrm{p}}$ divided by temperature,
$C_{\mathrm{p}}/T$, of [(C2H5)3NH]2Cu2(C2O4)3 as a function of $T$. Magnetic
fields are applied along the $b$ axis. The red solid line is a fit to the data
for $H$ = 0 T from 1 to 5 K using $C_{\mathrm{p}}/T=\gamma+\beta T^{2}+\delta
T^{4}$. The $\gamma$ value is obtained to be 98(1) mJ/mol K2. Inset:
Temperature dependence of $C_{\mathrm{p}}$ at $H=$ 0 T. Below about 2 K,
$C_{\mathrm{p}}(H=0~{}\mathrm{T})$ is approximately proportional to $T$.

##  II Experimental details

Single crystals of [(C2H5)3NH]2Cu2(C2O4)3 were synthesized by an aqueous-
solution method similar to that reported earlier [14]. Blue diamond-shaped
crystals were extracted by vacuum filtration [Fig. 1(c) inset]. We used
several crystals aligned along the $b$ axis. Specific-heat measurements were
conducted using the relaxation method with a homemade calorimeter. The
magnetic susceptibility $\chi$ was measured using a magnetic property
measurement system from Quantum Design, Inc. (1 T = 104 Oe). Torque
magnetometry was performed using a capacitive cantilever in a 35 T resistive
magnet at the National High Magnetic Field Laboratory, Tallahassee, Florida.

##  III Results and discussion

The specific heat, $C_{\mathrm{p}}(T)$, of [(C2H5)3NH]2Cu2(C2O4)3 shown in
Fig. 2 exhibits an unusual $T$ dependence at low $T$. We observe no magnetic
order in $C_{\mathrm{p}}$ down to $T=$ 180 mK, consistent with the previous
reports from $\mu$SR [14] and NMR [16], as shown in the inset to Fig. 2. The
zero-field specific heat is approximately proportional to $T$ at low $T$,
indicative of the presence of gapless charge-neutral magnetic excitations in
this Mott insulator.

Fitting $C_{\mathrm{p}}/T=\gamma+\beta T^{2}+\delta T^{4}$ to the data for $H$
= 0 T from 1 to 5 K (data are presented only below 4 K in Fig. 2), we extract
the linear-in-$T$ contribution to the zero-field specific heat, $\gamma T$,
associated with gapless charge-neutral magnetic excitations. The obtained
value of $\gamma$ per mole of formula units is 98(1) mJ/mol K2, somewhat
larger than that reported in the previous work ($\gamma=36$ mJ/Cu-mol K2 = 72
mJ/mol K2) [14]. We observe a deviation from the fitted curve below $T$ = 0.7
K. We will discuss this deviation later.

As shown in Fig. 2, the large linear-in-$T$ specific heat term due to gapless
magnetic excitations gradually diminishes with applied magnetic field parallel
to the $b$ axis, which is perpendicular to the quasi-one-dimensional chains
formed by Cu1 (along the $c$ axis) and Cu2 (along the $a$ axis) ions. At $H$ =
1 T, with decreasing temperature, $C_{\mathrm{p}}/T$ exhibits a broad peak at
$T\approx$ 0.8 K in Fig. 3 and then rapidly decreases. At larger magnetic
fields, the peak shifts to higher temperatures, and $C_{\mathrm{p}}/T$ becomes
zero at low temperatures, suggestive of the formation of an $H$-induced gap
for spin excitations. The gap can be seen more clearly in the temperature
dependence of the magnetic specific heat $C_{\mathrm{mag}}$ discussed below.

To investigate the magnetic contribution to the specific heat, we subtract the
phonon contribution, $C_{\mathrm{ph}}=\beta T^{3}+\delta T^{5}$, from the
total specific heat $C_{\mathrm{p}}$, and plot
$C_{\mathrm{mag}}(T)=C_{\mathrm{p}}(T)-C_{\mathrm{ph}}(T)$ in Fig. 3 for each
$H$. We clearly observe a crossover from linear-in-$T$ to gapped behavior with
increasing magnetic fields. The field-induced gap is also detected in the
field dependence of $C_{\mathrm{mag}}/T$, as shown in Fig. 4. At $T=$ 0.5 K,
$C_{\mathrm{mag}}/T$ due to gapless excitations at zero field decreases
rapidly with magnetic field, after showing a peak at about 0.5 T. Above 3 T,
$C_{\mathrm{mag}}(H)/T$ is nearly independent of magnetic field at $T=0.5$ K,
and its magnitude is about 5% of $C_{\mathrm{mag}}(H=0~{}\mathrm{T})/T$,
indicating that the ground state is completely gapped out in high magnetic
fields. Upon increasing the temperature, the peak moves to higher magnetic
fields and becomes broader.

Figure 3:  Temperature dependence of the magnetic contribution to the specific
heat of [(C2H5)3NH]2Cu2(C2O4)3 where
$C_{\mathrm{mag}}=C_{\mathrm{p}}-C_{\mathrm{ph}}$ and $C_{\mathrm{p}}$ is the
phonon contribution. The curves are successively shifted upward by 150 mJ/mol
K for clarity. The red dashed lines are fits to the sine-Gordon model (see the
main text). A clear development of an energy gap is observed upon applying
magnetic fields parallel to the $b$ axis. Figure 4:  Magnetic field $H$
dependence of the magnetic contribution to the specific heat
$C_{\mathrm{mag}}(T)/T$ of [(C2H5)3NH]2Cu2(C2O4)3at four temperatures. The
magnetic fields are applied along the $b$ axis. The peak position in
$C_{\mathrm{mag}}(T)/T$ shifts to a higher magnetic field with increasing
temperature.

The absence of magnetic order and the large linear-in-T contribution to
$C_{\mathrm{p}}$ at low temperatures observed in zero magnetic field
corroborate a charge-neutral gapless ground state, such as a quantum spin-
liquid state. To explain these observations, Zhang et al. proposed a
resonating-valence-bond (RVB) state in a hyper-honeycomb lattice of Cu2+ ions
associated with reduced dimensionality [14]. The proposed RVB state, however,
appears to be unlikely because in a spin-1/2 hyper-honeycomb lattice without
strong spin-orbit coupling, magnetic frustrations, which are key ingredients
to such spin-liquid states, are absent.

Instead, we attribute this behavior observed in [(C2H5)3NH]2Cu2(C2O4)3 to
weakly coupled one-dimensional antiferromagnetic Heisenberg chains (AFHCs)
stemming from a distorted hyper-honeycomb lattice elongated in the $b$-axis
direction due to the Jahn-Teller distortion of the CuO6 octahedra. Jacko and
Powell have pointed out that the Jahn-Teller distortion can cause the
formation of quasi-one-dimensional spin chains of Cu1 ions along the $c$ axis
and dimerization of Cu2 ions along the $a$ axis [15]. The values of the
exchange interactions extracted from the DFT and perturbation theory depend
somewhat on the assumed values of the on-site Coulomb (Hubbard) interaction,
$U$, and Hund’s rule coupling, $J_{\mathrm{H}}$. However, the relative
magnitudes of the interactions are $J_{4}\gg J_{1}\gg J_{2}\gg|J_{3}|$ for
reasonable choices of $U$ and $J_{\mathrm{H}}$, and the numerical values are
consistent with those estimated in Ref. [14] from orbital analysis. Given that
$J_{3}$ is negligible compared with $J_{4}$, the interchain coupling strength
$J_{\perp}$ can be estimated to be
$J_{\perp}/\mathrm{k_{B}}=J_{2}^{2}/(2J_{4}\mathrm{k_{B}})\simeq$ 0.1 K,
suggesting that [(C2H5)3NH]2Cu2(C2O4)3 can be treated as an AFHC.

In this situation, the zero-field linear-in-$T$ contribution, $\gamma T$, to
the specific heat can be attributed to the spinon contribution in one-
dimensional spin-1/2 AFHCs. The spinon specific-heat coefficient is given by
$\gamma=\frac{2R\mathrm{k_{B}}}{3J}$, where $R$ is the molar gas constant and
$J$ the exchange interaction between spins [18]. Extracted from the measured
$\gamma$ in the present work, $J/\mathrm{k_{B}}$ is estimated to be 56.7(8) K,
in excellent agreement with $J_{1}$ between Cu1 ions estimated from an orbital
analysis in Ref. [14].

The magnetic susceptibility $\chi(T)$ of this material also supports this
scenario. The overall $T$ dependence of $\chi$ at $H=$ 0.1 T ($H\parallel b$)
is plotted in Fig. 5. Associated with an order-disorder transition of
(C2H5)3NH+ cations [14], an anomaly indicated by the black arrow is observed
at $T\sim$ 170 K. Upon decreasing temperature, $\chi$ exhibits a broad maximum
at $T_{\mathrm{max}}\sim$ 35 K, indicating the presence of an AF exchange
interaction between the Cu ions in one-dimensional chains. At low $T$, $\chi$
shows a Curie-Weiss-type divergence, most likely due to impurity and/or defect
spins. This $T$ dependence of $\chi$ in [(C2H5)3NH]2Cu2(C2O4)3 is well
described by the sum of the contributions from one-dimensional AFHCs,
antiferromagnetic Heisenberg dimers (AFHDs), impurity/defect spins, and a
diamagnetic term. We analyze $\chi(T)$ quantitatively by using the following
equations,

$\displaystyle\chi_{\mathrm{fit}}$
$\displaystyle=\chi_{\mathrm{BF}}+\chi_{\mathrm{dimer}}+\chi_{\mathrm{CW}}+\chi_{0},$
(1) $\displaystyle\chi_{\mathrm{BF}}$ $\displaystyle=$
$\displaystyle\frac{N_{\mathrm{A}}g^{2}\mu_{\mathrm{B}}^{2}}{\mathrm{k_{B}}T}\frac{0.25+0.14995x+0.30094x^{2}}{1+1.9862x+0.68854x^{2}+6.0626x^{3}},$
(2) $\displaystyle\chi_{\mathrm{dimer}}$
$\displaystyle=\frac{N_{\mathrm{A}}g^{2}\mu_{\mathrm{B}}^{2}}{J_{\mathrm{dimer}}}\frac{e^{-1/t}}{t}\frac{1}{1+3e^{-1/t}},$
(3) $\displaystyle\chi_{\mathrm{CW}}$
$\displaystyle=\frac{C_{\mathrm{CW}}}{T-\theta},$ (4)

where $\chi_{\mathrm{BF}}$ is the Bonner–Fisher result [19], parameterized
with $x=J/(2\mathrm{k_{B}}T)$ by [20], $N_{\mathrm{A}}$ is Avogadro’s number,
$g$ is the spectroscopic-splitting factor, $\mu_{\mathrm{B}}$ is the Bohr
magneton, $\theta$ is the Weiss temperature, $\chi_{\mathrm{CW}}$ is a Curie-
Weiss contribution, $\chi_{\mathrm{dimer}}$ is the magnetic susceptibility of
the AFHDs [21, 18], $J_{\mathrm{dimer}}$ is the exchange interaction between
dimerized Cu2+ ions, $t=J_{\mathrm{dimer}}/\mathrm{k_{B}}T$, and $\chi_{0}$ is
the diamagnetic susceptibility ($-2.77\times 10^{-4}$ cm3/mol) estimated from
Pascal’s constants [22, 11]. As shown in Fig. 5, the fit is in excellent
agreement with our data from 2 to 140 K. The obtained fitting parameters are
$J/\mathrm{k_{B}}=54(1)$ K, $g=2.48(3)$, $C_{\mathrm{CW}}=1.78(2)\times
10^{-2}$ cm3K/mol, $\theta=0.98(1)$ K, and $J_{\mathrm{dimer}}=280(10)$ K. The
obtained $J=$ 54(1) K agrees with that found from our specific-heat
measurements. The obtained $J_{\mathrm{dimer}}$ is also consistent with
$J_{4}$ from theoretical calculations [14, 15]. Combining the obtained $g$
factor, $\chi_{\mathrm{BF}}(T\rightarrow 0~{}\mathrm{K})$, and $\gamma$ in the
present work, we evaluate the Wilson ratio,

$\displaystyle
R_{\mathrm{W}}=\frac{4}{3}\left(\frac{\pi\mathrm{k_{B}}}{g\mu_{\mathrm{B}}}\right)^{2}\frac{\chi}{\gamma}.$
(5)

We obtain $R_{\mathrm{W}}=2.1(1)$, which agrees with the expected value of
$R_{\mathrm{W}}=2$ for one-dimensional spin-1/2 AFHCs [18].

Figure 5:  Temperature dependence of magnetic susceptibility $\chi$ in
magnetic field $H$ = 0.1 T. The magnetic field is applied parallel to the $b$
axis. The red line is a fit to the data from 2 to 140 K using
$\chi_{\mathrm{fit}}=\chi_{\mathrm{BF}}+\chi_{\mathrm{CW}}+\chi_{\mathrm{dimer}}+\chi_{0}$,
where $\chi_{\mathrm{BF}}$ is the magnetic susceptibility for one-dimensional
AFHCs (the blue dashed line), $\chi_{\mathrm{CW}}$ is a contribution from
impurities/defects (the purple dashed line), $\chi_{\mathrm{dimer}}$ is a
contribution from AFHDs (the green dashed line), and $\chi_{0}$ is a
diamagnetic susceptibility estimated from the Pascal’s constants [22, 11]. The
arrow indicates an anomaly associated with an order-disorder transition of
(C2H5)3NH+ cations at $\sim$170 K [11]. Inset: Field dependence of magnetic
torque divided by magnetic field, $\tau/H=M_{\perp}V$, where $M_{\perp}$ is
the magnetization perpendicular to applied magnetic field and $V$ is the
sample volume, of [(C2H5)3NH]2Cu2(C2O4)3 at $T=0.4$ K. The magnetic field is
applied along the $b$ axis. No discernible anomaly associated with a magnetic
transition is observed.

However, the one-dimensional spin-1/2 AFHC remains gapless in magnetic fields
up to the saturation field. Our magnetic-torque data demonstrate no magnetic
transition nor saturation up to 35 T, as shown in the inset to Fig. 5. To
better understand this, we performed four-component relativistic DFT
calculations of the electronic structure in an all-electron full-potential
local orbital basis code (FPLO) [23], using the Perdew-Burke-Ernzerhof
generalized gradient approximation [24]. The density was converged on an
$(8\times 8\times 8)$ $k$ mesh. FPLO implements the full-potential non-
orthogonal local-orbital minimum basis band-structure scheme [23]. All
calculations were performed for the measured crystal structure [14]. These
calculations are thus the four-component analogue of the scalar relativistic
calculations reported in Ref. [15].

$n$ | Cu1 | Cu2 | 1 | 2 | 3 | 4
---|---|---|---|---|---|---
$t_{n}$ | $\left(\begin{array}[]{cc}-376&-149\\\ -149&-271\end{array}\right)$ | $\left(\begin{array}[]{cc}-417&-94\\\ -94&-237\end{array}\right)$ | $\left(\begin{array}[]{cc}-29&-90\\\ 119&169\end{array}\right)$ | $\left(\begin{array}[]{cc}245&-37\\\ -44&-36\end{array}\right)$ | $\left(\begin{array}[]{cc}36&124\\\ 124&141\end{array}\right)$ | $\left(\begin{array}[]{cc}36&-102\\\ -102&199\end{array}\right)$
$\lambda_{n}^{a*}$ | $\left(\begin{array}[]{cc}0&\pm 7\\\ \mp 7&0\end{array}\right)$ | $\left(\begin{array}[]{cc}0&7\\\ -7&0\end{array}\right)$ | $\left(\begin{array}[]{cc}\mp 2&\mp 12\\\ \mp 4&\pm 9\end{array}\right)$ | $\left(\begin{array}[]{cc}2&18\\\ -16&-2\end{array}\right)$ | $\left(\begin{array}[]{cc}0&-9\\\ 9&0\end{array}\right)$ | $\left(\begin{array}[]{cc}0&-8\\\ 8&0\end{array}\right)$
$\lambda_{n}^{b}$ | $\left(\begin{array}[]{cc}0&-25\\\ 25&0\end{array}\right)$ | $\left(\begin{array}[]{cc}0&-12\\\ 12&0\end{array}\right)$ | $\left(\begin{array}[]{cc}-3&-15\\\ -14&4\end{array}\right)$ | $\left(\begin{array}[]{cc}-4&2\\\ 1&2\end{array}\right)$ | $\left(\begin{array}[]{cc}0&14\\\ -14&0\end{array}\right)$ | $\left(\begin{array}[]{cc}0&-12\\\ 12&0\end{array}\right)$
$\lambda_{n}^{c}$ | $\left(\begin{array}[]{cc}0&\pm 13\\\ \mp 13&0\end{array}\right)$ | $\left(\begin{array}[]{cc}0&5\\\ -5&0\end{array}\right)$ | $\left(\begin{array}[]{cc}\pm 1&\pm 9\\\ \pm 9&\pm 9\end{array}\right)$ | $\left(\begin{array}[]{cc}-3&-10\\\ 17&-2\end{array}\right)$ | $\left(\begin{array}[]{cc}0&3\\\ -3&0\end{array}\right)$ | $\left(\begin{array}[]{cc}0&9\\\ -9&0\end{array}\right)$
Table 1: Tight-binding parameters calculated from DFT. Here $n=\text{Cu1}$ or
Cu2 indicates on-site parameters for the given sublattice, and $n=1-4$
indicates a bond labeled in Fig. 1(c) by the nearest-neighbor interaction
$J_{n}$. The symbols $\pm$ and $\mp$ in the Cu1 ($n=1$) columns indicate that
the spin-orbit-coupling parameter alternates on neighboring ions (bonds). All
values are in units of meV.

As expected for the light elements in [(C2H5)3NH]2Cu2(C2O4)3, we find only
minor difference in the band structure calculated from the four-component and
scalar relativistic theories. From the overlaps of the Wannier orbitals
extracted from the four-component calculation, we obtain a tight-binding model
of the form [25, 26]

$\displaystyle\hat{{\cal H}}_{\mathrm{tb}}=$
$\displaystyle\sum_{i,j}\hat{\bm{d}}_{i\alpha}^{\dagger}\left(t_{ij}\delta_{\alpha\beta}+i\bm{\lambda}_{ij}\cdot\bm{\sigma}_{\alpha\beta}\right)\hat{\bm{d}}_{j\beta},$
(6)

where
$\hat{\bm{d}}_{i\alpha}^{\dagger}=(\hat{d}_{i1\alpha}^{\dagger},\hat{d}_{i2\alpha}^{\dagger})$,
$\hat{d}_{i\mu\alpha}^{\dagger}$ $(\hat{d}_{i\mu\alpha})$ creates
(annihilates) an electron with spin $\alpha$ in the $\mu$th orbital on the
$i$th Cu ion, $\bm{\sigma}$ is the vector of Pauli matrices, $t_{ij}$ is the
scalar hopping integral, and $\bm{\lambda}_{ij}$ is the vectorial spin-orbit-
coupling hopping integral. The values of the hopping integrals calculated from
DFT are given in Table 1.

We model the electronic interactions on each atom by an $e_{g}$-symmetry
Kanamori interaction [27, 28]

$\displaystyle\hat{V}_{K}^{e_{g}}$ $\displaystyle=$ $\displaystyle
U\sum_{i,m}\hat{n}_{im\uparrow}\hat{n}_{im\downarrow}+U^{\prime}\sum_{i,m\neq
m^{\prime}}\hat{n}_{im\uparrow}\hat{n}_{im^{\prime}\downarrow}$ (7)
$\displaystyle+(U^{\prime}-J_{\mathrm{H}})\sum_{i,m<m^{\prime},\sigma}\hat{n}_{im\sigma}\hat{n}_{im^{\prime}\sigma}$
$\displaystyle-J_{\mathrm{H}}\sum_{i,m\neq
m^{\prime}}\hat{d}_{im\uparrow}^{\dagger}\hat{d}_{im\downarrow}\hat{d}_{im^{\prime}\downarrow}^{\dagger}\hat{d}_{im^{\prime}\uparrow}$
$\displaystyle+J_{\mathrm{H}}\sum_{i,m\neq
m^{\prime}}\hat{d}_{im\uparrow}^{\dagger}\hat{d}_{im\downarrow}^{\dagger}\hat{d}_{im^{\prime}\downarrow}\hat{d}_{im^{\prime}\uparrow},$

where $U$ ($U^{\prime}$) is the intra-(inter-)orbital Coulomb repulsion,
$J_{\mathrm{H}}$ is the Hund’s rule coupling, and
$\hat{n}_{im\sigma}=\hat{d}_{im\sigma}^{\dagger}\hat{d}_{im\sigma}$. For
simplicity, we approximate the on-site inter-orbital Coulomb repulsion as
$U^{\prime}=U-2J_{\mathrm{H}}$, which is exact in octahedral symmetry [28].

We calculate the magnetic interactions by perturbing around the large $U$
limit to second order in $t_{ij}$ and first order in $\bm{\lambda}_{ij}$ [29,
30]. The values of $U$ and $J_{\mathrm{H}}$ are not accurately known a priori
and have limited transferability between materials. Therefore, we estimate
them by comparing the calculated value of $J_{1}$ to that obtained
experimentally above. For example, for $U=6$ eV and $J_{\mathrm{H}}=0.3$ eV we
find $J_{1}=54$ K, consistent with the experimental estimate; therefore we fix
$U$ and $J_{\mathrm{H}}$ at these values.

Since we include spin-orbit coupling, we also find Dzyaloshinskii-Moriya
interactions (DMIs), Table 2. Below we focus on the DMI along the $J_{1}$
bonds, $\bm{D}_{1}=(D_{1}^{a^{*}},D_{1}^{b},D_{1}^{c})$, where the $a^{*}$
axis is perpendicular to both the $b$ and $c$ axes and $D_{1}^{n}$ is the
component of the DMI along the $n$ axis, because this is crucial for
understanding the experiments described above. We find, independent of the
values of $U$ and $J_{\mathrm{H}}$, that $D_{1}^{b}$ is the same for all
bonds, whereas $D_{1}^{a^{*}}$ and $D_{1}^{c}$ alternate in sign from bond to
bond. Thus the Hamiltonian for a single chain is

$\displaystyle\hat{{\cal
H}}_{1}=\sum_{i}\Big{\\{}J_{1}\hat{\bm{S}}_{i}\cdot\hat{\bm{S}}_{i+1}+\left[\bm{D}_{1}^{\|}+(-1)^{i}\bm{D}_{1}^{\perp}\right]\cdot\hat{\bm{S}}_{i}\times\hat{\bm{S}}_{i+1}$
$\displaystyle+g\mu_{\mathrm{B}}H\hat{S}_{i}^{b}\Big{\\}},$ (8)

where $\bm{D}_{1}^{\|}=(0,D_{1}^{b},0)$ is the (non-alternating) component of
the DMI parallel to the $b$ axis and
$\bm{D}_{1}^{\perp}=(D_{1}^{a^{*}},0,D_{1}^{c})$ is the (alternating)
component of the DMI perpendicular to the $b$ axis.

$n$ | 1 | 2 | 3 | 4
---|---|---|---|---
$J_{n}$ | 54 | 21 | 3 | 477
$\bm{D}_{n}=(D_{n}^{a^{*}},D_{n}^{b},D_{n}^{c})$ | $(\pm 13,3,\mp 16)$ | $(1.1,1.5,-0.6)$ | $\bm{0}$ | $\bm{0}$
Table 2: Calculated exchange interactions and DMIs for $U=6$ eV and
$J_{\mathrm{H}}=0.3$ eV. All values are in units of K.

In the presence of an external magnetic field, the Heisenberg chain with an
alternating DMI perpendicular to the field can be mapped to an XXZ model in an
alternating magnetic field [17] of magnitude

$\displaystyle
h=H\sin\left[\frac{1}{2}\arctan\left(\frac{|\bm{D}_{1}^{\perp}|}{J_{1}}\right)\right].$
(9)

This staggered field opens a gap, which can be described by the sine-Gordon
model for one-dimensional spin-1/2 AFHCs [31, 17]. In this model, the free
energy of the system is expressed in terms of the solution of a single
nonlinear integral equation for the complex quantity $\epsilon(\theta)$,

$\displaystyle\epsilon(\theta)$
$\displaystyle=-i\Delta\beta\sinh(\theta+i\eta^{\prime})$
$\displaystyle-\int_{-\infty}^{\infty}d\theta^{\prime}G_{0}(\theta-\theta^{\prime})\ln[1+\exp(-\epsilon(\theta^{\prime}))]$
(10)
$\displaystyle+\int_{-\infty}^{\infty}d\theta^{\prime}G_{0}(\theta-\theta^{\prime}+2i\eta^{\prime})\ln[1+\exp(-\bar{\epsilon}(\theta^{\prime}))],$

where $\beta=1/\mathrm{k_{B}}T$, $\Delta$ is a soliton mass, or gap,
$\bar{\epsilon}(\theta)$ is the complex conjugate of $\epsilon(\theta)$, and

$\displaystyle
G_{0}(\theta)=\int_{0}^{\infty}\frac{d\omega}{\pi^{2}}\frac{\cos(2\omega\theta/\pi)\sinh[\omega(\xi-1)]}{\sinh(\omega\xi)\cosh(\omega)}.$
(11)

The free energy density is given by

$\displaystyle f(\beta)=-\frac{\Delta N_{\mathrm{A}}}{\pi
v_{s}\beta}\textrm{Im}\int_{-\infty}^{\infty}d\theta\sinh(\theta+i\eta^{\prime})\ln[1+\exp(-\epsilon(\theta))],$
(12)

where $v_{s}$ is the spinon velocity. The free energy density is independent
of the value of $\eta^{\prime}$ if $0<\eta^{\prime}<\pi\xi/2$. Assuming that
the effect of the SU(2) symmetry breaking due to the applied magnetic field is
negligible, we fix $\xi=1/3$ and $v_{s}=\pi J_{1}/2$. Utilizing
$C_{\mathrm{p}}=-T\frac{\partial^{2}f(\beta)}{\partial T^{2}}$ and setting
$\Delta$ to be a fitting parameter, we fit the specific heat obtained from the
sine-Gordon model to the data in magnetic fields. The red dashed lines in Fig.
3 are the calculated specific heat curves, in excellent agreement with our
data.

Figure 6:  Field dependence of the soliton gap extracted from the magnetic
specific heat and NMR data [16]. The solid line is the soliton gap predicted
[17] from the sine-Gordon model, Eq. (13), for the DMI and exchange
interaction obtained from our DFT calculations. This curve contains no
adjustable parameters once $U$ and $J_{\mathrm{H}}$ are fixed and depends only
extremely weakly on the choice of these parameters.

The gap opened by the applied magnetic field is given by [17]

$\displaystyle\Delta(H)=1.85\left(\frac{g\mu_{\mathrm{B}}h}{2J_{1}}\right)^{2/3}J_{1}\left|\ln\left(\frac{g\mu_{\mathrm{B}}h}{2J_{1}}\right)\right|^{1/6}.$
(13)

The soliton gap extracted from the fits to the sine-Gordon model (Fig. 3) is
compared in Fig. 6 to this prediction. The agreement is excellent. Note that
once $U$ and $J_{\mathrm{H}}$ have been fixed to reproduce the experimental
estimate of $J_{1}$, there are no free parameters in Eq. (13). Furthermore,
the curve predicted by Eq. (13) depends extremely weakly on the choice of $U$
and $J_{\mathrm{H}}$. These results strongly indicate that the large linear-
in-$T$ contribution to the specific heat and magnetic-field-induced gapped
behavior stem from gapless spinons and gapped solitons in weakly coupled
spin-1/2 one-dimensional AFHCs, rather than the RVB state proposed in Ref.
[14].

Finally, we comment on the deviation from the fitting curve
$C_{\mathrm{p}}/T=\gamma+\beta T^{2}+\delta T^{4}$ in the specific heat of
[(C2H5)3NH]2Cu2(C2O4)3 in Fig. 2 at temperatures below 0.7 K $\sim
0.01J_{1}/\mathrm{k_{B}}$. In the similar temperature range below 1 K, NMR
reveals a clear signature of slowing down of Cu2+ spin dynamics [16]. The
deviation of the spinon contribution from linear-in-$T$ behavior in this
temperature range may be related to this slow spin dynamics. We note that a
similar deviation is also reported below $T=$ 1 K $\sim 0.01J/\mathrm{k_{B}}$
in the sine-Gordon spin system KCuGaF6 [32, 33]. The resemblance between the
two distinct materials may share a universal origin in one-dimensional
spin-1/2 AFHCs, requiring further experimental investigations.

##  IV Summary

We have measured the specific heat of [(C2H5)3NH]2Cu2(C2O4)3 single crystals
at very low temperatures in magnetic fields. We observe no anomaly in zero
field and a sizable linear-in-$T$ contribution to the specific heat associated
with spinon excitations. Upon applying magnetic fields, a soliton gap is
induced, varying as $\Delta(H)\propto H^{2/3}$. These results can be explained
by the sine-Gordon model for weakly coupled one-dimensional AFHCs as proposed
in [15], rather than an RVB state proposed in previous work [14]. Since
[(C2H5)3NH]2Cu2(C2O4)3 undergoes no magnetic transition down to 50 mK and has
weak interchain coupling, this compound is a canonical example of a one-
dimensional spin-1/2 AFHC. Interestingly, this weak interchain coupling in a
bipartite lattice is due to the strong dimerization of the Cu2 ions.
[(C2H5)3NH]2Cu2(C2O4)3 therefore represents a new route to extremely isolated
AFHCs complementing structural isolation and interchain geometric frustration.

###### Acknowledgements.

This work was supported by startup funding from the University of Central
Florida. C.D., H.S., and Y.N. were supported by an NSF Career DMR-1944975.
A.C.J. and B.J.P. were supported by the Australian Research Council,
DP160100060. The work at Ames National Laboratory was supported by the U.S.
Department of Energy, Division of Materials Sciences and Engineering. Ames
National Laboratory is operated for the U.S. Department of Energy by Iowa
State University under Contract No. DE-AC02-07CH11358. The National High
Magnetic Field Laboratory (NHMFL) is supported by the National Science
Foundation through DMR-1644779 and the State of Florida. S.Y., X.H., and Y.T.
were supported by the NHMFL UCGP Program.

## References

  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 EP (2010).
  * Savary and Balents [2017] L. Savary and L. Balents, Quantum spin liquids: a review, Rep. Prog. Phys. 80, 016502 (2017).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. 321, 2 (2006).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Rev. Mod. Phys. 89, 025003 (2017).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367, eaay0668 (2020).
  * Furukawa  _et al._ [2013] H. Furukawa, K. E. Cordova, M. O’Keeffe, and O. M. Yaghi, The Chemistry and Applications of Metal-Organic Frameworks, Science 341, 974 (2013).
  * Yamada  _et al._ [2017a] M. G. Yamada, H. Fujita, and M. Oshikawa, Designing Kitaev Spin Liquids in Metal-Organic Frameworks, Phys. Rev. Lett. 119, 057202 (2017a).
  * Yamada  _et al._ [2017b] M. G. Yamada, V. Dwivedi, and M. Hermanns, Crystalline Kitaev spin liquids, Phys. Rev. B 96, 155107 (2017b).
  * Wells [1977] A. Wells, Three Dimensional Nets and Polyhedra, Wiley-Interscience (Wiley, 1977).
  * O’Keeffe  _et al._ [2008] M. O’Keeffe, M. A. Peskov, S. J. Ramsden, and O. M. Yaghi, The Reticular Chemistry Structure Resource (RCSR) Database of, and Symbols for, Crystal Nets, Acc. Chem. Res. 41, 1782 (2008).
  * Zhang  _et al._ [2012] B. Zhang, Y. Zhang, and D. Zhu, $[(\mathrm{C}_{2}\mathrm{H}_{5})_{3}\mathrm{NH}]_{2}\mathrm{Cu}_{2}(\mathrm{C}_{2}\mathrm{O}_{4})_{3}$: a three-dimensional metal–oxalato framework showing structurally related dielectric and magnetic transitions at around 165 K, Dalton Trans. 41, 8509 (2012).
  * Momma and Izumi [2011] K. Momma and F. Izumi, VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Crystallogr. 44, 1272 (2011).
  * O’Brien  _et al._ [2016] K. O’Brien, M. Hermanns, and S. Trebst, Classification of gapless ${\mathbb{Z}}_{2}$ spin liquids in three-dimensional Kitaev models, Phys. Rev. B 93, 085101 (2016).
  * Zhang  _et al._ [2018] B. Zhang, P. J. Baker, Y. Zhang, D. Wang, Z. Wang, S. Su, D. Zhu, and F. L. Pratt, Quantum Spin Liquid from a Three-Dimensional Copper-Oxalate Framework, J. Am. Chem. Soc. 140, 122 (2018).
  * Jacko and Powell [2021] A. C. Jacko and B. J. Powell, Quasi-one dimensional magnetic interactions in the three-dimensional hyper-honeycomb framework $[({\mathrm{C}}_{2}{\mathrm{H}}_{5}{)}_{3}\mathrm{NH}{]}_{2}{\mathrm{Cu}}_{2}({\mathrm{C}}_{2}{\mathrm{O}}_{4}{)}_{3}$, Phys. Chem. Chem. Phys. 23, 5012 (2021).
  * Ding  _et al._ [2022] Q.-P. Ding, C. Dissanayake, S. Pakhira, W. J. Newsome, F. Uribe-Romo, D. C. Johnston, Y. Nakajima, and Y. Furukawa, Slow spin dynamics in the hyperhoneycomb lattice $[({\mathrm{C}}_{2}{\mathrm{H}}_{5}{)}_{3}\mathrm{NH}{]}_{2}{\mathrm{Cu}}_{2}({\mathrm{C}}_{2}{\mathrm{O}}_{4}{)}_{3}$ revealed by ${}^{1}\mathrm{H}$ NMR studies, Phys. Rev. B 105, L100405 (2022).
  * Oshikawa and Affleck [1997] M. Oshikawa and I. Affleck, Field-Induced Gap in $\mathit{S}\phantom{\rule{0.0pt}{0.0pt}}=\phantom{\rule{0.0pt}{0.0pt}}1/2$ Antiferromagnetic Chains, Phys. Rev. Lett. 79, 2883 (1997).
  * Johnston  _et al._ [2000] D. C. Johnston, R. K. Kremer, M. Troyer, X. Wang, A. Klümper, S. L. Bud’ko, A. F. Panchula, and P. C. Canfield, Thermodynamics of spin $S=1/2$ antiferromagnetic uniform and alternating-exchange Heisenberg chains, Phys. Rev. B 61, 9558 (2000).
  * Bonner and Fisher [1964] J. C. Bonner and M. E. Fisher, Linear Magnetic Chains with Anisotropic Coupling, Phys. Rev. 135, A640 (1964).
  * Estes  _et al._ [1978] W. E. Estes, D. P. Gavel, W. E. Hatfield, and D. J. Hodgson, Magnetic and structural characterization of dibromo- and dichlorobis(thiazole)copper(II), Inorg. Chem. 17, 1415 (1978).
  * Bonner  _et al._ [1979] J. C. Bonner, H. W. J. Blöte, J. W. Bray, and I. S. Jacobs, Susceptibility calculations for alternating antiferromagnetic chains, J. Appl. Phys. 50, 1810 (1979).
  * Bain and Berry [2008] G. A. Bain and J. F. Berry, Diamagnetic Corrections and Pascal’s Constants, J. Chem. Educ. 85, 532 (2008).
  * Koepernik and Eschrig [1999] K. Koepernik and H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme, Phys. Rev. B 59, 1743 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Jacko  _et al._ [2017] A. C. Jacko, A. L. Khosla, J. Merino, and B. J. Powell, Spin-orbit coupling in ${\mathrm{Mo}}_{3}{\mathrm{S}}_{7}{(\mathrm{dmit})}_{3}$, Phys. Rev. B 95, 155120 (2017).
  * Khosla  _et al._ [2017] A. L. Khosla, A. C. Jacko, J. Merino, and B. J. Powell, Spin-orbit coupling and strong electronic correlations in cyclic molecules, Phys. Rev. B 95, 115109 (2017).
  * Kanamori [1963] J. Kanamori, Electron Correlation and Ferromagnetism of Transition Metals, Prog. Theor. Phys. 30, 275 (1963).
  * Georges  _et al._ [2013] A. Georges, L. d. Medici, and J. Mravlje, Strong Correlations from Hund’s Coupling, Annu. Rev. Condens. Matter Phys. 4, 137 (2013).
  * Powell  _et al._ [2017] B. J. Powell, J. Merino, A. L. Khosla, and A. C. Jacko, Heisenberg and Dzyaloshinskii-Moriya interactions controlled by molecular packing in trinuclear organometallic clusters, Phys. Rev. B 95, 094432 (2017).
  * Jacko  _et al._ [2020] A. C. Jacko, E. P. Kenny, and B. J. Powell, Interplay of dipoles and spins in $\kappa-{(\mathrm{BEDT}-\mathrm{TTF})}_{2}X$, where $X=\mathrm{Hg}{(\mathrm{SCN})}_{2}\mathrm{Cl},\mathrm{Hg}{(\mathrm{SCN})}_{2}\mathrm{Br},\mathrm{Cu}[\mathrm{N}{(\mathrm{CN})}_{2}]\mathrm{Cl},\mathrm{Cu}[\mathrm{N}{(\mathrm{CN})}_{2}]\mathrm{Br}$, and ${\mathrm{Ag}}_{2}{(\mathrm{CN})}_{3}$, Phys. Rev. B 101, 125110 (2020).
  * Essler and Tsvelik [1998] F. H. L. Essler and A. M. Tsvelik, Dynamical magnetic susceptibilities in copper benzoate, Phys. Rev. B 57, 10592 (1998).
  * Umegaki  _et al._ [2012] I. Umegaki, H. Tanaka, T. Ono, M. Oshikawa, and K. Sakai, Thermodynamic properties of quantum sine-Gordon spin chain system $\mathrm{KCuGaF}_{6}$, Phys. Rev. B 85, 144423 (2012).
  * Umegaki  _et al._ [2015] I. Umegaki, H. Tanaka, N. Kurita, T. Ono, M. Laver, C. Niedermayer, C. Rüegg, S. Ohira-Kawamura, K. Nakajima, and K. Kakurai, Spinon, soliton, and breather in the spin-$\frac{1}{2}$ antiferromagnetic chain compound ${\mathrm{KCuGaF}}_{6}$, Phys. Rev. B 92, 174412 (2015).