# Temperature induced Neutral to Ionic phase Transition of the charge transfer
crystal Tetrathiafulvalene-Fluoranil

Elena Ferrari  Francesco Mezzadri  Matteo Masino<EMAIL_ADDRESS>Dipartimento di Scienze Chimiche, della Vita e della Sostenibilità Ambientale
(SCVSA) and INSTM-UdR Parma, Università di Parma, Parco Area delle Scienze,
IT-43124 Parma, Italy

###### Abstract

The Temperature Induced Neutral to Ionic phase Transition (TI-NIT) is a rare
phenomenon occurring in mixed stack charge transfer (CT) crystals made up of
alternating $\pi$-electron Donor (D) and Acceptor (A) molecules. We were able
to grow crystals of Tetrathiafulvalene-Fluoranil (TTF-FA), and to show that it
undergoes TI-NIT like the prototype CT crystal TTF-Chloranil. We characterized
both room and low T phases through IR and Raman spectroscopy and XRD,
demonstrating that while TTF-FA is quasi-neutral at room T, its ionicity jumps
from 0.15 to 0.7 at low T, therefore crossing the Neutral-Ionic borderline.
The transition, occuring around 150K, is first order, with large thermal
hysteresis and accompanied by crystal cracking. In the high T phase D and A
molecules lie on inversion center, i.e. the stacks are regular, whereas the
low T phase is characterized by the loss of the inversion symmetry along the
stack as the stacks are strongly dimerized and by the doubling of the unit
cell.

###### pacs:

##  I Introduction

Mixed stack (ms) charge transfer (CT) organic cocrystals may undergo a rare
phase transition named Neutral to Ionic Transition (NIT)[1]. These crystals
contain planar $\pi-$electron donor (D) and acceptor (A) molecules that
overlap their frontier orbitals forming DADADA stack chains, due to the CT
interaction. Depending on the ionicity ($\rho$), that is the average charge
located on D and A, ms CT crystals can be Neutral (N) if $\rho<0.5$ or Ionic
(I) otherwise. In few cases, low temperatures or high pressures drive the
system from a N ground state to a I one, originating a collective CT called
NIT. Due to Peierls instability, the NIT is always associated with stack
dimerization, that breaks the inversion symmetry of the chain.

Since dimerized stacks are polar, the I phase is potentially ferroelectric, as
in the case of the prototypical Tetrathiafulvalene-Chloranil (TTF-CA)[1, 2]
that undergoes NIT at 81 K. Interestingly, TTF-CA ferroelectricity has an
electronic origin, as the polarization is dominated by intermolecular charge
transfer[3]. Further studies on TTF-CA revealed that the same transition can
be also photoinduced [4, 5, 6]. Since the first discovery of the NIT, other
exotic phenomena such as dielectric anomaly [7, 8, 9], soft modes[10, 11],
quantum phase transitions[12, 13] have been observed in correspondence with
it.

For the above reasons, the mechanisms of the NIT are still being investigated,
both experimentally and theoretically. It is known that the NIT is
characterized by the competition between two instabilities: a valence one and
a structural one. While the former is favored by 3D Coulomb interactions, that
induce discontinuous NIT, the latter is due to electron-lattice phonon
coupling, that drives the stack to a continuous dimerization[14, 15]. The role
of these microscopic parameters has been demonstrated by comparing TTF-CA with
other ms CT crystals undergoing NIT[16].

In this context, the TTF-haloquinone series has been extensively studied by
chemical modification of the component molecules[17, 18]. However, the elusive
TTF-FA, first obtained in 1979 together with TTF-CA[19], has never been
experimentally reproduced and characterized. We were able to grow it as single
crystals and to observe that it also undergoes NIT at low temperature. Both
phases were characterized by vibrational spectroscopy and XRD.

##  II Experimental

Several attempts to grow TTF-FA crystals with slow methods always yielded a
red-brown amorphous product. Good quality single crystals were finally
obtained by mixing saturated toluene solutions of the two components in 1:1
stochiometric ratio and rapidly evaporating the solvent. The crystals were
unstable at room temperature due to FA sublimation and were stored at low
temperature.

Single crystal X-rays diffraction was collected at RT by using graphite
monochromatized Mo K$\alpha$ wavelength on a Bruker Smart diffractometer
equipped with an APEX II CCD detector. Low temperature data were collected
with a Bruker D8 Venture instrument, equipped with a Photon II CCD area
detector and micro-focused Mo K$\alpha$ radiation source. Temperature control
was achieved with an Oxford cryostream system working in nitrogen flux. Due to
crystal breaking at the NIT transition, the low temperature data collection
required the sample to be placed into a glass capillary filled with Fomblin®,
a perfluorinated polymer that increases its viscosity as temperature is
lowered. In such a way crystal cracking was avoided without inhibiting the
transition. The single crystals data reduction was carried out by using the
SADABS program [20] The software Shelxt was used for structures solution while
refinement was carried out full-matrix by using the Shelxl program [21].

The IR spectra were recorded with a Bruker IFS-66 FT-IR spectrometer coupled
to a Hyperion 1000 IR microscope, on the $ac$ and on the $ab$ face with the
electric field vector polarized parallel or perpendicular to the stack
direction $c$. The crystal faces were previously indexed by XRD.

The Raman spectra were recorded with a Renishaw 1000 spectrometer equipped
with the appropriate edge filter and coupled to a Leica M microscope using the
568.2 nm line of a Lexel Kr laser. Low frequency spectra have been measured
using the 633 nm HeNe laser and a Horiba LabRAM HR Evolution Raman
spectrometer equipped with a ULF Bragg filter. The laser power was always set
below 0.1 mW to avoid sample heating. Raman spectra were measured with the
exciting and scattered light polarized parallel or perpendicular to the stack
on both ac and ab crystal faces.

The crystals were cooled down to 80 K using a Linkam HFS 91 stage and were
fixed to it to prevent sample disintegration during the phase transition. In
the case of Raman measurements, the crystals were pasted to the stage with a
conductive silver paste, while they were covered with a thin and flat KBr
slide for transmission IR measurements.

Standard DFT computational methods (B3LYP, 6-31G(d)) were exploited for the
calculation of equilibrium geometry and vibrational frequencies of the
constituent molecules, in both the neutral and ionized state, using Gaussian
16 B.10[22]. The frequencies higher than 1000 cm${}^{-1}~{}$were scaled by the
factor 0.9613, as suggested in [23]. First principles calculations following
the approach described in [16] have been performed on the crystal structures
to estimate the tendency towards valence instability.

##  III Results and Discussion

###  III.1 Room temperature crystal structure

Single crystal XRD analysis confirmed the triclinic ($P-1$) structure found by
Torrance et al.[19], characterized by a mixed regular stack, with TTF and FA
molecules alternating along $a$ and laying on the inversion centers. The green
and needle-shaped TTF-FA crystals (Fig S1) are elongated in the stack
direction and usually show two different faces, parallel to the $ac$ and $ab$
planes (Fig 1). The stack direction can be easily recognized by optical
spectroscopy due to the presence of the CT excitation, occurring at 6000
cm${}^{-1}~{}$(Fig S4). The needle-like crystal shape reflects the strong and
directional CT interaction, that is the main intermolecular force and drives a
faster crystal growth along the stack direction. However, sizeable interstack
interactions are present also along the $c$-direction through C-H…F and S…O
contacts between TTF and FA units (Fig S2 and S3).

Figure 1: Crystal structure of TTF-FA at room temperature, viewed along the
[010] (left) and [001] (right) directions. The spectra have been recorded on
these two faces. The C=O $b_{1u}\nu_{10}$ mode can be excited in the ac-face
only.

TTF-FA crystals undergo a temperature induced first order NIT, associated with
sample cracking and a sudden colour change from green to brown. During the
first cooling the transition temperature is not well defined and varies in the
range 120-160 K depending on sample history and cooling rate. This is probably
due to sample defects and internal strains. Indeed, the spectra recorded on
different points of the same crystal sometimes reveal the coexistence of the
two phases[24]. The transition is reversible with a large thermal hysteresis,
with the ionic to neutral (I-N) transformation taking place at 200 K on
heating and the N-I one at 170 K on cooling again, in a reproducible way.
While the optical axes and the polarization direction of the CT excitation are
unchanged by the transition, the strongly discontinuous evolution of the IR
and Raman spectra (Fig S6) indicates a large ionicity jump, stack dimerization
and doubling of the unit cell, as described in the following sections.

###  III.2 IR spectra: Ionicity

The NIT is characterized by two order parameters: the ionicity and the stack
dimerization. Both can be conveniently studied by vibrational spectroscopy:
the ionicity affects the vibrational frequencies of D and A while the
dimerization reduces the lattice symmetry.

The degree of charge transfer can be obtained from the frequency shifts of
some selected “charge sensitive modes” [25] of the constituent molecules.
Since the frequencies of the totally symmetric vibrations can be perturbed by
e-mv interaction [26, 27], the most reliable modes are in-plane, antisymmetric
with respect to the inversion center and can be observed in the IR spectra
polarized perpendicular to the stack. Since in a first approximation the
$D_{2h}$ point group of the TTF and FA molecules is retained in the crystal,
these vibrations, polarized along the two in-plane molecular axes, have
$b_{1u}$ and $b_{2u}$ symmetry. While the spectra recorded on the $ab$ face
are dominated by the FA $b_{2u}$ bands (see Fig. 1 and Fig S5, upper panel),
the spectra recorded on the $ac$ face, reported in Fig 2, show mainly the
$b_{1u}$ vibrations and a small projection of $b_{2u}$ ones. The complete
assignment of the IR spectra is reported in Table S1.

Figure 2: IR spectra of the N- and I- phases of TTF-FA, recorded on $ac$ face with the electric field polarized perpendicular to the stack. The spectra are upshifted for clarity. Inset: Ionicity changes of TTF-FA estimated from the FA $b_{1u}\nu_{10}$ frequency: first cooling (blue) and heating (red).  Table 1:  Assignment of some selected modes of FA in TTF-FA N and I phases. The frequencies are in cm${}^{-1}~{}$and the calculated values have been scaled.  Mode | FA | TTF-FA 293K | TTF-FA 80K | FA-
---|---|---|---|---
| expt (calc) |  |  | (calc)
$a_{g}\nu_{1}$ | 1704 (1709) | 1642 | 1547 | (1545)
$a_{g}\nu_{3}$ | 1251 (1290) | 1233 | 1170 | (1269)
$b_{1u}\nu_{10}$ | 1700 (1704) | 1670 | 1572 | (1520)
$b_{2u}\nu_{18}$ | 1669 (1649) | 1651 | 1616 | (1532)
Table 2:  Assignment of some selected modes of TTF in TTF-FA N and I phases. The frequencies are in cm${}^{-1}~{}$and the calculated values have been scaled.  Mode | TTF | TTF-FA 293K | TTF-FA 80K | TTF+
---|---|---|---|---
| expt(calc) |  |  | expt(calc)
$a_{g}\nu_{2}$ | 1555(1582) | 1540 | 1522 | 1505(1511)
$a_{g}\nu_{3}$ | 1518(1538) | 1475 | 1384 | 1420(1393)
$b_{1u}\nu_{14}$ | 1530(1560) | 1527 | 1478 | 1478(1492)
$b_{1u}\nu_{16}$ | 781 (782) | 785 | 803 | 836 (831)
$b_{3u}\nu_{34}$ | 639 (654) | 657 | 680 | 705 (709)

