# Coupling of structure and magnetism to spin splitting in hybrid organic-
inorganic perovskites

Ravi Kashikar1<EMAIL_ADDRESS>D. DeTellem1 P. S. Ghosh2 Yixuan Xu3 S. Ma3 S.
Witanachchi1 Manh-Huong Phan1 S. Lisenkov1 I. Ponomareva1<EMAIL_ADDRESS>1\.
Department of Physics, University of South Florida, Tampa, Florida 33620, USA
2\. Glass & Advanced Materials Division, Bhabha Atomic Research Centre, Mumbai
400 085, India  3\. Department of Chemistry, University of North Texas, CHEM
305D, 1508 W Mulberry Street, Denton, Texas 76201, USA

###### Abstract

Hybrid organic-inorganic perovskites are famous for the diversity of their
chemical compositions, phases and phase transitions, and associated physical
properties. We use a combination of experimental and computational techniques
to reveal strong coupling between structure, magnetism, and spin splitting in
a representative of the largest family of hybrid organic-inorganic
perovskites: the formates. With the help of first-principles simulations, we
find spin splitting in both conduction and valence bands of [NH2NH3]Co(HCOO)3
induced by spin-orbit interactions, which can reach up to 14 meV. Our magnetic
measurements reveal that this material exhibits canted antiferromagnetism
below 15.5 K. The direction of the associated antiferromagnetic order
parameter is strongly coupled with the spin splitting already in the
centrosymmetric phase, allowing for the creation and annihilation of spin
splitting through the application of a magnetic field. Furthermore, the
structural phase transition into experimentally observed polar Pna21 phase
completely changes the aforementioned spin splitting and its coupling to
magnetic degrees of freedom. This reveals that in [NH2NH3]Co(HCOO)3 the
structure and magnetism are strongly coupled to spin splitting in a way that
allows for its manipulation through both magnetic and electric fields. As an
example, for a given point inside the Brillouin zone of centrosymmetric Pnma
phase of [NH2NH3]Co(HCOO)3, spin splitting can be turned on/off by aligning
the antiferromagnetic vector along certain crystallographic directions or
through inducing a polar phase by the application of an electric field. In
addition, the spin textures in [NH2NH3]Co(HCOO)3 are fully coupled to these
order parameter vectors. We believe that our findings offer an important step
toward fundamental understanding and practical applications of materials with
coupled properties.

##  I Introduction

The ability to couple structural and magnetic degrees of freedom to materials
properties is of tremendous scientific and technological importance. For
example, coupling between structural and electron spin degrees of freedom is
promising for applications in spintronicsBaltz  _et al._ (2018); Jungwirth
_et al._ (2016). One of the ways to access spin degrees of freedom is through
spin-orbit interaction, which breaks the Hamiltonian symmetry and splits spin-
degenerate bandsDresselhaus (1955); Even  _et al._ (2013). Spin-orbit
interactions in solids give origin to the so-called Rashba and Dresselhaus
effects, which originate from structural inversion asymmetry and bulk
inversion asymmetry, respectivelyRashba (1960); Dresselhaus (1955);
Soumyanarayanan  _et al._ (2016); Bihlmayer  _et al._ (2022). These effects
give origin to spin splitting in the electronic bands and, consequently, the
opportunity to access a particular spin channel and also lead to the emergence
of spin texturing in either real or reciprocal space. Certain types of spin
textures are known to enhance spin states’ lifetimes and, therefore, are
critical for applications in spintronicsSchliemann (2017); Tao and Tsymbal
(2018, 2021). Structural degrees of freedom are typically associated with
ionic positions and lattice vectors and can be manipulated through temperature
and external fieldsLines and Glass (2001). Ferroics usually offer great
opportunities for structural manipulations as they have a tendency to undergo
phase transitions and often couple to electric fields. Magnetic degrees of
freedom are usually associated with localized spins and may be either
disordered or ordered. Examples are paramagnetism and (anti) ferromagnetism,
respectively. While ferromagnetism results in spin dependent electronic band
structure, more subtle effect of spin splitting has been predicted from group
theoretical analysis in antiferromagnetic (AFM) materialsYuan  _et al._
(2020); Egorov  _et al._ (2021); Yuan  _et al._ (2021a, b, 2022); Hayami  _et
al._ (2019, 2020); Zhao  _et al._ (2022). Remarkably, some AFM materials can
exhibit spin splitting even in the absence of spin-orbit interactionsYuan  _et
al._ (2020). It appears that materials in the ferroic, or more precisely
multiferroic, family may have significant potential for coupling structure and
magnetic ordering to spin splitting. This motivated us to look into the hybrid
organic-inorganic perovskites class of materials, that share the chemical
formula ABX3, where A typically is an organic cation, while B is an inorganic
cation, while X could be either organic or inorganicLi  _et al._ (2020). In
particular, we will look into one of the largest families in this class:
hybrid formate perovskites with chemical formula AB(HCOO)3Popoola  _et al._
(2023). What makes them attractive for our search of materials with coupled
degrees of freedom, is that many of them exhibit AFM ordering and many undergo
phase transitions into polar space groups above room temperatureMa and Sun
(2020). For example, [NH2NH3][Co(HCOO)3] and [CH3NH2NH2][Mn(HCOO)3] exhibit
transitions into polar phases below 363 K and 309 K, respectively Chen  _et
al._ (2014); Maczka  _et al._ (2017). Ferroelectricity in [NH2NH3][Co(HCOO)3]
was reported in Ref.Ghosh  _et al._ (2022), while ferroelectric ordering in
[CH3NH2NH2][Mn(HCOO)3] was suggested from pyroelectric measurements Maczka
_et al._ (2017). It is important to realize that in such materials magnetic
ordering is mostly induced by the magnetic moments of the transition metal
atoms on the B site, while ferroelectric polarization originates from the
dipole moment of the A site, which may allow researchers to overcome the well-
known limitations of inorganic materials for coexistence of magnetic and
electric orderings Hill (2000). Indeed, we can find multiferroics with desired
features of AFM in [(CH3)2NH2]M(HCOO)3 (M=Mn, Fe, Co, Ni ) and NH4M(HCOO)3
(M=Mn, Co) family Ma and Sun (2020). Among them, one promising candidate for
coupling between multiple degrees of freedom to spin splitting is
[NH2NH3]Co(HCOO)3. It exhibits a high-temperature Pnma phase and undergoes a
phase transition to the polar Pna21 phase at 363 K Gunatilleke  _et al._
(2017). As a result, it is predicted to be ferroelectric at room temperature
with a remnant polarization of 2.6 $\mu$C/cm2 Ghosh  _et al._ (2022). The
structural phase transition in this compound is driven by the hydrogen bond
stabilization Ghosh  _et al._ (2022). [NH2NH3]Co(HCOO)3 also exhibits
competing, possibly antipolar, P212121 phase, which show spin-canted AFM long-
range-ordering, with Néel temperatures of 13.9 K Chen  _et al._ (2014).
Therefore, it may exhibit the desired coupling of structure and magnetism to
spin splitting.

Previously spin splitting coupled to the direction of polarization and AFM
order parameters, tunable by both magnetic ordering and polarization
direction, has been predicted in hybrid organic-inorganic perovskite TMCM-
MnCl3Lou  _et al._ (2020). Recently, Yananose et al., have examined
C(NH2)3]M(HCOO)3 (M = Cr, Cu) and predicted the persistent and irregular spin
texture Yananose  _et al._ (2023). Among nonmagnetic hybrid perovskites TMCM-
CdCl3 and MPSnBr3 exhibit persistent spin texture in the valence band and
conduction band respectivelyKashikar  _et al._ (2021, 2022). However, what is
presently missing is the understanding of the interplay between structural and
magnetic degrees of freedom in the establishment of spin splitting in hybrid
organic-inorganic perovskites.

The goal of this study is to utilize a combination of computational and
experimental tools to: (i) predict the existence of spin splitting in hybrid
formate perovskites, which is strongly coupled to both structural and magnetic
degrees of freedom; (ii) reveal the possibility to create, annihilate and
manipulate spin splitting through application of electric and magnetic fields;
(iii) rationalize our findings through decoupling contributions from different
degrees of freedom; (iv) propose practical ways to achieve spin splitting
manipulation in such materials.

##  II  Experimental and Computational Details

[NH2NH3]Co(HCOO)3 was synthesized using the mild solution method. The solvent
diffusion was carried out using Co(No3)2.6H2O (1.0 mmol, 0.291 g) in 5.0 mL
methanol solution. Later, it was layered onto 5.0 mL methanol with 95% (78
mmol, 3.2 mL) formic acid and 98% (6.2 mmol, 0.6 mL) hydrazine monohydrate.
The magenta-colored rectangular crystals are obtained after 24 hours which are
separated from the bulk phase and washed with reagent grade ethanol after the
mother solvent was removed. The perovskite structures are confirmed with
single-crystal X-ray diffraction measurements. Magnetic measurements were
carried out in a Quantum Design Dynacool Physical Property Measurement System
(PPMS) using a vibrating sample magnetometer (VSM) with the magnetic field
aligned parallel to the a-b plane. Magnetization versus applied field M(H) was
measured at 10 K, with an applied field of 9 T. Magnetization versus
temperature, $M(T)$, was measured between 5 K and 100 K, using the zero-field-
cooled warming protocol and under 0.1 T applied field.

