# High-Temperature Decomposition of Diisopropyl Methylphosphonate (DIMP) on
Alumina: Mechanistic Predictions from Ab Initio Molecular Dynamics

Sohag Biswas  Bryan M. Wong<EMAIL_ADDRESS>Website: http://www.bmwong-
group.com [

###### Abstract

The enhanced degradation of organophosphorous-based chemical warfare agents
(CWAs) on metal-oxide surfaces holds immense promise for neutralization
efforts; however, the underlying mechanisms in this process remain poorly
understood. We utilize large-scale quantum calculations for the first time to
probe the high-temperature degradation of diisopropyl methylphosphonate
(DIMP), a nerve agent simulant. Our Born-Oppenheimer molecular dynamics (BOMD)
calculations show that the $\gamma$-Al2O3 surface shows immense promise for
quickly adsorbing and destroying CWAs. We find that the alumina surface
quickly adsorbs DIMP at all temperatures, and subsequent decomposition of DIMP
proceeds via a propene elimination. Our BOMD calculations are complemented
with metadynamics simulations to produce free energy paths, which show that
the activation barrier decreases with temperature and DIMP readily decomposes
on $\gamma$-Al2O3. Our first-principle BOMD and metadynamics simulations
provide crucial diagnostics for sarin decomposition models and mechanistic
information for examining CWA decomposition reactions on other candidate metal
oxide surfaces.

Unknown University] Department of Chemical & Environmental Engineering,
Materials Science & Engineering Program, Department of Physics & Astronomy,
and Department of Chemistry, University of California-Riverside, Riverside,
California 92521, USA

##  1 Introduction

The neutralization of chemical warfare agents (CWAs) continues to be a
pressing area of interest for the safe and effective removal of these
hazardous compounds. Among the various CWAs, the most nefarious are
organophosphate nerve agents (such as sarin), which contain P=O, P–O–C, and
P–C bonds that enable lethal phosphorylating mechanisms.1 Over the past few
decades, a variety of destruction and neutralization methods have been used to
safely eliminate these hazardous compounds. For example, destruction-based
methods (typically pyrolysis) allow a one-step approach for the complete
disposal of CWAs at the expense of using specialized equipment under extreme
conditions1, 2. On the other hand, neutralization methods offer potentially
reversible chemical treatments, leading to possible CWA precursors under less
severe conditions.1

Recent studies have shown that metal oxides can effectively destroy CWAs due
to their high surface area, a large number of highly reactive edges, corner
defect sites, unusual lattice planes, and high surface-to-volume ratios. In
particular, metal oxides, such as CaO,3, 4 MgO,5, 6 ZnO,7, 8, 9 TiO2,10, 11,
12, 13, 14, 15, 16 Al2O3,17, 18, 19, 20 Fe3O4,21 and CuO22 are candidates as
adsorbents for enhancing the decomposition of CWAs. In $\gamma$-Al2O3, Al
atoms in the bulk exhibit either a tetrahedral or octahedral coordination.
However, depending on the exposed crystallographic surface, Al atoms on the
surface can display penta-, tetra, and tri-coordination and exhibit Lewis
acidity.23, 24, 25, 26 Because of its high degree of surface heterogeneity,
$\gamma$-Al2O3 offers a high catalytic activity and is a promising candidate
for the decomposition of various CWAs.27

In this work, we present the first ab initio molecular dynamics study for
probing high-temperature decomposition dynamics of diisopropyl
methylphosphonate (DIMP) on the $\gamma$-Al2O3 surface. Due to its structural
similarity with sarin, DIMP has been used in experiments to mimic the
decomposition reaction mechanism of CWAs. DIMP is the only surrogate with the
same isopropyl (–O–C3H7) group found in sarin gas (the only structural
difference between sarin and DIMP is the substitution of a fluorine for an
isopropyl group in the former). Several techniques have also been used to
probe DIMP decomposition, including microwave discharge approaches,28
pyrolysis on porous substrates,29 laser heating,30 and thermal decomposition
at 700 – 800 K.31 Despite these experimental studies, theoretical analyses of
DIMP decomposition on metal oxides are scarce, except for a previous study on
mechanism and rates of thermal decomposition of DIMP.32 In this work, we
utilize large-scale ab initio molecular dynamics simulations to probe the
adsorption dynamics and time scales of DIMP decomposition at various
temperatures. In addition, we present new ab initio-based metadynamics
simulations to calculate free-energy barriers for various bond breaking
decomposition reactions of DIMP. These computational techniques allow us to
(1) predict activation energies and detailed mechanistic pathways at various
temperatures and (2) establish accurate sarin decomposition models on metal-
oxide surfaces to guide experimental efforts for neutralizing DIMP.

##  2 Simulation Details

###  2.1 Molecular Dynamics Simulations

Density functional theory (DFT) based Born-Oppenheimer molecular dynamics
(BOMD) simulations were carried out using the CP2K33 software suite. We have
specifically chosen to use this software package since the implementation of
linear-scaling Kohn-Sham approaches in CP2K allows robust and efficient
electronic structure calculations for large systems. The Perdew-Burke-
Ernzerhof (PBE)34 was used for the DFT calculations with Grimmes’s D3 method
to account for dispersion forces.35 To obtain reasonable accuracy, we utilized
a DZVP (double zeta valence polarized) basis set for Al in the DFT
calculations and the TZV2P basis set for C, O, H, and P atoms with the
Goedecker, Teter, and Hutter (GTH) pseudopotentials36, 37 for atomic core
electrons. A similar basis set for Al atoms was also used to accurately
calculate methane activation38 and alkane dehydrogenation39 on the
$\gamma$-Al2O3 surface. The orbital transformation method with an electronic
gradient tolerance value of $1\times 10^{-5}$ atomic units was adopted as the
convergence criteria for the SCF cycle.40 A kinetic energy cut-off of 400 Ry
for the auxiliary plane-wave basis with a 0.5 fs timestep was employed to
integrate the equations of motion. The initial guess was furnished by the
stable predictor-corrector extrapolation method at each molecular dynamics
step.41, 42 We have carried out several molecular dynamics simulations in the
200 – 1000∘ C temperature range in increments of 100∘ C (i.e., a total of nine
independent trajectories were considered).

Figure 1: (a) Bulk structure of $\gamma$-Al2O3. (b) Molecular structure of
diisopropyl-methylphosphonate (DIMP). (c) Top view of DIMP adsorbed on the
$\gamma$-Al2O3 surface. (d) Side view of DIMP adsorbed on the $\gamma$-Al2O3
surface. The blue, red, orange, brown, and light pink spheres represent Al, O,
P, C, and H atoms, respectively.

All simulations were carried out using a Nose-Hoover chain thermostat with the
canonical ensemble (NVT).43, 44

The nonspinel model of the bulk $\gamma$-alumina unit cell (shown in Figure
1a) was built using the crystallographic model by Dinge et al.23, 45, which
has been shown to match well with experimental structural parameters.46 Our
calculated cell parameters ($a=7.90$ Å, $b=7.93$ Å, and $c=8.07$ Å)45, 23 for
bulk $\gamma$-alumina are in excellent agreement and within 2$\%$ of the
experimental cell parameters ($a=b=7.96$ Å and $c=7.81$ Å).46 In our
simulations, we used the (100) facet of this model, which is the lowest-energy
facet reported for this material.23 It is worth mentioning that our
calculations represent a simplified model of the metal-oxide surface, and
impurities (such as H2O, OH, SOx, etc.) may be present and play a crucial role
in its reactivity. Nevertheless, our large-scale BOMD calculations still
provide critical atomistic insight into the reactivity of the original
pristine material, which serves as a baseline for comparing its catalytic
activity against other mixed metal oxide surfaces (and their associated
impurities), which we save for future studies.

For our NVT simulations, a $3\times 1\times 2$ supercell of $\gamma$-alumina
and a single DIMP molecule (Figure 1b) containing a total of 268 atoms were
used. We introduced a single DIMP molecule 5 – 6 Å above the center of the
alumina slab along the $y$-direction, as shown in Figure 1c. Periodic boundary
conditions were applied in the $x$ and $z$ directions. Thus, the $xz$ plane of
the slab is parallel to the surface, and the $y$-axis forms the surface normal
where DIMP interacts with the alumina surface. We introduced a vacuum layer of
15 Å on top of the surface to avoid any significant overlap between the
electronic density of periodically translated cells. Our BOMD simulations show
that the DIMP molecule moves extremely fast and explores a large space of
configurations. To limit this exploration to regions where DIMP dissociation
on alumina might occur, it is essential to restrict the movement of the DIMP
molecule. Specifically, an external spherical potential, which only acts on
the DIMP molecule, was placed at the center of the system. Due to the
computationally expensive nature of these simulations, we performed 12 ps NVT
simulations for each trajectory. We also calculated adsorption energies using
conventional geometry optimizations using the Broyden-Fletcher-Gold-farb-
Shanno (BFGS) minimization algorithm until the forces converged to $4.0\times
10^{-4}$ Bohr with an SCF convergence criteria of $1\times 10^{-5}$ au. For
the geometry optimization, we used a relatively small supercell ($1\times
2\times 2$) compared to our BOMD simulations.

###  2.2 Metadynamics Simulations

Throughout our BOMD simulations, we did not observe any decomposition of DIMP
on alumina within the 200 – 600∘ C temperature range (decomposition did occur
at higher temperatures, which is discussed later in this paper). At lower
temperatures, the high computational cost of these simulations only permits
explorations of short time intervals; therefore, BOMD simulations alone do not
permit a routine exploration of the complete reaction dynamics. One can
sidestep this limitation by using advanced computational techniques such as
metadynamics simulations (MetaD) 47, which we describe further below. In
short, MetaD bypasses the sampling limitations of traditional molecular
dynamics by applying a history-dependent biasing potential as a function of
time to enable efficient sampling of the free energy surface. The free energy
surface itself is defined over a set of collective variables (CVs), which are
carefully chosen to provide a complete description of the system’s slow
degrees of freedom. The decomposition process is characterized by the breaking
of various bonds within the DIMP molecule on the alumina surface. For this
reason, we utilized two CVs, shown in Figure 2, for our accelerated sampling
simulations.

Figure 2: Collective variables (CVs) used to describe the adsorption dynamics
of DIMP on the alumina surface.

The CV specifies the coordination number (CN) during the MetaD simulations,
which is expressed as a function of the distance between two atoms:

$\text{CV}\text{ or
}\text{CN}={{1-\left({{d_{\text{AB}}}\over{d_{0}}}\right)^{p}}\over{1-\left({{d_{\text{AB}}}\over{d_{0}}}\right)^{p+q}}},$
(1)

where $d_{\text{AB}}$ is the distance between atoms A and B, and $d_{0}$ is
the reference distance or fixed cut-off parameter (this parameter
characterizes the standard bond distance between atoms A and B). The variables
$p$ and $q$ in Equation 1 are constants, which were fixed to $p=q=6$.

In this work, we chose CV1 to be the distance between the C and O atoms within
DIMP, and CV2 as the distance between the O atom in DIMP and an Al atom on the
alumina surface, as shown in Figure 2. MetaD simulations were carried out by
depositing Gaussians with heights of 0.02 and 0.001 Hartree for the 200 – 500∘
C and 600 – 1000∘ C temperature ranges, respectively. The widths of the
Gaussians were set to 0.1 for both the CV1 and CV2 simulations. In this work,
well-tempered MetaD (wt-MetaD) simulations were used 48, 49 with the
deposition rate of the Gaussian hills set to 20 steps. The well-tempering was
implemented using a Gaussian height damping factor of $\Delta$T such that the
ratio $\frac{\Delta\text{T}+\text{T}}{\text{T}}$ was equal to 6. We performed
two sets of wt-MetaD simulations for each temperature from different initial
conditions, which were extracted from the pre-equilibrated BOMD simulations.
We confirmed a representative sampling and convergence of the reactant,
product, transition state, and free energy differences, particularly for the
dissociation of the C–O bond in DIMP on the alumina surface.

##  3 Results and Discussions

###  3.1 Adsorption of DIMP on the Alumina Surface

Figure 3: Time-dependent fluctuations of various Al–O distances, which confirm
DIMP adsorption throughout the BOMD simulations for all temperatures.

Metal oxide surfaces are typically catalytic and can adsorb and subsequently
decompose CWAs into benign products. Prior experiments50 and theoretical27
calculations have shown that sarin simulants can adsorb on metal oxides by
forming a bond with metal atoms via a phosphoryl oxygen. Initially, we
observed the DIMP molecule interacting with an Al atom via the O atom of the
P–OC3H7 moiety. To further explore these dynamical effects with our BOMD
simulations, Figure 3 plots the distance between one of the surface Al atoms
and the oxygen of the P–OC3H7 moiety in DIMP at various temperatures. These
calculations show that DIMP adsorption occurs within two picoseconds for all
temperatures, and the oxygen of the P–OC3H7 moiety of DIMP interacts with the
tetra-coordinated Al center. We did not observe any desorption of DIMP on the
alumina surface, as indicated by the small Al48–O3 bond distances in Figure 3
(as a side note, the 600∘ C plot in the center of Figure 3 does show
dissociation of the Al48–O3 bond, but a new Al11–O2 bond quickly forms
thereafter, and the DIMP molecule still remains on the alumina surface).
Figure 4 shows snapshots from our simulations depicting various DIMP
adsorption configurations on the alumina surface. As suggested by prior
theoretical calculations,51, 52, 53 the adsorption via the O atom of the
P–OC3H7 moiety is crucial for propene elimination. As the MD simulations
progresses, we also observe interactions between the O (in the –P=O group) and
Al atoms. In short, the DIMP molecule remained adsorbed on the alumina surface
in our simulations for all temperatures.