In the spectrum of the high temperature phase the C=O and C=C stretching bands
of FA, named $b_{1u}\nu_{10}$ and $b_{2u}\nu_{18}$ respectively, are
redshifted compared to the same bands of the pure compound, as expected [28].
In fact, in passing from haloquinone neutral molecules to their radical anions
these double bonds acquire partially single bond character and for this reason
their stretching frequencies are strongly lowered [29]. Furthermore, their
frequency shifts are known to show a linear behavior, making both modes good
ionicity probes: $b_{1u}\nu_{10}$ in the whole ionicity range and
$b_{2u}\nu_{18}$ for neutral ($\rho<0.5$) compounds only [25]. In the case of
FA the ionization frequency shifts were estimated through DFT calculations due
to the lack of experimental data on FA anion. According to the calculations,
$b_{1u}\nu_{10}$ and $b_{2u}\nu_{18}$ are redshifted 184 and 117
cm${}^{-1}~{}$, respectively (Table 1). The $b_{1u}\nu_{10}$ frequency full
ionicity shift is consistent with that of other haloquinones, which is 160
cm${}^{-1}~{}$for both chloranil and bromanil [29]. At room temperature
$b_{1u}\nu_{10}$ and $b_{2u}\nu_{18}$ occur at 1670 and 1651 cm${}^{-1}~{}$,
indicating that TTF-FA is on the neutral side with a degree of charge transfer
of 0.15. Also the frequency shifts of other selected TTF charge sensitive
modes [30], $b_{1u}\nu_{14}$, $b_{1u}\nu_{16}$ and $b_{3u}\nu_{34}$ (Table 2),
yielding values in the range 0.1-0.3, confirm the above estimate.

After the phase transition all the charge sensitive bands of TTF and FA shift,
showing a large ionicity increase (Fig 2). In particular, the two strong bands
observed at 1572 and 1616 cm${}^{-1}~{}$can be assigned to FA $b_{1u}\nu_{10}$
and $b_{2u}\nu_{18}$, respectively. In fact, the former band is present only
in the spectrum recorded on the ac face, while the latter is found in the
spectra of both $ac$ and $ab$ faces, as for the high temperature phase (Fig
S5, lower panel). The $b_{1u}\nu_{10}$ frequency corresponds to a ionicity
value of 0.7, showing that TTF-FA crosses the neutral-ionic borderline at the
transition. The temperature evolution of $\rho$ estimated from the
$b_{1u}\nu_{10}$ mode, shown in the inset of Fig 2, is strongly discontinuous,
with very little changes on cooling above and below the transition. In the
case of TTF, the antisymmetric C=C stretching band is shifted to 1478
cm${}^{-1}~{}$, as in the completely ionized molecules ($\rho=1$). However,
the other charge sensitive modes indicate a lower degree of charge transfer
(0.6 for $b_{3u}\nu_{34}$ and 0.4 for $b_{1u}\nu_{16}$). Also the frequency
shifts of two Raman active TTF totally symmetric modes, $a_{g}\nu_{2}$ and
$\nu_{3}$, yield a ionicity value of 0.6-0.7, if the effects of e-mv coupling
are considered[27, 31], assuming the CT transition frequency to be 6000
cm${}^{-1}~{}$as in the neutral phase. The small discrepancies between the TTF
bands are probably due to the deformation of the molecules in the ionic phase,
as described in the following.

###  III.3 Vibronic effects: stack dimerization

While the IR spectrum polarized perpendicular to the stack gives information
on the ionicity, the parallel polarized one, being affected by vibronic
effects, is very sensitive to stack symmetry [26, 27]. In the IR spectrum of
the neutral phase polarized parallel to the stack (Fig 3, upper panel) the
absence of vibronic activation of the totally symmetric modes demonstrates
that the stack is regular, with both molecules laying on inversion centers.
However, some strong and broad IR absorptions are found at the same distance,
about 65 cm${}^{-1}~{}$, from the most strongly coupled totally symmetric
modes, i.e. TTF $a_{g}\nu_{3}$ and FA $a_{g}\nu_{1}$ and $\nu_{2}$. Similar
features were observed in the IR spectra of the centrosymmetric neutral phase
of TTF-CA [10]. They are the so-called sidebands, two phonon excitations due
to sum and difference combinations between the totally symmetric
intramolecular vibrations and an effective Peierls phonon, the lattice mode
that drives the stack dimerization. As in the case of TTF-CA, all three pairs
of sidebands shift towards the corresponding totally symmetric bands on
lowering temperature (Fig S7). This evolution reflects the softening of the
Peierls mode that remains incomplete, interrupted by the discontinuous first
order transition.

Figure 3: Comparison of TTF-FA parallel polarized IR spectrum and Raman
spectrum recorded on $ac$ face with both exciting and scattered light
polarized perpendicular to the stack, at 293 K (upper panel) and 80 K (lower
panel). $\lambda_{exc}$= 568.2 nm

Differently, the parallel polarized IR spectrum of the low temperature phase
shows saturated absorptions in correspondence with the three most strongly
coupled totally symmetric vibrations of TTF and FA, observed in the Raman
spectrum at nearly the same frequencies (Fig 3, lower panel). The IR
activation of these vibrations is a proof of the stack dimerization [27].
Unluckily, it was not possible to follow the temperature evolution of these
bands as they saturate the absorption spectrum and the reflectivity
measurements were not possible due to the damage of the crystal surface during
the transition.

###  III.4 Lattice phonons: doubling of the unit cell

Some information on the stack symmetry and on specific lattice vibrations,
such as the effective Peierls mode, have been obtained from the parallel
polarized IR spectra in the range of the intramolecular vibrations thanks to
vibronic effects. A deeper insight on lattice symmetry can be achieved by
direct Raman observation of the lattice vibrations in the range 5-200
cm${}^{-1}~{}$. These spectra are very sensitive to crystal packing and the
number of Raman active intermolecular vibrations can provide information on
the symmetry and multiplicity of the unit cell.

In the case of the neutral phase, the $P-1$ unit cell contains one D-A pair
and both molecules occupy sites with $C_{i}$ symmetry. Consequently, 9
optically active lattice phonons are expected. Since the mutual exclusion rule
holds, the six rotations should be Raman active while the three translations
IR active. In fact, the low frequency Raman spectra of the neutral phase show
six bands (Fig 4, upper panel). Both the polarizations reported here (spectra
recorded on $ac$ face) share the same bands with different relative
intensities. The other polarizations do not add relevant information.

Figure 4: Low frequency Raman spectra of TTF-FA N- (upper panel) and I- phases
(lower panel), recorded on $ac$ face. The polarization of exciting and
scattered light is the same: perpendicular (red continuous line) or parallel
(blue dashed line) to the stack. The relative magnitudes of all the spectra
reflect the intensity of the light scattered by the sample (the parallel
polarized $(a,a)$ spectrum of the I-phase has been scaled by a factor of 3 for
clarity). $\lambda_{exc}$= 633 nm

The phase transition is revealed by the appearance of many new bands and by
the intensity enhancement of the $(a;a)$ polarized spectrum, in particular in
the range 80-170 cm${}^{-1}~{}$(Fig 4, lower panel), due to the activation of
translational modes coupled to the CT excitations. This spectrum shows 12
bands, most of which are in common with the $(c^{*},c^{*})$ polarization.
Overall, 15 different bands can be distinguished in the spectral range 0-180
cm${}^{-1}~{}$. Differently, in the $(b^{*},b^{*})$ polarization the majority
of these bands are very weak (Fig S8), suggesting the occurrence of cell
doubling along the $c$-axis. If a triclinic unit cell is doubled, the
arrangement of the two adjacent dimerized stacks can be either ferrielectric,
if the stacks are inequivalent, or antiferroelectric if the stacks are
equivalent and related by inversion symmetry. In both cases 21 optically
active phonons are expected, but the number of Raman active ones is different
in each situation. In the former case all the phonons would be both IR and
Raman active due to the lack of inversion symmetry. Differently, in the latter
case, due to the presence of inversion center situated between the chains, the
mutual exclusion rule would be still valid and only 12 Raman active phonons (6
rotations and 6 translations) are predicted.

Considering the fact that at 80 K the increased resolution due to band
narrowing makes the weakest bands detectable, the ferrielectric arrangement is
unlikely. In fact, both TTF and FA have low frequency intramolecular
vibrations [28, 30]. All of them might be Raman active due to the symmetry
lowering of the molecule in the crystal. In these circumstances, in absence of
a symmetry center more than 21 Raman bands would be expected, in disagreement
with the 15 ones observed.

###  III.5 Crystal structure of the ionic phase

Below the transition temperature, single crystal diffraction clearly indicates
a change of the lattice parameters being now $a$=7.2023(2), $b$=7.8625(3),
$c$=12.6272(4) Å, $\alpha$=83.040(2), $\beta$= 73.717(1), $\gamma$=77.224(2)°.
Structure refinement in $P-1$ gave R1 agreement factor of 3.30 (Table 3), no
relevant improvement was observed in the $P1$ symmetry, despite doubling of
the parameters, as a consequence the centrosymmetric solution was chosen.
Respect to the spatial arrangement of the molecules, the $a$ and $b$ axes of
the neutral phase respectively correspond to the $b$ and $a$ axes of the ionic
phase, due to a crystallographic convention. In the following, all the axis
labels are referred to the neutral phase.

The ionic phase is characterized by the doubling of the unit cell along the
$c$-axis, with two equivalent stacks, strongly dimerized and related by
inversion symmetry in an antiferroelectric arrangement (Fig 5, right). The
intradimer and interdimer distances (plane-centroid) are 3.17 and 3.29 Å,
respectively. The intradimer relative orientation of TTF and FA is similar to
that observed in the neutral phase. On the contrary, the interdimer
interaction is reduced as the dimers are shifted respect to each other
approximately along the $c$-axis. This increases the angle between the stack
axis and the FA molecular planes from 18.3° to 35.3°. For the same reason the
$\beta$ angle changes from 93.1° to 83.0° and the a-axis increases from 7.054
to 7.863 Å(Table 3). Furthermore, the antiferroelectric arrangement of the
ionic phase might be related to the interchain interactions that are also
present in the neutral phase (Fig S2 and S3). This motif is maintained in the
ionic phase. Finally, we noticed that TTF molecules, no longer laying on the
inversion centers, are not planar, being twisted around the central C=C bond.
This is probably a consequence of the double bond weakening due to ionization.

Figure 5: Crystal structures of TTF-FA at 293 K (left) and 100 K (right), viewed along the [010] direction. The inversion centers are marked in green. Intra- and inter-dimer D-A distances (molecular plane to centroid) are reported.  Table 3:  Structural parameters of N and I phases of TTF-FA. *Respect to the spatial arrangement of the molecules, the a and b axes of the N phase respectively correspond to the b and a axes of the I phase, due to a crystallographic convention. | Neutral phase (293K) | Ionic phase (100K)∗
---|---|---
Crystal system | Triclinic | Triclinic
Space group | $P-1$ | $P-1$
$a$ (Å) | 7.054(5) | 7.2023(2)
$b$ (Å) | 7.304(6) | 7.8625(3)
$c$ (Å) | 7.359(5) | 12.6272(4)
$\alpha$(deg) | 106.150(12) | 83.040(2)
$\beta$(deg) | 93.101(14) | 73.7170(10)
$\gamma$(deg) | 102.419(14) | 77.224(2)
Z | 1 | 2
V (Å3) | 353.0(5) | 668.04
R | 4.03 | 3.30

##  IV Conclusions

In this study the temperature induced Neutral to Ionic phase transition of
TTF-FA has been investigated by vibrational spectroscopy and single crystal
XRD. The transition is first order and associated with a large ionicity jump
($\Delta\rho\simeq 0.5$), stack dimerization and doubling of the unit cell,
taking place suddenly and simultaneously, together with an explosive sample
cracking. The antiferroelectric arrangement of the chains in the ionic phase
is the most likely from both spectroscopic and structural data. Also the
absence of splitting of all the charge sensitive bands suggests the
equivalence of the two stacks.

In each phase the temperature evolution of the IR and Raman spectra is
strongly discontinuous and does not anticipate the corresponding spectra of
the other phase approaching the transition. The same behavior is common with
the first order NIT observed in TTF-CA, but with a stronger discontinuity
between the spectra of the two phases and a much larger thermal hysteresis.
However, the space groups of the two CT systems are different: the TTF-CA
N-phase is monoclinic $P_{21}/n$, with already two stacks in the unit cell
that dimerize at the transition in the polar $P_{n}$ space group[32]. On the
contrary, TTF-FA N-phase is triclinic with one D-A pair in the unit cell. At
the transition the stacks dimerize in antiphase resulting in an
antiferroelectric order, associated with the doubling of the unit cell. From
this perspective TTF-FA broken symmetry is analogous to that of dimethyl-TTF-
CA [33]. However, in the latter system the phase transition is close to second
order and probably dominated by 1-d stack dimerization over that of charge
instability. In fact, the softening of the effective Peierls mode is almost
complete and the ionicity of the low temperature phase remains below 0.5 [33].