For computations we use spin-polarized density functional theory as
implemented in the Vienna Abinitio Simulation Package (VASP) with a Perdew-
Burke-Ernzerhof (PBE) exchange-correlation functional Kresse and Furthmüller
(1996); Blöchl (1994); Perdew  _et al._ (1996). The energy cutoff is set to
700 eV for the plane-wave basis, and the Brillouin zone integration is carried
out using Gamma-centered k-mesh of 5 $\times$ 5 $\times$ 3\. We use Ueff = 2.0
eV correction for Co-d states as proposed by Dudarev et al. Dudarev  _et al._
(1998). The zero-damping DFT-D3 approach, proposed by Grimme et al., is
employed to incorporate the van der Waals interactionsGrimme  _et al._ (2010).
Structural relaxations are performed using a conjugate gradient algorithm
until the forces on each atom are less than 1 meV/Å, while stresses are less
than 0.001 GPa. The electronic structure calculations are performed with both
spin-orbit coupling (SOC) on and off. The electric polarization is computed
using the modern theory of polarization developed by King-Smith and
VanderbiltKing-Smith and Vanderbilt (1993).

Figure 1:  (a) Temperature dependence of the magnetization (blue) and its
derivative dM/dT (red); (b) Inverse susceptibility as a function of
temperature $\chi^{-1}$(T) and a fit to the Curie-Weiss law (red line); (c)
Magnetic field dependence of magnetization curve taken at 10 K, showing a
ferromagnetic hysteresis at low fields and a spin-flop transition around 6.5
T.

##  III Results and Discussions

Although AFM ordering has been reported for [NH2NH3]Co(HCOO)3 in P212121 phase
Chen  _et al._ (2014), no reports are available for the polar Pna21 phase. Our
magnetic data analysis for the Pna21 phase are given in Fig. 1 which suggest
canted AFM ordering below 15.5 K. Figure 1(a) shows the temperature dependence
of magnetization along with its derivative, $dM/dT$, revealing a magnetic
phase transition at 15.5 K. Inverse susceptibility $\chi^{-1}$ was calculated
from this data $H/M(T)$ and fit with the Curie-Weiss law, $\chi=C/(T-\theta)$
(Fig. 1(b)), where a Curie-Weiss temperature of $\theta$ = -14.49 K is
obtained. It is clear that at high temperatures in the paramagnetic regime,
the $\chi^{-1}$(T) data obeys the Curie-Weiss law. The negative value of
$\theta$ confirms AFM couplings between adjacent Co2+ ions and that the phase
transition at TN $\sim$ 15.5 K is AFM in nature. From the Curie constant
obtained through the fit (C = 1.707 cm3 K mol-1), an effective magnetic moment
of $\mu_{eff}\sim$3.296 $\mu$B/Co is obtained, using the relation
$C=N_{A}\mu_{B}^{2}\mu_{eff}^{2}/3k_{B}$. It is also worth noting from Fig.
1(b) that there is a sudden decline in $\chi^{-1}$ at $\sim$15 K, signaling
the characteristic of a spin-canted AFM system Wang  _et al._ (2021); Kuo  _et
al._ (2023); Cao  _et al._ (2017). To confirm the occurrence of spin canting
in [NH2NH3]Co(HCOO)3 the magnetic field dependence of magnetization M(H) was
measured at temperatures below the TN. Fig. 1(c) displays such a M(H) loop
taken at 10 K, confirming that there is a small canted AFM moment at low
fields of $\sim$0.4 emu/g, with a coercive field of 0.25 T. When the applied
magnetic field exceeds HSF $\sim$ 6.5 T, there is an increase in the
magnetization characteristic of a spin-flop transition (Fig. 1(c)) Wang  _et
al._ (2021). This provides further evidence for the spin-canted AFM order in
[NH2NH3]Co(HCOO)3 Chen  _et al._ (2014). The spin-canted AFM behavior was also
observed in Cu3(TeO4)(SO4).H2O Wang  _et al._ (2021) and
{Co(N3)(bpmb)(H2O)2.H2O}n Kuo  _et al._ (2023). Thus, at low-temperature
[NH2NH3]Co(HCOO)3 in Pna21 phase is a multiferroic as it exhibits the
coexistence of electric (spontaneous polarization) and magnetic orderings. As
a result, it could potentially exhibit spin splitting due to AFM ordering and
noncentrosymmetric nature.

To find out whether these two orderings are indeed capable of inducing spin
splitting in this material, we turn to simulations. To decouple the
contributions from magnetic and electric orderings, both centrosymmetric and
noncentrosymmetric structures are required. To obtain the noncentrosymmetric
phase, we carried out the structural relaxation of the Pna21 phase reported in
Ref. Ghosh  _et al._ (2022). The structure is visualized in Fig. 2(a). The
ground state was shown to be G-type AFM in Ref. Ghosh  _et al._ (2022).
Therefore, we use G-type AFM arrangements of magnetic moments on Co atoms but
simulate three different orientations of the magnetic moments: along the -x,
-y, and -z directions, as shown in Fig. 2(c). Note that the axes of our
Cartesian coordinate system are aligned along the crystallographic a, b, and c
axis. It should be noted that we do not find canted magnetic moment in
computations, which could be both due to our computational resolution and the
size of the supercell. The centrosymmetric structure for [NH2NH3]Co(HCOO)3 is
Pnma which is also a high-temperature phase reported in Ref.Ghosh  _et al._
(2022). However, it is a disordered phase with partial occupancies and,
therefore, cannot be used in simulations. To construct a prototype of such a
phase that does not have partial occupancy but still has experimental
relevance, we create a polarization reversal path for the Pna21 phase
following Ref.Kingsland  _et al._ (2021). The structure that corresponds to
zero polarization is the centrosymmetric Pnma phase. This structure is next
subjected to full structural relaxation, which preserved the space group. The
structure is visualized in Fig. 2(b).

Experimentally, such a phase could occur during polarization reversal at the
fields close to the coercive field. At zero Kelvin, the Pnma structure is 280
meV/f.u. higher in energy than the ground state. The structural relaxation has
been carried out for different types of magnetic orderings: ferromagnetic,
A-type, C-type, and G-type AFM orderings, and it was found that G-type AFM
ordering is the lowest in energy. Furthermore, we considered G-type AFM
orderings with different directions of magnetic moments, resulting in the
alignment of the AFM vector along x, y, or z-directions as shown in Fig. 2(c).
We found that AFM-Gx (Pnma) has the lowest energy followed by AFM-Gz (Pn′ma′)
and AFM-Gy (Pn′m′a), with 0.002 and 0.02 meV/f.u. higher energy, respectively.
The structures are provided in Ref.our

Figure 2:  Crystal structure of [NH2NH3]Co(HCOO)3 in Pna21 and Pnma phases in
different cross sections (a)-(b). Magnetic moment arrangements in G-AFM
ordering and corresponding magnetic space group (c).

The magnetic space groups Pn′m′a, Pn′ma′, and Pnma do not possess $\Theta I$
symmetry and, therefore, may exhibit SOC-induced spin splitting Yuan  _et al._
(2020). Here, $\Theta$ and $I$ are time reversal and spatial inversion
symmetries. Moreover, these magnetic space groups belong to types III, III,
and I, respectively, which allows for AFM-induced spin splitting, as proposed
in Ref. Yuan  _et al._ (2020). Following the classification of Ref. Yuan  _et
al._ (2021b) they belong to spin split type-4A and allow for spin splitting
with or without SOC, based on the symmetry considerations. To explore that, we
carry out DFT computations of electronic structure with SOC turned off and on.
Of course, when SOC is turned off, we use spin-polarized computations. Figure
3 presents our data for both cases. We find no spin splitting in the absence
of SOC. The partial density of states shown in the inset to Fig. 3(a) infers
that bands near the Fermi level are from Co-d states. In addition, the O-p
states are populated at the valence band region. Among the Co-d states, dxy
and d${}_{z^{2}}$ have double occupancy, and hence, they peak in the valence
band spectrum. The remaining d orbitals have single occupancy, thus offering a
magnetic moment of 2.68 $\mu_{B}$ at each Co site.

Figure 3: Band structure of [NH2NH3]Co(HCOO)3 in Pnma phase without (a) and
with (b)-(d) SOC.

Turning SOC on results in spin splitting (see Fig. 3(b)-(d)), which shows a
strong dependence on the direction of the AFM order parameter. To quantify
such dependence we compute spin splitting along the different directions of
the Brillouin zone and present it in Fig. 4.

Figure 4 confirms spin splitting dependence on the direction of AFM order
parameter throughout the Brillouin zone. For example, the spin splitting along
Y-$\Gamma$-Z direction in CB can be turned on or off by switching the AFM
vector between the y or the z direction. Similar trends are observed in VB.
Experimentally, the change of the AFM vector direction can be achieved through
the spin-flop mechanism, which opens the possibility to create/erase or tune
spin splitting by the application of a magnetic field. For instance, the field
at which the spin-flop transition occurs is 6.5 T for T = 10 K (Fig. 1(c)). We
note that the magnitude of spin splitting also shows a strong dependence on
the direction of the AFM vector. Thus, our simulations on the centrosymmetric
phase prototype for AFM [NH2NH3]Co(HCOO)3 revealed: (i) the presence of SOC-
induced spin splitting up to 4 meV in both CB and VB; (ii) strong dependence
of spin splitting magnitude on the direction of AFM vector, which allows for
both creation and erasing of spin splitting, as well as its manipulations by
the magnetic field.