Figure 4: Snapshots illustrating the adsorption of DIMP on the alumina
surface. From 200 – 500∘ C, adsorption occurs via the formation of the Al48–O3
bond. At 600∘ C, adsorption proceeds via the formation of the Al48–O3 and
Al11–O2 bonds. At 700 and 1000∘ C, the Al48–O3 bond is formed during the
adsorption process, and from 800 – 900∘ C, an Al43–O2 bond is formed. Figure
5: Various adsorption configurations of DIMP on the $\gamma$-Al2O3 (100)
surface. Bond distances are shown in angstroms.

To complement our MetaD simulations, we also calculated adsorption energies,
$E_{\text{ads}}$, using the following expression:

$E_{\text{ads}}=E_{\text{DIMP}+\text{surface}}-E_{\text{surface}}-E_{\text{DIMP}},$
(2)

where $E_{\text{surface}}$ is the energy of the alumina surface,
$E_{\text{DIMP}}$ is the energy of an isolated gas-phase DIMP molecule, and
$E_{\text{DIMP + surface}}$ represents the energy of the adsorbed molecule on
the surface. A negative value of $E_{\text{ads}}$ corresponds to an exothermic
process and a stable adsorption configuration. The tri-coordinated Al surface
atoms are known to be strong Lewis acid-type catalytic sites 54 and give large
adsorption energies when the DIMP molecule binds to them. Figure 5 illustrates
selected optimized structures of DIMP on the alumina surface. In panel 5(a),
the DIMP molecule is bonded to two tetra-coordinated Al atoms via the
phosphoryl oxygen, resulting in a bridging adsorption with
$E_{\text{ads}}=-48.63$ kcal/mol. In panel 5(b), the alkoxy oxygen is bonded
to a tri-coordinated Al atom, giving an adsorption energy of -101.40 kcal/mol.
The largest adsorption energy (-153.55 kcal/mol) was obtained for the
configuration shown in panel 5(c) where the alkoxy and phosphoryl oxygen atoms
are bonded to two different tri-coordinated Al centers. In panel 5(d), the
phosphoryl oxygen is bonded to a tri-coordinated Al center, giving an
adsorption energy of -93.12 kcal/mol. The adsorption energy for configuration
in panel 5(e) is -117.97 kcal/mol in which the DIMP molecule binds with two
tetra-coordinated Al atoms via two alkoxy oxygen atoms. Our calculated
adsorption energy for the configuration shown in panel 5(a) is in good
agreement with a previously studied similar configuration of sarin adsorption
on $\gamma$-Al2O3 (-49.2 kcal/mol).27

###  3.2 Decomposition of DIMP on the Alumina Surface

Figure 6: Evolution of C–O distances as a function of time on the alumina
surface at various temperatures, which describes DIMP decomposition. The
corresponding atom labeling is shown in the inset panel of (a).

We next examine DIMP decomposition mechanisms on the alumina surface. Previous
studies have shown that one of the decomposition pathways of DIMP proceeds
through the breaking of a C–O bond.28 Figure 6a displays time-dependent
variations of C–O bond distances from 200 – 600∘ C, which fluctuate around
1.50 Å (the covalent C–O bond distance), indicating that the C–O bond remains
intact. On the other hand, Figure 6b shows that when the temperature is raised
to 700 – 1000∘ C, the C–O bond lengthens significantly, signifying
decomposition. It is worth noting that the timescales for DIMP decomposition
varies with temperature, and we observe a very rapid decomposition at 1000∘ C.

Figure 7 depicts snapshots at crucial points for the decomposition of DIMP on
the alumina surface for a few representative MD trajectories. As mentioned in
the previous section, adsorption initially occurs via the oxygen atom of the
P–OC3H7 moiety in DIMP and an Al atom on the alumina surface (which occurs at
all temperatures).

Figure 7: Snapshots of DIMP decomposition from representative MD trajectories
at various temperatures.

At 700∘ C, the 3 panels at the top of Figure 7 depict a possible decomposition
pathway in which the C3–O3 bond breaks after 2 ps and an Al–O bond forms with
the alumina surface. During this process, a surface-bound n-alkyl (–C3H7)
species is formed. At $\sim$3 ps, a single terminal C–H bond from the n-alkyl
molecule is broken when a proton is abstracted by an oxygen atom on the
alumina surface. This forms propene, an unsaturated compound, as a by-product.
As the simulation further proceeds, the P–O bond in DIMP dissociates at 7 ps
and forms a new Al–OCHCH3CH3 adsorbate on the alumina surface. In addition,
once the P–O bond in DIMP is cleaved, a new P–O bond is formed between the
phosphorus atom and a di-coordinated surface oxygen atom (see the 3 panels at
the top of Figure 7). At 800 – 1000∘ C (depicted in the bottom 9 panels of
Figure 7), DIMP decomposition involves only the propene elimination via a two-
step process. The first step of the propene elimination is associated with the
dissociation of the C2–O2 (800 – 900∘ C) or C3–O3 (1000∘ C) bond. The second
step comprises the migration of the hydrogen atom from the methyl group of the
–C3H7 fragment to one of the surface oxygen atoms. Within these simulated
timescales, we did not observe any further decomposition of DIMP within the
800 – 1000∘ C temperature range.

###  3.3 Free energy profiles

The propene end-products predicted by our BOMD simulations corroborates
previous experiments on DIMP, which include thermal decomposition studies, 55,
56 pyrolysis and combustion in nitrogen/oxygen-rich environments, 57 as well
as microwave 28 and laser-induced 30 decomposition under inert environments.
Collectively, all of these prior experimental studies detected propene as one
of the main by-products of DIMP decomposition, which further supports our BOMD
predictions.

Figure 8: Reconstructed free energy surface for the decomposition of DIMP on
the alumina surface at (a) 1000∘ C, (b) 900∘ C, (c) 800∘ C, and (d) 700∘ C.
The R, TS, I, and P labels in each free energy surface correspond to the
reactant, transition state, intermediate, and product, respectively.

To further investigate the decomposition of DIMP on alumina, we utilized well-
tempered metadynamics simulations by adopting two collective variables. We
performed two sets of metadynamics simulations (for each temperature within
the 200 – 1000∘ C range) using different initial conditions. Obtaining
accurate free energy profiles from metadynamics simulations requires (1)
longer simulation times until all accessible regions of the potential are
explored (with trajectories spanning forward and backward many times between
reactant and product states) and (2) a careful selection of collective
variables. Our well-tempered metadynamics simulations at various temperatures
suggests that when the C–O bond is broken, subsequent reforming of this bond
is prohibited (Figure S19), which indicates that the reaction is irreversible.
The Supporting Information (Figures S20 and S21) provides further details on
the convergence of our free energy profiles between the transition state and
reactant basin as a function of time, which demonstrates that our metadynamics
simulations are fully converged.

Figure 9: Snapshots from the well-tempered metadynamics simulations within the
700 – 1000∘ C temperature range. The R, TS, and P labels illustrate the
reactant, transition state, and product, respectively.

The free energy surface at 1000∘ C is shown in Figure 8 (a). The two CVs that
capture the energetics of DIMP decomposition are presented in Figure S1. In
this mechanism, the reactant (R) goes to a product (P) via a transition state
(TS) and a small tiny intermediate (I) in between the transition state and
product. The reactant state is defined by CV$1=0.75$ (C–O $=1.66$ Å) and
CV$2=0.11$ (Al–O $=2.17$ Å ), where DIMP is bonded to the surface Al atom. The
decomposed DIMP product, P, is characterized by CV$1=0.03$ (C–O $=3.12$ Å) and
CV$2=0.40$ (Al–O $=1.77$ Å), in which the C–O bond of DIMP is cleaved. At this
stage, the H atom from the –C3H7 moiety is transferred to a surface oxygen
atom, and propene is formed. The C–O bond length stretches to 1.87 Å at the
transition state (which has an activation barrier of only 4.11 kJ/mol) and
subsequently dissociates. The reactant, transition state, and product from our
well-tempered metadynamics simulations at 1000∘ C are shown in Figure 9(a). In
the Supplementary Material, we provide an additional free energy profile
(Figure S2) at 1000∘ C that utilizes different initial conditions. For this
separate case, the reactant proceeds to the product state via a stable
transition state with a tiny minimum, and the activation barrier is $\sim$6.36
kJ/mol. Collectively, the average free energy value from these independent
trajectories is 5.24 kJ/mol.

The reaction mechanism at 900∘ C is very similar to that of 1000∘ C, which
also shows a tiny intermediate between the product and transition state, as
shown in Figure 8(b). The free energy barrier value is $\sim$11.85 kJ/mol, and
Figure S3 in the Supplementary Material shows the corresponding evolution of
CV1 and CV2. The geometries of the reactant, transition state, and product
along the metadynamics trajectory are reported in Figure 9(b). We have also
performed additional metadynamics calculations using different initial
conditions to test the reproducibility of these results. In these additional
calculations (depicted in Figure S4 in the Supplementary Information), the
reactant proceeds to the product through a single transition state with a free
energy barrier value of $\sim$15.73 kJ/mol. The average free energy barrier
value obtained from these different initial conditions is 13.79 kJ/mol.

The reconstructed free energy surfaces at 800 and 700∘ C are shown in Figures
8(c) and 8(d), and fluctuations of the corresponding CV values are presented
in Figures S5 and S6, respectively. In contrast to the higher temperatures
discussed previously, the decomposition of DIMP at 800 and 700∘ C proceeds via
a single transition state. At 800∘ C, the collective variables CV1 (0.71) and
CV2 (0.12) define the reactant state R, at which DIMP forms a bond with the
surface Al atom (typical C–O and Al–O bond distances at this geometry are 1.53
and 2.25 Å, respectively). The C–O bond in DIMP stretches to 1.91 Å to form a
transition state, and the reaction progresses to the product where the C–O
bond subsequently dissociates, and propene is formed. The net free energy
barrier for this activation process is $\sim$22.53 kJ/mol. The additional free
energy profile at 800∘ C from the different initial conditions is shown in
Figure S7. A very similar reaction mechanism is also observed at 700∘ C with a
slightly higher activation barrier of $\sim$26.65 kJ/mol. Figure S8 in the
Supplementary Material represents the topology of another free energy surface
at 700∘ C from different initial conditions. The reactant, transition state,
and product geometries from these metadynamics simulations are shown in
Figures 9(c) and 9(d) for 800 and 700∘ C, respectively.

Figure 10: Free energy activation barrier as a function of temperature for the
decomposition of DIMP on the alumina surface. The free energy activation
barrier value is calculated by averaging over two metadynamics simulations at
each temperature.

A detailed analysis of free energy profiles for DIMP decomposition on alumina
at various temperatures (200 – 600∘ C) with different initial conditions is
given in the Supplementary Material (Figures S9-S18). Figure 10 summarizes the
free energy activation barrier as a function of temperature, which shows that
the free energy activation barrier decreases with temperature. Static DFT
calculations have obtained free energy activation barrier values of 113 kJ/mol
(ZnO surface),51 108.0 kJ/mol (rutile surface),53 122.6 kJ/mol (anatase
surface),53 and 53.7 kJ/mol (MoO2 surface).52 However, the free energy
activation barriers for C–O bond breaking from our simulations are much lower
than these previously reported values. We also note that the free energy
activation barrier value at 600∘ C is very similar to that of 700∘ C.

Figure 11: Panel (a) depicts collective variables used in studying the
formation of an Al–OCHCH3CH3 adsorbate at 700∘ C, and panel (b) shows the 3D
reconstructed free energy surface, where R, TS1, I, TS2, and P correspond to
the reactant, transition state 1, intermediate, transition state 2, and
product, respectively. Panel (c) shows snapshots of R, TS1, I, TS2, and P
along the metadynamics trajectory.

At 700∘ C, we observed formation of an Al–OCHCH3CH3 adsorbate in addition to
propene. To further analyze these energetics with metadynamics simulations, we
utilized two collective variables: CV3 [P1–O2] and CV4 [P1–064], which denote
the coordination number of the P1 phosphorus with respect to the O2 and O64
oxygen atoms, respectively, in the DIMP molecule. The CVs are shown in Figure
11(a), and the corresponding free energy surface is depicted in Figure 11(b).
The reaction mechanism takes place in two substeps, as shown in Figure 11(c).
The reactant is described by CV$3=0.80$ and CV$4=0.05$, in which the P1–O2 and
P1–O64 bond distances are 1.69 and 3.02 Å, respectively. The P1–O2 bond first
stretches to 1.95 Å and forms the first transition state (TS1). The P1–O2 bond
then dissociates to 2.46 and creates a stable intermediate, I (CV$3=0.15$,
CV$4=0.04$), giving a free energy activation of 18.51 kJ/mol. The reaction
proceeds over a second transition state (TS2) in which the P1–O2 bond is
further stretched and the resulting –OCHCH3CH3 fragment binds to the surface
Al atom. In the product state, the distance from the P1 atom to the surface
oxygen (O64) shortens (1.74 Å), indicating the formation of a P1–O64 covalent
bond. In contrast, the P1–O2 bond length further increases, showing the
complete decomposition of the DIMP molecule.