Calculations along the lines of Ref. [16] based on the modified Hubbard model,
performed on the crystal structures, predict instability in a large ionicity
range for TTF-FA (Fig S9 and Table S2). This agrees with the strong ionicity
jump observed. Moreover, the computed model parameters allow the comparison
between TTF-FA and other ms CT systems. In this way TTF-FA is close to TTF-CA:
despite the differences in crystal packing and D-A overlap, the two systems
share similar microscopic parameters leading to a discontinuous behavior.

Since many simultaneous processes usually play a role in the NIT, the
mechanism is often difficult to identify in detail. However, a quite plausible
scenario can be put forward on the basis of the general understanding gained
through the analysis of the extended Peierls-Hubbard model [14, 34]. NIT is in
fact characterized by the interplay between the Peierls electron-phonon
coupling, yielding the dimerization of the stack, and the Madelung energy
which causes the $\rho$ increase. In this case it is clear that the Madelung
energy increase following thermal lattice contraction prevails, inducing the
observed discontinuous ionicity variation. Furthermore, the low steric
hindrance of FA, the smallest haloquinone, reduces the D-A intermolecular
distances increasing 3-d interactions and the Madelung energy. This is the
most likely explanation for the relatively high transition temperature of TTF-
FA, compared to other TTF-haloquinone complexes [17, 18]. Also, the sizeable
interstack interactions between TTF and FA certainly play an important role in
the doubling of the unit cell, leading to an antiferroelectric ordering.
Finally,the presence of metastable phases and the large thermal hysteresis
suggest a high activation energy, probably due to the drastic structural
rearrangement during the transition.

Interestingly, systems showing discontinuous NIT with a high activation
barrier might be subject to photoinduced NIT as well. This means that a pulsed
excitation can potentially generate a long-lived metastable phase, as in the
case of TTF-CA [5]. In addition, it is well known that the few mixed stack CT
showing temperature induced NIT undergo the transition at a lower pressure
($\simeq 1$GPa), compared to other systems undergoing pressure induced NIT
only [18]. Thus, TTF-FA is expected to show NIT upon the application of a
modest pressure, as suggested in Ref. [1]. Further studies are needed to
confirm this hypothesis and to compare the results with those of the widely
studied TTF-CA [35], thus achieving a deeper understanding of the transition.

##  V Acknowledgements

This project benefited from the equipment and support of the COMP-HUB
Initiative, funded by the “Departments of Excellence” program of the Italian
Ministry for Education, University and Research (MIUR, 2018-2022). The authors
gratefully acknowledge Prof A. Girlando for the enlightening discussions.

## References

  * Torrance  _et al._ [1981a] J. B. Torrance, J. E. Vazquez, J. J. Mayerle, and V. Y. Lee, Phys. Rev. Lett. 46, 253 (1981a).
  * Torrance  _et al._ [1981b] J. Torrance, A. Girlando, J. Mayerle, J. Crowley, V. Lee, P. Batail, and S. LaPlaca, Phys. Rev. Lett. 47, 1747 (1981b).
  * Kobayashi  _et al._ [2012] K. Kobayashi, S. Horiuchi, R. Kumai, F. Kagawa, Y. Murakami, and Y. Tokura, Physical Review Letters 108, 1 (2012).
  * Koshihara  _et al._ [1990] S. Koshihara, Y. Tokura, T. Mitani, G. Saito, and T. Koda, Phys. Rev. B 42, 6853 (1990).
  * Koshihara  _et al._ [1999] S.-y. Koshihara, Y. Takahashi, H. Sakai, Y. Tokura, and T. Luty, The Journal of Physical Chemistry B 103, 2592 (1999), https://doi.org/10.1021/jp984172i .
  * Collet  _et al._ [2003] E. Collet, M.-H. Lemée-Cailleau, M. B.-L. Cointe, H. Cailleau, M. Wulff, T. Luty, S.-Y. Koshihara, M. Meyer, L. Toupet, P. Rabiller, and S. Techert, Science 300, 612 (2003), https://www.science.org/doi/pdf/10.1126/science.1082001 .
  * Mitani  _et al._ [1984] T. Mitani, G. Saito, Y. Tokura, and T. Koda, Phys. Rev. Lett. 53, 842 (1984).
  * Mitani  _et al._ [1987] T. Mitani, Y. Kaneko, S. Tanuma, Y. Tokura, T. Koda, and G. Saito, Phys. Rev. B 35, 427 (1987).
  * Okamoto  _et al._ [1991] H. Okamoto, T. Mitani, Y. Tokura, S. Koshihara, T. Komatsu, Y. Iwasa, T. Koda, and G. Saito, Phys. Rev. B 43, 8224 (1991).
  * Masino  _et al._ [2003] M. Masino, A. Girlando, and Z. G. Soos, Chemical Physics Letters 369, 428 (2003).
  * Girlando  _et al._ [2008] A. Girlando, M. Masino, A. Painelli, N. Drichko, M. Dressel, A. Brillante, R. G. Della Valle, and E. Venuti, Phys. Rev. B 78, 045103 (2008).
  * Horiuchi  _et al._ [2003] S. Horiuchi, Y. Okimoto, R. Kumai, and Y. Tokura, Science 299, 229 (2003), https://www.science.org/doi/pdf/10.1126/science.1076129 .
  * Kagawa  _et al._ [2017] F. Kagawa, S. Horiuchi, and Y. Tokura, Crystals 7, 10.3390/cryst7040106 (2017).
  * Painelli and Girlando [1988] A. Painelli and A. Girlando, Phys. Rev. B 37, 5748 (1988).
  * Girlando and Painelli [1986] A. Girlando and A. Painelli, Phys. Rev. B 34, 2131 (1986).
  * Delchiaro  _et al._ [2017] F. Delchiaro, A. Girlando, A. Painelli, A. Bandyopadhyay, S. K. Pati, and G. D’Avino, Phys. Rev. B 95, 155125 (2017).
  * Horiuchi  _et al._ [2001] S. Horiuchi, Y. Okimoto, R. Kumai, and Y. Tokura, Journal of the American Chemical Society 123, 665 (2001), pMID: 11456578, https://doi.org/10.1021/ja0016654 .
  * Horiuchi  _et al._ [2006] S. Horiuchi, R. Kumai, Y. Okimoto, and Y. Tokura, Chemical Physics 325, 78 (2006).
  * Mayerle  _et al._ [1979] J. J. Mayerle, J. B. Torrance, and J. I. Crowley, Acta Crystallographica Section B 35, 2988 (1979).
  * Sheldrick [2015] G. M. Sheldrick, Acta Crystallographica Section C 71, 3 (2015).
  * Bruker [2008] Bruker,  _APEX2, SAINT and SADABS_ (Bruker AXS Inc., Madison, Wisconsin, USA, 2008).
  * Frisch  _et al._ [2016] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian˜16 Revision C.01 (2016), gaussian Inc. Wallingford CT.
  * Merrick  _et al._ [2007] J. P. Merrick, D. Moran, and L. Radom, The Journal of Physical Chemistry A 111, 11683 (2007), pMID: 17948971, https://doi.org/10.1021/jp073974n .
  * Cointe  _et al._ [2003] M. B.-L. Cointe, M. H. Lemée-Cailleau, H. Cailleau, B. Toudic, A. Moréac, F. Moussa, C. Ayache, and N. Karl, Phys. Rev. B 68, 064103 (2003).
  * Ranzieri  _et al._ [2007a] P. Ranzieri, M. Masino, and A. Girlando, The Journal of Physical Chemistry B 111, 12844 (2007a), pMID: 17941667, https://doi.org/10.1021/jp075510v .
  * Girlando  _et al._ [1982] A. Girlando, R. Bozio, C. Pecile, and J. B. Torrance, Phys. Rev. B 26, 2306 (1982).
  * Painelli and Girlando [1986] A. Painelli and A. Girlando, The Journal of Chemical Physics 84, 5655 (1986), https://doi.org/10.1063/1.449926 .
  * Girlando and Pecile [1975] A. Girlando and C. Pecile, J. Chem. Soc., Faraday Trans. 2 71, 689 (1975).
  * Girlando  _et al._ [1978] A. Girlando, I. Zanon, R. Bozio, and C. Pecile, The Journal of Chemical Physics 68, 22 (1978), https://aip.scitation.org/doi/pdf/10.1063/1.435487 .
  * Bozio  _et al._ [1979] R. Bozio, I. Zanon, A. Girlando, and C. Pecile, The Journal of Chemical Physics 71, 2282 (1979), https://aip.scitation.org/doi/pdf/10.1063/1.438564 .
  * Masino  _et al._ [2017] M. Masino, N. Castagnetti, and A. Girlando, Crystals 7, 10.3390/cryst7040108 (2017).
  * Le Cointe  _et al._ [1995] M. Le Cointe, M. H. Lemée-Cailleau, H. Cailleau, B. Toudic, L. Toupet, G. Heger, F. Moussa, P. Schweiss, K. H. Kraft, and N. Karl, Phys. Rev. B 51, 3374 (1995).
  * Ranzieri  _et al._ [2007b] P. Ranzieri, M. Masino, A. Girlando, and M.-H. Lemée-Cailleau, Phys. Rev. B 76, 134115 (2007b).
  * D’Avino  _et al._ [2017] G. D’Avino, A. Painelli, and Z. G. Soos, Crystals 7, 10.3390/cryst7050144 (2017).
  * Hanfland  _et al._ [1988] M. Hanfland, A. Brillante, A. Girlando, and K. Syassen, Phys. Rev. B 38, 1456 (1988).


# Scattering transforms on the sphere, application to large scale structure
modelling

L. Mousset1 E. Allys1 M. A. Price2 J. Aumont3 J.M. Delouis4 L. Montier3 J. D.
McEwen2 1Laboratoire de Physique de l’École Normale Supérieure, ENS, Univ.
PSL, CNRS, Sorbonne Univ., Univ. Paris Cité, 75005 Paris, France
2Mullard Space Science Laboratory, UCL, Holmbury St Mary, Dorking, Surrey RH5
6NT, UK
3Institut de Recherches en Astrophysique et Planétologie, Univ. de Toulouse,
CNRS, CNES, UPS
4Laboratoire d’Océanographie Physique et Spatiale, Univ. Brest, CNRS, Ifremer,
IRD, Brest, France

###### Abstract

Scattering transforms are a new type of summary statistics recently developed
for the study of highly non-Gaussian processes, which have been shown to be
very promising for astrophysical studies. In particular, they allow one to
build generative models of complex non-linear fields from a limited amount of
data. In the context of upcoming cosmological surveys, the extension of these
tools to spherical data is necessary. We develop scattering transforms on the
sphere and focus on the construction of maximum-entropy generative models of
astrophysical fields. The quality of the generative models, both statistically
and visually, is very satisfying, which therefore open up a wide range of new
applications for future cosmological studies.

##  1 Introduction

Scattering transforms (ST) are a recently developed class of summary
statistics for the study of non-Gaussian processes [1]. They are inspired from
neural networks but do not require any training step to be computed.
Introduced recently in astrophysics [2] [3], ST have since demonstrated their
ability to characterize highly non-Gaussian processes, for instance for
parameter estimation and classification tasks.

Another feature of ST is that they allow one to build very efficient
generative models of physical fields, in a maximum entropy framework [4]. This
allows one to sample new approximate realisations of a given process relying
only on its ST statistics, that can be estimated even from a small amount of
data, sometime even a single example image [3] [5] [6] [7].

While these promising ST generative models have mainly been developed for 2D
planar data, the adaptation of these tools to spherical data is necessary for
cosmological analysis, especially for the next generation of full sky surveys
such as LiteBIRD for the cosmic microwave background polarization, or Rubin-
LSST and Euclid for study of the large scale structures of the Universe. The
extension of ST to spherical data however raises some difficulties: the
definition of a directional convolution with oriented filters [8] [9], as well
as the transposition of the planar translations which appear in certain ST
representations. In this paper, we propose an adaptation of state-of-the-art
ST to spherical fields. As a first step, we restrict ourselves to homogeneous
fields with properties that do not depend on the position on the sphere. This
naturally leads us to cosmological fields. In this proceeding, we only show
the results for a weak lensing field from the Large Scale Structures (LSS) of
the Universe using the CosmoGrid data set [10]. For this field, we construct
and validate a ST generative model from one single example image. Generative
models built from other fields, such as the thermal Sunyaev-Zeldovitch effect
are presented in Mousset et al. (2024) [11].