Figure 4:  Spin splitting between the top two VBs (a) and (c) and bottom two
CBs (b) and (d) along high symmetry k-path in different phases of
[NH2NH3]Co(HCOO)3.  Figure 5:  Band structure of Pna21 [NH2NH3]Co(HCOO)3
without SOC (a), and with SOC in (b)-(d) for x, y, and z spin orientations of
AFM-G magnetic ordering.

Next, we turn to an investigation of the role that spatial asymmetry plays in
spin splitting and the aforementioned findings. Removal of the mirror plane
perpendicular to the c-axis results in the ground state structure of
[NH2NH3]Co(HCOO)3 with the Pna21 space group. Introducing G-type AFM ordering
with different directions of $\mathbf{G}$ vector followed by full structural
relaxation resulted in the structures with magnetic space groups Pna21,
Pn′a′21, and Pn′a2${}_{1}^{\prime}$ for AFM-Gx, AFM-Gy, and AFM-Gz,
respectively. Energetically, AFM-Gx has the lowest energy, followed by AFM-Gz
(0.18 meV/f.u.) and AFM-Gy (0.47 meV/f.u.). All these structures belong to the
spin splitting type-4B as per classification of Ref.Yuan  _et al._ (2021b),
and by symmetry, allow for spin splitting with and without SOC. Figure 5
presents the electronic structure calculations for both of these cases.

Again, we find that just AFM on its own is not capable of inducing spin
splitting in this material. However, the introduction of SOC results in such
spin splitting for all directions of the AFM vector considered here. Moreover,
just like before, we find that both magnitudes of spin splitting in its
location in the Brillouin zone exhibit strong dependence on the direction of
AFM vector. This is quantified in Fig. 4. The magnitude of spin splitting is
up to 7 meV in CB and up to 4 meV in VB, which offers enhancement with respect
to centrosymmetric phases. We find a direct band gap of 2.1 eV, which is
likely to underestimate the true band gap, as known for DFT. Remarkably, the
bands are spin-splitted at VBM and CBM (located at $\Gamma$ point) despite
being in time-reversal invariant momenta point. This is the consequence of the
AFM ordering, discussed in Ref.Yuan  _et al._ (2020, 2021b). It is also a
rare, but technologically significant property. Comparison between our data
for centrosymmetric and noncentrosymmetric phases with different magnetic
order parameter directions (see Fig. 4) reveals that [NH2NH3]Co(HCOO)3
exhibits strong coupling between multiple degrees of freedom: structural,
magnetic, and spin splitting. For example, in $\Gamma$ point of VB spin
splitting is zero for AFM-Gx. It can be turned on by aligning the AFM vector
along the y-axis, via a spin-flop mechanism induced by a magnetic field.
Inducing the centrosymmetric phase through electric field or temperature
manipulations will result in near disappearance of spin splitting. Figure 4
also reveals greater overlap in the spin splitting peaks between space groups
containing time-reversed symmetry elements as compared to those that don’t,
suggesting strong contributions to spin splitting from AFM.

Figure 6:  Spin splitting along various k-paths in the Brillouin of the polar
phases of [NH2NH3]Co(HCOO)3 for the lowermost conduction bands (a)-(c) and
uppermost valence bands (d)-(f) for x, y, and z orientations for the AFM order
parameter. Comparison of uppermost VB spin splitting for centrosymmetric and
polar phases (g)-(i).

To find out whether even larger values of spin splitting can be achieved in
the vicinity of VBM and CBM as well as to assess anisotropy of the spin
splitting, we explore several directions in the vicinity of CBM and VBM.
Figure 6 gives spin splitting along the directions investigated. The largest
value can reach 14 meV in VB and occurs between $\Gamma$-T. To estimate the
contribution to spin splitting from antiferromagnetism, we overlap the data
for centrosymmetric and noncentrosymmetric structures in Figs. 6(g)-(i). For
this analysis, the path in the Brillouin zone associated with the largest spin
splitting is chosen. From the comparison, we conclude that AFM makes a
significant contribution to spin splitting in all cases. Further analysis can
be carried out to decouple structural contribution to spin splitting into
contributions from the [NH2NH3]-1 molecule and the Co(HCOO)3 framework. To
compute the contribution due to the molecule, we keep the framework in its
Pnma structure while molecules are in their Pna21 configuration. To compute
contributions due to the framework, we keep molecules in their Pnma positions
while the framework is in Pna21 structure. The decomposition is shown in Fig.
6(g). We can see that the largest contribution to the spin splitting comes
from the framework. For a tentative explanation of such a finding, we
decompose the polarization of the Pna21 phase into contributions from the
molecule and the framework. The framework contributes 73% to the total
polarization, while the molecule contributes 27%. Therefore, it may be that
structural distortions of the framework make the largest contribution to the
local electric field, responsible for the spin-orbit interaction and
associated spin splitting.

Figure 7:  Spin textures around the CBM (a)-(c) and VBM (d)-(f) in various
k-planes for the polar phases of [NH2NH3]Co(HCOO)3.

Figure 7 shows spin textures in CB and VB for Pna21 type structures. We find
that the type of spin texture is fully coupled to the direction of the AFM
vector. Interestingly, for the Pna21 structure with AFM-Gx orientation, the
spin textures are mostly persistent in a large portion of the Brillouin zone
which is highly desirable for the application of spintronics. Switching the
direction of the AFM order parameter changes the type of spin textures and,
therefore, will cause relaxation.

In order to find out whether the discussed findings are specific to
[NH2NH3]Co(HCOO)3 or likely to occur in other materials from the same family
we have repeated calculations for [NH2NH3]Mn(HCOO)3. This material exhibits a
canted AFM phase below 7.9 K and undergoes a polar (Pna21) to nonpolar (Pnma)
phase transition at 355 KChen  _et al._ (2014), similar to the case of
[NH2NH3]Co(HCOO)3. All computational approaches remained the same except that
we used Ueff=4.0 eV for Mn. The structures are provided in Ref. our . The
simulated electronic structure for Pnma and Pna21 phases are shown in Fig. S1
and S2 of the Supplementary Material for -x, -y and -z orientation of the AFM
order parameter. In both cases, we observe spin splitting in the conduction
and valence bands. The magnitude of the spin splitting is 4 meV, near $\Gamma$
for VB for both phases. The magnitude of the spin splitting varies with
magnetic moment directions. The dependence of spin splitting in VB and CB on
the direction of the AFM vector of both phases is shown in Fig. S3. The spin
texture around the VBM for [NH2NH3]Mn(HCOO)3 is persistent in nature (see Fig.
S4 of the supplementary material). Thus, [NH2NH3]Mn(HCOO)3 also exhibits the
main features established for [NH2NH3]Co(HCOO)3. These include strong coupling
between structural, magnetic, and spin degrees of freedom; the presence of
spin splitting even in the centrosymmetric phase of the material; and
persistent spin textures.

##  IV Summary

In summary, we have synthesized [NH2NH3]Co(HCOO)3 in the polar Pna21 space
group and found it to be a canted antiferromagnet below 15.5 K. DFT
computations have then been utilized to predict spin splitting in this
compound. Interestingly, SOC-induced spin splitting exists already in the
prototype of the centrosymmetric phase of the material and even in the time-
reversal invariant momenta points of the Brillouin zone. These findings
originate from the antiferromagnetic ordering. The spin splitting couples
strongly to the direction of the AFM vector and can be turned on or off
through AFM vector reorientation, which opens a way to spin splitting
tunability by the application of a magnetic field, for example, via a spin-
flop mechanism. Phase transition from the centrosymmetric to the polar phase
of the material results in a dramatic change in the spin splitting landscape
in the Brillouin zone, thus offering additional tunability by the electric
field or temperature. For example, a macroscopically centrosymmetric phase can
be achieved at low temperatures during electric-field-induced polarization
reversal. Just like in the prototypical centrosymmetric phase of the material,
the polar phase exhibits a strong dependence of the spin splitting on the
direction of the AFM vector. Thus, our work reveals that the low-temperature
polar phase of [NH2NH3]Co(HCOO)3 (and [NH2NH3]Mn(HCOO)3) is magnetically
ordered and exhibits strong coupling between structure, magnetic ordering,
spin splitting and spin textures. These findings enrich our fundamental
understanding of property manipulation and materials functionality and are
likely to stimulate further research. From a practical viewpoint, both
structure and magnetic ordering can be manipulated by external fields,
promising novel applications in spintronics.

## Acknowledgment

This work was supported by the National Science Foundation under Grant No.
EPMD-2029800.