Collectively, our AIMD calculations show that the decomposition of DIMP most
likely progresses via a propene elimination on the alumina surface. The final
step of the propene elimination occurs via the abstraction of a hydrogen atom
by the surface oxygen atom of gamma-alumina. We also performed metadynamics
simulations to evaluate the free energy activation of this process. The C-H
coordination (CV5) and H–O coordination (CV6) were selected as collective
variables, as shown in Figure S22 in the Supporting Information. A Gaussian
with a height of 0.0001 Hartree was used for these simulations. Figure S23 in
the Supporting Information illustrates free energy profiles for the C–H
activation within the 200 – 500∘ C temperature range. In Figure S24, we show
free energy profiles in the 600 – 900∘ C temperature range. The free energy
profile at 1000∘ C is shown in Figure S25a, and Figure S25b summarizes the
free energy activation barrier as a function of temperature. The free energy
activation barriers range from 3.91 – 1.29 kJ/mol for the 200 – 1000∘ C
temperature range. We find that the free energy activation barrier value
decreases with increasing the temperature. The small free energy barrier value
in the 600 – 1000∘ C temperature range lies between 5 – 28 kJ/mol, which can
be easily accessible by experiment. Our theoretical findings are also
consistent with experimental studies, which also identified propene as the
main gas-phase product of the DIMP decomposition. In summary, our AIMD
simulations show that the $\gamma$-Al2O3 surface can trap and subsequently
decompose DIMP due to strong electrostatic attractions between the phosphoryl
oxygen and surface Al atoms.

##  4 Conclusions

In this work, we have harnessed large-scale ab initio molecular dynamics
calculations to investigate the adsorption and decomposition of DIMP on the
$\gamma$-alumina surface over a wide range of temperatures. Our DFT-based
molecular dynamics calculations predict a spontaneous decomposition of DIMP
with propene as the main by-product within the 700 – 1000∘ C temperature range
(the decomposition reaction leads to propene and an Al–OCHCH3CH3 adsorbate at
700∘ C). Due to the short-time scales inherent to BOMD simulations, it is
likely that a similar decomposition would also occur at lower temperatures but
would take longer to happen. Well-tempered metadynamics AIMD simulations were
performed for temperatures ranging from 200 – 1000∘ C to provide atomistic-
level details of the reaction path and associated energetics of DIMP
decomposition. Our metadynamics calculations also reveal that the free energy
barrier value decreases with temperature. The low free energy barrier values
at higher temperatures suggest that that reaction is extremely fast at higher
temperatures and is likely to occur at a lower rate at lower temperatures.

In our study, we obtained free energy values of 62.75 to 5.24 kJ/mol within
the 200 – 1000∘ C temperature range. Recent experiments on CWAs have also been
carried out under similar temperatures, including vapor phase decomposition of
DMMP (a sarin surrogate) at 500 – 800 K,58 sarin decomposition on TiO2
nanoparticles at 1000 K,53 and decomposition of CWAs at 2000 K.31, 59, 60
Collectively, all of these experiments showed that the temperature ranges
studied in this work can also be achieved under operational conditions, and
DIMP would decompose on the $\gamma$-Al2O3 surface. Due to the structural
similarity between DIMP and sarin, our calculations provide additional insight
into decomposition mechanisms of both these molecules and elucidate atomic
details of sarin decomposition on candidate metal-oxide surfaces. As a final
remark, this study serves as a convincing demonstration of the use of DFT-
based molecular dynamics simulations for investigating the interactions of
CWAs with existing metal-oxides, which can be used to guide experimental
efforts on these hazardous compounds.

The project or effort depicted was or is sponsored by the Department of the
Defense, Defense Threat Reduction Agency under the Materials Science in
Extreme Environments University Research Alliance, HDTRA1-20-2-0001. The
content of the information does not necessarily reflect the position or the
policy of the federal government, and no official endorsement should be
inferred.

Additional materials on free energy profiles, fluctuations of collective
variables, convergence tests of metadynamics calculations, and fluctuations of
the C–O bond distance in DIMP from metadynamics calculations.

## References

  * Kim et al. 2011 Kim, K.; Tsay, O. G.; Atwood, D. A.; Churchill, D. G. Destruction and Detection of Chemical Warfare Agents. _Chem. Rev._ 2011, _111_ , 5345–5403
  * Picard et al. 2019 Picard, B.; Chataigner, I.; Maddaluno, J.; Legros, J. Introduction to chemical warfare agents, relevant simulants and modern neutralisation methods. _Org. Biomol. Chem._ 2019, _17_ , 6528–6537
  * Michalkova et al. 2007 Michalkova, A.; Paukku, Y.; Majumdar, D.; Leszczynski, J. Theoretical study of adsorption of tabun on calcium oxide clusters. _Chem. Phys. Lett._ 2007, _438_ , 72–77
  * Wagner et al. 2000 Wagner, G. W.; Koper, O. B.; Lucas, E.; Decker, S.; Klabunde, K. J. Reactions of VX, GD, and HD with Nanosize CaO: Autocatalytic Dehydrohalogenation of HD. _J. Phys. Chem. B_ 2000, _104_ , 5118–5123
  * Michalkova et al. 2004 Michalkova, A.; Ilchenko, M.; Gorb, L.; Leszczynski, J. Theoretical Study of the Adsorption and Decomposition of Sarin on Magnesium Oxide. _J. Phys. Chem. B_ 2004, _108_ , 5294–5303
  * Wagner et al. 1999 Wagner, G. W.; Bartram, P. W.; Koper, O.; Klabunde, K. J. Reactions of VX, GD, and HD with Nanosize MgO. _J. Phys. Chem. B_ 1999, _103_ , 3225–3228
  * Mahato et al. 2009 Mahato, T.; Prasad, G.; Singh, B.; Acharya, J.; Srivastava, A.; Vijayaraghavan, R. Nanocrystalline zinc oxide for the decontamination of sarin. _J. Hazard. Mater._ 2009, _165_ , 928–932
  * Prasad et al. 2007 Prasad, G.; Mahato, T.; Singh, B.; Ganesan, K.; Pandey, P.; Sekhar, K. Detoxification reactions of sulphur mustard on the surface of zinc oxide nanosized rods. _J. Hazard. Mater._ 2007, _149_ , 460–464
  * Prasad et al. 2011 Prasad, G.; Ramacharyulu, P.; Singh, B.; Batra, K.; Srivastava, A. R.; Ganesan, K.; Vijayaraghavan, R. Sun light assisted photocatalytic decontamination of sulfur mustard using ZnO nanoparticles. _J. Mol. Catal. A: Chem._ 2011, _349_ , 55–62
  * Hirakawa et al. 2010 Hirakawa, T.; Sato, K.; Komano, A.; Kishi, S.; Nishimoto, C. K.; Mera, N.; Kugishima, M.; Sano, T.; Ichinose, H.; Negishi, N.; Seto, Y.; Takeuchi, K. Experimental Study on Adsorption and Photocatalytic Decomposition of Isopropyl Methylphosphonofluoridate at Surface of TiO2 Photocatalyst. _J. Phys. Chem. C_ 2010, _114_ , 2305–2314
  * Prasad et al. 2008 Prasad, G. K.; Mahato, T. H.; Singh, B.; Ganesan, K.; Srivastava, A. R.; Kaushik, M. P.; Vijayraghavan, R. Decontamination of sulfur mustard and sarin on titania nanotubes. _AIChE Journal_ 2008, _54_ , 2957–2963
  * Prasad et al. 2009 Prasad, G.; Singh, B.; Ganesan, K.; Batra, A.; Kumeria, T.; Gutch, P.; Vijayaraghavan, R. Modified titania nanotubes for decontamination of sulphur mustard. _J. Hazard. Mater._ 2009, _167_ , 1192–1197
  * Ramacharyulu et al. 2012 Ramacharyulu, P.; Prasad, G.; Ganesan, K.; Singh, B. Photocatalytic decontamination of sulfur mustard using titania nanomaterials. _J. Mol. Catal. A: Chem._ 2012, _353-354_ , 132–137
  * Sato et al. 2011 Sato, K.; Hirakawa, T.; Komano, A.; Kishi, S.; Nishimoto, C. K.; Mera, N.; Kugishima, M.; Sano, T.; Ichinose, H.; Negishi, N.; Seto, Y.; Takeuchi, K. Titanium dioxide photocatalysis to decompose isopropyl methylphosphonofluoridate (GB) in gas phase. _Appl. Catal. B: Environ._ 2011, _106_ , 316–322
  * Štengl et al. 2012 Štengl, V.; Grygar, T. M.; Opluštil, F.; Němec, T. Ge4+ doped TiO2 for stoichiometric degradation of warfare agents. _J. Hazard. Mater._ 2012, _227-228_ , 62–67
  * Wagner et al. 2008 Wagner, G. W.; Chen, Q.; Wu, Y. Reactions of VX, GD, and HD with Nanotubular Titania. _J. Phys. Chem. C_ 2008, _112_ , 11901–11906
  * Kuiper et al. 1976 Kuiper, A.; van Bokhoven, J.; Medema, J. The role of heterogeneity in the kinetics of a surface reaction: I. Infrared characterization of the adsorption structures of organophosphonates and their decomposition. _J. Catal._ 1976, _43_ , 154–167
  * Saxena et al. 2009 Saxena, A.; Sharma, A.; Srivastava, A. K.; Singh, B.; Gutch, P. K.; Semwal, R. P. Kinetics of adsorption of sulfur mustard on Al2O3 nanoparticles with and without impregnants. _J. Chem. Technol. Biotechnol._ 2009, _84_ , 1860–1872
  * Saxena et al. 2010 Saxena, A.; Srivastava, A. K.; Singh, B.; Gupta, A. K.; Suryanarayana, M. V.; Pandey, P. Kinetics of adsorptive removal of DEClP and GB on impregnated Al2O3 nanoparticles. _J. Hazard. Mater._ 2010, _175_ , 795–801
  * Wagner et al. 2001 Wagner, G. W.; Procell, L. R.; O’Connor, R. J.; Munavalli, S.; Carnes, C. L.; Kapoor, P. N.; Klabunde, K. J. Reactions of VX, GB, GD, and HD with Nanosize Al2O3. Formation of Aluminophosphonates. _J. Am. Chem. Soc._ 2001, _123_ , 1636–1644
  * Walenta et al. 2020 Walenta, C. A.; Xu, F.; Tesvara, C.; O’Connor, C. R.; Sautet, P.; Friend, C. M. Facile Decomposition of Organophosphonates by Dual Lewis Sites on a Fe3O4(111) Film. _J. Phys. Chem. C_ 2020, _124_ , 12432–12441
  * Trotochaud et al. 2017 Trotochaud, L.; Tsyshevsky, R.; Holdren, S.; Fears, K.; Head, A. R.; Yu, Y.; Karslıoğlu, O.; Pletincx, S.; Eichhorn, B.; Owrutsky, J.; Long, J.; Zachariah, M.; Kuklja, M. M.; Bluhm, H. Spectroscopic and Computational Investigation of Room-Temperature Decomposition of a Chemical Warfare Agent Simulant on Polycrystalline Cupric Oxide. _Chem. Mater._ 2017, _29_ , 7483–7496
  * Digne et al. 2004 Digne, M.; Sautet, P.; Raybaud, P.; Euzen, P.; Toulhoat, H. Use of DFT to achieve a rational understanding of acid–basic properties of $\gamma$-alumina surfaces. _J. Catal._ 2004, _226_ , 54–68
  * Gu et al. 2018 Gu, J.; Wang, J.; Leszczynski, J. Structure and Energetics of (111) Surface of $\gamma$-Al2O3: Insights from DFT Including Periodic Boundary Approach. _ACS Omega_ 2018, _3_ , 1881–1888
  * Roy et al. 2012 Roy, S.; Mpourmpakis, G.; Hong, D.-Y.; Vlachos, D. G.; Bhan, A.; Gorte, R. J. Mechanistic Study of Alcohol Dehydration on $\gamma$-Al2O3. _ACS Catal._ 2012, _2_ , 1846–1853
  * Christiansen et al. 2013 Christiansen, M. A.; Mpourmpakis, G.; Vlachos, D. G. Density Functional Theory-Computed Mechanisms of Ethylene and Diethyl Ether Formation from Ethanol on $\gamma$-Al2O3(100). _ACS Catal._ 2013, _3_ , 1965–1975
  * Bermudez 2007 Bermudez, V. M. Quantum-Chemical Study of the Adsorption of DMMP and Sarin on $\gamma$-Al2O3. _J. Phys. Chem. C_ 2007, _111_ , 3719–3728
  * Bailin et al. 1975 Bailin, L. J.; Sibert, M. E.; Jonas, L. A.; Bell, A. T. Microwave decomposition of toxic vapor simulants. _Environ. Sci. Technol._ 1975, _9_ , 254–258
  * Gibson and Sibener 2018 Gibson, K. D.; Sibener, S. J. Fate of Some Chemical Warfare Simulants Adsorbed on an Inert Surface when Exposed to Rapid Laser Initiated Heating. _J. Phys. Chem. C_ 2018, _122_ , 24684–24689
  * Thompson et al. 2019 Thompson, R. S.; Brann, M. R.; Purdy, E. H.; Graham, J. D.; McMillan, A. A.; Sibener, S. J. Rapid Laser-Induced Temperature Jump Decomposition of the Nerve Agent Simulant Diisopropyl Methylphosphonate under Atmospheric Conditions. _J. Phys. Chem. C_ 2019, _123_ , 21564–21570
  * Senyurt et al. 2021 Senyurt, E. I.; Schoenitz, M.; Dreizin, E. L. Rapid destruction of sarin surrogates by gas phase reactions with focus on diisopropyl methylphosphonate (DIMP). _Def. Technol._ 2021, _17_ , 703–714
  * Glaude et al. 2002 Glaude, P.; Melius, C.; Pitz, W.; Westbrook, C. Detailed chemical kinetic reaction mechanisms for incineration of organophosphorus and fluoroorganophosphorus compounds. _Proc. Combust. Inst._ 2002, _29_ , 2469–2476
  * VandeVondele et al. 2005 VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. _Comput. Phys. Commun._ 2005, _167_ , 103–128
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Grimme 2006 Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. _J. Comput. Chem._ 2006, _27_ , 1787–1799
  * Goedecker et al. 1996 Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. _Phys. Rev. B_ 1996, _54_ , 1703–1710
  * Hartwigsen et al. 1998 Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. _Phys. Rev. B_ 1998, _58_ , 3641–3662
  * Cholewinski et al. 2018 Cholewinski, M. C.; Dixit, M.; Mpourmpakis, G. Computational Study of Methane Activation on $\gamma$-Al2O3. _ACS Omega_ 2018, _3_ , 18242–18250
  * Dixit et al. 2018 Dixit, M.; Kostetskyy, P.; Mpourmpakis, G. Structure–Activity Relationships in Alkane Dehydrogenation on $\gamma$-Al2O3: Site-Dependent Reactions. _ACS Catal._ 2018, _8_ , 11570–11578
  * VandeVondele and Hutter 2003 VandeVondele, J.; Hutter, J. An efficient orbital transformation method for electronic structure calculations. _J. Chem. Phys._ 2003, _118_ , 4365–4369
  * Kolafa 2004 Kolafa, J. Time-reversible always stable predictor–corrector method for molecular dynamics of polarizable molecules. _J. Comput. Chem._ 2004, _25_ , 335–342
  * Kühne et al. 2007 Kühne, T. D.; Krack, M.; Mohamed, F. R.; Parrinello, M. Efficient and Accurate Car-Parrinello-like Approach to Born-Oppenheimer Molecular Dynamics. _Phys. Rev. Lett._ 2007, _98_ , 066401
  * Nosé 1984 Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. _J. Chem. Phys._ 1984, _81_ , 511–519
  * Hoover 1985 Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. _Phys. Rev. A_ 1985, _31_ , 1695–1697
  * Digne et al. 2002 Digne, M.; Sautet, P.; Raybaud, P.; Euzen, P.; Toulhoat, H. Hydroxyl Groups on $\gamma$-Alumina Surfaces: A DFT Study. _J. Catal._ 2002, _211_ , 1–5
  * Wilson 1979 Wilson, S. The dehydration of boehmite, $\gamma$-AlOOH, to $\gamma$-Al2O3. _J. Solid State Chem._ 1979, _30_ , 247–255
  * Biswas et al. 2020 Biswas, S.; Kwon, H.; Barsanti, K. C.; Myllys, N.; Smith, J. N.; Wong, B. M. Ab initio metadynamics calculations of dimethylamine for probing pKb variations in bulk vs. surface environments. _Phys. Chem. Chem. Phys._ 2020, _22_ , 26265–26277
  * Barducci et al. 2008 Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. _Phys. Rev. Lett._ 2008, _100_ , 020603
  * Biswas and Wong 2021 Biswas, S.; Wong, B. M. Ab initio metadynamics calculations reveal complex interfacial effects in acetic acid deprotonation dynamics. _J. Mol. Liq._ 2021, _330_ , 115624
  * Mitchell et al. 1997 Mitchell, M. B.; Sheinker, V. N.; Mintz, E. A. Adsorption and Decomposition of Dimethyl Methylphosphonate on Metal Oxides. _J. Phys. Chem. B_ 1997, _101_ , 11192–11203
  * Tsyshevsky et al. 2019 Tsyshevsky, R.; Holdren, S.; Eichhorn, B. W.; Zachariah, M. R.; Kuklja, M. M. Sarin Decomposition on Pristine and Hydroxylated ZnO: Quantum-Chemical Modeling. _J. Phys. Chem. C_ 2019, _123_ , 26432–26441
  * Tsyshevsky et al. 2020 Tsyshevsky, R.; Head, A. R.; Trotochaud, L.; Bluhm, H.; Kuklja, M. M. Mechanisms of Degradation of Toxic Nerve Agents: Quantum-chemical Insight into Interactions of Sarin and Soman with Molybdenum Dioxide. _Surf. Sci._ 2020, _700_ , 121639
  * Tsyshevsky et al. 2021 Tsyshevsky, R.; McEntee, M.; Durke, E. M.; Karwacki, C.; Kuklja, M. M. Degradation of Fatal Toxic Nerve Agents on Dry TiO2. _ACS Appl. Mater. Interfaces_ 2021, _13_ , 696–705
  * Wischert et al. 2011 Wischert, R.; Copéret, C.; Delbecq, F.; Sautet, P. Dinitrogen: a selective probe for tri-coordinate Al “defect” sites on alumina. _Chem. Commun._ 2011, _47_ , 4890–4892
  * Zegers and Fisher 1998 Zegers, E.; Fisher, E. Gas-Phase Pyrolysis of Diisopropyl Methylphosphonate. _Combust. Flame_ 1998, _115_ , 230–240
  * Senyurt et al. 2020 Senyurt, E. I.; Schoenitz, M.; Dreizin, E. L. Rapid destruction of sarin surrogates by gas phase reactions with focus on diisopropyl methylphosphonate (DIMP). _Def. Technol._ 2020,
  * Yuan and Eilers 2019 Yuan, B.; Eilers, H. T-jump pyrolysis and combustion of diisopropyl methylphosphonate. _Combust. Flame_ 2019, _199_ , 69–84
  * Mukhopadhyay et al. 2021 Mukhopadhyay, S.; Schoenitz, M.; Dreizin, E. L. Vapor-phase decomposition of dimethyl methylphosphonate (DMMP), a sarin surrogate, in presence of metal oxides. _Def. Technol._ 2021, _17_ , 1095–1114
  * Neupane et al. 2018 Neupane, S.; Barnes, F.; Barak, S.; Ninnemann, E.; Loparo, Z.; Masunov, A. E.; Vasu, S. S. Shock Tube/Laser Absorption and Kinetic Modeling Study of Triethyl Phosphate Combustion. _J. Phys. Chem. A_ 2018, _122_ , 3829–3836
  * Mathieu et al. 2019 Mathieu, O.; Kulatilaka, W. D.; Petersen, E. L. Shock-tube studies of Sarin surrogates. _Shock Waves_ 2019, _29_ , 441–449