##  2 Scattering covariance on the sphere

ST refer to a family of summary statistics which includes in particular the
Wavelet Scattering Transforms (WST) [2] or the Wavelet Phase Harmonics (WPH)
[3]. In this work, we have considered the scattering covariances (SC), or
scattering spectra [5]. We chose the SC, because they only rely on
convolutions, and not on translations as the Wavelet Phase Harmonics [3],
which are difficult to univocally define on spherical maps.

SC statistics are computed from wavelet transforms, which are obtained by
convolving an initial map with a set of wavelet filters, where each filter
extracts the local information at a particular scale (labelled by $j$) and
orientation (labelled by $\gamma$). Wavelet filters need to be localized both
in pixel and harmonic space.

To compute the wavelet transforms, various convolutions on the sphere can be
considered. In this work we follow the standard directional convolution
formalism presented in, e.g., McEwen et al. (2015) [8]. The directional
convolution $I\star\Psi^{j}$ of a field $I$ with a wavelet $\Psi^{j}$ consists
in applying a rotation by a set $\rho=(\alpha,\beta,\gamma)$ of Euler angles
of the wavelet $\Psi^{j}$ initially located at the north pole, before
computing an inner product between the wavelet and the field $I$:

$(I\star\Psi^{j})(\rho)=\int_{\Omega}I(\bm{\omega})[R_{\rho}\Psi^{j}(\bm{\omega})]^{*}\text{d}\Omega,$
(1)

where $R_{\rho}$ is the rotation by Euler angles $\rho$, and $*$ stands for
complex conjugation. From $(I\star\Psi^{j})(\rho)$, we can identify
($\beta,\alpha$) with the spherical coordinates $\bm{\omega}=(\theta,\varphi)$
and $\gamma$ to the orientation which is probed in the convolution. In this
way, we obtain oriented wavelet coefficients, that we describe with the
following shorthand notation

$(I\star\Psi^{j,\gamma})(\bm{\omega})\equiv(I\star\Psi^{j})(\alpha=\varphi,\beta=\theta,\gamma).$
(2)

SC statistics characterize the power and sparsity at each scale, as well as
interaction between different scales. They are built from successive
applications of wavelet transforms and modulus operators, followed by average
and covariance computations [5]. We consider two coefficients at a single
oriented scale $\lambda_{1}=(j_{1},\gamma_{1})$:

$S_{1}^{\lambda_{1}}=\langle|I\star\Psi^{\lambda_{1}}|\rangle\,,\quad
S_{2}^{\lambda_{1}}=\langle|I\star\Psi^{\lambda_{1}}|^{2}\rangle\,,$ (3)

and two coefficients that characterize the couplings between two and three
oriented scales:

$S_{3}^{\lambda_{1},\lambda_{2}}=\text{Cov}\left[I\star\Psi^{\lambda_{1}},|I\star\Psi^{\lambda_{2}}|\star\Psi^{\lambda_{1}}\right],\quad
S_{4}^{\lambda_{1},\lambda_{2},\lambda_{3}}=\text{Cov}\left[|I\star\Psi^{\lambda_{3}}|\star\Psi^{\lambda_{1}},|I\star\Psi^{\lambda_{2}}|\star\Psi^{\lambda_{1}}\right],$
(4)

where $\langle\cdot\rangle$ corresponds to the mean over the sphere, and where
covariances are defined as $\text{Cov}[XY]=\langle XY^{*}\rangle-\langle
X\rangle\langle Y^{*}\rangle$ for two complex fields $X$ and $Y$.

##  3 Generative model of the large scale structures

###  3.1 Maximum entropy generative model

We build generative models under scattering covariance constraints. These are
maximum entropy microcanonical models, which are approximately sampled by
gradient descent [4]. Such models are constructed from statistics $\Phi$
estimated from a target field $x_{t}$; in this paper the target field is a
single full-sky map. The associated microcanonical set $\Omega_{\varepsilon}$
of width $\varepsilon$ is

$\Omega_{\varepsilon}=\left\\{x:||\Phi(x)-\Phi(x_{t})||^{2}<\varepsilon\right\\},$
(5)

where $||.||$ is the Euclidean norm. The microcanonical maximum entropy model
is the model of maximal entropy defined over $\Omega_{\varepsilon}$, which has
an uniform distribution over this set.

In this paper, we approximate this sampling with a gradient descent approach,
which consists in transporting a higher entropy Gaussian white noise
distribution into a distribution supported in $\Omega_{\varepsilon}$. In
practice, each new sample is obtained by first drawing a white noise
realization, and then performing a gradient descent in pixel space using a
loss

$\mathcal{L}(x)=||\Phi(x)-\Phi(x_{t})||^{2}.$ (6)

The typical width $\varepsilon$ of the microcanonical ensemble is then fixed
by the number of iterations used in the gradient descent. In practice we find
that $\sim 100$ iterations is typically sufficient, in which case an image at
$\texttt{nside}=128$ may be generated in $\sim 4$ seconds.

In our case, the summary statistics $\Phi(x)$ that we consider are the mean
over pixels $\langle x\rangle$, its variance $\text{Var}(x)$ and the SC
statistics.

###  3.2 Validation of the generative models

We compare the generated fields with the target one. We first do a visual
comparison of the maps to asses the quality of the spatial texture
reproduction. We then compute various standard statistics, such as the
probability density function (PDF) and the angular power spectrum, in order to
quantitatively evaluate our generative model. As we have drawn 50 samples of
the microcanonical ensemble, we compute the mean and the standard deviation of
these statistics over those 50 realisations. We also propose a comparison with
samples from a Gaussian model built from the power spectrum of the target
field. This allows us to quantify the contribution of SC statistics compared
to purely Gaussian statistics.

In Fig. 1 we show the full-sky maps. We plot the logarithm of the maps in
order to better see the textures. The generated map appears to be visually
very similar to the target one, which clearly shows that the SC statistics
capture an important part of the non-Gaussian texture of the field. On the
contrary, the structures are not reproduced in the Gaussian realisation.

Figure 1: The original target field, a generated map and a Gaussian
realisation. We plot the logarithm of the field and color bars are identical
for all maps.

The PDF is shown in Fig. 2, both with a linear (left) and a logarithmic
(middle) y-axis scaling in order to better exhibit the tails of the
distributions on several orders of magnitude. The target field is shown in
blue, the generated ones in red and the Gaussian realisations in yellow. The
comparison with the Gaussian PDF allows us to better see the non symmetric
shape, which is characteristic of non-Gaussian features. As we can see, the
PDF is well reproduced up to five orders of magnitude.

The angular power spectrum is shown in Fig. 2 (right) and it is well
reproduced over all scales. However, small oscillations around the target can
be seen in the generated power spectra. These are residual features related to
the frequency bands of the wavelets, which illustrate the trade-off between
the quality of reproduction we want to achieve and the number of filters we
use; i.e. the computational efficiency of our generative model. In Mousset et
al. (2024) [11], another validation using the Minkowski functionals, which are
standard non-Gaussian statistics, is made, showing that they are also very
well reproduced.

Figure 2: PDF, both with a linear (left) and a logarithmic (middle) y-axis
scaling, and angular power spectrum for the LSS field. The target is shown in
blue, the generated fields in red, and the Gaussian realisations in yellow. We
plot the mean (solid line) and the standard deviation (shadow envelope) over
50 realisations.

##  4 Conclusions

The main result of this work is the extension of state-of-the-art ST for
generative modelling to spherical fields. We have worked with the last
generation of ST statistics, named scattering covariances, which were
previously introduced for one dimension and two dimensions planar fields. They
have the advantage of relying only on successive wavelet transforms and
modulus, as well as on covariances, and do not require any translations. We
have also used state-of-the-art directional convolutions on the sphere [12],
computed in spherical harmonic space.

These developments allow us to build generative models of full sky spherical
fields without the need for large training datasets. In fact, our method holds
even in the limit of a single data realisation. In this proceeding, we have
shown the result for the LSS weak lensing field for which they performed
extremely well. In fact, the performance of those generative models were
validated quantitatively on different fields and the diversity in terms of
structures between the maps shows the impressive ability of SC to
comprehensively characterize very different non-Gaussian textures [11].

This work introduces a new powerful innovative approach for spherical data,
and it opens interesting perspectives for astrophysical applications. In
particular, we plan to use it for the study and the modeling of CMB
astrophysical foregrounds. The first goal will be to have a tool to produce
multiple realisations of the different astrophysical components. Then, ST
could play a role in component separation, relying both on recently developed
ST-based statistical component separations approaches, e.g., Auclair et al.
(2024) [13], as well as investigating how classical component separation
methods could benefit from ST, using the non-Gaussianities as an additional
lever arm to disentangle different components.

## References

## References

  * [1] J. Bruna and S. Mallat, IEEE Trans. Pattern Anal. Mach. Intell. 35, 1872 (2013)
  * [2] E. Allys et al, A&A 629, A115 (2019)
  * [3] E. Allys et al, Phys. Rev. D 102, 103506 (2020)
  * [4] J. Bruna and S. Mallat, Mathematical Statistics and Learning 1, 257 (2019)
  * [5] S. Cheng et al, PNAS Nexus 3, 103 (2024)
  * [6] B. Régaldo-Saint Blancard et al, The Astrophysical Journal 943, 9 (2023)
  * [7] M. A. Price et al, The Open Journal of Astrophysics 6, 35 (2023)
  * [8] J. D. McEwen et al, IEEE Trans. Sig. Proc. , arXiv:1509.06749 (2015)
  * [9] J. D. McEwen, C. Durastanti and Y. Wiaux, Applied Comput. Harm. Anal. 44, 59 (2018)
  * [10] T. Kacprzak et al, JCAP 2023, 050 (2023)
  * [11] L. Mousset et al, submitted to A&A , (2024)
  * [12] M. A. Price and J. D. McEwen, Journal of Computational Physics 510, 113109 (2024)
  * [13] C. Auclair et al, A&A 681, A1 (2024)


#  First-principles based study of magnetic states and high-pressure enthalpy
landscape of manganese sulfide polymorphs

Artem Chmeruk<EMAIL_ADDRESS>Maribel Núñez-Valdez  Deutsches
GeoForschungsZentrum GFZ, Telegrafenberg, 14473 Potsdam, Germany  Institut für
Geowissenschaften, Goethe-Universität Frankfurt, D-60438 Frankfurt am Main,
Germany

(August 27, 2024)

###### Abstract

Using first-principles calculations in combination with special quasirandom
structure and occupation control matrix methods, we study the magnetic
ordering and the effect of pressure on manganese sulfide polymorphs. At
ambient conditions, MnS is commonly observed in paramagnetic rock-salt
structure, but as temperature decreases at constant pressure it becomes
antiferromagnetic. On the other hand, at room temperature MnS has shown to
undergo structural transformations as pressure increases. Here, we show that
our approach involving the ordering/disordering of the local magnetic moments
in addition to the explicit control of the localization of the Mn
$d$-electrons produces energy band gaps and local magnetic moments in
excellent agreement with those observed experimentally, particularly for
paramagnetic MnS. Finally, we focus on how MnS evolves under pressure and from
its enthalpy landscape we identify at about 21 GPa, the structural
transformation from rock-salt to orthorhombic MnP-type. This structural
transformation resembles closely experimental results in which a new stable
but unidentified MnS phase was previously reported.

††preprint: APS/123-QED

##  I INTRODUCTION

Manganese sulfide has had a renew draw of attention due to a number of its
possible applications. For example, MnS, similarly to Cd and Zn chalcogenides,
is known to form part of diluted magnetic semiconductors that possess
attractive magneto-optical properties Tappero  _et al._ (1997). Recently,
various nanostructured polymorphs of MnS have attracted wide interest as they
have shown a combination of electrochemical properties that make them
promising materials for applications in Li-ion and Li-S batteries as well as
in supercapacitors He  _et al._ (2014); Zhang  _et al._ (2004, 2012); Tang
_et al._ (2015).