## References

  * Baltz  _et al._ (2018) V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Rev. Mod. Phys. 90, 015005 (2018).
  * Jungwirth  _et al._ (2016) T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich, Nature Nanotechnology 11, 231 (2016).
  * Dresselhaus (1955) G. Dresselhaus, Physical Review 100, 580 (1955).
  * Even  _et al._ (2013) J. Even, L. Pedesseau, J.-M. Jancu, and C. Katan, The Journal of Physical Chemistry Letters 4, 2999 (2013).
  * Rashba (1960) E. Rashba, Sov. Phys.-Solid State 2, 1109 (1960).
  * Soumyanarayanan  _et al._ (2016) A. Soumyanarayanan, N. Reyren, A. Fert, and C. Panagopoulos, Nature 539, 509 (2016).
  * Bihlmayer  _et al._ (2022) G. Bihlmayer, P. , D. V. Vyalikh, E. V. Chulkov, and A. Manchon, Nature Reviews Physics 4, 642 (2022).
  * Schliemann (2017) J. Schliemann, Reviews of Modern Physics 89, 011001 (2017).
  * Tao and Tsymbal (2018) L. Tao and E. Y. Tsymbal, Nature communications 9, 2763 (2018).
  * Tao and Tsymbal (2021) L. Tao and E. Y. Tsymbal, Journal of Physics D: Applied Physics 54, 113001 (2021).
  * Lines and Glass (2001) M. E. Lines and A. M. Glass,  _Principles and Applications of Ferroelectrics and Related Materials_ (Oxford University Press, 2001).
  * Yuan  _et al._ (2020) L.-D. Yuan, Z. Wang, J.-W. Luo, E. I. Rashba, and A. Zunger, Phys. Rev. B 102, 014422 (2020).
  * Egorov  _et al._ (2021) S. A. Egorov, D. B. Litvin, and R. A. Evarestov, The Journal of Physical Chemistry C 125, 16147 (2021).
  * Yuan  _et al._ (2021a) L.-D. Yuan, Z. Wang, J.-W. Luo, and A. Zunger, Phys. Rev. B 103, 224410 (2021a).
  * Yuan  _et al._ (2021b) L.-D. Yuan, Z. Wang, J.-W. Luo, and A. Zunger, Phys. Rev. Materials 5, 014409 (2021b).
  * Yuan  _et al._ (2022) L.-D. Yuan, X. Zhang, C. Mera, and A. Zunger, “Uncovering hidden spin polarization of energy bands in antiferromagnets,” (2022).
  * Hayami  _et al._ (2019) S. Hayami, Y. Yanagi, and H. Kusunose, Journal of the Physical Society of Japan 88, 123702 (2019).
  * Hayami  _et al._ (2020) S. Hayami, Y. Yanagi, and H. Kusunose, Phys. Rev. B 102, 144441 (2020).
  * Zhao  _et al._ (2022) H. J. Zhao, X. Liu, Y. Wang, Y. Yang, L. Bellaiche, and Y. Ma, Phys. Rev. Lett. 129, 187602 (2022).
  * Li  _et al._ (2020) W. Li, A. Stroppa, Z. Wang, and S. Gao,  _Hybrid Organic‐Inorganic Perovskites_ (2020).
  * Popoola  _et al._ (2023) A. Popoola, P. S. Ghosh, M. Kingsland, R. Kashikar, D. DeTellem, Y. Xu, S. Ma, S. Witanachchi, S. Lisenkov, and I. Ponomareva, The Journal of Chemical Physics 159, 074702 (2023).
  * Ma and Sun (2020) Y. Ma and Y. Sun, Journal of Applied Physics 127, 080901 (2020).
  * Chen  _et al._ (2014) S. Chen, R. Shang, K.-L. Hu, Z.-M. Wang, and S. Gao, Inorg. Chem. Front. 1, 83 (2014).
  * Maczka  _et al._ (2017) M. Maczka, A. Gagor, M. Ptak, W. Paraguassu, T. A. da Silva, A. Sieradzki, and A. Pikul, Chemistry of Materials 29, 2264 (2017), https://doi.org/10.1021/acs.chemmater.6b05249 .
  * Ghosh  _et al._ (2022) P. S. Ghosh, D. DeTellem, J. Ren, S. Witanachchi, S. Ma, S. Lisenkov, and I. Ponomareva, Phys. Rev. Lett. 128, 077601 (2022).
  * Hill (2000) N. A. Hill, The Journal of Physical Chemistry B 104, 6694 (2000).
  * Gunatilleke  _et al._ (2017) W. D. C. B. Gunatilleke, K. Wei, Z. Niu, L. Wojtas, G. Nolas, and S. Ma, Dalton Trans. 46, 13342 (2017).
  * Lou  _et al._ (2020) F. Lou, T. Gu, J. Ji, J. Feng, H. Xiang, and A. Stroppa, npj Computational Materials 6, 114 (2020).
  * Yananose  _et al._ (2023) K. Yananose, E. R. Clark, P. J. Saines, P. Barone, A. Stroppa, and J. Yu, Inorganic Chemistry 62, 17299 (2023).
  * Kashikar  _et al._ (2021) R. Kashikar, P. S. Ghosh, S. Lisenkov, B. R. K. Nanda, and I. Ponomareva, Phys. Rev. B 104, 235132 (2021).
  * Kashikar  _et al._ (2022) R. Kashikar, P. S. Ghosh, S. Lisenkov, A. Stroppa, and I. Ponomareva, Phys. Rev. Materials 6, 104603 (2022).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * 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).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, The Journal of Chemical Physics 132, 154104 (2010).
  * King-Smith and Vanderbilt (1993) R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  * Wang  _et al._ (2021) Z.-C. Wang, K. Thanabalasingam, J. P. Scheifers, A. Streeter, G. T. McCandless, J. Gaudet, C. M. Brown, C. U. Segre, J. Y. Chan, and F. Tafti, Inorganic Chemistry 60, 10565 (2021), pMID: 34176270.
  * Kuo  _et al._ (2023) Y.-L. Kuo, H.-K. Liu, and C.-I. Yang, Inorganics 11 (2023).
  * Cao  _et al._ (2017) C. Cao, S.-J. Liu, S.-L. Yao, T.-F. Zheng, Y.-Q. Chen, J.-L. Chen, and H.-R. Wen, Crystal Growth & Design 17, 4757 (2017).
  * Kingsland  _et al._ (2021) M. Kingsland, P. S. Ghosh, S. Lisenkov, and I. Ponomareva, The Journal of Physical Chemistry C 125, 8794 (2021).
  * (42) https://github.com/USFmatscilab/Co-Mn-Formate-AFM_G.

Supplementary Material

Figure S1: Band structure of [NH2NH3]Mn(HCOO)3 in Pnma phase with SOC.  Figure
S2: Band structure of [NH2NH3]Mn(HCOO)3 in Pna21 phase with SOC.  Figure S3:
Spin-splitting between the top two valence bands and bottom two conduction
bands along high symmetry k-path in different phases of [NH2NH3]Mn(HCOO)3.
Figure S4:  Spin textures of [NH2NH3]Mn(HCOO)3around the valence band maximum
in various k-planes.


#  Simulating $\mathbb{Z}_{2}$ Lattice Gauge Theory with the Variational
Quantum Thermalizer

Michael Fromm<EMAIL_ADDRESS>Institut für Theoretische Physik,
Goethe-Universität Frankfurt
Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany  Owe Philipsen
<EMAIL_ADDRESS>Institut für Theoretische Physik, Goethe-
Universität Frankfurt
Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany  Michael Spannowsky
<EMAIL_ADDRESS>University of Durham, Science Laboratories
South Rd, Durham DH1 3LE, UK  Christopher Winterowd<EMAIL_ADDRESS>frankfurt.de Institut für Theoretische Physik, Goethe-Universität Frankfurt
Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany

###### Abstract

The properties of strongly-coupled lattice gauge theories at finite density as
well as in real time have largely eluded first-principles studies on the
lattice. This is due to the failure of importance sampling for systems with a
complex action. An alternative to evade the sign problem is quantum
simulation. Although still in its infancy, a lot of progress has been made in
devising algorithms to address these problems. In particular, recent efforts
have addressed the question of how to produce thermal Gibbs states on a
quantum computer. In this study, we apply a variational quantum algorithm to a
low-dimensional model which has a local abelian gauge symmetry. We demonstrate
how this approach can be applied to obtain information regarding the phase
diagram as well as unequal-time correlation functions at non-zero temperature.

Hamiltonian lattice gauge theory, Quantum computing

††preprint: IPPP/23/25

##  I Introduction

Much interest has been generated recently regarding the potential of quantum
computing and quantum simulation to address long-standing problems in high-
energy physics [1, 2]. Traditionally, lattice gauge theory has been applied in
the Euclidean formulation, where advances in computing power have led to great
success in studying properties of strongly-interacting matter in thermal
equilibrium. However, due to the notorious sign problem, there has been a
growing interest in the Hamiltonian formulation and its application to address
non-zero chemical potential and real-time dynamics. For digital quantum
computers in the medium term, one problem that has arisen is how to produce
thermal states.