## TOC graphic


for maximal functions  associated to hypersurfaces in $\bR^3$ with
height $h<2:$ Part I, Transactions of the American Mathematical
Society, Volume 372, 2019,  No. 2, pages 1363-1406.

[4]
Greenblatt, M.,
Maximal averages over hypersurfaces  and the Newton polyhedron.
J. Funct. Anal.,  262,  (2012), 2314–2348.

[5]
Greenblatt, M.,
$L^p$ boundedness of maximal averages over hypersurfaces in $\RR^3$.
Trans. Amer. Math. Soc.,  365,  (2013), 1875–1900.

[6] Ikromov, I. A., Kempe, M.,  Müller, D.,
Estimates for maximal functions  associated to hypersurfaces in $\bR^3$  and related problems of harmonic analysis.  Acta Math. 204 (2010), 151–271.

[7]
Ikromov, I. A.,  Müller, D.,
On adapted coordinate systems.
Trans. Amer. Math. Soc., 363 (2011), no. 6, 2821–2848.

[8]
Ikromov, I. A.,  Müller, D.,
Fourier restriction for hypersurfaces in three dimensions and Newton polyhedra;
Annals of Mathematics Studies 194, Princeton University Press, Princeton and Oxford 2016; 260 pp.

[9]  Ikromov, I. A.; Usmanov, S. È. On the boundedness of maximal operators associated with hypersurfaces.
(Russian) Sovrem. Mat. Fundam. Napravl. 64 (2018), no. 4, 650–681. 42B25

[10]
Iosevich, A. and Sawyer, E.,
Oscillatory integrals and maximal averages over homogeneous surfaces.
Duke Math. J., 82 (1996), no. 1, 103–141.

[11]
Iosevich, A.,  Sawyer, E. and Seeger, A.,
On averaging operators associated with convex hypersurfaces of finite
J. Anal. Math., 79 (1999), 159–187.

[12]
Iosevich, A., Sawyer, E.,
Maximal averages over surfaces.
Adv. Math., 132 (1997), 46–119.

[13]
Nagel, A.,  Seeger, A.  and Wainger, S.,
Averages over convex hypersurfaces.
Amer. J. Math., 115 (1993), no. 4, 903–927.

[14]
Phong, D.H. and  Stein, E.M.,
The Newton polyhedron and oscillatory integral operators.
Acta Math., 179 (1997), no. 1, 105–152.

[15] Phong, D. H., Stein, E. M.,  Sturm, J. A.,
On the growth and stability of real-analytic functions.
Amer. J. Math., 121 (1999), no. 3, 519-554.

[16]
Stein, E.M.,
Maximal functions. I. Spherical means.
Proc. Nat. Acad. Sci. U.S.A., 73 (1976), no. 7, 2174–2175.

[17]
Stein, E. M.,
Harmonic analysis: Real-variable methods, orthogonality, and
oscillatory integrals. Princeton Mathematical Series 43.
Princeton University Press, Princeton, NJ, 1993.

[18]
Varchenko, A. N.,
Newton polyhedra and estimates of oscillating integrals.
Funkcional. Anal. i Priložen, 10 (1976), 13–38 (Russian); English translation in
Funktional Anal. Appl., 18 (1976), 175–196.

[19]
Zimmermann, E.,
On $L^p$-estimates for maximal averages over hypersurfaces
not satisfying the transversality condition.
Doctoral thesis, Kiel 2014.


# Complex magnetic phase diagram with a small phase pocket in a three-
dimensional frustrated magnet CuInCr4S8

M. Gen<EMAIL_ADDRESS>Institute for Solid State Physics, The
University of Tokyo, Kashiwa 277-8581, Japan  Department of Advanced Materials
Science, University of Tokyo, Kashiwa 277-8561, Japan  H. Ishikawa  Institute
for Solid State Physics, The University of Tokyo, Kashiwa 277-8581, Japan  A.
Ikeda  Institute for Solid State Physics, The University of Tokyo, Kashiwa
277-8581, Japan  A. Miyake  Institute for Solid State Physics, The University
of Tokyo, Kashiwa 277-8581, Japan  Z. Yang  Institute for Solid State Physics,
The University of Tokyo, Kashiwa 277-8581, Japan  Y. Okamoto  Department of
Applied Physics, Nagoya University, Nagoya 464-8603, Japan  M. Mori
Department of Applied Physics, Nagoya University, Nagoya 464-8603, Japan  K.
Takenaka  Department of Applied Physics, Nagoya University, Nagoya 464-8603,
Japan  T. Kurumaji  Department of Advanced Materials Science, University of
Tokyo, Kashiwa 277-8561, Japan  Y. Tokunaga  Department of Advanced Materials
Science, University of Tokyo, Kashiwa 277-8561, Japan  T. Arima  Department of
Advanced Materials Science, University of Tokyo, Kashiwa 277-8561, Japan  M.
Tokunaga  Institute for Solid State Physics, The University of Tokyo, Kashiwa
277-8581, Japan  K. Kindo  Institute for Solid State Physics, The University
of Tokyo, Kashiwa 277-8581, Japan  Y. H. Matsuda  Institute for Solid State
Physics, The University of Tokyo, Kashiwa 277-8581, Japan  Y. Kohama
Institute for Solid State Physics, The University of Tokyo, Kashiwa 277-8581,
Japan

###### Abstract

Frustrated magnets with a strong spin-lattice coupling can show rich magnetic
phases and the associated fascinating phenomena. A promising platform is the
breathing pyrochlore magnet CuInCr4S8 with localized $S=3/2$ Cr3+ ions, which
is proposed to be effectively viewed as an $S=6$ Heisenberg antiferromagnet on
the face-centered-cubic lattice. Here, we unveil that CuInCr4S8 exhibits a
complex magnetic phase diagram with a small phase pocket (A-phase) by means of
magnetization, magnetostriction, magnetocapacitance, and magnetocaloric-effect
measurements in pulsed high magnetic fields of up to 60 T. Remarkably, the
appearance of A-phase is accompanied by anomalous magnetostrictive and
magnetocapacitive responses, suggesting the emergence of helimagnetism in
contrast to the neighboring commensurate magnetic phases. Besides, the high-
entropy nature is confirmed in the high-temperature side of A-phase. These
features are potentially related to a thermal fluctuation-driven multiple-$q$
state caused by the magnetic frustration, which has been theoretically
predicted but yet experimentally undiscovered.