At room temperature ($T_{R}$) and ambient pressure (0 GPa), MnS is a
paramagnetic (PM) insulator and it can be found crystallized into face-
centered cubic rock-salt (RS) or hexagonal wurtzite (WZ) structures, with the
former being the most common in nature known as the mineral alabandite. RS-
MnS, below its Néel temperature ($T_{N}=150$ K), exhibits antiferromagnetic
(AFM) ordering of the second kind (AFM-II), i.e., sheets of ferromagnetically
coupled magnetic moments stacked antiferromagnetically in the [111] direction
Corliss  _et al._ (1956); Lines and Jones (1966). This particular magnetic
configuration in RS-MnS leads to a slight distortion from the ideal cubic
structure towards a trigonal symmetry, making RS-MnS a magnetostrictive
material Morosin (1970). On the other hand, WZ-MnS has a $T_{N}$ of 80 K
Danielian and Stevens (1961) and displays AFM ordering of the third kind (AFM-
III) in which two-thirds of the nearest neighbors are antiparallel and one-
third parallel while next-nearest neighbors are arranged in the opposite way
Corliss  _et al._ (1956).

As a transition metal (TM) compound, MnS is also a source of active research
among theorists and modelers in materials science. The electronic structure of
TM systems is notoriously challenging to describe and often requires resorting
to specialized theories that go beyond the single-particle description. Thus,
MnS, as its oxide counterpart Trimarchi  _et al._ (2018), is a fruitful and
compelling TM system that can be used to test whether or not a single-particle
theory such as density functional theory (DFT) Hohenberg and Kohn (1964); Kohn
and Sham (1965) can correctly determine the system’s ground-state and other
physical properties. Several works have been carried out in that regard using
various implementations of DFT Oguchi  _et al._ (1983); Raybaud  _et al._
(1997a, b); Tappero  _et al._ (1997); Tappero and Lichanot (1998); Hobbs and
Hafner (1999); Rohrbach  _et al._ (2003). Although, some studies correctly
predicted the AFM-II configuration as the ground-state of RS-MnS, they showed
some discrepancy between their predicted local magnetic moment ($m_{loc}$)
values of 4.92 $\mu$B Tappero  _et al._ (1997) and 4.082 $\mu$B Hobbs and
Hafner (1999), and the neutron diffraction experimental measurement of
$m_{loc}=4.54$ $\mu$B at $4.2$ K Fender  _et al._ (1968). On the other hand,
in order to study the PM phase of MnS, simulations fall into two categories:
I) MnS is modeled as a _non-magnetic_ system, i.e., $m_{loc}=0$; or II) MnS is
modeled with finite but randomly oriented local magnetic moments.
Investigations employing the non-magnetic approach have the shortcoming that
within the single-particle band theory, MnS is always expected to be metallic
Raybaud  _et al._ (1997b); Hobbs and Hafner (1999), since the electron density
is uniformly distributed over the orbital states, thus, leading to partially
unfilled $d$-bands. The second approach, more compatible with the actual
behavior of the Mn local magnetic moments, makes use of the one-electron
energy spectrum determined from multiple-scattering theory, namely the
Korringa-Kohn-Rostoker coherent potential approximation (KKR-CPA) method
Oguchi  _et al._ (1983). Within this picture, non-collinear paramagnetism is
achieved when the electron is treated as moving through a medium of completely
randomly oriented magnetic moments, then the energy spectrum is obtained by
the scattering events that the electron experiences. This scheme, usually
called _disordered local moments_ (DLM), is inherently better for modeling PM
phases than the _non-magnetic_ approach described above, as it considers an
assembly of finite magnetic moments, which corresponds to the experimentally
observed PM state of MnS. However, as the CPA effective potential is
constructed in such a way as to represent an average effect of an alloy of
randomly oriented moments, the possibility for individual magnetic moments to
develop their own local magnetic environments is excluded. Hence, this factor
would constitute a critical drawback in the modeling of PM MnS. An alternative
way to modeling magnetic disorder, and the one employed in this work, is the
so-called method of _special quasirandom structures_ (SQS) Zunger  _et al._
(1990). Although collinear by construction, it offers an upper hand over the
CPA in that, a SQS is built to represent an _average quantity_ , in this case
the net magnetization. This subtle difference in the treatment of finite
magnetic moments with respect to the CPA has been shown to be a source of
gapping in TM oxides and perovskites Trimarchi  _et al._ (2018); Varignon  _et
al._ (2019), and we expect it to play an essential role in TM sulphides as
well. For a more detailed discussion on the difference between DLM and SQS
techniques, see e.g. Trimarchi  _et al._ (2018) and references therein.

Besides the significance of MnS at ambient conditions, there have been
experimental and theoretical attempts to study the effect of pressure on it.
Using X-ray diffraction and diamond-anvil-cell (DAC) experiments, a phase
transition from RS (B1) to orthorhombic GeS-type (B16) structure was reported
at $\sim$7.2 GPa Kraft and Greuling (1988). However, a subsequent independent
DAC experiment up to $\sim$21 GPa did not confirm it McCammon (1991), while
another one also did not observe the B1$\rightarrow$B16 transformation and
instead found a transition from RS to an unknown phase at about 26 GPa that
remained stable up to at least 46 GPa Sweeney and Heinz (1993). More recent
DAC experiments reported the fabrication of a MnS quenchable high-pressure
nanostructure, which was indexed as orthorhombic MnP-type (B31) Xiao  _et al._
(2015). This B1$\rightarrow$B31 transformation was observed at the transition
pressure ($P_{T}$) of approximately 22 GPa, which is in good agreement with
the previous findings in bulk MnS, thus shedding light on the nature of the
unknown structure. However, the question whether or not there is another
structural transition above 46 GPa still remained unanswered. Finally, static
DFT simulations of the B1- and B31-MnS phases under pressure suggested that
the B31 phase was preferable at all pressures, including ambient conditions
Xiao  _et al._ (2015). Though, this DFT study did not specify what magnetic
configuration was considered or under what scheme (non-magnetic or disorder
PM) the MnS phases were simulated.

Therefore, the purpose of our investigations is twofold. Firstly, we extend
the work of Trimarchi  _et al._ (2018); Varignon  _et al._ (2019) for TM
oxides and perovskites to the case of TM sulphides. Our goal here is to show
that without resorting to specialized theories of highly correlated materials
beyond on site corrections (DFT+$U$), the insulating PM phase of MnS can be
achieved by introducing magnetic disorder using SQS’s and occupation matrix
control methodologies to handle localization of Mn $d$-states. In this way, we
are able to produce comparable physical properties such as lattice parameter,
energy gap, and local magnetic moment to those observed experimentally.
Secondly, using different PM MnS polymorphs’ ground states as candidates to
high-pressure phases, we calculate their enthalpies as a function of pressure
in order to pinpoint possible structural transformations of RS-MnS. Here, we
emphasize that the PM MnS polymorphs’ structures are constructed and optimized
following the SQS formalism previously described. Thus, in Section II we
define the various methods and theoretical considerations used in our work,
then we present our results in Section III-A on the AFM phases of MnS in order
to compare to previous studies employing different approaches and to available
experimental data. Similarly, in Section III-B we report our findings for
selected PM MnS polymorphs, and then, in Section III-C we analyze the high-
pressure enthalpy landscape of these polymorphs to evaluate the feasibility of
structural transformations as a function of pressure. Lastly, we present our
summary and conclusions in Section IV.

Figure 1: AFM arrangements in the MnS polymorphs investigated in this work.
(a) B1 with AFM-II ordering and (b) B4 with AFM-III ordering were observed
experimentally Corliss  _et al._ (1956); AFM configurations for the (c) B16
and (d) B31 polymorphs were generated in this study for comparison. S atoms
are not shown for simplification. Structures are visualized using VESTA Momma
and Izumi (2008)

.

Figure 2: DFT$+U$ partial density of states (DOS) for a Mn ion projected on
the majority spin (positive) and minority spin (negative) channels with $U=3$
eV within AFM configurations for (a) B1, (b) B4, (c) B16 and (d) B31 MnS
polymorphs. The Fermi energy ($E_{F}$) level is represented by the dashed line
at zero. The characters of the $d$-states are interpreted through group theory
analysis of the corresponding crystal-field surroundings before the onset of
the AFM ordering. This is done for the purpose of comparison with the PM case.
AFM ordering lowers the symmetry as follows: in B1 $O_{h}\rightarrow D_{3d}$,
in B4 $C_{6v}\rightarrow C_{1h}$, in B16 $D_{2h}\rightarrow C_{2h}$, and in
B31 $D_{2h}\rightarrow C_{2h}$.

##  II COMPUTATIONAL METHODS

All our calculations were performed within DFT+$U$, using the projector-
augmented plane wave method (PAW) Blöchl (1994) as implemented in the VASP
(version 5.4.4) code Kresse and Furthmüller (1996a, b). The valence
configurations were 4 _s 1_3 _d 6_ for Mn, and 3 _s 2_3 _p 4_ for S,
respectively. The exchange-correlation (XC) term in the effective Kohn-Sham
potential was approximated according to the Perdew-Burke-Ernzerhof
parameterization for solids (PBEsol) of the generalized gradient approximation
(GGA) Csonka  _et al._ (2009). To treat the Coulomb repulsion of the Mn
$d$-electrons, we added the Hubbard-$U$ correction Liechtenstein  _et al._
(1995) within the rotationally invariant Dudarev prescription Dudarev  _et
al._ (1998):

$E_{U}=\frac{U}{2}\sum_{I,\sigma}\sum_{i}\lambda^{I,\sigma}_{i}(1-\lambda^{I,\sigma}_{i}).$
(1)

As it is well known, $E_{U}$ represents a penalty energy proportional to $U$
(here $U$ represents the effective difference between the on-site Coulomb and
exchange interactions) for atom $I$ and spin channel $\sigma$.
$\lambda^{I,\sigma}_{i}$ are the eigenvalues (with values between 0 and 1) of
the occupation matrix (OM) $n^{I,\sigma}_{m,m^{\prime}}$ for an orthogonal set
of localized orbitals $i$, which, in general, are the linear combinations of
the atomic $d$-orbitals $m$. For all our simulations we chose $U=3$ eV, as it
was already tested for the AFM B1-MnS phase at 0 GPa Rohrbach  _et al._
(2003). We also assessed the local density approximation (LDA) Ceperley and
Alder (1980) and the standard PBE Perdew  _et al._ (1996) XC functionals +$U$,
but we found that PBEsol+$U$ produced more accurate structural parameters when
compared to available experimental values for both AFM and PM phases of MnS.
Therefore, we only present results within this scheme in the remaining
sections. Integration in the Brillouin zone was done on a $\Gamma$-centered
grid of uniformly distributed k-points with a spacing of $2\pi\times 0.3$ Å-1.
The selected plane-wave kinetic energy cutoff was 500 eV and convergence of
our structural optimizations was assumed when the total energy changes were
less than $10^{-8}$ eV and the forces on each atom smaller than $10^{-3}$
eV/Å. To simulate magnetic ordering (AFM below $T_{N}$) and disordering (in PM
phases at $T_{R}$) large enough supercells have to be built in order to
accomodate the appropriate AFM orderings as well as to allow for multiple
relaxation patterns. These two features have been shown to lead to gapping in
TM oxides and perovskites Trimarchi  _et al._ (2018); Varignon  _et al._
(2019). Overall, we considered four MnS polymorphs, namely, RS (B1,
$Fm\overline{3}m$), wurtzite (WZ or B4, $P6_{3}mc$), GeS-type (B16, $Pnma$)
and MnP-type (B31, $Pnma$).