In recent years, there has already been an explosion of progress in tackling
this problem. Each approach to producing a thermal state has had to address
the issue of how to produce a mixed state. This can typically be done by
enlarging the system size and then entangling the ancillary system with the
original system. One idea which takes this approach to create a Gibbs state
via purification is the so-called thermofield double states method [3].
Another interesting method is known as active cooling [4]. This takes a
general initial state, couples it to a thermal heat bath, and uses the concept
of the “Maxwell demon" to obtain a thermal state at the desired temperature.
Alternatively, as the Boltzmann operator is non-unitary, one can also attempt
to directly approximate its action on quantum states. Using judiciously chosen
initial states produced by Haar-random circuits, one can use this
approximation of the Boltzmann operator to obtain Gibbs states for spin models
[5] as well as gauge theories [6].

Figure 1:  Schematic depiction of the hybrid approach used in the VQT
algorithm. The action of the variational circuits $VQC_{1}$ and $VQC_{2}$ and
the measurements are performed on the quantum device. The results of the
measurements are then sent to the classical device, which is depicted in grey.
This information allows one to construct the cost function and propose a new
value for the variational parameters using one’s optimizer of choice.

One can also study the problem of producing thermal states using variational
methods. One such method, the variational quantum thermalizer (VQT),
introduced in [7], serves as a natural extension of the variational quantum
eigensolver (VQE) [8, 9] to the case of non-zero temperature. The latter uses
the variational principle at $T=0$ to approximate the ground-state of a
Hamiltonian $\hat{H}$ by a parametrized ansatz
$\ket{\psi_{\boldsymbol{\xi}}}$, subject to

$\displaystyle
E_{0}\leq\bra{\psi_{\mathbf{\boldsymbol{\xi}}}}\hat{H}\ket{\psi_{\mathbf{\boldsymbol{\xi}}}}.$
(1)

Similarily, for the case of $T\neq 0$, VQT addresses the problem of finding a
variational approximation $\hat{\rho}_{var}$ to the thermal state,
$\hat{\rho}_{Gibbs}\sim e^{-\beta\hat{H}}$ for inverse temperature
$\beta=1/T$. Similar to its zero temperature analogue in Eq.(1), this problem
can also be reformulated as an optimization problem for the free energy with
respect to variational parameters $\boldsymbol{\xi}$, such that

$\displaystyle
F\equiv\min_{\mathbf{\boldsymbol{\xi}}}F(\mathbf{\boldsymbol{\xi}})=\min_{\mathbf{\boldsymbol{\xi}}}\left(\mathrm{tr}(\hat{H}\hat{\rho}_{\boldsymbol{\xi}})+T\mathrm{tr}(\hat{\rho}_{\boldsymbol{\xi}}\log{\hat{\rho}_{\boldsymbol{\xi}}})\right)\,,$
(2)

where the second term on the right-hand side involves the von Neumann entropy
$S=-\mathrm{tr}(\hat{\rho}\log{\hat{\rho}})$. As the density operator appears
non-linearly in the expression for the entropy and thus is not simply
determined by a quantum-mechanical expectation value, several avenues have
been taken in the literature: The original formulation of VQT [7] employs
classical sampling of the distribution corresponding to
$\rho_{\boldsymbol{\xi}}$ which necessitates a model ansatz. Through this
ansatz, the entropy is then given _classically_ by a closed-form expression.
To prepare the latent distribution corresponding to $\rho_{var}$ in a
_quantum_ circuit, the noise-assisted variational quantum thermalizer [10]
(NAVQT) uses (simulated) parameterized depolarization gates that allow for the
preparation of a mixed state, again with a closed-form, approximate expression
for the entropy in terms of the noise level $\lambda$, which thus becomes a
variational parameter. Of particular interest to us, due to both its
simplicity and versatility, is a variant of the VQT introduced in [11], see
Fig.1. While the first two examples used a classical sampling of input states
[7] or stochastic mixtures of unitary circuits [10], respectively, to create
mixed quantum states, the approach in [11] utilizes intermediate projective
measurements. As measurements are non-unitary operations, they act as a source
of entropy, ultimately leading to a variational approximation of the thermal
state111For a version of this algorithm using entanglement of the system with
$n_{q}$ ancillary qubits followed by a projective measurement on the latter,
see [12]. .

The subject of our study is the Hamiltonian formulation of the
$\mathbb{Z}_{2}$ lattice gauge theory (LGT) in $1+1d$
($\mathbb{Z}_{2}^{1+1}$), in which we investigate the application of the VQT
to this theory. This is performed in its formulation with a gauge-redundant
Hilbert space [13] and its recent, resource-efficient form derived
systematically in [14, 15]. In addition, we investigate several observables of
interest at both non-zero temperature $T$ and non-zero chemical potential
$\mu$.

Figure 2:  Depiction of a gauge-invariant state for $Z_{2}^{1+1}$ LGT with
staggered fermions on a lattice with linear extent $L$. The links are all in
the eigenstate of $X$ with eigenvalue $+1$, the even sites are empty
(fermions), and the odd sites are filled (anti-fermions). We note here that
the links at the boundary are set to unity.

@C=1.0em @R=0.2em @!R
q_0 : & q_0 : X [0em]2 R_Z ($\Delta_{0}$) [0em]2 2 H R_Z ($\Delta_{3}$) R_Z
($\Delta_{4}$) H 2 [0em]2
q_1 : q_1 : H R_X ($\Delta_{1}$) -1 -1
q_2 : q_2 : R_Z ($\Delta_{2}$) -2 -2

Figure 3:  Depiction of a single-layer of a variational quantum circuit in the
gauge-redundant formulation. The layer is composed of one sublayer of single-
qubit rotation gates and a sublayer containing the decomposed three-qubit gate
which is formed from the gauge-matter term in the gauge-redundant formulation.
All relevant gates contain variational parameters $\Delta_{i}$.

##  II Background

###  II.1 Variational Quantum Thermalizer

As displayed in Fig.1, starting from an $n_{q}$-qubit initial state
$|\psi_{0}\rangle$ with density matrix
$\rho_{0}=|\psi_{0}\rangle\langle\psi_{0}|$, we set up a concatenation of two
separate variational circuits $\mathrm{VQC}_{1}$ and $\mathrm{VQC}_{2}$. The
first variational circuit, $\mathrm{VQC}_{1}$, creates the state
$|\psi_{\mathrm{VQC}_{1}}\rangle$, which is given by

$\displaystyle\ket{\psi_{\mathrm{VQC}_{1}}}\equiv
U_{1}(\boldsymbol{\theta})\ket{\psi_{0}}=\sum_{i=0}^{2^{n_{q}}-1}a_{i}(\boldsymbol{\theta})\ket{i}\,.$
(3)

An intermediate measurement in the computational basis collapses this state
and statistically yields the latent density matrix

$\displaystyle\rho_{1}(\boldsymbol{\theta})=\sum_{i=0}^{2^{n_{q}}-1}|a_{i}(\boldsymbol{\theta})|^{2}\ket{i}\bra{i}\,.$
(4)

The subsequent application of the variational circuit $\mathrm{VQC}_{2}$
transforms this into

$\displaystyle\rho_{2}(\boldsymbol{\theta},\boldsymbol{\phi})$
$\displaystyle=$ $\displaystyle
U_{2}(\boldsymbol{\phi})\rho_{1}(\boldsymbol{\theta})U^{\dagger}_{2}(\boldsymbol{\phi}),$
(5) $\displaystyle\equiv$
$\displaystyle\sum_{i=0}^{2^{n_{q}}-1}p_{i}(\boldsymbol{\theta})\ket{\psi_{i}(\boldsymbol{\phi})}\bra{\psi_{i}(\boldsymbol{\phi})},$

after which a measurement of the energy is performed. The von Neumann entropy
is then calculated as
$S=-\sum_{i}p_{i}\log{p_{i}}$, with $p_{i}\approx n_{i}/n_{s}$, where $n_{i}$
represents the counts for state $i$ in the computational basis and $n_{s}$
represents the number of shots. By varying the parameter set
$\boldsymbol{\xi}=(\boldsymbol{\theta},\boldsymbol{\phi})$ in an optimal way,
one obtains a variational approximation to the Gibbs state
$\hat{\rho}_{Gibbs}\sim e^{-\beta\hat{H}}$. The fact that the circuit is
broken up into two parts is intended for noisy intermediate-scale quantum
(NISQ) devices, as the depth of each variational circuit can be kept
relatively shallow. This approach, however, naturally leads to a truncation in
the number of states in Eq.(5), as both the number of shots, $n_{s}$, and the
memory necessary to store the histogram describing the counts scales with the
dimension of the Hilbert space of the underlying problem. We pursue this
discussion in further detail in Sect.III.2.4. In [11], the algorithm was
demonstrated to reach convergence for the transverse field Heisenberg-Model in
$d=1,2$. Applying qVQT to gauge theories will introduce the need for a gauge-
invariant algorithm implementation. We show below for the case of
$\mathbb{Z}_{2}^{1+1}$ how this can be achieved.

###  II.2 $\mathbb{Z_{2}}$ LGT in $1+1d$

This section briefly outlines our two approaches to $\mathbb{Z}_{2}$ lattice
gauge theory in the Hamiltonian formulation. In the first approach, introduced
in [13] and referred to as the gauge-redundant formulation, local gauge
invariance must be enforced on the states to define the physical Hilbert
space. This can be achieved, for example, by introducing a penalty term as in
[6]. In the second approach, used in [14] and referred to as the resource-
efficient formulation, the matter states are decoupled, and the constraint
imposed by Gauss’ law is built into the Hamiltonian.