##  I Introduction

Magnetic frustration is a key ingredient for the emergence of exotic magnetic
states in condensed matters. One promising route to induce strong frustration
is to arrange localized spins on a crystal lattice involving triangular
motifs, as represented by the Heisenberg antiferromagnet on the pyrochlore
lattice, a three-dimensional network of corner-sharing tetrahedra. In the
classical-spin picture, competing nearest-neighbor (NN) exchange interactions
arising from geometrical constraints prevent the development of a magnetic
long-range order even at low temperatures 1998_Moe . In real systems, however,
the lattice degrees of freedom often come into play to lift the macroscopic
degeneracy 2002_Tch ; 2004_Pen , yielding a rich magnetic phase diagram with
spin-lattice-coupled ordered states in the parameter space of temperature,
magnetic field, and pressure.

Chromium spinels ACr2X4 (A=Mg, Zn, Cd, Hg; X=O, S, Se) are known as typical
insulating frustrated magnets to realize versatile magnetic phenomena caused
by the interplay between the magnetic frustration and the spin-lattice
coupling (SLC): e.g., a 1/2-magnetization plateau 2006_Ued ; 2008_Koj ;
2011_Miy ; 2014_Miy , negative thermal expansion 2007_Hem ; 2013_Kit ;
2019_Ros , and spin-driven electric polarization 2006_Web ; 2008_Mur ;
2020_Ros . Here, tetrahedrally coordinated nonmagnetic A2+ ions and
octahedrally coordinated $S=3/2$ Cr3+ ions form a diamond and pyrochlore
network, respectively. The NN Cr-Cr exchange interaction arises from the
antiferromagnetic (AFM) direct exchange due to the overlapping of t2g orbitals
and the ferromagnetic (FM) superexchange via Cr-X-Cr path. In oxides, the
former contribution is predominant because of the short Cr-Cr length 2008_Yar
, so that the system can be approximated to the NN Heisenberg antiferromagnet
with geometrical frustration. In sulfides and selenides, by contrast, the
superexchange contribution exceeds the direct one, and further-neighbor
exchange interactions via Cr-X-X-Cr and Cr-X-A-X-Cr paths also become strong
2008_Yar , so that the magnetic properties are governed by bond frustration
2011_Tsu .

Recently, A-site ordered Cr spinels AA’Cr4X8 (A=Li, Cu; A’=Ga, In), where
nonmagnetic monovalent A+ and trivalent A’3+ cations are regularly arranged,
have attracted new attention in terms of the ground state control 2013_Oka ;
2014_Tan ; 2015_Nil ; 2016_Sah ; 2019_Gen ; 2020_Kan ; 2018_Oka ; 2020_Gen ;
2021_Gao ; 2020_Pok ; 2019_Gho . The difference in their chemical pressure
renders the Cr4 tetrahetra expanding and contracting alternately, resulting in
the so-called breathing pyrochlore lattice with nonequivalent NN exchange
interactions $J$ and $J^{\prime}$ in the small and large tetrahedra,
respectively [Fig. 1(a)]. A large majority of previous studies have focused on
Li(Ga, In)Cr4O8 with $J>J^{\prime}>0$ (AFM) 2013_Oka ; 2014_Tan ; 2015_Nil ;
2016_Sah ; 2019_Gen ; 2020_Kan , where the system can be viewed as an
intermediate state between the regular pyrochlore lattice ($J^{\prime}/J=1$)
and the isolated tetrahedral clusters ($J^{\prime}/J=0$). For small values of
$J^{\prime}/J$, the spin-gapped behavior associated with the tetramer singlet
formation rather than the classical spin-liquid behavior can emerge, as
observed in LiInCr4O8 2013_Oka ; 2014_Tan ; 2015_Nil . Furthermore, it has
been theoretically suggested that a variety of magnetic long-range orders can
appear depending on the value of $J^{\prime}/J$ and the strength of the SLC
2019_Aoy ; 2021_Aoy , which seems relevant to the magnetic structures on
Li(Ga, In)Cr4O8 at low temperatures 2015_Nil ; 2016_Sah .

In the breathing pyrochlore Cr spinels, not only the magnitude but the sign of
$J$ and $J^{\prime}$ can be flexibly tuned by the chemical substitution, which
can bring about a dramatic modification in the magnetic model rather than a
perturbative change. Of particular interest is CuInCr4S8 2018_Oka ; 2020_Gen ;
2021_Gao , where the special condition of $J>0$ (AFM) and $J^{\prime}<0$ (FM)
is satisfied due to the delicate balance between the direct and superexchange
interactions [Fig. 1(b)]; the exchange parameters at room temperature are
estimated to $J=14.7$ K, $J^{\prime}=-26.0$ K, $J_{2}=1.1$ K, $J_{\rm 3a}=6.4$
K, and $J_{\rm 3b}=4.5$ K by the density-functional-theory calculation [Fig.
1(c)] 2019_Gho . Since the FM $J^{\prime}$ does not compete with any other
exchange interactions listed above, magnetic moments on each large Cr4
tetrahedron are aligned ferromagnetically at low temperatures, as evidenced by
the observation of the cluster excitation using inelastic neutron scattering
2021_Gao . The magnetism of CuInCr4S8 can hence be mapped onto $S=6$ spins on
the face-centered-cubic (FCC) lattice with the NN AFM interaction $J_{\rm
t1}\equiv J+4J_{2}+2(J_{\rm 3a}+J_{\rm 3b})$ [see Fig. 1(c)]. Since the
geometrical frustration remains in this effective FCC lattice, a robust
1/2-magnetization plateau with a 3-up-1-down magnetic structure (C-phase)
appears via the exchange-striction mechanism 2020_Gen , as in Cr spinel oxides
2006_Ued ; 2008_Koj ; 2011_Miy ; 2014_Miy ; 2019_Gen . On the other hands, the
presence of an intermediate-field phase (Y-phase) in a wide field range
between the low-field AFM phase (X-phase) and C-phase is a distinctive feature
of this compound, which seems to owe to the peculiar combination of the
emergent FCC lattice and the SLC 2020_Gen . The previous focus on the in-field
properties of CuInCr4S8 was limited to the ground state 2018_Oka ; 2020_Gen ,
while the effect of thermal fluctuations on the stability of C- and Y-phases
have remained unexplored.

Figure 1: (a) Breathing pyrochlore lattice comprised of magnetic ions with two
kinds of the NN exchange interactions, $J$ and $J^{\prime}$, in the small and
large tetrahedra, respectively, and further-neighbor exchange interactions up
to the 3rd-NN. (b) Circular representation of the magnetic models of several
breathing pyrochlore Cr spinels as a function of $J$ and $J^{\prime}$
normalized to ${\overline{J}}\equiv\sqrt{J^{2}+J^{\prime 2}}$, which are
obtained from the previous DFT calculation 2019_Gho . (c) Effective magnetic
model of CuInCr4S8 mapped on the face-centered-cubic lattice with the NN AFM
exchange interaction $J_{\rm t1}\equiv J+4J_{2}+2(J_{\rm 3a}+J_{\rm 3b})$ and
the 2nd- and 3rd-NN interactions, $J_{\rm t2}$ and $J_{\rm t3}$, respectively.
The illustrations in (a) and (c) are drawn by VESTA software 2011_Mom .

Here, we report on a fascinating magnetic-field-versus-temperature ($H$–$T$)
phase diagram of CuInCr4S8 [Fig. 6(a)], revealed by means of magnetization,
magnetostriction, magnetocapacitance, and magnetocaloric-effect (MCE)
measurements in pulsed high magnetic fields of up to 60 T. The most striking
finding is the existence of a small phase-pocket (A-phase), which extends from
the phase boundary with the paramagnetic phase towards lower temperatures
instead of Y-phase. Remarkably, the transition from X- to A-phase is
accompanied by anomalous magnetostrictive and magnetocapacitive responses.
These observations suggest the appearance of helimagnetism in A-phase through
a commensurate-to-incommensurate phase transition. Moreover, the high-entropy
nature is confirmed in the high-temperature side of A-phase. Given the
similarity of the identified $H$–$T$ phase diagram to that of a number of
skyrmion hosts 2009_Muh ; 2012_Sek ; 2021_Tok , the emergence of a
multiple-$q$ vortex-like magnetic structure is also possible. We propose that
the unique magnetic properties of CuInCr4S8 could be ascribed to the
cooperative effect of the strong SLC and the breathing bond-alternation
introducing the FCC-lattice magnetism.

##  II Experimental methods