To model the structures below $T_{N}$, we imposed the AFM-II and AFM-III
orderings observed experimentally in the B1 and B4 polymorphs, respectively
Corliss  _et al._ (1956). For the B16- and B31-MnS structures there is no
experimental or computational data on their precise magnetic ordering, except
for one report stating that the B31 polymorph is PM above 5 K Xiao  _et al._
(2015). Consequently, we tested several possible AFM arrangements and adopted
the lowest energy configurations found. (All AFM input structures are provided
in the Supplemental Material See Supplemental Material at URL
http://will.be.inserted.by.publisher .) We optimized all MnS polymorphs using
64-atom $2\times 2\times 2$ supercells, except for B4-MnS, for which a 36-atom
$3\times 3\times 1$ supercell was sufficient to realize the experimentally
observed AFM-III ordering.

On the other hand, to model the PM MnS polymorphs, we constructed SQS
supercells using the Alloy Theoretic Automated Toolkit (ATAT) software package
Van De Walle  _et al._ (2002). Under this scheme, the PM state is created as a
disordered alloy of up $\uparrow$ and down $\downarrow$ moments located at
different sites. As the construction of a SQS is based on the computation of
the correlation function between the species that constitute the alloy
($\uparrow$ and $\downarrow$ moments in this case), the size of the supercell
(number of atoms) used is vitally important to obtain magnetic configurations
whose components are not spatially correlated among themselves. The SQS degree
of randomness is improved by the number of atomic figures (pairs, triplets,
quadruplets, etc.) included in the calculation of the correlation function and
by the interaction distance between the atoms in a given figure. For example,
one could start by only considering nearest neighbors for atomic pairs, then
gradually add more pairs (increasing the interaction distance), and/or higher
order figures into consideration. The larger the interaction distance and the
more figures are considered, the larger the supercell becomes to achieve total
randomness, and the SQS generation quickly becomes computationally demanding.
Because the creation of a SQS is purely configurational, in order to save
resources, for this study we used a previously produced SQS supercell for the
B1 structure of TM oxides Trimarchi  _et al._ (2018) with 64 atoms. For the
B16- and B31-MnS polymorphs, pairs and triplets were included to obtain
64-atom $2\times 2\times 2$ SQS’s supercells. Similarly, for the B4 structure,
we proved 2$\times$2$\times$2, 3$\times$3$\times$2 and 4$\times$4$\times$1
supercells with 32, 72 and 64 atoms, respectively, resulting in ground state
energy differences among them of less than $10^{-3}$ eV per formula unit
(f.u.), thus, we took the 32-atom supercell. (All constructed PM structures
are given in the Supplemental Material See Supplemental Material at URL
http://will.be.inserted.by.publisher .)

Finally, under pressure ($P$) we could not converge some PM MnS structures to
stable ground states. In order to deal with this obstacle, it has been noticed
that non-integer occupancies in Eq. 1 may lead to local minima and not allow a
system to achieve its true ground state. Thus, by controlling orbital filling
explicitly in $n^{I,\sigma}_{m,m^{\prime}}$, one can help a non-converging
trapped system circumvent a metastable OM ill-setup. This method was shown to
find stable states, previously unaccessible, for $d$ and $f$ oxides Allen and
Watson (2014). Therefore, in our collinear magnetism approach, to keep the
magnetic ordering as determined at 0 GPa, and at the same time isolate energy
changes due only to structural transformations for $P>0$ GPa, we specified for
a given Mn ion in all polymorphs the same diagonal unitary occupation of its
$d_{m}$ orbitals ($m=-2,-1,\dots,2\Rightarrow 5\times 5$ matrices) using the
“occupation-matrix-control-in-VASP” algorithm Allen and Watson (2014). In this
way, our formerly problematic cases were able to converge.

Figure 3: Generated SQS supercells representing the PM state of MnS polymorphs
as an alloy of randomly distributed Mn $\uparrow$ and $\downarrow$ magnetic
moments for (a) cubic B1, (b) hexagonal B4, and orthorhombic (c) B16 and (d)
B31 structures. S atoms are not shown for simplification. Structures are
visualized using VESTA Momma and Izumi (2008).

##  III RESULTS AND DISCUSSION

###  III.1 AFM ordering in MnS

Firstly, our static DFT+$U$ optimized MnS polymorphs with the AFM ordering
considered for each structure are shown in Fig. 1. Our predicted lattice
parameter $a_{0}=5.172$ Å for B1-MnS (Fig. 1(a)) is only $\sim$0.8% smaller
than its experimental counterpart of 5.212 Å Corliss  _et al._ (1956). As it
was mentioned in Section I, the AFM-II onset in the B1-MnS structure leads to
a trigonal distortion, the degree of this distortion can be estimated by the
magnitude of the cube corner angle given by $\frac{\pi}{2}+\Delta$, where
$\Delta$ measures the deviation from the ideal cubic symmetry. In our relaxed
B1-MnS polymorph we find $\Delta\approx 0.095^{\circ}$, which is in good
agreement with the observed deviation of $0.099^{\circ}\pm 0.015^{\circ}$
Morosin (1970). For the hexagonal B4-MnS polymorph we obtain the lattice
parameters $a_{0}=3.963$ Å and $c_{0}=6.437$ Å, which are in excellent
agreement with experimental values of $a_{0}=3.987$ Å and $c_{0}=6.438$ Å
Corliss  _et al._ (1956). The optimized volumes per f.u. for each MnS
polymorph are listed in Table 1 for comparison, and show that B4-MnS is the
least dense.

Figure 4:  DFT$+U$ partial density of states (DOS) for a Mn ion projected on the majority spin (positive) and minority spin (negative) channels with $U=3$ eV within PM SQS’s for (a) B1, (b) B4, (c) B16 and (d) B31 MnS polymorphs. The Fermi energy ($E_{F}$) level is represented by the dashed line at zero. Table 1: Computed volume $V$ per formula unit (f.u.), energy gap $E_{g}$, and local magnetic moment $m_{loc}$ of MnS polymorphs with AFM configurations. Magnetic | Polymorph |  $V$/f.u. | $E_{g}$ | $m_{loc}$
---|---|---|---|---
Ordering |  |  | (Å3) | (eV) | ($\mu$B)
AFM-II | RS | B1 | 34.58 | 2.0 | 4.45
AFM-III | WZ | B4 | 43.76 | 2.4 | 4.43
AFM | GeS | B16 | 34.66 | 2.1 | 4.45
AFM | MnP | B31 | 34.38 | 1.6 | 4.42

As it is well known, computed DFT band gap energies ($Eg$) are ordinarily
underestimated, but from our simulations we find all MnS polymorphs clearly
insulating, as can be seen from $E_{g}$ values in Table 1. For B1-MnS we
obtain $E_{g}=2.0$ eV, which is to some extent, in better agreement with the
experimental value of about 3.1 eV Huffman and Wild (1967), than reports of
$E_{g}=1.36$ eV as predicted from KKR-CPA Oguchi  _et al._ (1983), $\sim$1 eV
from Perdew-Wang (PW) GGA without $U$ Hobbs and Hafner (1999), and $\sim$1.5
eV from GGA-PW$+U$, with $U=3$ eV Rohrbach  _et al._ (2003). Overall, our
findings imply that the magnitude of $E_{g}$ increases with $V/$f.u., with the
hexagonal B4-MnS polymorph having the largest energy band gap, while the
orthorhombic B31-MnS structure the most narrow (Table 1).

The projected density of states of Mn $s$ and $d$ orbitals for all MnS
polymorphs are shown in Fig. 2. The band characters are assigned roughly by
looking at the $m$, $l$ orbital quantum numbers and considering the crystal-
field splitting observed for the point group symmetries of the respective
structures. For simplicity in this symmetry analysis, we do not consider spin
and magnetic anti-unitary operator. Of course, one could do a complete formal
analysis using magnetic point groups and magnetic irreducible co-
representations Bradley and Cracknell (2009) if one were after a thorough
understanding of the bands, but for the purpose of our study it is not
necessary. It should be noted, however, that the onset of the AFM ordering
lowers the site symmetry, leading to additional lifting of orbital
degeneracies. For example, the $d$-states of octahedrally coordinated Mn2+ in
B1-MnS split into two levels, a high energy doublet ($e_{g}^{2}$) and a low
energy triplet ($t_{2g}^{3}$). Here, the band gap is opened between
$e_{g}$–majority and $t_{2g}$–minority bands, Fig. 2(a). However, the onset of
AFM-II ordering lowers the symmetry from octahedral $O_{h}$ to rhombohedral
$D_{3d}$, resulting in additional splitting of the $t_{2g}$ states into a
singlet and a doublet. In B4-MnS the hexagonal crystal field splits the
$d$-shell into a singlet $A_{1}$ and two doublets $E_{1}$ and $E_{2}$.
According to our calculation, in this structure the band gap is opened between
the $A_{1}$–majority singlet and the $E_{2}$–minority doublet, Fig. 2(b),
while the AFM-III ordering lowers the original $C_{6v}$ site symmetry to
$C_{1h}$. Lastly, both B16- and B31-MnS structures belong to the same $D_{2h}$
point group. The orthorhombic crystal field of this symmetry lifts the five-
fold degeneracy of the spherical $d$-shell leading to five singlets. Fig. 2(c)
and (d) show that the band gap is opened between $B_{1g}$–majority and
$A_{g}$–minority bands in both phases.
Our calculated local magnetic moments of the four AFM-MnS polymorphs are very
similar in magnitude as one can see from Table 1. For RS-MnS we find a
$m_{loc}=4.45$ $\mu$B, which is in much better correspondence with
experimental observations of 4.54 $\mu$B Fender  _et al._ (1968) than earlier
calculations of 4.92 $\mu$B under a Hartree-Fock scheme Tappero  _et al._
(1997), 4.39 $\mu$B from KKR-CPA Oguchi  _et al._ (1983), and 4.082 $\mu$B
within GGA-PW Hobbs and Hafner (1999).

###  III.2 Magnetically disordered PM MnS

Our converged MnS polymorphs in PM states as modeled by the construction of
SQS’s are shown in Fig. 3. Here we emphasize the fact that in our
calculations, all degrees of freedom (lattice parameters and ionic positions)
were allowed to relax, unlike other studies using SQS structures for TM oxides
in which the symmetry and volume were kept fixed Trimarchi  _et al._ (2018).
At first sight, these constrictions may appear reasonable to apply, as
experimentally, there is no observable distortion from the cubic symmetry of
the PM B1-MnS polymorph. Yet, allowing full relaxations in ionic coordinates
and lattice parameters of our SQS structures and comparing to their respective
ideal symmetry could be a criterion to judge how well the created SQS
represents its polymorph. For example, for PM B1-MnS (Fig. 3(a)) we find, as
in the case when AFM ordering is achieved, a small trigonal distortion. The
computed deviation from the ideal cubic angle for our 64-atom SQS supercell is
$\Delta\approx 0.024^{\circ}$, whereas for a 216-atom SQS supercell we obtain
$\Delta\approx 0.003^{\circ}$, one order of magnitude smaller than for the
64-atom SQS supercell. Using the B1 64-atom and 216-atom SQS’s we find lattice
parameters of $a_{0}=5.181$ Å and 5.1807 Å , respectively. These values are
somewhat underestimated, but still in good agreement, with respect to
measurements of 5.225 Å McCammon (1991), 5.225 Å Sweeney and Heinz (1993), and
5.29 Å Xiao  _et al._ (2015). Similarly to the AFM-ordered MnS phases, we find
all PM MnS polymorphs to be insulators, although with smaller energy band gaps
(Table 2) than their AFM counterparts (Table 1). . Our PM $E_{g}$ values
(Table 2) suggest, however, that magnetic disorder affects the least the
energy band gap of the hexagonal B4-MnS polymorph in comparison to its AFM
analog, while B16-MnS shows the largest $E_{g}$ reduction. We particularly
highlight the fact that, even though for B1-MnS our predicted $E_{g}=1.4$ eV
is fairly underestimated as compared to the observed values of 2.7 eV Sato
_et al._ (1997) and 2.8 eV Huffman and Wild (1967), we obtain an insulating PM
B1-MnS phase, unlike other calculations in which it was found to be metallic
Oguchi  _et al._ (1983); Raybaud  _et al._ (1997a).
Finally, examining our calculated local magnetic moments in the PM phases
listed in Table 2, we observe a slight increase in magnitude when compared to
the results in the AFM structures (Table 1). For the PM B1-MnS our predicted
value is 4.50 $\mu$B.

Table 2: Computed volume $V$ per formula unit (f.u.), energy gap $E_{g}$, and local magnetic moment $m_{loc}$ of MnS polymorphs in PM state. Magnetic | Polymorph |  $V$/f.u. | $E_{g}$ | $m_{loc}$
---|---|---|---|---
Ordering |  |  | (Å3) | (eV) | ($\mu$B)
PM | RS | B1 | 34.76 | 1.4 | 4.50
| WZ | B4 | 44.06 | 2.3 | 4.44
| GeS | B16 | 36.94 | 1.2 | 4.50
| MnP | B31 | 34.94 | 0.9 | 4.46

To understand better the decrease in the magnitude of the PM energy band gaps
with respect to the AFM cases, we show in Fig. 4 the projected density of
states of a Mn ion in the PM state of the four MnS polymorphs. Our
calculations show that the introduction of magnetic disorder in the MnS
polymorphs spreads and produces new spin-majority and spin-minority states in
the $E_{g}$ range of the AFM ordered phases (Fig. 2). Consequently, magnetic
disorder exclusively (as the occupation matrix is constrained) forces a shift
of the Fermi energy level in reference to the ordered AFM cases, but in such a
way that the resulting PM MnS phases remain insulating and hence with smaller
energy band gaps. In principle, different magnetic disorder given by distinct
SQS’s would give rise to different energy band gap openings, this effect was
reported in the context of NbMnSb by disordering Mn and Ni with respect to the
sites they occupy in the ordered phase Orgassa  _et al._ (1999). Additionally,
in Fig. 4(a) and (c) we can see that for PM B1- and B31-MnS phases, the energy
band gap opens, as in their AFM analogs, between the $e_{g}$–majority and the
$t_{2g}$–minority bands, and between $B_{1g}$–majority and $A_{g}$–minority,
respectively. However, for PM B4- and B16-MnS phases, magnetic disorder
changes the character of the band gap openings in reference to their AFM
counterparts, to be between $A_{1}$–majority and $A_{1}$–minority (Fig. 4(b)),
and $A_{g}$–majority and $B_{1}$–minority (Fig. 4(d)), respectively.

###  III.3 AFM to PM evolution of the band gap

Figure 5: Representation of slight deviations from the AFM-II ordering. Blue-
up and red-down arrows are Mn atoms with $\uparrow$\- or $\downarrow$-magnetic
moment, yellow balls are S atoms. The violet lines indicate two cases of
introducing “magnetic disorder” by swapping Mn atoms of opposite magnetic
moments. In the text, a structure with D1-disorder corresponds to having only
one pair of Mn-$\uparrow$ and Mn-$\downarrow$ swapped (S1), whereas a
structure with D2-disorder has swaps S1 and S2 as indicated. Half-solid-half-
dashed loops depict $p$-$d$ hybridized orbitals. When two Mn atoms with
parallel magnetic moments are coupled through a 180o
Mn$\uparrow$-S-Mn$\uparrow$ bond, an additional coupling between $d$-derived
orbitals also occurs $\tilde{d}_{i\uparrow}^{\dagger}\tilde{d}_{j\uparrow}$.
Figure 6: Orbital-projected density of states on the majority spin (positive)
and minority spin (negative) channels with $U=3$ eV of four different
configurations of the B1-MnS with increasing level of disorder. The magnitude
of the $p$-states is magnified by a factor of three for visibility. (a) The
DOS of the AFM state is largely determined by the $p$-$d$ and $s$-$d$
hybridizations between S and Mn centered orbitals. (b) and (c) The shrinkage
of the band gap in the spin-majority channel is mainly caused by the formation
of a 180o Mn-$\uparrow$-S-Mn-$\uparrow$ chain that is absent in the AFM case.
Additionally, there is also splitting in the $t_{2g}$-derived states. This is
mostly determined by the overall number of Mn atoms with parallel magnetic
moment in a given plane. For example, the structure D2 has the same number of
Mn-$\uparrow$ atoms in $xy$ and $yz$ planes which is why we observe degeneracy
in these states, whilst this number is different in all planes in D1. (d) The
fully random SQS structure exhibits the features of D1 in the $e_{g}$
splitting whilst retaining the $xy+yz$ degeneracy shown in D2.

In order to gain more insights into the narrowing of the band gap of MnS in
going from an AFM ordered phase to its magnetically disordered PM one above
its $T_{N}$, we gradually introduce disorder in steps. We accomplish this task
by studying in more detail the projected Mn DOS for a couple of configurations
slightly departing from the original AFM ordering. First, we recall that the
most general non-interacting Hamiltonian within the DFT+$U$ scheme is given by
Hubbard (1963); Bruus and Flensberg (2004):

$\begin{split}\hat{H}=\sum_{i,j,m,m^{\prime},\sigma}\big{(}t_{i,j,\sigma}^{m,m^{\prime}}c^{\dagger}_{im,\sigma}c_{jm^{\prime},\sigma}+h.c.\big{)}\\\
+U\sum_{i,m^{\prime\prime}}n_{im^{\prime\prime}\sigma}<n_{im^{\prime\prime}-\sigma}>,\end{split}$
(2)

where $c^{\dagger}_{im,\sigma}$ and $c_{jm^{\prime},\sigma}$ create and
destroy an electron with spin $\sigma$ on site $i$, orbital $m$, and site $j$,
orbital $m^{\prime}$, respectively. $t_{i,j,\sigma}^{m,m^{\prime}}$ are the
hybridization integrals, which one can assume to follow the two-center Slater-
Koster (SK) approximation Slater and Koster (1954) as
$t_{i,j,\sigma}^{m,m^{\prime}}\approx E_{i,j,\sigma}^{m,m^{\prime}}$, i.e.,
hybridization between orbital $m$ of character $p_{i}$ and orbital
$m^{\prime}$ of character $d_{j}$. The last term in Eq. (2) represents the
mean-field decoupling of the Hubbard interaction term and only acts on the
$d$-states of the same site ($m^{\prime\prime}$ runs only through
$d$-orbitals), with $n_{im^{\prime\prime}\sigma}$ being the number operator
($c^{\dagger}_{i}c_{i}$) and $<n_{im^{\prime\prime}-\sigma}>$ the average
occupation of the respective orbital and spin (in our work, these quantities
are fixed in the DFT cycle using the OMC).

Below we present our results focused solely on the B1-MnS structure, although
for the other three polymorphs a similar analysis would also apply. Thus, we
slightly disorder the AFM-II 64-atom supercell of B1-MnS and investigate two
cases before achieving a fully disordered PM structure. Firstly, we consider a
structure in which in a 180o chain, we swap a couple of Mn atoms with
antiparallel magnetic moments (D1-disorder), and then a structure with two of
such swaps (D2-disorder), as shown in Fig. 5. For simplification, we exclude
long range interaction effects between neighboring unit-cells. To isolate the
system’s reaction to disorder, we also suppress all relaxation effects and we
only use the $\Gamma$-point in the self-consistent calculation cycle. Taking
the converged ground state for the AFM-II B1-MnS phase as the reference state,
in which the largest non-zero hybridization integrals in Eq. 2 arise between
Mn $d$-orbitals, and $s$\- and $p$-orbitals of S (Fig. 6(a)), we observe that
upon introducing D1-disorder (i.e., after the creation of a FM coupling in one
of the 180o chains, Fig. 5), there is now a possibility of additional
combination between $p$-$d$ hybridized orbitals of neighboring Mn sites with
parallel magnetic moments. This re-hybridization seems to greatly contribute
to the lift of the degeneracy of the $t_{2g}$\- and $e_{g}$-derived states and
the spreading and pushing of the latter at the top of the valence band,
effectively shrinking the energy gap as shown in Fig. 6(b). The effect of this
re-hybridization in the shift of $e_{g}$ states is more pronounced when
D2-disorder (with S1 and S2 swaps) is considered (Fig. 6(c)). However, once
that a number of swaps occurs randomly and the PM state is reached in the SQS
of B1-MnS, the energy band gap reaches its final magnitude. Another possible
mechanism of re-hybridization, also partially responsible for the lift of the
$t_{2g}$ degeneracy, can take place within planes. In the AFM-II Mn
$d$-orbitals hybridize with $p$-orbitals of S ions located at 90o. That is,
according to SK rules, in octahedral geometry $d_{xy}$ orbitals hybridize only
with $p_{x}$ and $p_{y}$ in the (010) and (100) directions, as shown in Fig.
7(a). However, when the nearest Mn neighbors have parallel magnetic moments,
there is an additional coupling between the $p$-$d$ hybridized orbitals as
depicted in Fig. 7(b). Such re-hybridization is made possible through the
admixture of $p$-states that provide the necessary spatial extension of the
otherwise highly localized $d$-states. These re-hybridized
$\tilde{d}$-orbitals centered on the Mn atoms can further mix with each other
provided there is enough spatial overlap between them. In such a case, the
$t_{2g}$ degeneracy is predominantly determined by the number of Mn atoms with
parallel magnetic moments in $xy$, $yz$ and $xz$ planes, i.e., if this number
is different for two given planes the degeneracy between the respective
orbitals is lifted, Fig. 6(b)-(c).

Figure 7: Schematic coupling between $d$-derived orbitals of neighboring Mn
sites. (a) In AFM-II B1-MnS, $d$-orbitals hybridize with $p$-orbitals of
sulfur ions Sa and Sb in the (010) and (100) directions, respectively. (b)
Nearest Mn neighbors with parallel magnetic moments allow for an additional
coupling between $p$-$d$ hybridized orbitals
($\tilde{d}_{i\uparrow}\tilde{d}_{j\uparrow}$). When the number of Mn atoms
with parallel magnetic moment in the $xy$-, $yz$\- or $xz$-plane is changed
through disorder, the respective orbital degeneracy is lifted.

###  III.4 High-pressure landscape of MnS polymorphs

As discussed in the Introduction, room temperature experimental studies of the
stability of the B1-MnS polymorph under pressure, have reported a structural
change to a new phase. However, the results were not conclusive as the
transition pressure and the new phase were not uniquely determined Kraft and
Greuling (1988); McCammon (1991); Sweeney and Heinz (1993). Therefore, in
order to examine the likelihood of a pressure induced structural phase
transformation of B1-MnS, for each MnS polymorph we computed the enthalpy
($H$) as a function of pressure ($P$), volume ($V$), and internal energy
($E$), i.e.,

$H(P)=E[V(P)]+PV(P).$ (3)

In Fig. 8 we show the relative enthalpy $\Delta H$ per f.u. of the MnS
polymorphs studied here with respect to the B1-MnS structure, i.e., $\Delta
H=H(\text{MnS-phase})-H(\text{B1})$, as a function of pressure between 0 and
60 GPa. Although high-pressure experiments at $T_{R}$ deal with the PM phases
of MnS, we present results for both AFM and PM polymorphs.

Figure 8: Static PBE+$U$ relative enthalpy $\Delta H$ per f.u. between
different (a) AFM and (b) PM states of MnS polymorphs. All results are given
in reference to the B1 structure. Figure 9: Pressure dependence of volume per
f.u. for all PM MnS polymorphs as predicted by our computations. We compare
our PM MnS results [$\star$] to experimental values: EXP91 McCammon (1991),
EXP93 Sweeney and Heinz (1993), and EXP15 Xiao  _et al._ (2015). The dashed
line indicates our predicted transition pressure ($P_{T}=21$ GPa) at which
B1$\rightarrow$B31. Solid lines are fits of the 3rd order Birch-Murnaghan’s
equation of state to our calculated volumes.

On one hand, our static $\Delta H$ values suggest that if we were at a
sufficiently low temperature to achieve AFM ordering for all MnS polymorphs,
B31-MnS would be the most stable phase at 0 GPa, and it would remain so up to
60 GPa, Fig. 8(a). On the other hand, from our static calculations of the
enthalpy in PM phases of MnS, Fig. 8(b), we conclude that the RS structure
B1-MnS is the most stable phase at ambient conditions, but as pressure
increases, it undergoes a transformation to the orthorhombic B31-MnS polymorph
at $P_{T}\approx 21$ GPa, as illustrated for clarity by the inset in Fig.
8(b). Our predicted $P_{T}$ in bulk B1-MnS is in close agreement with the
experimentally observed structural transformation of B1-MnS to an unidentified
phase with lower symmetry than hexagonal (B4-MnS) at $\sim$26 GPa Sweeney and
Heinz (1993). Our determined B1$\rightarrow$B31 transformation was also
established at about 22.3 GPa in experiments synthesizing high-pressure MnS
nanorods Xiao  _et al._ (2015), in which these B31-MnS nanorods were
quenchable to 0 GPa. Furthermore, static LDA+$U$ DFT calculations (with an
effective $U=5.13$ eV) performed alongside the experimental study found that,
although the B1- and B31-MnS polymorphs are energetically very close below 8
GPa, B31 is the most stable at all pressures between 0 and 40 GPa. However,
unlike in our static PBEsol calculations using the constructed magnetically
disorder SQS PM states of MnS, the former LDA+$U$ modeling Xiao  _et al._
(2015) does not specify how the PM was simulated. We notice as well, that our
four PM polymorphs considered here, at the level of our calculations, are
remarkably close in energy at 0 GPa, but the B1 possesses the absolute lowest
energy per f.u., and as pressure increases, B16 and B31 enthalpies increase
slightly and then decrease to start competing for stability against the B1-MnS
structure, with B31-MnS eventually becoming more stable at $P_{T}\approx 21$
GPa. We also do not observe any other further transformation up to 60 GPa. The
predicted behavior from our computations for the AFM MnS polymorphs’ enthalpy
trends are surprisingly in qualitative agreement with the LDA+$U$ modeling of
the B1- and B31-MnS phases Xiao  _et al._ (2015).

Figure 10: Lattice parameters as a function of pressure for cubic B1- and
orthorhombic B31-MnS structures. Static PBE+$U$ results [$\star$] of our
magnetically disordered PM polymorphs are compared to experimental values:
EXP91 McCammon (1991), EXP93 Sweeney and Heinz (1993), and EXP15 Xiao  _et
al._ (2015). The dashed line indicates our predicted transition pressure
($P_{T}=21$ GPa) at which B1$\rightarrow$B31.

In Fig. 9 we show our predicted change in volume per f.u. as pressure
increases for the PM MnS structures modeled in this study. We can see that
experimental volumes for cubic B1-MnS McCammon (1991); Sweeney and Heinz
(1993) are exceptionally consistent with our findings. The largest deviation
from our computed $V/$f.u. occurs for B31-MnS with respect to measurements
provided from high-pressure experiments on MnS nanorods Xiao  _et al._ (2015).
As a consequence of the B1$\rightarrow$B31 structural transformation, from our
calculated trends we obtain a decrease in volume of MnS of only approximately
2% at $P_{T}$. This volume reduction, however, is one order of magnitude
smaller in contrast to the one observed in the B31 nanorods Xiao  _et al._
(2015).

A closer look at the pressure dependance of our predicted lattice parameters
for the B1-MnS shows that our calculated values are in excellent agreement
with available experimental data McCammon (1991); Sweeney and Heinz (1993);
Xiao  _et al._ (2015), as can be seen in Fig. 10. In the case of B31-MnS, our
calculated $a,b,$ and $c$ lattice constants are systematically overestimated
by about $\sim$10%, $\sim$5% and $\sim$4%, respectively, in comparison to
experimental values at high pressures Xiao  _et al._ (2015). At 0 GPa, the
agreement between our results and the experimental values becomes much better
for $a_{0}$ and $b_{0}$, but curiously our $c_{0}$ magnitude is $\sim$8%
larger than the only measurement reported up to date Xiao  _et al._ (2015).

Lastly, we obtain the bulk modulus at 0 GPa ($K_{0}$) from Birch-Murnaghan’s
third order equation of state (EOS) fittings to our $P-V$ data between 0 and
60 GPa. Our predicted $K_{0}\approx 92$ GPa (with its pressure derivative
$K^{\prime}_{0}=3$) for PM B1-MnS is in reasonable agreement with the
experimental value reported of 88$\pm$6 GPa as fit with a variable
$K^{\prime}_{0}$ Sweeney and Heinz (1993). Our result is, however,
overestimated by almost 20%, when the experimental data is fit using a
constant $K^{\prime}_{0}=4$. A comparison between our computed bulk moduli
indicates that B1-, B16-, and B31-MnS polymorphs oppose almost
indistinguishable resistance against volume compression under hydrostatic
pressure, while B4 is the easiest to compress with a $K_{0}\approx 57$ GPa.

##  IV CONCLUSIONS

In summary, we carried out static first-principles calculations to model the
AFM and PM states of four MnS polymorphs. We demonstrated that the combination
of PBE+$U$ with the construction of SQS supercells and localization of the Mn
$d$-electrons through occupation matrix control methods allowed us to achieve
not only convergence and accuracy of structural optimizations but was vitally
crucial to obtain finite energy band gaps and local magnetic moments in the PM
phases. This result was particularly important in the case of the PM rock-salt
B1-MnS polymorph, which experimentally has been observed to be an insulator,
but was predicted to be metallic by multiple simulations. In our study we also
showed that with our approach, we were able to isolate energy changes as a
function of pressure due purely to ionic and lattice parameters relaxation
under hydrostatic compression. In this manner, we computed the enthalpies of
the PM MnS polymorphs in order to explore their high-pressure landscape to
detect structural transformations, and indeed, we determined the cubic B1-MnS
structure to be the most stable at ambient pressure and up to approximately 21
GPa, pressure at which B1 undergoes a structural transformation to the
orthorhombic B31-MnS phase. The overall trends in the electronic structure as
well as phase transitions are expected to hold for different $U$ values. The
increase (decrease) of the $U$ value will result in the increase (decrease) of
the structural parameters and, therefore, in the change of magnitude of the
transition pressure. The chosen value in our work was supported by the close
agreement to experimental structural parameters. Our predicted
B1$\rightarrow$B31 transformation, in the context of all our modeling
considerations, is rather meaningful as it closely resembles a structural
phase transition observed from X-ray diffraction and high-pressure experiments
performed on B1-MnS, in which the new but unidentified phase was reported at
about 26 GPa. Overall, our approach proved to be accurate in the modeling of
manganese sulphide polymorphs, thus we plan to extend it in the future to the
investigation of other TM compounds.

###### Acknowledgements.

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V.
(www.gauss-centre.eu) for funding this project by providing computing time
through the John von Neumann Institute for Computing (NIC) on the GCS
Supercomputer JUWELS at Jülich Supercomputing Centre (JSC) under project
abinitiomodmatsgeo. AC and MNV were supported by the Helmholtz Association
through funding of first-time professorial appointments of excellent women
scientists (W2/W3)

## author declarations

The authors have no conflicts to disclose.

## References

  * Tappero  _et al._ (1997) R. Tappero, P. D’Arco, and A. Lichanot, Chem. Phys. Lett. 273, 83 (1997).
  * He  _et al._ (2014) G. He, C. J. Hart, X. Liang, A. Garsuch, and L. F. Nazar, ACS Appl. Mater. Interfaces 6, 10917 (2014).
  * Zhang  _et al._ (2004) X. V. Zhang, S. T. Martin, C. M. Friend, M. A. Schoonen, and H. D. Holland, J. Am. Chem. Soc. 126, 11247 (2004).
  * Zhang  _et al._ (2012) L. Zhang, L. Zhou, H. B. Wu, R. Xu, and X. W. Lou, Angew. Chem. Int. Ed. 51, 7267 (2012).
  * Tang  _et al._ (2015) Y. Tang, T. Chen, and S. Yu, Chem. Commun. 51, 9018 (2015).
  * Corliss  _et al._ (1956) L. Corliss, N. Elliott, and J. Hastings, Phys. Rev. 104, 924 (1956).
  * Lines and Jones (1966) M. Lines and E. Jones, Phys. Rev. 141, 525 (1966).
  * Morosin (1970) B. Morosin, Phys. Rev. B 1, 236 (1970).
  * Danielian and Stevens (1961) A. Danielian and K. Stevens, Proc. Phys. Soc. (1958-1967) 77, 124 (1961).
  * Trimarchi  _et al._ (2018) G. Trimarchi, Z. Wang, and A. Zunger, Phys. Rev. B 97, 035107 (2018).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Oguchi  _et al._ (1983) T. Oguchi, K. Terakura, and A. R. Williams, Phys. Rev. B 28, 6443 (1983).
  * Raybaud  _et al._ (1997a) P. Raybaud, G. Kresse, J. Hafner, and H. Toulhoat, J. Phys. Condens. Matter 9, 11085 (1997a).
  * Raybaud  _et al._ (1997b) P. Raybaud, J. Hafner, G. Kresse, and H. Toulhoat, J. Phys. Condens. Matter 9, 11107 (1997b).
  * Tappero and Lichanot (1998) R. Tappero and A. Lichanot, Chem. Phys. 236, 97 (1998).
  * Hobbs and Hafner (1999) D. Hobbs and J. Hafner, J. Phys. Condens. Matter 11, 8197 (1999).
  * Rohrbach  _et al._ (2003) A. Rohrbach, J. Hafner, and G. Kresse, J. Phys. Condens. Matter 15, 979 (2003).
  * Fender  _et al._ (1968) B. Fender, A. Jacobson, and F. Wedgwood, J. Chem. Phys 48, 990 (1968).
  * Zunger  _et al._ (1990) A. Zunger, S.-H. Wei, L. G. Ferreira, and J. E. Bernard, Phys. Rev. Lett. 65, 353 (1990).
  * Varignon  _et al._ (2019) J. Varignon, M. Bibes, and A. Zunger, Nat. Commun. 10, 1 (2019).
  * Kraft and Greuling (1988) A. Kraft and B. Greuling, Cryst. Res. Technol. 23, 605 (1988).
  * McCammon (1991) C. McCammon, Phys. Chem. Miner. 17, 636 (1991).
  * Sweeney and Heinz (1993) J. S. Sweeney and D. L. Heinz, Phys. Chem. Miner. 20, 63 (1993).
  * Xiao  _et al._ (2015) G. Xiao, X. Yang, X. Zhang, K. Wang, X. Huang, Z. Ding, Y. Ma, G. Zou, and B. Zou, J. Am. Chem. Soc. 137, 10297 (2015).
  * Momma and Izumi (2008) K. Momma and F. Izumi, J. Appl. Crystallogr. 41, 653 (2008).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Physical review B 54, 11169 (1996b).
  * Csonka  _et al._ (2009) G. I. Csonka, J. P. Perdew, A. Ruzsinszky, P. H. T. Philipsen, S. Lebègue, J. Paier, O. A. Vydrov, and J. G. Ángyán, Phys. Rev. B 79, 155107 (2009).
  * Liechtenstein  _et al._ (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
  * Dudarev  _et al._ (1998) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
  * Ceperley and Alder (1980) D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * (35) See Supplemental Material at URL http://will.be.inserted.by.publisher, for the files of the generated AFM and PM MnS polymorphs’ structures used in this work.
  * Van De Walle  _et al._ (2002) A. Van De Walle, M. Asta, and G. Ceder, CALPHAD 26, 539 (2002).
  * Allen and Watson (2014) J. P. Allen and G. W. Watson, Phys. Chem. Chem. Phys. 16, 21016 (2014).
  * Huffman and Wild (1967) D. R. Huffman and R. L. Wild, Physical Review 156, 989 (1967).
  * Bradley and Cracknell (2009) C. Bradley and A. Cracknell,  _The mathematical theory of symmetry in solids: representation theory for point groups and space groups_ (Oxford University Press, 2009).
  * Sato  _et al._ (1997) H. Sato, T. Mihara, A. Furuta, M. Tamura, K. Mimura, N. Happo, M. Taniguchi, and Y. Ueda, Physical Review B 56, 7222 (1997).
  * Orgassa  _et al._ (1999) D. Orgassa, H. Fujiwara, T. C. Schulthess, and W. H. Butler, Phys. Rev. B 60, 13237 (1999).
  * Hubbard (1963) J. Hubbard, Proc. Math. Phys. Eng. Sci. 276, 238 (1963).
  * Bruus and Flensberg (2004) H. Bruus and K. Flensberg,  _Many-body quantum theory in condensed matter physics: an introduction_ (OUP Oxford, 2004).
  * Slater and Koster (1954) J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).