####  II.2.1 Gauge-Redundant Formulation

In the so-called group basis, $\mathbb{Z}_{2}^{1+1}$ LGT with a single flavour
of staggered fermion is formally given by the Hamiltonian
$H=H_{E}+H_{GM}+H_{M}+H_{\mu}$ with

$\displaystyle H_{E}$ $\displaystyle=$
$\displaystyle-\epsilon\sum_{n=0}^{L-2}X_{n},$ (6) $\displaystyle H_{GM}$
$\displaystyle=$
$\displaystyle\frac{J}{2}\sum_{n}^{L-2}\left(\sigma^{x}_{n}Z_{n}\sigma^{x}_{n+1}+\sigma^{y}_{n}Z_{n}\sigma^{y}_{n+1}\right),$
(7) $\displaystyle H_{M}$ $\displaystyle=$
$\displaystyle\frac{m}{2}\sum_{n=0}^{L-1}(-1)^{n}\sigma_{n}^{z},$ (8)
$\displaystyle H_{\mu}$ $\displaystyle=$
$\displaystyle-\frac{\mu}{2}\sum_{n=0}^{L-1}\sigma_{n}^{z}\,.$ (9)

In the above equations, the fermionic creation and annihilation operators have
been replaced by Pauli matrices living at the sites of the lattice via the
usual Jordan-Wigner transformation. Here we have adopted the familiar
convention of denoting Pauli operators acting on the link Hilbert space of the
gauge bosons by $\\{X,Y,Z\\}$, and those acting on the fermionic matter space
by $\\{\sigma^{x},\sigma^{y},\sigma^{z}\\}$. The occurrence of both gauge-
bosonic and fermionic degrees of freedom implies a gauge redundant Hilbert
space. With our choice of basis, a gauge transformation at site $n$ is defined
by [16]

$\displaystyle{\Theta}_{g}(n)$ $\displaystyle=$ $\displaystyle
X_{n}X_{n-1}e^{i\pi\left[\psi^{\dagger}_{n}\psi_{n}+\frac{1}{2}((-1)^{n}-1)\right]},$
(10) $\displaystyle=$ $\displaystyle(-1)^{n+1}X_{n}X_{n-1}\sigma_{n}^{z}\,.$

We illustrate a gauge-invariant state $|\psi\rangle$ of the lattice system in
Fig.2, where the fermionic states are in the computational basis, and the
gauge bosons are in the $X$-eigenbasis, $X\ket{\pm}=\pm\ket{\pm}$. This
subspace of the total Hilbert space is needed to compute physical observables.

####  II.2.2 Resource-Efficient Formulation

It turns out that one can eliminate the fermionic matter fields while
simultaneously removing the Hilbert space redundancy of the theory. This is
accomplished in a two-step process whereby the fermions are first mapped to
hard-core bosons by a unitary transformation and are then decoupled from the
gauge bosons by a second unitary transformation. This method is, in fact, a
very powerful formalism which can also be applied in arbitrary dimensions to
continuous gauge groups such as $SU(N)$ and $U(N)$ [17]. For $\mathbb{Z}_{2}$
LGT in $1+1$ dimensions, the Hamiltonian in this formulation is as follows

$\displaystyle H^{(2)}_{E}$ $\displaystyle=$
$\displaystyle-\sum_{n=0}^{L-2}\left(\epsilon Z_{n}+\frac{J}{2}Y_{n}\right),$
(11) $\displaystyle H^{(2)}_{GM}$ $\displaystyle=$
$\displaystyle-\frac{J}{2}\sum_{n=1}^{L-3}Z_{n-1}Y_{n}Z_{n+1}$ (12)
$\displaystyle-\frac{J}{2}\left(Y_{0}Z_{1}+Z_{L-3}Y_{L-2}\right),$
$\displaystyle H^{(2)}_{M}$ $\displaystyle=$
$\displaystyle-\frac{m}{2}\sum_{n=1}^{L-1}Z_{n-1}Z_{n}-\frac{m}{2}\left(Z_{0}+Z_{L-2}\right),$
(13) $\displaystyle H^{(2)}_{\mu}$ $\displaystyle=$
$\displaystyle\frac{\mu}{2}\sum_{n=1}^{L-2}\varepsilon_{n}Z_{n-1}Z_{n}$ (14)
$\displaystyle+\frac{\mu}{2}\left(Z_{0}+\varepsilon_{L-1}Z_{L-2}\right)\,.$

One notices that the Hamiltonian obtained by performing the above-mentioned
transformations only vaguely resembles the original $1+1$-dimensional Kogut-
Susskind Hamiltonian. For example, the gauge-matter interaction in eq. (12)
couples sets of three adjacent links. Similarly, for the term originating from
the staggered mass in Eq. (13), one obtains a coupling between two adjacent
links. With this resource-efficient formulation, one can simulate larger
systems without worrying about the elements of the variational circuit
preserving gauge invariance.

Figure 4: _(Left)_ The fidelity of the approximation to the Gibbs state
produced by VQT on the classical simulator and $\Delta F$ as a function of
$T/m$ at $\mu=0$ on a system with $N=4$ matter sites, corresponding to
$n_{q}=7$ qubits in the gauge-redundant formulation. _(Right)_ The chiral
condensate $\langle\bar{\psi}\psi\rangle$ calculated on the same data sets
where for clarity the exact result has also been included.

##  III Implementation and Results

The $\mathbb{Z}_{2}$ lattice gauge theory has a chiral symmetry which is
spontaneously broken at low $T$. At large temperatures as well as at large
values of the chemical potential, this symmetry is restored. To benchmark our
approximation to the Gibbs state, we calculate the chiral condensate,
$\langle\,\bar{\psi}\psi\,\rangle$, which is sensitive to restoring chiral
symmetry. As this is a local, time-independent observable, one would like a
second, time-dependent observable that measures correlations in the system. We
chose to study the unequal-time density-density correlator, which provides
information about charge dynamics in the system.

###  III.1 Gauge redundant formulation

In this formulation of the theory, gauge invariance has to be considered
explicitly. Apart from the gauge-invariant initial state (c.f. Fig. 2), this
not only applies to the variational circuits $\mathrm{VQC_{1/2}}$ but also to
the intermediate measurement. A schematic depiction of the circuit
illustrating the above-mentioned components is shown in Fig. 1. For
$\mathrm{VQC_{1/2}}$, a gauge-invariant gate set can be constructed by
following the general considerations laid out in [18]. This consists of simply
taking the Pauli strings in our Hamiltonian Eq.(6) as elementary building
blocks. Due to the simplicity of the model, the gate set will hence decompose
into parameter-dependent single qubit $R_{x}$ and $R_{z}$ gates, together with
a 3-qubit gate $\sim e^{i\sigma^{x}Z\sigma^{x}+i\sigma^{y}Z\sigma^{y}}$
stemming from the hopping term in Eq.(6). Following the procedure presented in
[19], the latter can be easily decomposed into single-qubit and entangling
gates. In Fig. 3, one observes the resulting circuit layer consisting of a
sublayer of single-qubit gates along with a sublayer of the decomposed multi-
qubit gate for a three-qubit system where the $\Delta_{i}$ represents generic
variational parameters. Several such layers then yield the $\mathrm{VQC}_{i}$,
with parameter sets $\boldsymbol{\theta}$ and $\boldsymbol{\phi}$ for $i=1,2$,
respectively (Eq.(5)). With the composition of the variational quantum
circuits at hand, we still need to define the figures of merit quantifying the
quality of the optimization step depicted in the classical block of Fig. 1. In
addition to the relative difference in free energy, $\Delta
F=F_{\rho_{var}}/F_{\rho_{\text{Gibbs}}}-1$, we take the fidelity [20]

$\displaystyle
F_{2}(\rho_{1},\rho_{2})=\left(\operatorname{Tr}\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho_{2}}}\right)^{2},$
(15)

as a distance measure between two quantum states represented by $\rho_{1}$ and
$\rho_{2}$, respectively. In Fig. 4, one sees the results of the classically
simulated variational optimization process for our gauge-redundant theory for
$N=4$ matter sites corresponding to $n_{q}=7$ qubits using open boundary
conditions with a two-layer ansatz in the $\mathrm{VQC}_{i}$. Simulating at
$a\epsilon=am=0.5,a\mu=0$ and $n_{s}=10^{4}$ for varying $T/m$, chiral
symmetry is explicitly broken in our model. On the left, we plot the
performance measures defined above, where each point corresponds to the
minimum in free energy over up to $10^{2}$ samples, resulting in optimal
parameter sets $(\boldsymbol{\theta}^{*},\boldsymbol{\phi}^{*})$. The quality
of the approximation by the variational ansatz decreases slightly in the
crossover region until, with increasing $T/m$, the entropic term in Eq.(2)
becomes dominant, yielding a high-fidelity solution. On the right, the chiral
condensate $\langle\bar{\psi}\psi\rangle$ evaluated using
$\rho_{2}(\boldsymbol{\theta}^{*},\boldsymbol{\phi}^{*})$ is shown as a
function of $T/m$. The displayed errors were obtained from the variation of
the observable over solutions corresponding to the $20$th percentile of the
free energy distribution for each set of samples.