The magnetic ordering temperatures of CuInCr4S8 reported so far range from 26
K to 35 K 2021_Gao ; 1970_Pin ; 1980_Plu ; 2016_Ami . Our previous
polycrystalline samples (Sample #2) exhibit a magnetic phase transition at
$T_{\rm N}\approx 32$ K 2018_Oka , which is slightly different from any values
reported by other groups. We prepared new polycrystalline batches (Sample #1)
by the solid state reaction as in Ref. 2018_Oka , and found a higher
transition temperature as $T_{\rm N}\approx 35$ K [Figs. 2(a) and 2(b)], in
accordance with that reported in Refs. 2021_Gao ; 1980_Plu .

To compare the sample quality between Samples #1 and #2, powder X-ray
diffraction (XRD) patterns were obtained using a commercial X-ray
diffractometer (SmartLab, Rigaku), confirming the cubic $F{\overline{4}}3m$
structure with no detectable impurity for both the samples. Since no obvious
difference was observed in their XRD patterns, it was impossible to evaluate
which sample contains less crystallographic disorder and/or defects in the
nonmagnetic Cu+ and In3+ sites and the ligand S2- sites. A series of high-
field experiments revealed that field-induced phase transitions are much
sharper in Sample #1 than in Sample #2, suggesting that Sample #1 is of higher
quality, i.e., the Cr-based pyrochlore network may be better breathing. The
results for Sample #1 are shown in the main text, while those for Sample #2
are shown in Supplemental Material Supple .

Magnetization up to 7 T was measured using a commercial magnetometer (MPMS,
Quantum Design). Magnetization up to $\sim$60 T was measured by the
conventional induction method using a coaxial pickup coil in a non-destructive
(ND) pulsed magnet ($\sim$4 ms duration) at the Institute for Solid State
Physics (ISSP), University of Tokyo, Japan. Longitudinal magnetostriction up
to $\sim$60 T was measured by the optical fiber-Bragg-grating (FBG) technique
in a ND pulsed magnet ($\sim$36 ms duration) at ISSP 2018_Ike . The fiber was
adhered to a rod-shaped sintered sample by epoxy Stycast1266. The distortion
$\Delta L/L$ was detected by the optical filter method with a resolution of
$\sim$10-6. The thermal expansion at zero field was measured using the same
setup during a natural temperature rise from 2 to 60 K. The absolute value of
$\Delta L/L$ is not calibrated, so that the displayed value of $\Delta L/L$
should be smaller than actual lattice constant change. Dielectric constant at
zero field was measured at a frequency of 10 kHz by using an LCR meter
(E4980A, Agilent) in a commercial cryostat equipped with a superconducting
magnet (PPMS, Quantum Design). Magnetocapacitance along the field direction
($E\parallel B$) up to $\sim$60 T was measured at a frequency of 50 kHz by
using a capacitance bridge (1615-A, General Radio) in a ND pulsed magnet
($\sim$36 ms duration) at ISSP 2020_Miy . Silver paste was painted on the top
and bottom surfaces to form electrodes for a disk-shaped sintered sample. Heat
capacity at zero field was measured by the thermal relaxation method using a
PPMS. The magnetocaloric effect up to $\sim$60 T was measured under both the
nonadiabatic and adiabatic conditions in a ND pulsed magnet ($\sim$36 ms
duration) at ISSP. A sensitive Au16Ge84 film thermometer was sputtered on the
surface of a disk-shaped sintered sample, which was calibrated by a commercial
Cernox 1030 thermometer 2013_Kih .

##  III Results and discussions

###  III.1  Magnetic properties of the low-field phase (X-phase)

We first mention notable aspects of the low-field phase (X-phase) of
CuInCr4S8. Cr spinels typically undergo a first-order phase transition into an
AFM 2-up-2-down collinear structure 2015_Nil ; 2016_Sah ; 2007_Mat ; 2009_Ji
or an incommensurate helical structure 2008_Mur ; 2011_Tsu  accompanied by the
crystal symmetry lowering due to the SLC. For CuInCr4S8, the appearance of
X-phase is associated with a drop in the magnetic susceptibility $M/H$ at
$T_{\rm N}=35$ K [Fig. 2(a)] and a peak in the specific heat $C$ at $T_{\rm
p}=34$ K [Fig. 2(b)]. A recent neutron diffraction study has revealed that the
crystal structure remains the cubic $F{\overline{4}}3m$ symmetry down to 2 K
2021_Gao . Our further experiments also support the uniqueness of X-phase;
neither the normalized thermal expansion $\Delta L/L_{\rm 2K}$ nor the
dielectric constant $\varepsilon^{\prime}$ exhibit any sharp anomalies at
$T_{\rm N}$ [Figs. 2(c) and 2(d)], unlike in other Cr spinels 2007_Hem ;
2013_Kit ; 2019_Ros ; 2016_Sah ; 2020_Kan . Furthermore, no hysteretic
behavior is observed in the $M/H$–$T$ curve (not shown), indicating that the
phase transition at $T_{\rm N}$ is of the second order, and the SLC does not
seem to be responsible for the appearance of X-phase.

Figure 2: (a) Magnetization measured at 7 T after zero-field cooling. (b)
Specific heat divided by temperature at 0 T. (c) Thermal expansion compared to
the 2-K data at 0 T. (d) Dielectric constant at a frequency of 10 kHz at 0 T.
All the measurements were performed in a warming process.

In order to understand the characteristic of X-phase, it is useful to consider
the $J_{\rm t1}$-$J_{\rm t2}$-$J_{\rm t3}$ model on the effective FCC lattice,
where $J_{\rm t2}$ and $J_{\rm t3}$ stand for the 2nd- and 3rd-NN interactions
in the effective FCC model [Fig. 1(a)]. Theoretically, two kinds of
commensurate AFM order, Type-I with q=(1 0 0) and Type-III with q=(1 1/2 0),
appear at zero field for $J_{\rm t2}<4J_{\rm t3}$ and $J_{\rm t2}>4J_{\rm t3}$
respectively, if $|J_{\rm t2}|,|J_{\rm t3}|\ll J_{\rm t1}$ 2020_Bal . Note
that $J_{\rm t2}$ and $J_{\rm t3}$ arise from the 4th or higher order NN
exchange interactions in the original breathing pyrochlore lattice, and should
be much weaker than $J_{\rm t1}$ in CuInCr4S8. Indeed, Type-I order has been
confirmed by powder neutron diffraction, and the observed magnon dispersion is
accountable by incorporating FM $J_{\rm t2}$ ($<0$) 2021_Gao ; 1971_Plu .
Although the local spin configuration of X-phase has not been fully
identified, an “all-in-all-out”-like one 1971_Plu  and a coplanar one in which
the two AFM spin pairs are canted to each other 2021_Gao  have been proposed
as candidates rather than the 2-up-2-down collinear one. This is curious
because the SLC and thermal fluctuations tend to favor a collinear spin
configuration out of the continuous degeneracy regarding the directions of
four inequivalent spins in Type-I order 1987_Hen ; 2015_Ben . Moreover,
neither of the proposed magnetic structures can be explained just by
considering the Dzyaloshinskii-Moriya (DM) interactions which should exist to
a greater or lesser extent in the breathing pyrochlore lattice 1957_Dzy ;
1960_Mor ; 2005_Elh . We note that a small singularity at $T_{\rm p}$ and a
hump structure around 10 K in the $C/T$–$T$ curve [Fig. 2(b)] suggest the
existence of a large magnetic entropy even below $T_{\rm N}$ Supple . These
features are unique to CuInCr4S8 among Cr spinels 2013_Kit ; 2013_Oka ;
2018_Oka ; 2013_Kem  and would stem from the magnetic frustration inherent in
the effective FCC-lattice Heisenberg antiferromagnet.

Figure 3: (a) Magnetization, (b) its field-derivative, (c) longitudinal
magnetostriction, and (d) its field-derivative as a function of magnetic field
measured at various initial temperatures $T_{\rm ini}$. The thick (thin) lines
correspond to the data in field-elevating (descending) processes. All the data
except for the lowest-temperature one are shifted upward for clarity. The
inset of (b) is a magnetization-versus-temperature plot at 35 T. The inset of
(c) shows the enlarged view of the magnetostriction curves in field-elevating
processes for 10 K $\leq T_{\rm ini}\leq$ 30 K, where lattice shrinkage occurs
as the precursor of the phase transition from X- to A- phase at $H_{\rm A1}$.
The filled (open) triangles in (a) and (c) denote the phase boundaries,
$H_{\rm A1}$ and $H_{\rm A2}$ ($H_{\rm c1}$, $H_{\rm c2}$, and $H_{\rm P}$),
which are determined from the $dM/dH$ and $dL/dH$ anomalies, respectively. The
dashed gray lines in (b) and (d) trace the boundaries of A-phase.

###  III.2  Temperature dependence of magnetization, magnetostriction, and
magnetocapacitance up to 60 T

In order to clarify the in-field properties of CuInCr4S8, we measured
magnetization, longitudinal magnetostriction, and magnetocapacitance in pulsed
high magnetic fields of up to 60 T at various initial temperatures ($T_{\rm
ini}$). All the magnetization and magnetostriction data are summarized in Fig.
3, and some of them are arranged along with the magnetocapacitance data for
several selected $T_{\rm ini}$ in Fig. 4. Note that the observed results are
irreversible and hysteretic especially for low $T_{\rm ini}$ because of
significant sample-temperature changes in the field-descending process, as
discussed in Sec. III.4. We will hence mainly focus on the data in the field-
elevating process, which are useful for constructing the $H$–$T$ phase diagram
(the discussions on the field-descending process are given in Supplementary
Materials Supple ).

For the lowest $T_{\rm ini}$ of 1.4 K, two-step metamagnetic transitions take
place at $\mu_{0}H_{\rm c1}=32$ T and $\mu_{0}H_{\rm c2}=56$ T [Fig. 3(a)],
followed by a 1/2-magnetization plateau which continues up to $\mu_{0}H_{\rm
c3}\approx 110$ T 2020_Gen . Correspondingly, a substantial lattice expansion
is observed at $H_{\rm c1}$ and $H_{\rm c2}$ [Fig. 3(c)], suggesting the
importance of the SLC on these phase transitions. As discussed in Ref.
2020_Gen , the overall magnetization curve can be roughly reproduced by the
minimal Hamiltonian incorporating the SLC, which effectively produces the
biquadratic exchange term $-({\mathbf{S}}_{i}\cdot{\mathbf{S}}_{j})^{2}$
between NN spins 2004_Pen ; 2021_Aoy . According to this, the spin
configurations within each small tetrahedron in the low-field X-phase, the
intermediate-field Y-phase between $H_{\rm c1}$ and $H_{\rm c2}$, and C-phase
with the 1/2-magnetization plateau can be assigned to the canted 2:2, canted
2:1:1, and collinear 3-up-1-down state, respectively [Fig. 6(b)]. As $T_{\rm
ini}$ increases, each critical field is gradually decreased. The hump
structures at $H_{\rm c2}$ in $dM/dH$ and $dL/dH$ become blunt for $T_{\rm
ini}\geq 20$ K [Figs. 3(b) and 3(d)], making the phase boundary indiscernible.
For $T_{\rm ini}\geq 35$ K (=$T_{\rm N}$), the transition from the
paramagnetic to a field-induced ordered phase is observed. The critical field
$H_{\rm P}$ becomes higher on increasing $T_{\rm ini}$ up to 50 K, suggesting
that C-phase is stabilized by thermal fluctuations and may extend to higher
temperatures. A moderate magnetocapacitance change is also seen at the
transition into C-phase for $T_{\rm ini}\geq 10$ K [Figs. 4(f), 4(i), 4(l),
and 4(o)], which may be attributed to the local lattice distortion to
stabilize the 3-up-1-down magnetic structure via the exchange-striction
mechanism.

Figure 4: Top, middle, and bottom panels show the magnetic-field dependence of
magnetization, longitudinal magnetostriction, and magnetocapacitance,
respectively, for (a)-(c) $T_{\rm ini}$ = 4.2 K, (d)-(f) 10 K, (g)-(i) 15 K,
(j)-(l) 25 K, and (m)-(o) 40 K. The field-derivatives of magnetization and
magnetostriction are also shown by gray lines in the right axes. The
magnetocapacitance was measured along the field direction at a frequency of 50
kHz. The thick (thin) lines correspond to the data in field-elevating
(descending) processes. The phase boundaries are roughly drawn by dashed black
lines, which are determined from the data in the field-elevating process.

###  III.3  Magnetostrictive and magnetocapacitive effects in a newly
discovered phase (A-phase)

In addition to the previously reported X-, Y-, and C-phases, a novel field-
induced phase (A-phase) is found in the series of high-field experiments. Let
us here focus on the magnetization and magnetostriction curves for an
intermediate $T_{\rm ini}$ range between 10 and 35 K. Taking a closer look at
Fig. 3(a), one can notice that the lowest-field metamagnetic transition
becomes exceptionally steep for $T_{\rm ini}=10$ K, which is clearly shown as
the peak sharpness in the $dM/dH$ curve [Fig. 3(b)]. This is reflected in the
plot of magnetization-versus-temperature at 35 T showing a maximum at 10 K
[inset of Fig. 3(b)]. Besides, a dip structure of $\Delta L/L_{\rm 0T}$
appears in a narrow field range for 10 K $\leq T_{\rm ini}\leq$ 30 K [inset of
Fig. 3(c)]. These features suggest that the nature of the lowest-field
metamagnetic transition changes across $T_{\rm ini}=10$ K. Accordingly, the
critical field above 10 K should be distinguished from $H_{\rm c1}$ and
instead termed $H_{\rm A1}$. Importantly, another tiny $dM/dH$ peak appears at
around the midpoint between $H_{\rm A1}$ and $H_{\rm c2}$ (termed $H_{\rm
A2}$) for $T_{\rm ini}=10$ K, developing toward higher $T_{\rm ini}$ up to
32.5 K [Fig. 3(b)]. An anomaly is also visible at $H_{\rm A2}$ in the $\Delta
L/L_{\rm 0T}$ data for 20 K $\leq T_{\rm ini}\leq$ 30 K [Fig. 3(c)]. These
observations indicate the existence of another phase intervened between X- and
Y-phases.

According to the microscopic magnetoelastic theory 2019_Ros , the classical
Heisenberg antiferromagnet on the regular pyrochlore lattice exhibits a volume
expansion, equivalently a positive magnetostriction, proportional to $M^{2}$
during the successive field-induced phase transitions until saturation, when
assuming a constant 16-sublattice magnetic unit cell. This tendency is roughly
seen in the $\Delta L/L_{\rm 0T}$ data of CuInCr4S8 below 10 K [Figs. 3(a) and
3(c)], where $\Delta L/L_{0T}$ follows a parabolic field-dependence within
X-phase and smoothly increases even across the phase transitions at $H_{\rm
c1}$ and $H_{\rm c2}$. This suggests that the magnetic unit cells in Y- and
C-phases remain unchanged as in X-phase. From this viewpoint, the observed dip
structure of $\Delta L/L_{\rm 0T}$ for 12 K $\leq T_{\rm ini}\leq$ 30 K is
unusual. The negative magnetostriction is not seen within A-phase, but occurs
as a precursor to the metamagnetic transition from X- to A-phase at $H_{\rm
A1}$ [Figs. 4(g), 4(h), 4(j), and 4(k)]. We infer that such a nonmonotonic
$\Delta L/L_{\rm 0T}$ behavior would arise from the modification of the
magnetic unit cell, i.e., the reconstruction of $q$-vector.

Interestingly, it is uncovered that A-phase is associated with a substantial
magnetodielectric coupling in addition to the anomalous magnetostrictive
behavior. No signature of a magnetocapacitance change is observed at $H_{\rm
c1}$ within our experimental resolution for $T_{\rm ini}=4.2$ K [Fig. 4(c)],
while the magnetocapacitance is suddenly enhanced by $\sim$0.3 % at $H_{\rm
A1}$ for $T_{\rm ini}=10$ K [Fig. 4(f)], which is exactly the characteristic
temperature of the appearance of A-phase. Such a magnetocapacitive response is
visible also for $T_{\rm ini}=15$ K and 25 K [Figs. 4(i) and 4(l)], signaling
the emergence of multiferroicity in A-phase originating from the helical spin
modulation 2005_Kat ; 2007_Ari .

###  III.4  Sample-temperature change in our experimental condition with
millisecond pulsed-field duration

All the aforementioned measurements were performed under the nonadiabatic
condition by immersing the sample in 4He liquid or gas. However, the heating
or cooling effects of the sample are usually non-negligible in millisecond
pulsed-field experiments 2021_Miy . Accordingly, we measured the sample
temperature $T(H)$ in pulsed magnetic fields of up to 60 T under the
nonadiabatic condition at $T_{\rm ini}=4.2$ K. Surprisingly, $T(H)$ was found
to drastically increase by more than 10 K at $H_{\rm c2}$ in the field-
elevating process [Fig. 5(a)]. This is possibly due to the hysteresis loss
accompanied by the first-order phase transition into C-phase. Since the
thermal equilibrium was not sufficiently achieved around the maximum field,
the sample was not cooled down below 10 K in the field-descending process.
Here, successive phase transitions more complicated than those from C-, Y-,
A-, to X-phase would take place, judging from at least three anomalies in the
$T(H)$ curve in the field-descending process Supple . Note that such a strong
heating effect does not seem to occur in the magnetization measurements at
$T_{\rm ini}=1.4$ and 4.2 K because the magnetization curves in the field-
descending process are less hysteretic and no double-hump structure is seen in
the $dM/dH$ curves around 25 T unlike for $T_{\rm ini}\geq 7.5$ K [Figs. 3(a),
(b), and 4(a)]. The difference in the degree of the heating effect might be
due to the sample form as well as the thermal bath environment; powder samples
were used for the magnetization measurements, whereas sintered samples for
other measurements. We stress that a nearly isothermal condition is confirmed
below $H_{\rm c2}$ for all the measurements at $T_{\rm ini}=4.2$ K. This would
be applicable also for $T_{\rm ini}>4.2$ K because the heating effect becomes
less pronounced for higher $T_{\rm ini}$ as shown next. Thus, it would be
reasonable to use all the data in the field-elevating process for mapping out
the $H$–$T$ phase diagram.

Figure 5: $T(H)$ curves obtained (a) under the nonadiabatic condition at
$T_{\rm ini}=4.2$ K and (b)(c) under the adiabatic condition at various
$T_{\rm ini}$. The thick (thin) lines correspond to the data in the field-
elevating (descending) process. In (b), the $T(H)$ curves in field-descending
processes for $T_{\rm ini}=5$ and 8 K completely overlap with each other.
Dashed gray lines are guides to the eye to make it easier to see the possible
phase transition from Y- to A-phase at $H_{\rm A0}$ (denoted by red
triangles).

###  III.5  Magnetic entropy changes through successive phase transitions

In order to extract information about the magnetic entropy ($S_{\rm M}$) in
each phase, we further investigated the magnetocaloric effect (MCE) under the
adiabatic condition. The MCE is the sample-temperature change induced by the
evolution of $S_{\rm M}$ upon the application or removal of an external
magnetic field. When there is no irreversible first-order transition, the
obtained $T(H)$ curve can be regarded as the isentropic curve. As a
consequence, an $S_{\rm M}$ increase (decrease) is detected as a $T(H)$
decrease (increase) 2013_Kih .

For $T_{\rm ini}\leq T_{\rm N}$, no apparent $T(H)$ change is observed below
$H_{\rm c1}$ or $H_{\rm A1}$ [Figs. 5(b) and 5(c)], indicating that there is
almost no $S_{\rm M}$ change in X-phase upon the application of a magnetic
field. Considering the magnetocaloric relation fulfilled under the isentropic
condition

$\left(\frac{\partial T}{\partial
H}\right)_{S}=-\frac{T}{C_{H}}\left(\frac{\partial M}{\partial T}\right)_{H},$
(1)

where $C_{H}$ is the specific heat at the constant magnetic field, this
observation is compatible with the $M/H$–$T$ data at a low field showing a
slight temperature dependence below $T_{\rm N}$ [inset of Fig. 3(b)]. For
$T_{\rm ini}=5$ and 8 K, however, the $T(H)$ curves show a dramatic increase
above $H_{\rm c1}$ and eventually overlap with each other around the maximum
field [Fig. 5(b)]. Those in field-descending processes are completely
irreversible, indicating that the contribution of the hysteresis loss is
larger than that of an intrinsic MCE in this temperature range. Importantly,
an additional $T(H)$ kink is observed at $H_{\rm A0}$ just above $H_{\rm c1}$,
presumably indicating the phase transition from Y- to A-phase [see dashed gray
lines in Fig. 5(b)]. By drawing the $T(H)$ curves in field-elevating processes
on the $H$–$T$ phase diagram, it can be seen that $T(H)$ significantly
increases within A-phase [Fig. 6(c)]. Bearing $(\partial M/\partial T)_{H}<0$
just above 10 K at around 35 T [inset of Fig. 3(b)] and Eq. (1) in mind, this
would reflect the intrinsic entropy reduction in the low-temperature region of
A-phase.

As $T_{\rm ini}$ increases, the $T(H)$ change at $H_{\rm A1}$ and the
accompanying hysteresis become smaller [Fig. 5(c)], assuring the achievement
of the quasi-isentropic condition. Notably, the $T(H)$ behaviors for 20 K
$\leq T_{\rm ini}\leq$ 32 K are qualitatively different from those for $T_{\rm
ini}=12$ and 16 K; the inflection point at $H_{\rm A1}$ disappears for the
former, indicating that $S_{\rm M}$ in A-phase is as high as that in X-phase
above 20 K (we note again that a large $S_{\rm M}$ remains in X-phase as
mentioned in Sec. III.1). This implies that A-phase would be further
subdivided into two phases, as shown in Fig. 6(c). At the transitions from Y-
and the paramagnetic phases to C-phase, $T(H)$ commonly starts rising for all
the $T_{\rm ini}$ because of the lowest-$S_{\rm M}$ nature in the 3-up-1-down
collinear magnetic structure [Figs. 5(b) and 5(c)]. These observations are
helpful for determining the phase boundary, especially between Y- and
C-phases.

Figure 6: (a) Magnetic-field-versus-temperature ($H$–$T$) phase diagram of
CuInCr4S8 derived from a series of experiments in the present work. (b)
Schematic of the magnetic structures in X-, Y-, and C-phases predicted by the
previous magnetoelastic theory 2020_Gen . (c) Adiabatic $T(H)$ curves in the
field-elevating process shown in Figs. 5(b) and 5(c) are plotted on the phase
diagram, focusing on the $H$–$T$ region where multiple phase transitions take
place. Judging from the magnetoentropic signatures, A-phase is possibly
further subdivided into two distinct phases in the low- (pink) and high-
temperature regions (red).

###  III.6  H–T phase diagram

Based on the above experimental results, we construct the $H$–$T$ phase
diagram of CuInCr4S8 for Sample #1 ($T_{\rm N}=35$ K), as shown in Fig. 6(a).
The newly discovered A-phase is found to appear in a closed $H$–$T$ regime
around 25-40 T and 10-35 K. It should be noted that a similar phase diagram
was identified for Sample #2 ($T_{\rm N}=32$ K), but A-phase seems to further
extend down to below 1.4 K Supple . These facts indicate that the appearance
of A-phase is an intrinsic property of CuInCr4S8, and its metastability is
quite sensitive against the sample quality.

Although the perfect theoretical description which can reproduce the $H$–$T$
phase diagram of CuInCr4S8 is elusive at this stage, we propose one possible
origin for the appearance of A-phase. In all the observed field-induced phase
transitions, local lattice distortions as well as the lattice constant change
should take place. This induces a modulation of the strengths of $J$ and
$J^{\prime}$ from bond to bond; e.g., the trigonal distortion produces three
$J(1+\delta_{1})$ and three $J(1-\delta_{2})$ per one small tetrahedron in the
3-up-1-down collinear state 2011_Kim . It has been theoretically shown that
the introduction of tetragonal or trigonal distortions into the FCC-lattice
Heisenberg antiferromagnet can induce incommensurate helical orders at zero
field 1972_Yam . This may be related to the observed in-field properties of
CuInCr4S8; the lattice distortion would lead to the competition between
commensurate Y- and incommensurate A-phases, and the latter is stabilized on
the higher temperature side by the entropy effect.

Finally, we note that the present assignment of A-phase is tentative and it
may be further subdivided into the low-$S_{\rm M}$ and high-$S_{\rm M}$
phases, as implied from the magnetoentropic signatures [Fig. 6(c)]. If this is
the case, the emergence of a multiple-$q$ state with a superposition of
multiple spin modulations, such as a skyrmion lattice (SkL), is possible in
the high-temperature side of A-phase. The SkL phase has been so far found in a
number of materials with a noncentrosymmetric crystal lattice, i.e., MnSi and
Cu2OSeO3 2009_Muh ; 2012_Sek ; 2021_Tok . There, the spin helix is induced by
the competition between the FM exchange and DM interactions, and the
triple-$q$ state is stabilized relative to the corresponding single-$q$ state
by thermal fluctuations. Indeed, the high-$S_{\rm M}$ nature of the SkL phase
has been experimentally demonstrated by analyzing the magnetization data based
on the magnetocaloric relation Eq. (1) 2018_Boc . Accordingly, the SkL phase
usually appears in a closed $H$–$T$ regime (often called “A-phase pocket”).
There also exist theoretical predictions on the realization of thermal
fluctuation-driven multiple-$q$ states caused by the magnetic frustration,
such as SkL (triple-$q$) in the $J_{1}$-$J_{3}$ model on the triangular
lattice 2012_Oku  and the meron/antimeron-like lattice (double-$q$) in the
$J_{1}$-$J_{2}$ model on the honeycomb lattice 2019_Shi_PRL ; 2019_Shi_PRB .
The proposed $H$–$T$ phase diagrams are rather complicated, similar to those
of the DM-induced skyrmion hosts accommodating the A-phase pocket. Although
several compounds such as NiGa2S4 2005_Nak  and Bi3Mn4O12(NO3) 2009_Smi  have
been raised as candidates for realizing these predictions, the requirement of
high-field experiments has challenged their verifications. Only recently,
MnSc2S4, characterized by the $J_{1}$-$J_{2}$-$J_{3}$ model on the
centrosymmetric diamond lattice, was found to exhibit a fractional
antiferromagnetic SkL, which seems to be attributed to the anisotropic
exchange couplings as well as the magnetic frustration and is stabilized even
at zero temperature 2020_Gao . Furthermore, it has been theoretically proposed
that the $J$-$J^{\prime}$-$J_{3}$ model on the breathing pyrochlore lattice
can host a quadrupole-$q$ hedgehog-lattice spin texture at zero field
2021_Aoy_HG . In light of these backgrounds, it would be an intriguing subject
whether or not the multiple-$q$ state is realized within A-phase of CuInCr4S8.
It is necessary to perform NMR or neutron scattering experiments on single
crystals, which are left for future work. Further theoretical investigations
are also desired to verify the proposed scenario.

##  IV Summary

In summary, we have revealed that the breathing pyrochlore magnet CuInCr4S8,
whose magnetism is governed by the magnetic frustration and the SLC, exhibits
an exotic $H$–$T$ phase diagram with a A-phase pocket. The appearance of
A-phase is found to be associated with anomalous magnetostrictive and
magnetocapacitive responses. This suggests the manifestation of helimagnetism
in contrast to the neighboring commensurate magnetic phases. Besides, the
high-entropy nature is confirmed in the high-temperature side of A-phase.
These features are reminiscent of the emergence of a thermal-fluctuations-
driven multiple-$q$ phase such as the SkL. To the best of our knowledge, there
have been no reports on the $H$–$T$ phase diagram dressed with the A-phase
pocket in insulating frustrated magnets. We believe that the present discovery
will attract more interest in breathing pyrochlore systems as well as open up
new avenues for exploring multiple-$q$ states.

## Acknowledgements

The authors are grateful for helpful discussions with S. Gao. The authors
thank S. Kitani for providing the data of the specific heat of CdCr2O4 shown
in Ref. 2013_Kit , which is useful for discussions Supple . This work was
financially supported by the JSPS KAKENHI Grants-In-Aid for Scientific
Research (No. 19H05823, No. 20J10988, and No. 20K20892) and Basic Science
Program No. 18-001 of Tokyo Electric Power Company (TEPCO) memorial
foundation. M.G. was supported by the JSPS through a Grant-in-Aid for JSPS
Fellows.

## References

  * (1) R. Moessner and J. T. Chalker, Low-temperature properties of classical geometrically frustrated antiferromagnets, Phys. Rev. B 58, 12049 (1998).
  * (2) O. Tchernyshyov, R. Moessner, and S. L. Sondhi, Spin-Peierls phases in pyrochlore antiferromagnets, Phys. Rev. B 66, 064403 (2002).
  * (3) K. Penc, N. Shannon, and H. Shiba, Half-Magnetization Plateau Stabilized by Structural Distortion in the Antiferromagnetic Heisenberg Model on a Pyrochlore Lattice, Phys. Rev. Lett. 93, 197203 (2004).
  * (4) H. Ueda, H. Mitamura, T. Goto, and Y. Ueda, Successive field-induced transitions in a frustrated antiferromagnet HgCr2O4, Phys. Rev. B 73, 094415 (2006).
  * (5) E. Kojima, A. Miyata, S. Miyabe, S. Takeyama, H. Ueda, and Y. Ueda, Full-magnetization of geometrically frustrated CdCr2O4 determined by Faraday rotation measurements at magnetic fields up to 140 T, Phys. Rev. B 77, 212408 (2008).
  * (6) A. Miyata, H. Ueda, Y. Ueda, H. Sawabe, and S. Takeyama, Magnetic Phases of a Highly Frustrated Magnet ZnCr2O4, up to an Ultrahigh Magnetic Field of 600 T, Phys. Rev. Lett. 107, 207203 (2011).
  * (7) A. Miyata, H. Ueda, and S. Takeyama, Canted 2:1:1 Magnetic Supersolid Phase in a Frustrated Magnet MgCr2O4 as a Small Limit of the Biquadratic Spin Interaction, J. Phys. Soc. Jpn. 83, 063702 (2014).
  * (8) J. Hemberger, H.-A. Krug von Nidda, V. Tsurkan, and A. Loidl, Large Magnetostriction and Negative Thermal Expansion in the Frustrated Antiferromagnet ZnCr2Se4, Phys. Rev. Lett. 98, 147203 (2007).
  * (9) S. Kitani, M. Tachibana, N. Taira, and H. Kawaji, Thermal study of the interplay between spin and lattice in CoCr2O4 and CdCr2O4, Phys. Rev. B 87, 064402 (2013).
  * (10) L. Rossi, A. Bobel, S. Wiedmann, R. Küchler, Y. Motome, K. Penc, N. Shannon, H. Ueda, and B. Bryant, Negative Thermal Expansion in the Plateau State of a Magnetically Frustrated Spinel, Phys. Rev. Lett. 123, 027205 (2019).
  * (11) S. Weber, P. Lunkenheimer, R. Fichtl, J. Hemberger, V. Tsurkan, and A. Loidl, Colossal Magnetocapacitance and Colossal Magnetoresistance in HgCr2S4, Phys. Rev. Lett. 96, 157202 (2006).
  * (12) H. Murakawa, Y. Onose, K. Ohgushi, S. Ishiwata, and Y. Tokura, Generation of Electric Polarization with Rotating Magnetic Field in Helimagnet ZnCr2Se4, J. Phys. Soc. Jpn. 77, 043709 (2008).
  * (13) L. Rossi, D. Brüning, H. Ueda, Y, Skourski, T. Lorenz, and B. Bryant, Magnetoelectric coupling in a frustrated spinel studied using high-field scanning probe microscopy, Appl. Phys. Lett. 116, 262901 (2020).
  * (14) A. N. Yaresko, Electronic band structure and exchange coupling constants in ACr2X4 spinels (A=Zn, Cd, Hg; X=O, S, Se), Phys. Rev. B 77, 115106 (2008).
  * (15) V. Tsurkan, S. Zherlitsyn, V. Felea, S. Yasin, Yu. Skourski, J. Deisenhofer, H.-A. Krug von Nidda, P. Lemmens, J. Wosnitza, and A. Loidl, Magnetostructural Transitions in a Frustrated Magnet at High Fields, Phys. Rev. Lett. 106, 247202 (2011).
  * (16) Y. Okamoto, G. J. Nilsen, J. P. Attfield, and Z. Hiroi, Breathing Pyrochlore Lattice Realized in A-Site Ordered Spinel Oxides LiGaCr4O8 and LiInCr4O8, Phys. Rev. Lett. 110, 097203 (2013).
  * (17) Y. Tanaka, M. Yoshida, M. Takigawa, Y. Okamoto, and Z. Hiroi, Novel Phase Transitions in the Breathing Pyrochlore Lattice: 7Li-NMR on LiInCr4O8 and LiGaCr4O8, Phys. Rev. Lett. 113, 227204 (2014).
  * (18) G. J. Nilsen, Y. Okamoto, T. Masuda, J. Rodriguez-Carvajal, H. Mutka, T. Hansen, and Z. Hiroi, Complex magnetostructural order in the frustrated spinel LiInCr4O8, Phys. Rev. B 91, 174435 (2015).
  * (19) R. Saha, F. Fauth, M. Avdeev, P. Kayser, B. J. Kennedy, and A. Sundaresan, Magnetodielectric effects in A-site cation-ordered chromate spinels LiMCr4O8 (M=Ga and In), Phys. Rev. B 94, 064420 (2016).
  * (20) M. Gen, D. Nakamura, Y. Okamoto, and S. Takeyama, Ultra-high magnetic field magnetic phases up to 130 T in a breathing pyrochlore antiferromagnet LiInCr4O8, J. Magn. Magn. Mater. 473, 387 (2019).
  * (21) T. Kanematsu, M. Mori, Y. Okamoto, T. Yajima, and K. Takenaka, Thermal Expansion and Volume Magnetostriction in Breathing Pyrochlore Magnets LiACr4X8 (A=Ga, In, X=O, S), J. Phys. Soc. Jpn. 89, 073708 (2020).
  * (22) Y. Okamoto, M. Mori, N. Katayama, A. Miyake, M. Tokunaga, A. Matsuo, K. Kindo, and K. Takenaka, Magnetic and Structural Properties of A-Site Ordered Chromium Spinel Sulfides: Alternating Antiferromagnetic and Ferromagnetic Interactions in the Breathing Pyrochlore Lattice, J. Phys. Soc. Jpn. 87, 034709 (2018).
  * (23) M. Gen, Y. Okamoto, M. Mori, K. Takenaka, and Y. Kohama, Magnetization process of the breathing pyrochlore magnet CuInCr4S8 in ultrahigh magnetic fields up to 150 T, Phys. Rev. B 101, 054434 (2020).
  * (24) S. Gao, A. F. May, M.-H. Du, J. A. M. Paddison, H. S. Arachchige, G. Pokharel, C. dela Cruz, Q. Zhang, G. Ehlers, D. S. Parker, D. G. Mandrus, M. B. Stone, and A. D. Christianson, Hierarchical excitations from correlated spin tetrahedra on the breathing pyrochlore lattice, Phys. Rev. B 103, 214418 (2021).
  * (25) G. Pokharel, H. S. Arachchige, T. J. Williams, A. F. May, R. S. Fishman, G. Sala, S. Calder, G. Ehlers, D. S. Parker, T. Hong, A. Wildes, D. Mandrus, J. A. M. Paddison, and A. D. Christianson, Cluster Frustration in the Breathing Pyrochlore Magnet LiGaCr4S8, Phys. Rev. Lett. 125, 167201 (2020).
  * (26) P. Ghosh, Y. Iqbal, T. Müller, R. Thomale, J. Reuther, M. J. P. Gingras, and H. O. Jeschke, Breathing chromium spinels: a showcase for a variety of pyrochlore Heisenberg Hamiltonians, npj Quantum Mater. 4, 63 (2019).
  * (27) K. Aoyama and H. Kawamura, Spin ordering induced by lattice distortions in classical Heisenberg antiferromagnets on the breathing pyrochlore lattice, Phys. Rev. B 99, 144406 (2019).
  * (28) K. Aoyama, M. Gen, and H. Kawamura, Effects of spin-lattice coupling and a magnetic field in classical Heisenberg antiferromagnets on the breathing pyrochlore lattice, Phys. Rev. B 104, 184411 (2021).
  * (29) S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Skyrmion Lattice in a Chiral Magnet, Science 323, 915 (2009).
  * (30) S. Seki, X. Z. Yu, S. Ishiwata, and Y. Tokura, Observation of Skyrmions in a Multiferroic Material, Science 336, 198 (2012).
  * (31) Y. Tokura and N. Kanazawa, Magnetic Skyrmion Materials, Chem. Rev. 121, 2857 (2021).
  * (32) K. Momma and F. Izumi, VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Crystallogr. 44, 1272 (2011).
  * (33) H. L. Pinch, M. J. Woods, and E. Lopatin, Some new mixed A-site chromium chalcogenide spinels, Mat. Res. Bull. 5, 425 (1970).
  * (34) R. Plumier, M. Sougi, M. Lecomte, and A. Miedan-Gros, High Magnetic Field Study of Sulphur Spinel Cu1/2In1/2Cr2S4. Z. Phys. B 40, 227 (1980).
  * (35) T. G. Aminov, G. G. Shabunina, E. V. Busheva, and V. M. Novotortsev, Magnetic Properties of Cox(Cu0.5In0.5)1-xCr2S4 Solid Solutions, Russ. J. Inorg. Chem. 61, 461 (2016).
  * (36) See Supplemental Material at http://xxxx for additional information of experiments and analyses.
  * (37) A. Ikeda, Y. H. Matsuda, and H. Tsuda, Note: Optical filter method for high-resolution magnetostriction measurement using fiber Bragg grating under millisecond-pulsed high magnetic fields at cryogenic temperatures, Rev. Sci. Instrum. 89, 096103 (2018).
  * (38) A. Miyake, H. Mitamura, S. Kawachi, K. Kimura, T. Kimura, T. Kihara, M. Tachibana, and M. Tokunaga, Capacitive detection of magnetostriction, dielectric constant, and magneto-caloric effects in pulsed magnetic fields, Rev. Sci. Instrum. 91, 105103 (2020).
  * (39) T. Kihara, Y. Kohama, Y. Hashimoto. S. Katsumoto, and M. Tokunaga, Adiabatic measurements of magneto-caloric effects in pulsed magnetic fields up to 55 T, Rev. Sci. Instrum. 84, 074901 (2013).
  * (40) M. Matsuda, H. Ueda, A. Kikkawa, Y. Tanaka, K. Katsumata, Y. Narumi, T. Inami, Y. Ueda, and S.-H. Lee, Spin-lattice instability to a fractional magnetization state in the spinel HgCr2O4. Nat. Phys. 3, 397 (2007).
  * (41) S. Ji, S.-H. Lee, C. Broholm, T. Y. Koo, W. Ratcliff, S.-W. Cheong, and P. Zschack, Spin-Lattice Order in Frustrated ZnCr2O4. Phys. Rev. Lett. 103, 037201 (2009).
  * (42) P. Balla, Y. Iqbal, and K. Penc, Degenerate manifolds, helimagnets, and multi-Q chiral phases in the classical Heisenberg antiferromagnet on the face-centered-cubic lattice, Phys. Rev. Res. 2, 043278 (2020).
  * (43) R. Plumier, F. K. Lotgering, and R. P. van Stapele, Magnetic properties of Cu1/2In1/2Cr2S4 and some related compounds, J. Phys. C 1, 324 (1971).
  * (44) C. L. Henley, Ordering by disorder: Ground-state selection in fcc vector antiferromagnets, J. Appl. Phys. 61, 3962 (1987).
  * (45) O. Benton and N. Shannon, Ground State Selection and Spin-Liquid Behavior in the Classical Heisenberg Model on the Breathing Pyrochlore Lattice, J. Phys. Soc. Jpn. 84, 104710 (2015).
  * (46) I. Dzyaloshinskii, A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics. J. Phys. Chem. Solids 4, 241 (1958).
  * (47) T. Moriya, New Mechanism of Anisotropic Superexchange Interaction. Phys. Rev. Lett. 4, 228 (1960).
  * (48) M. Elhajal, B. Canals, R. Sunyer, and C. Lacroix, Ordering in the pyrochlore antiferromagnet due to Dzyaloshinsky-Moriya interactions. Phys. Rev. B 71, 094420 (2005).
  * (49) M. C. Kemei, P. T. Barton, S. L. Moffitt, M. W. Gaultois, J. A. Kurzman, R. Seshadri, M. R. Suchomel, and Y. Kim, Crystal structures of spin Jahn-Teller-ordered MgCr2O4 and ZnCr2O4, J. Phys.: Condens. Matter 25, 326001 (2013).
  * (50) H. Katsura, N. Nagaosa, and A. V. Balatsky, Spin Current and Magnetoelectric Effect in Noncollinear Magnets, Phys. Rev. Lett. 95, 057205 (2005).
  * (51) T. Arima, Ferroelectricity Induced by Proper-Screw Type Magnetic Order, J. Phys. Soc. Jpn. 76, 073702 (2007).
  * (52) A. Miyake, Y. Shimizu, Y. J. Sato, D. Li, A. Nakamura, Y. Homma, F. Honda, J. Flouquet, M. Tokunaga, and D. Aoki, Enhancement and Discontinuity of Effective Mass through the First-Order Metamagnetic Transition in UTe2, J. Phys. Soc. Jpn. 90, 103702 (2021).
  * (53) S. Kimura, M. Hagiwara, T. Takeuchi, H. Yamaguchi, H. Ueda, Y. Ueda, and K. Kindo, Large change in the exchange interactions of HgCr2O4 under very high magnetic fields, Phys. Rev. B 83, 214401 (2011).
  * (54) Y. Yamamoto and T. Nagamiya, Spin Arrangements in Magnetic Compounds of the Rocksalt Crystal Structure, J. Phys. Soc. Jpn. 32, 1248 (1972).
  * (55) J. D. Bocarsly, R. F. Need, R. Seshadri, and S. D. Wilson, Magnetoentropic signatures of skyrmionic phase behavior in FeGe, Phys. Rev. B 97, 100404(R) (2018).
  * (56) T. Okubo, S. Chung, and H. Kawamura, Multiple-$q$ States and the Skyrmion Lattice of the Triangular- Lattice Heisenberg Antiferromagnet under Magnetic Fields, Phys. Rev. Lett. 108, 017206 (2012).
  * (57) T. Shimokawa and H. Kawamura, Ripple State in the Frustrated Honeycomb-Lattice Antiferromagnet, Phys. Rev. Lett. 123, 057202 (2019).
  * (58) T. Shimokawa, T. Okubo, and H. Kawamura, Multiple-$q$ states of the $J_{1}$-$J_{2}$ classical honeycomb-lattice Heisenberg antiferromagnet under a magnetic field, Phys. Rev. B 100, 224404 (2019).
  * (59) S. Nakatsuji, Y. Nambu, H. Tonomura, O. Sakai, S. Jonas, C. Broholm, H. Tsunetsugu, Y. Qiu, and Y. Maeno, Spin Disorder on a Triangular Lattice, Science 309, 1697 (2005).
  * (60) O. Smirnova, M. Azuma, N. Kumada, Y. Kusano, M. Matsuda, Y. Shimakawa, T. Takei, Y. Yonesaki, and N. Kinomura, Synthesis, Crystal Structure, and Magnetic Properties of Bi3Mn4O12(NO3) Oxynitrate Comprising $S=3/2$ Honeycomb Lattice, J. Am. Chem. Soc. 131, 8313 (2009).
  * (61) S. Gao, H. D. Rosales, F. A. G. Albarracín, V. Tsurkan, G. Kaur, T. Fennell, P. Steffens, M. Boehm, P. Čermák, A. Schneidewind, E. Ressouche, D. C. Cabra, C. Rüegg, and O. Zaharko, Fractional antiferromagnetic skyrmion lattice induced by anisotropic couplings, Nature 586, 37 (2020).
  * (62) K. Aoyama and H. Kawamura, Hedgehog-lattice spin texture in classical Heisenberg antiferromagnets on the breathing pyrochlore lattice, Phys. Rev. B 103, 014406 (2021).