@C=1.0em @R=0.2em @!R
q_0 : & q_0 : R_Z ($\Delta$) [0em]2
q_1 : q_1 : S^† H -1 -1 H S
q_2 : q_2 : -1 -1

Figure 5:  The decomposed three-qubit gate corresponding to the three boson
coupling in Eq.(11).

###  III.2 Resource efficient formulation

Turning now to the Hamiltonian Eqs. (11) with matter degrees of freedom
eliminated and gauge invariance _built-in_ , the choice of variational ansatz
circuits for the $\mathrm{VQC}_{i}$ is in principle unrestricted, leaving us
with the freedom to use, e.g. hardware efficient circuits such as EfficientSU2
from qiskit’s standard circuit library [21]. As for the set of physical gates,
by inspecting the Hamiltonian one sees that the single and multi-qubit terms
can be employed as elementary ansätze. As an exaample, the decomposed three-
qubit term is depicted in Fig. 5. An alternation of a single-qubit sublayer
containing $R_{x}$ and $R_{y}$ gates, followed by a staggered multi-qubit
sublayer, will hence serve as our variational ansatz for both
$\mathrm{VQC}_{1/2}$.

####  III.2.1 Observables at $T,\mu>0$

As was done in the gauge-redundant formulation, we quantitatively investigate
the quality of our variational approximation to the Gibbs state and calculate
the chiral condensate. This has been done as a function of temperature $T/m$
and varying chemical potential $\mu/m$ at $\epsilon/m=0.5$ with
$n_{s}=10^{4}$. Owing to the resource efficiency of the formulation, $n_{q}=7$
classically simulated qubits now correspond to twice the system size ($N=8$)
when compared to the gauge-redundant formulation. While simulations for
varying $T/m$ across the cross-over region, Fig. 6, show good agreement
between exact and simulated results, entering the transition region for fixed
temperature and increasing $\mu/m$, we see that our algorithm has difficulty
to match the exact result. One notes that as we lower the temperature, the
transition becomes more and more discontinuous. In Fig. 8 and Fig 7, this
observation is confirmed by our estimation of the statistical error which
measures the variation of $\langle\bar{\psi}\psi\rangle$ over the top $20$th
percentile of variational solutions. We expect an increase in performance can
be achieved here by changing from the gradient-free optimization algorithm
(COBYLA) used in this study to gradient-based methods. This is an approach
which we have not pursued further.

Figure 6: _(Left)_ The fidelity of the approximation to the Gibbs state
produced by VQT on the classical simulator and $\Delta f$ as a function of
$T/m$ for $\mu=0$ on a system with $N=8$ matter sites, corresponding to
$n_{q}=7$ qubits in the resource-efficient formulation. _(Right)_ The chiral
condensate $\langle\bar{\psi}\psi\rangle$ calculated on the same data sets.

Figure 7: _(Left)_ The fidelity of the approximation to the Gibbs state
produced by VQT on the classical simulator and $\Delta f$ as a function of
$\mu/m$ for $T/m=0.8$ on a system with $N=8$ matter sites, corresponding to
$n_{q}=7$ qubits in the resource-efficient formulation._(Right)_ The chiral
condensate $\langle\bar{\psi}\psi\rangle$ calculated on the same data sets.

Figure 8: _(Left)_ The fidelity of the approximation to the Gibbs state
produced by VQT on the classical simulator and $\Delta f$ as a function of
$\mu/m$ for $T/m=0.4$ on a system with $N=8$ matter sites, corresponding to
$n_{q}=7$ qubits in the resource-efficient formulation. _(Right)_ The chiral
condensate $\langle\bar{\psi}\psi\rangle$ calculated on the same data sets.

####  III.2.2 Thermal Unequal-Time Correlators

In addition to the study of quantum systems at finite fermion density, quantum
simulation in the NISQ era and beyond will allow the investigation of real-
time observables. These include thermal unequal-time correlation functions. In
general, an unequal-time correlator of two operators takes the form

$\displaystyle
C_{AB}(t)\equiv\langle\,e^{iHt}Ae^{-iHt}B\,\rangle=\langle\,A(t)B\,\rangle,$
(16)

where $H$ is the Hamiltonian of the system and $A(t)$ is the operator in the
Heisenberg picture. In this study, the expectation value in (16) will be taken
with respect to our variational density matrix,
$\rho_{2}(\boldsymbol{\theta}^{*},\boldsymbol{\phi}^{*})$.

The question of how one measures such a correlation function naturally arises
when studying lattice gauge theories. Fortunately, a general procedure has
already been devised to measure the correlation function of two arbitrary
unitary operators $A$ and $B$ on a universal quantum computer [22, 23]. This
procedure is commonly referred to as Ramsey interferometry and only involves
adding a single ancillary qubit entangled with our system’s original $n_{q}$
qubits. The time evolution can be performed with a Trotterization of the time
evolution operator whose circuit depth increases linearly with the number of
steps. A variational approach also exists to obtain the time evolution of an
arbitrary quantum state with a fixed circuit depth [15]. In this work, a
simple Trotterization has been used.

In Fig. 9 we show the unequal-time thermal density-density correlation
function $\langle O_{0}(0)O_{i}(t)\rangle$ at $T/m=1$ for $N=8$, where
$O_{i}\equiv\varepsilon_{i}Z_{i-1}\otimes Z_{i}$. The different data sets
represent different spatial separations $i=1,\ldots,5$ of the density operator
on the lattice plotted as a function of time with a remarkable agreement
between simulated and exact results for the time range displayed. In these
classical simulations of our variational quantum circuit, both
$\mathrm{VQC}_{1}$ and $\mathrm{VQC}_{2}$ contain two layers of alternating
single and multi-qubit sublayers. For simplicity’s sake, we have assumed
infinite statistics in the computation of
$\rho_{2}(\boldsymbol{\theta}^{*},\boldsymbol{\phi}^{*})$, hence excluding
shot noise. One is ultimately interested in obtaining transport coefficients
from real-time dynamics. As the latter involves time integrals of unequal-time
correlators, it remains to be seen if such simulations yield sufficient
accuracy.

Figure 9:  (Left) The unequal-time thermal correlation function $\langle
O_{0}(0)O_{i}(t)\rangle$ for varying distance $i=1,\ldots,5$ (shifted in
magnitude for displaying purposes) at $T/m=1$ for $N=8$, where
$O_{i}=\varepsilon_{i}Z_{i-1}\otimes Z_{i}$, i.e. the relevant term of the
number operator. The solid lines are the exact solutions while the crosses
represent the values calculated with the Trotterized time-evolution operator
($a\delta_{t}=0.2$) for the variatonal Gibbs state $\rho_{var}$ produced with
infinite statistics. (Right) The absolute value of the meson string operator
as a function of the string length (lattice units). The results for two
temperatures are plotted using both the variational Gibbs state and the exact
density matrix. One notices that at larger temperatures the current
variational ansätze give states which deviate from the exact result.

####  III.2.3 $\bar{q}q$ meson screening

It is useful to check that our algorithm correctly reproduces physical
observables, which are more easily accessible by standard Eulcidean methods.
As an example, we consider the thermal expectation value of a gauge-invariant
$\bar{q}q$ meson operator, i.e. a quark anti-quark pair connected by a flux
tube [14, 24],

$\displaystyle
M_{i}\equiv\psi^{\dagger}_{0}\left(\prod^{i-1}_{k=0}X_{k}\right)\psi_{i}.$
(17)

We have written the above operator in the gauge-redundant version for
simplicity but have used the corresponding resource-efficient expression in
our calculations. As demonstrated recently with the help of DMRG [25], at
$T=0$, this equal-time two-point function shows an exponential decay for any
$\epsilon>0$. The exponents of a spatial decay are screening masses, which in
this case give the energy of the flux tube and its excitations, whose linear
increase with separation indicates confinement. We expect this behaviour to
persist for small temperatures until the onset of string breaking, when the
flux gets increasingly screened and its energy no longer grows with
separation.

Using the same variational setup as described in Sect.III.2.2, Fig. 9 (right)
shows $|\langle M_{i}\rangle|$ at $T/m=0.5$ and $T/m=1.25$. Whereas the low-
temperature correlator exhibits marked similarities to the zero-temperature
decay (c.f. [25]), which is accurately captured by the variational Gibbs
state, at higher temperatures, the algorithm fails to reproduce the exact
correlator accurately, possibly due to the low complexity of the ansatz
circuits in our simulations (two layers of variational gates for the
$\mathrm{VQC}_{i}$). While this precludes a detailed picture of the long
distance behavior with increasing temperature, our variational calculation
does exhibit two qualitative features known from other techniques: i) the
screening masses generally grow with temperature as expected (faster decay);
ii) the difference between high and low temperature is small at short
distances, but pronounced at large distances, because of the flux screening.

####  III.2.4 Resources and Entropy Estimation

We now turn to the question of resources in runtime and classical memory.
Considering the limited system sizes currently studied in the NISQ era,
estimating the von Neumann entropy $S=-\sum_{i}p_{i}\log{p_{i}}$ from the
intermediate measurements via $p_{i}\approx n_{i}/n_{s}$ is a viable approach.
Applying the algorithm with the information flow depicted in Fig. 1 for
realistic system sizes will necessarily involve a certain level of
approximation. This stems from the fact that the number of $p_{i}$’s grows
with the dimension of the Hilbert space, $2^{n_{q}}$. One possible way
proposed in [11], is the partitioning of the full $n_{q}$ qubit system into
$n_{q}/n_{ss}$ independent subsystems of size $n_{ss}$ in the entropic part of
the algorithm ($\mathrm{VQC}_{1}$) that prepares the latent distribution
$\rho_{1}(\boldsymbol{\theta})$. In the version of the algorithm which employs
a classical estimation of the entropy [7] this corresponds to a particular
choice of the latent distribution known as factorized latent space models. It
remains to be seen if such an approximation is a valid approach for gauge
theories where thermalization and subsystem entanglement are non-trivially
connected [26].

There are, however, alternative methods that can be used to estimate the
entropy efficiently. For example, it turns out that an evaluation of the
Taylor series approximation to the entropy can be useful [27]. This method
requires the addition of $n_{q}$ ancillary bits to perform a qubit-efficient,
repeated SWAP-test evaluation of the trace of various powers of the latent
density matrix, $\mathrm{tr}(\rho^{n}_{1}(\theta))$. The results of this
approach for the simple gauge theory used in this study are displayed in Fig.
10. For temperatures up to $T/m\sim 1$, only a few orders are needed to give
reasonable values of the fidelity.

Figure 10: The fidelity of the Gibbs state produced using the Taylor-series
approximated von Neumann entropy as a function of $T/m$ at $a\mu=0$ for
various truncation orders $K$ on a system with $N=8$ matter sites. One sees
that for low temperatures this approximation works well.

##  IV Conclusion

We have explored the use of variational quantum algorithms in producing
thermal states. These approaches were applied to $\mathbb{Z}_{2}$ lattice
gauge theory in $1+1$ dimensions using two independent formulations, and
several observables of interest were computed. We find this variant of the VQT
a suitable algorithm to approximate the thermal states of quantum systems in
the NISQ era as estimating the entropy of larger systems will necessitate the
introduction of truncations in the estimation. A particular aspect we have
left for future work is the performance of intermediate measurements for
mixed-state creation on actual quantum devices compared to the usual approach
of entangling ancillary systems followed by a terminal measurement.
A further possible extension of our work is the formulation of the VQT to
study non-Abelian gauge theories at finite temperatures. As the latter will
generally come with a gauge-redundant Hilbert space, even after reformulating
in a resource-efficient manner, this amounts to respecting gauge invariance or
equivalent local constraints throughout the variational quantum circuit.

###### Acknowledgements.

The work is supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through the grant CRC-TR 211 “Strong-interaction matter
under extreme conditions” – project number 315477589 – TRR 211 and by the
State of Hesse within the Research Cluster ELEMENTS (Project ID 500/10.006).
M.S. and M.F. acknowledge support by the Munich Institute for Astro-, Particle
and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy –
EXC-2094 – 390783311.

## References

  * Bañuls  _et al._ [2020] M. C. Bañuls _et al._ , Simulating Lattice Gauge Theories within Quantum Technologies, Eur. Phys. J. D 74, 165 (2020), arXiv:1911.00003 [quant-ph] .
  * Bauer  _et al._ [2023] C. W. Bauer _et al._ , Quantum Simulation for High-Energy Physics, PRX Quantum 4, 027001 (2023), arXiv:2204.03381 [quant-ph] .
  * Wu and Hsieh [2018] J. Wu and T. H. Hsieh, Variational thermal quantum simulation via thermofield double states, Phys. Rev. Lett. 123, 220502 (2019) 123, 220502 (2018), arXiv:1811.11756 [cond-mat.str-el] .
  * Ball and Cohen [2022] C. Ball and T. D. Cohen, Boltzmann distributions on a quantum computer via active cooling 10.48550/ARXIV.2212.06730 (2022), arXiv:2212.06730 [quant-ph] .
  * Powers  _et al._ [2021] C. Powers, L. B. Oftelie, D. Camps, and W. A. de Jong, Exploring finite temperature properties of materials with quantum computers, (2021), arXiv:2109.01619 [quant-ph] .
  * Davoudi  _et al._ [2022] Z. Davoudi, N. Mueller, and C. Powers, Toward quantum computing phase diagrams of gauge theories with thermal pure quantum states, (2022), arXiv:2208.13112 [hep-lat] .
  * Verdon  _et al._ [2019] G. Verdon, J. Marks, S. Nanda, S. Leichenauer, and J. Hidary, Quantum hamiltonian-based models and the variational quantum thermalizer algorithm, (2019), arXiv:1910.02071 [quant-ph] .
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5, 10.1038/ncomms5213 (2014).
  * McClean  _et al._ [2015] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New J. Phys. 18 (2016) 023023 10.1088/1367-2630/18/2/023023 (2015), arXiv:1509.04279 [quant-ph] .
  * Foldager  _et al._ [2021] J. Foldager, A. Pesah, and L. K. Hansen, Noise-assisted variational quantum thermalization, (2021), arXiv:2111.03935 [quant-ph] .
  * Selisko  _et al._ [2022] J. Selisko, M. Amsler, T. Hammerschmidt, R. Drautz, and T. Eckl, Extending the variational quantum eigensolver to finite temperatures, (2022), arXiv:2208.07621 [quant-ph] .
  * Consiglio  _et al._ [2023] M. Consiglio, J. Settino, A. Giordano, C. Mastroianni, F. Plastina, S. Lorenzo, S. Maniscalco, J. Goold, and T. J. G. Apollaro, Variational gibbs state preparation on nisq devices 10.48550/ARXIV.2303.11276 (2023), arXiv:2303.11276 [quant-ph] .
  * Horn  _et al._ [1979] D. Horn, M. Weinstein, and S. Yankielowicz, Hamiltonian approach to $z(n)$ lattice gauge theories, Phys. Rev. D 19, 3715 (1979).
  * Greenberg  _et al._ [2022] T. Greenberg, G. Pardo, A. Fortinsky, and E. Zohar, Resource-efficient quantum simulation of lattice gauge theories in arbitrary dimensions: Solving for gauss’ law and fermion elimination, (2022), arXiv:2206.00685 [quant-ph] .
  * Irmejs  _et al._ [2022] R. Irmejs, M. C. Banuls, and J. I. Cirac, Quantum simulation of z2 lattice gauge theory with minimal requirements, (2022), arXiv:2206.08909 [quant-ph] .
  * Zohar and Burrello [2014] E. Zohar and M. Burrello, A formulation of lattice gauge theories for quantum simulations, Phys. Rev. D 91, 054506 (2015) 10.1103/PhysRevD.91.054506 (2014), arXiv:1409.3085 [quant-ph] .
  * Zohar and Cirac [2018] E. Zohar and J. I. Cirac, Eliminating fermionic matter fields in lattice gauge theories, Phys. Rev. B 98, 075119 (2018) 10.1103/PhysRevB.98.075119 (2018), arXiv:1805.05347 [quant-ph] .
  * Mazzola  _et al._ [2021] G. Mazzola, S. V. Mathis, G. Mazzola, and I. Tavernelli, Gauge invariant quantum circuits for $u(1)$ and yang-mills lattice gauge theories, Phys. Rev. Research, 2021 10.1103/PhysRevResearch.3.043209 (2021), arXiv:2105.05870 [quant-ph] .
  * Murairi  _et al._ [2022] E. M. Murairi, M. J. Cervia, H. Kumar, P. F. Bedaque, and A. Alexandru, How many quantum gates do gauge theories require?, (2022), arXiv:2208.11789 [hep-lat] .
  * Michael A. Nielsen [2010] I. L. C. Michael A. Nielsen,  _Quantum Computation and Quantum Information_ (2010).
  * Qiskit contributors [2023] Qiskit contributors, Qiskit: An open-source framework for quantum computing (2023).
  * Ortiz  _et al._ [2000] G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Quantum algorithms for fermionic simulations, Phys.Rev.A64:022319,2001 64, 022319 (2000), arXiv:cond-mat/0012334 [cond-mat] .
  * Somma  _et al._ [2002] R. Somma, G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Simulating physical phenomena by quantum networks, Phys. Rev. A 65, 042323 (2002).
  * Zohar [2019] E. Zohar, Local manipulation and measurement of nonlocal many-body operators in lattice gauge theory quantum simulators, Phys. Rev. D 101, 034518 (2020) 101, 034518 (2019), arXiv:1911.11156 [quant-ph] .
  * Borla  _et al._ [2019] U. Borla, R. Verresen, F. Grusdt, and S. Moroz, Confined phases of one-dimensional spinless fermions coupled to $z_{2}$ gauge theory, Phys. Rev. Lett. 124, 120503 (2020) 124, 120503 (2019), arXiv:1909.07399 [cond-mat.str-el] .
  * Mueller  _et al._ [2022] N. Mueller, T. V. Zache, and R. Ott, Thermalization of gauge theories from their entanglement spectrum, Physical Review Letters 129, 011601 (2022).
  * Wang  _et al._ [2021] Y. Wang, G. Li, and X. Wang, Variational quantum gibbs state preparation with a truncated taylor series, Phys. Rev. Appl. 16, 054035 (2021).