# Solar rotation elements i, $\mathbf{\Omega}$ and period determined using
sunspot observations by Ruđer Bošković in 1777111Contributed paper (oral
presentation) at the International Advent Workshop 17th - 18th December 2020,
Zagreb - Graz, Hvar Observatory, Faculty of Geodesy, University of Zagreb,
Croatia and IGAM, Institute of Physics, Karl-Franzens University Graz,
Austria.

M. HUSAK1 R. BRAJŠA2 and D. ŠPOLJARIĆ2

1Trakošćanska 20 42000 Varaždin  Croatia
2Faculty of Geodesy University of Zagreb  Kačićeva 26
10000 Zagreb Croatia

###### Abstract

Ruđer Bošković developed methods for determination of solar rotation elements:
the solar equator inclination $i$, the longitude of the node $\Omega$ and the
period of solar rotation. In his last work Opera pertinentia ad opticam et
astronomiam, published in 1785, in the chapter Opuscule II he described his
methods, the formulae with figure descriptions and an example for calculation
of the solar rotation elements with detailed numerical explanation using his
own observations performed in September 1777. The original numerical procedure
was performed using logarithmic formulae. In present work we give a
description of the original results of R. Bošković and compare them with our
recalculated values.

###### keywords:

Ruđer Bošković - solar rotation elements - methods

##  1 Introduction

Ruđer Bošković (1711–1787) was a Croatian Jesuit priest and his main fields of
interest were physics, mathematics, geodesy, cartography, and astronomy,
including both theoretical and practical aspects. Bošković performed research
on determination of the solar rotation elements which are called Carrington’s
elements: the longitude of the ascending node $\Omega$ and the inclination of
the solar equator to the ecliptic $i$ (Figure 1).

Figure 1: The Carrington solar rotation elements: the longitude of the
ascending node, $\Omega$, and the inclination of the solar equator regarding
ecliptic, $i$ (according to Figure 4-2 in Brajša (1989)).

In his first dissertation De maculis solaribus (Boscovich, 1736) he described
his methods for solar rotation elements determination using sunspot
observations. In September 1777 he observed sunspots in Noslon near Sens 120
km south from Paris where he was a guest of the cardinal Paul d’Albert de
Luynes (1703–1788), an amateur astronomer. The detailed research report about
the solar observations was published in Opuscule II in fourteen chapters and
appendix Appendice (Boscovich, 1785, 75-178). The Opuscule II gives (in the
chapters §. I. – §. XIV.) as follows: the formulae with figure descriptions
(Préface) and observations (§. I.), the sunspot positions determination (§.
II. – §. III.), the methods (§. IV. – §. VII.), the calculation instructions
(§. IX. – §. XIII.), and additionally the reflexions (Réflexions in §. VIII.
and §. XIV.).

##  2 Methods

In September 1777 Bošković observed four sunspots using the
telescope222Télescope dans Lalande T2 - Noslon - Wikipédia (wikipedia.org).
https://fr.wikipedia.org/wiki/Noslon#/media/Fichier:Télescope_dans_Lalande_T2.png,
accessed 11. 05. 2021. (Figure 2) with optical micrometre and pendulum for
time measurement. As an example, he used observations of one sunspot. Cardinal
de Luynes had contemporary and professional instruments, which his guests used
for research.

Figure 2: Telescope that was described by Lalande in Astronomie T II.

For his methods Bošković developed step-by-step procedures in a form of twelve
tables which are convenient for calculations with logarithmic tables
(Boscovich, 1785, Tab. I. to Tab. XII., 166-169). He performed measurements
and data reduction in two steps using his own formulae (Boscovich, 1785,
Préface, №12 – №18, 81-84). First, he measured the sunspot position at the
solar disk and the true solar time, as well as the position of the sunspot in
ecliptic coordinate system, the longitude and the latitude (T.M., lon.t,
lat.B.t) (Boscovich, 1785, Tab. I. and Tab. II., 166-167), in the moment of
the sunspot observation (Boscovich, 1785, §. I., 87-89). Then he calculated
the solar rotation elements: $i$, $\Omega$ and the sidereal and the synodic
rotational period. The Bošković’s methods use the mean solar time and sunspot
position combinations (T.M., lon.t, lat.B.t) for solar rotation elements
determination: the longitude of the node $N=\Omega$ using two positions of the
same sunspot in Tab. III. and Tab. IV. and the solar equator inclination $i$
using two positions of the same sunspot and the longitude of the node
$N=\Omega$ in Tab. V. and Tab. VI. The periods of the solar rotation, the
sidereal and synodic ones, are given in Tab. VII. to Tab. XI. The
determination of the longitude of the node $N=\Omega$ and the solar equator
inclination $i$ using three positions of the same sunspot are presented in
Tab. XII. These tables are reproduced with detailed descriptions in the paper
by Husak, Brajša, and Špoljarić (2021).

##  3 Results

Bošković used observational data, astronomic almanac data and logarithmic
tables for the determination of the mean solar time and positions of a
sunspot. Observational data are given in Opuscule II and other input data for
calculation were taken from astronomic almanac. The values from the
astronomical almanac Connoissance des temps333URL 1: Connoissance des temps.
Connoissance des temps, Pour l’ Année commune 1777. publiée Par l’ordre de l’
Académie Royale des Sciences, el calculée Par M. Jeaurat, de la même Académie.
A Paris de l’imprimerie royale. M. DCCLXXVI. (Observatoire de Paris),
https://gallica.bnf.fr/ark:/12148/cb327469896/date are: the correction of the
local meridian to Paris meridian $\Delta t_{\rm Sens}=3^{m}48^{s}$ at the page
268, the correction of the true solar time to the mean solar time $\Delta
t_{\rm MST}$ at the page 103 (T.M., French Temps moyen au Midi vrai – Average
time at true noon), the longitude $\lambda_{\odot}$=$lon._{\odot}$ and the
declination $\delta_{\odot}$=dec.⊙ of the Sun in Paris interpolated to the
true solar time in minutes for a day of observation at the pages 102-103, the
apparent diameter of the Sun $D_{\odot}$ interpolated to the day of
observation at the page 108 and the inclination of the ecliptic rounded to
whole minutes $\varepsilon=23^{\circ}28^{\prime}$ at the page 4. Ruđer
Bošković determined constant of the micrometre C=1237/1915 (relative units of
micrometre per arc sec) experimentally by multiple observations of the
vertical size of the solar disc on 11th September 1777 using interpolated
solar diameter $D_{\odot}$. In present work we recalculated results using
original Bošković’s and corrected mean solar time and positions (T.M., lon.t,
lat.B.t). In that way we can compare the results of the example of Bošković
and the present work results. The input data for the solar rotation elements
are data for the sunspot (T.M., lon.t, lat.B.t) and, depending on the method,
the solar rotation elements determined earlier ($\Omega=N$, $i$, $T$).
Bošković determined the longitude of the node $\Omega=N$ and the solar equator
inclination $i$ combining different pairs of observational data. So, he
obtained the following values: as the arithmetic mean of six pairs
$\Omega_{6}=70^{\circ}21^{\prime}$, of eight pairs
$\Omega_{8}=71^{\circ}32^{\prime}$ and of ten pairs
$\Omega_{10}=73^{\circ}09^{\prime}$ of the positions of the same sunspot which
Bošković used in the Tab. IV. Further, the solar equator inclination
$i=7^{\circ}44^{\prime}$ was calculated as the arithmetic mean of five pairs
of the positions of the same sunspot which Bošković used in the Tab. VI. The
sidereal period $T^{\prime}=26.77$ days was determined as the arithmetic mean
of six pairs of times and sunspot positions of the same sunspot which Bošković
used in the Tab. X. and synodic period $T^{\prime\prime}=28.89$ days in the
Tab. XI. The longitude of the node $N=\Omega_{136}=74^{\circ}03^{\prime}$ and
the solar equator inclination $i=6^{\circ}49^{\prime}$ using three positions
(1, 3, 6) of the same sunspot are given in the Tab. XII. We recalculated the
results using the methodology of Bošković as follows. The longitude of the
node $\Omega=N$ was calculated as the arithmetic mean of six pairs
$\Omega_{6}=70^{\circ}21^{\prime}$, eight pairs
$\Omega_{8}=71^{\circ}32^{\prime}$ and ten pairs
$\Omega_{10}=73^{\circ}12^{\prime}$ of the positions of the same sunspot which
Bošković used in the Tab. IV. The solar equator inclination
$i=7^{\circ}45^{\prime}$ was calculated as the arithmetic mean of five pairs
of the positions of the same sunspot which Bošković used in the Tab. VI. The
sidereal period $T^{\prime}=26.76$ days was calculated as the arithmetic mean
of six pairs of corrected times and sunspot positions of the same sunspot
which Bošković used in the Tab. X. The corresponding synodic period is
$T^{\prime\prime}=28.87$ days and with Bošković’s times $T^{\prime}=26.77$
days and $T^{\prime\prime}=28.89$ days. The longitude of the node
$N=\Omega_{136}=74^{\circ}03^{\prime}$ and the solar equator inclination
$i=6^{\circ}48^{\prime}$ using three positions (1, 3, 6) of the same sunspot
are given in the Tab. XII.

##  4 Analysis and discussion

In modern research we do not use logarithmic tables and step-by-step
calculations which Bošković used in his example (Boscovich, 1785, Tab. I. to
Tab. XII., 166-169). Modern computers enable fast complex calculations in
double precision. This research is performed in spreadsheet (Microsoft
Excel®). It could be performed at any other computer platform which supports
elemental mathematical functions. The main problem in reproduction of
Bošković’s example was discovering calculation chains at the places where
formulae are missing. The formulae were successfully reproduced using
Bošković’s results and descriptions of calculations in Opuscule II. We
determined the mean solar time T.M. using correction for the mean solar time
which resulted with T.M. values different from these which Bošković published.
We assume that Bošković used a wrong table from the astronomical almanac, but
we do not have the final confirmation for this hypothesis yet.

##  5 Conclusion

In present work we successfully reproduced the sunspot positions and the solar
rotation elements using the original formulae and original observational data
of R. Bošković. The reproduced values of the mean solar time and the positions
(T.M., lon.t, lat.B.t) are almost the same. We then repeated determination of
the solar rotation elements with Bošković’s original coordinates and time and
additionally periods with corrected time. The reproduction of the results uses
the methodology which Bošković used in his example with the same input data
for the time and sunspot positions. The present work results confirm that
Bošković did extraordinary work, from the beginning (his own methods, with his
own observations and calculations) to the final results. The example confirms
and verifies his methods for solar rotation elements using his own sunspot
observations performed in the late 18th century with very similar results
compared with present work results. The most important achievement of the
present work is a confirmation of the methodology and reproduction of the
results with original formulae for the methods for the solar rotation elements
determination using sunspot observations of Ruđer Bošković. In September 1777
Bošković observed four sunspots which we can now further reduce and apply
modern formulae convenient for performing the methods of Ruđer Bošković. This
research is in progress and in the next step we will analyse the precision of
input parameters (solar ephemeris, apparent solar diameter, etc.) on the
precision of output data in solar rotation determination. Further, we will
develop formulae applicable for modern computer programs and try to fix
numerical calculation problems and difficulties. These results will enable
application of the methods for solar rotation determination to modern data and
will give a possibility to compare the results of R. Bošković with modern
results.

## Acknowledgements

This work has been supported by the Croatian Science Foundation under the
project 7549 "Millimeter and submillimeter observations of the solar
chromosphere with ALMA".

## References

  * Boscovich, R. J.: 1736, A. M. G. D. De Maculis Solaribusn exercitatio astronomica habita (About sunspots), Colegio Romano Societatis Jesu, Typographia Komarek, Romae, 3-10, (in Latin).

  * Boscovich, R. J.: 1785, Opera pertinentia as opticam, et astronomiam (… optics and astronomy), Maxima ex parte nova, & omnia hucusque inedita, In quinque Tomos distributa Ludovico XVI. Galliarum Regi Potentissimo Dicata. Tomus quintus, Opuscule II, Bassani, 75-178. (in French).

  * Brajša, R.: 1989, Determination of the solar differential rotation at high latitudes using polar crown filaments, Diploma work, University of Zagreb, Faculty of Science and Mathematics, Zagreb, (in Croatian).

  * Husak, M., Brajša, R., Špoljarić, D.: 2021, RGN Zbornik 36 No. 3, in press.


# Orbitally driven spin reorientation in Mn doped YBaCuFeO5

Mukesh Sharma  T. Maitra<EMAIL_ADDRESS>Department of Physics,
Indian Institute of Technology Roorkee, Roorkee-247667, Uttarakhand,India

###### Abstract

Oxygen-deficient layered perovskite YBaCuFeO5 (YBCFO) is one rare type-II
multiferroic material where ferroelectricity, driven by incommensurate spiral
magnetic order, is believed to be achievable up to temperatures higher than
room temperature. A cycloidal spiral rather than helical spiral order is
essential ingredient for the existence of ferroelectricity in this material.
Motivated by a recent experimental work on Mn-doped YBCFO where the spiral
plane is observed to cant more towards the crystallographic $c$ axis upon Mn
doping at Fe sites compared to that in the parent compound, we performed a
detailed theoretical investigation using the density functional theory
calculations to understand the mechanism behind such spin reorientation. Our
total energy calculations, within GGA+U+SO approximation, reveal that Fe/Cu
spin moments indeed align more towards $c$ axis in the Mn-doped compounds than
they were in parent compound YBCFO. Largest exchange interaction (Cu-Cu in the
$ab$-plane) is observed to decrease systematically with Mn doping
concentration reflecting the lowering of transition temperature seen in
experiment. Further, the inter-bilayer Cu-Mn exchange interaction becomes
ferromagnetic in the doped compound whereas the corresponding Cu-Fe exchange
was antiferromagnetic in the parent compound giving rise to frustration in the
commensurate magnetic order. Most importantly, our electronic structure
calculations reveal that in the doped compounds because of hybridization with
the d${}_{z^{2}}$ orbital of Mn, the highest occupied Cu orbital becomes
d${}_{z^{2}}$ as opposed to the parent compound where the highest occupied Cu
orbital is d${}_{x^{2}-y^{2}}$. Therefore, we believe that the occupancy of
out of plane oriented Cu d${}_{z^{2}}$ orbital in place of planar
d${}_{x^{2}-y^{2}}$ orbital along with the frustrating exchange interaction
drive the spins to align along $c$ in doped compound.

##  I Introduction

Multiferroic materials with high transition temperature as well as strong
coupling between magnetization and polarization are highly sought after
because of their tremendous potential in device applications such as in
spintronics, storage devices, sensors etc.Yoshinori_Tokura2006 ;
W_Eerenstein2006 . Type-I multiferroics, where magnetization and polarization
have different origins, are found to have high transition temperatures but due
to weak coupling between magnetization and polarization these materials are
less useful. In contrast, the type-II multiferroics where the ferroelectricity
is driven by the magnetic ordering, the coupling between magnetization and
polarization is strong and hence are more useful for device
applicationsBuurma2016 ; Khomskii2006 . Unfortunately, most type-II
multiferroics have much lower transition temperatures compared to room
temperatureKhomskii2009 .

Cupric oxide (CuO), and bilayer perovskite, YBaCuFeO5 (YBCFO), are the two so
far discovered type-II multiferroics, which show ferroelectric ordering at
higher temperatures ($>$ 200 K). In CuO, a spiral magnetic order drives the
ferroelectricity which survives up to higher temperatures but in a narrow
temperature window (213 K - 230 K)T_kimura2008 . YBCFOL_Er-Rakho1988 ;
B_Kundys2009 ; Yuji_Kawamura2010 , on the other hand, has spiral
(incommensurate) magnetic ordering up to a much higher temperature and in a
wide range (T $<$ 230 K). YBCFO displays two magnetic transitions,
paramagnetic to commensurate (CM) antiferromagnetic (T${}_{N1}\sim$ 440 K) and
CM to incommensurate (ICM) (T${}_{N2}\sim$ 230 K) spiral magnetic ordering,
respectivelyB_Kundys2009 . This ICM or spiral magnetic order is believed to be
accompanied by ferroelectricity as some earlier experimental work
reportedB_Kundys2009 ; Yuji_Kawamura2010 ; morin2015 .

Of late, there has been a debate on the nature of spiral magnetic state
(whether helical or cycloidal) and existence of ferroelectricity in this
materialLai2017 ; Dey2018 . Cycloidal spiral state (where the spin moments are
directed out of $ab$-plane) is necessary to have finite polarization in this
systemDey2018 . In non geometrically frustrated system which is the case in
YBCFO, competing nearest (JNN) and next nearest neighbour (JNNN) exchange
interaction along $c$ direction can stabilize a cycloidal spiral state.
However, it has been found from density functional theory (DFT) calculations
that neither (JNNN) nor the Dzyaloshinskii-Moriya interaction (DMI) is strong
enough to give rise to such a state morin2015 ; Dey2018 . An alternative
mechanism involving Fe/Cu chemical disorder giving rise to randomly
distributed Fe-Fe antiferromagnetic bonds was subsequently proposed for
causing the desired frustration in the commensurate magnetic structure and
stabilizing the spiral order therebymorin2015 ; Scaramucci2018 . Later on,
experimentally it was shown that by introducing Fe/Cu chemical disorder during
the preparation of sample indeed gives rise to stable spiral order up to a
much higher temperaturemorin2016 . A combination of both chemical disorder and
substitution at the Ba sites by other rare-earth ions has been shown recently
to raise the spiral ordering temperature to way beyond room
temperatureshang2018 ; Lal2017 . However, these authors did not report any
observation of ferroelectricity in their samples.

Figure 1: A $2\times 2\times 2$ YBCFO supercell is shown with (a) ABAB and (b)
AABB Fe/Cu arrangements respectively. (c) One unit cell of YBCFO is shown
where inter-bipyramidal distance (d1) and bipyramidal thickness (d2) are
marked.

Previous experimental reportsRuiz1998 ; Mombru1998 ; Caignaert1995 ; morin2015
; Lai2017  suggest that YBCFO crystallizes in to two possible space groups
P4/mmm and P4mm depending on Fe/Cu ion distribution where the former is
centrosymmetric and the latter is non-centrosymmetric. In Fig. 1 we present
the crystal structure of YBCFO where bipyramidal layers of MO5 (M = Fe, Cu)
separated by Y+3 layer can be seen. Due to very similar sizes of Fe3+ and Cu2+
ions both can go to the same site. However, it has been found that the
experimentally observed ferromagnetic (FM) interaction within bipyramidal
layers along $c$ direction is only possible if they are occupied
preferentially by the Fe-Cu pairsmorin2015 . Keeping this in view and
considering experimentally observed antiferromagnetic (AFM) interaction in the
$ab$-plane, we have considered two possible Fe/Cu arrangements as shown in
Fig. 1: (a) ABAB arrangement where FeO5-CuO5 bipyramidal layers are stacked
along $c$ direction in such a way that Fe and Cu layers alternate along $c$
and (b) AABB arrangement where FeO5-CuO5 bipyramidal layers are stacked along
$c$ direction in such a way that two layers of Fe and two layers of Cu
alternate along $c$. Fig. 1(c) depicts one unit cell of YBCFO, the bipyramidal
thickness (d2) and inter-bipyramidal separation (d1).

In a very recent experimental work, Zhang et al.zhang2021  have looked into
the effect of doping of Mn at Fe sites, (i.e. YBaCuFe1-xMnxO5 (YBCFMO) with
x=0, 0.01, 0.05, 0.1, 0.15 and 0.2) on the nature the spiral state in this
system. The authors have observed that on increasing Mn doping the spiral
plane tilts more towards $c$-axis from $ab$-plane even though the transition
temperatures of spiral state as well as commensurate magnetic state (TN2 and
TN1 respectively) decrease with doping. The tilting of spiral plane towards
$c$-axis is a desirable property for having finite electric polarization in
the systemKatsura2005 . In this work we investigate the underlying mechanism
behind such a change in spin moments’ orientation through a detail electronic
structure calculation using density functional theory.

##  II Methods

We used density functional theory (DFT), a first principles approach, to
investigate the magnetic order, electronic structure, orbital state in Mn
doped YBCFO. Vienna Abinitio Simulation Package (VASP) is used for DFT
calculations which uses pseudopotentials and plane-wave basis setsvasp ; vasp1
. We considered Perdew-Burke-Ernzerof generalized gradient approximation (PBE-
GGA)gga  as exchange correlation functional for our calculations. We set the
plane wave cut off to 600 eV and sampled the Brillouin zone by $8\times
8\times 4$ Monkhorst-pack k-grid points. YBCFO and YBCFMO are strongly
correlated systems and so the effect of Coulomb correlation (U) for localized
d-states cannot be ignored. We incorporated the on-site Coulomb correlation
(U) and Hund’s coupling strength (J) within the GGA+U approximation. For this
purpose, we have set U (J) values at 5 (1), 5 (1) and 8 (0) eV for d-states of
Fe, Mn and Cu respectively as used in previous literaturemorin2015 . We have
considered YBCFO supercell as shown in Fig.1(a-b) which is the parent compound
where we substituted Fe by Mn to create the doped compounds YBCFMO with x =
0.125, 0.25. We used experimentalzhang2021  lattice parameters for YBCFO (x =
0) and to find the lattice parameters for YBCFMO with x=0.125 and 0.25 we
interpolated and extrapolated the experimental lattice parameters available
for x=0.01, 0.05, 0.1, 0.15 and 0.2zhang2021 . For x=0.25 we have also
performed c/a optimization. These lattice parameters for x=0.0, 0.125 and 0.25
are listed in Table 1. We then performed full ionic relaxation until the
Hellman-Feynman force reached 0.01 $eV/\AA{}$ for parent as well as doped
compounds. Non-collinear calculations were performed within GGA+U+SO
approximation where SO is the spin-orbit (SO)so  coupling. Magnetic exchange
interactions are calculated using the method prescribed by Xiang et al.xiang .

##  III Results and Discussion

###  III.1 Crystal Structure Analysis

We have investigated in detail the effect of Mn doping on the YBCFO crystal
structure using density functional theory calculations as discussed below. As
it is observed experimentally that Mn ions have strong preference for Fe
siteszhang2021 , we have performed calculations considering Mn ions substitute
Fe rather than Cu in YBCFO. In Table 1 we list the lattice parameters and the
total energies of ABAB and AABB arrangements (see Fig.1) for undoped and Mn-
doped cases considered in our calculations. Looking at Table 1 we observe that
the lattice constants $a$ and $b$ decrease slowly but monotonically with Mn
doping whereas the lattice constant $c$ increases with doping exactly similar
to what is observed in the experimentzhang2021 . We further observe from Table
1 that ABAB arrangement has lower energy than AABB arrangement both in undoped
and doped cases considered which is consistent with the P4mm symmetry observed
in experiment. Therefore, in the discussion below we will mostly focus on ABAB
structure.

Table 1: Lattice parameters (in $\AA$) for $2\times 2\times 2$ supercell for YBCFO and YBCFMO are listed for different $x$ values. $\Delta$E($x$) = Ex(AABB) - Ex(ABAB) is the energy difference, in $meV$, per formula unit between the two Fe/Cu ordering considered within GGA+U+SO approximation. x | 2a | 2b | 2c |  $\Delta$E($x$)
---|---|---|---|---
0.0 | 7.7493 | 7.7493 | 15.3253 | 72.01
0.125 | 7.7447 | 7.7447 | 15.3386 | 4.06
0.25 | 7.7405 | 7.7405 | 15.3503 |
Figure 2: (a)-(c) show the calculated average basal and apical M-O bond
lengths for M=Fe, Cu and Mn respectively as a function of doping concentration
(x) in YBCFMO. (d), (e), and (f) represent the inter-bilayer separation (d1),
bilayer thickness (d2), and their ratios (d1/d2), respectively as a function
of doping.

To explore how local Jahn-Teller splitting is affected by doping we look at
the average M-O bonds (M = Fe/Mn/Cu) in the $ab$-plane (basal) and along
$c$-direction (apical). We present in Fig.2 (a)-(c) these bond lengths as a
function of doping. We would like to note here that for x = 0.25 we have
considered four different Mn distributions as shown in Fig. 3. The data used
in Fig.2 for x = 0.25 are extracted from the lowest energy configuration out
of the four. From Fig.2 we observe that in case of Fe-O, the basal bond is
longer than the corresponding apical bond whereas it is the opposite in Cu-O
case. This is consistent with the experimental observations by Zhang et
al.zhang2021 . In case of Mn-O, the apical bond is longer than the basal ones.
Thus the doping of Mn at Fe site would decrease the average basal M-O bond
length and increase the average apical M-O bond length. The difference between
apical and basal M-O bond lengths is maximum for Cu-O and minimum for Mn-O.
Hence local Jahn-Teller effect is more in CuO5 than in MnO5. In Fig.2 (d)-(f)
we present the variation of inter-bilayer separation (d1) (see Fig.1(c)) ,
bilayer thickness (d2) and their ratio (d1/d2) respectively. Comparing with
experimentzhang2021 , we find the behaviour of these three parameters are
quite similar to that observe in experiment up to x = 0.125. However, their
behaviour differs from experiment when we go from x = 0.125 to 0.25.

Figure 3: Ground state spin orientation in YBCFMO calculated within GGA+U+SO
for (a) x=0, (b) x=0.125, (c) x=0.25 (confg-I), (d) x=0.25 (confg-II), (e)
x=0.25 (confg-III) and (f) x=0.25 (confg-IV).

###  III.2 Magnetic state and exchange interaction

In the experimental measurements by Zhang et al.zhang2021 , it has been
observed that Mn doping drives the system more towards a cycloidal spiral
state where moments tend to move away from $ab$-plane unlike the helical
spiral state of the parent compound observed previously by Lai et al.Lai2017 .
To explore the effect of Mn doping on the spin orientation, we performed total
energy calculations for various non-collinear spin orientations such as [100],
[110],[001] and [111] within GGA+U+SO approximation considering the
commensurate magnetic state. The total energy comparison is presented in Table
2. Here for x = 0.25 case we have considered four different Mn distributions
(see Fig. 3, denoted as confg-I, confg-II, confg-III and confg-IV).

Table 2: Total energy comparison among various orientations of spin moments for different Mn doping values calculated within GGA+U+SO approximation. | Orientation of spin moments
---|---
Doping (x) | [100] | [110] | [001] | [111]
0.0 | 0.44 | 0.47 | 1.11 | 0.0
0.125 | 1.84 | 1.94 | 0.0 | 1.42
0.25 (confg-I) | 3.46 | 3.69 | 0.0 | 2.59
0.25 (confg-III) | 0.0 | 0.05 | 0.55 | 0.50

The ground state spin orientation is shown for various undoped and doped cases
considered in our calculations in Fig. 3. We observe that in the parent
compound the spins prefer to lie in the $ab$-plane if we consider AABB
stacking as reported previouslyDey2018 . However, if we consider ABAB stacking
spins prefer to lie along [111] direction (see Fig 3(a)) which is consistent
with data reported by Zhang et al.zhang2021 . With Mn doping of 12.5$\%$ (x =
0.125), the spin moments are seen to reorient towards $c$-axis (see Fig.
3(b)). For x = 0.25, the ground state spin orientations for four different Mn
configurations (confg-I to confg-IV) are shown in Fig. 3(c)-(f) respectively.
We observe that except confg-III (i.e. Fig. 3(e)) which is not the lowest
energy state, all the remaining three configurations have their spin easy axis
oriented towards $c$-axis. In confg-III the spin easy axis lies in the
$ab$-plane. Amongst these four configurations, confg-I (Fig. 3(c)) has the
lowest energy. Therefore, we also conclude from our calculations that with Mn
doping the spin easy axis is driven towards $c$ axis as observed in the
experiment. To understand the spin reorientation with Mn doping, we evaluated
various exchange interactions as listed in Table III using the method
prescribed by Xiang et al.xiang .

Table 3: Calculated nearest neighbour magnetic exchange interaction (J), in meV, within GGA+U. Here Jab, J${}_{c}^{BP}$, and J${}_{c}^{IP}$ denote the interaction between magnetic moments within ab-plane, along $c$ within bi-pyramids (BP), and along $c$ inter-bilayer (IP) respectively. Positive sign indicate antiferromagnetic interaction and negative sign indicate ferromagnetic interaction. J |  x$\xrightarrow{}$ 0.0 | 0.125 | 0.25 (confg-I) | 0.25 (confg-III)
---|---|---|---|---
$J_{ab}$ |  $J_{Cu-Cu}$ = 200.09 |  $J_{Cu-Cu}$ = 183.08 |  $J_{Cu-Cu}$ = 165.26 |  $J_{Cu-Cu}$ = 151.65
|  $J_{Fe-Fe}$ = 14.54 |  $J_{Fe-Fe}$ = 12.92 |  $J_{Fe-Fe}$ = 15.03 |  $J_{Fe-Fe}$ = 7.19
|  |  $J_{Fe-Mn}$ = 4.18 |  $J_{Mn-Mn}$ = -0.15 |  $J_{Fe-Mn}$ = 3.71
|  |  |  $J_{Mn-Fe}$ = 4.33 |
$J_{c}^{BP}$ |  $J_{Fe-Cu}$ = -1.79 |  $J_{Fe-Cu}$ = -1.94 |  $J_{Fe-Cu}$ = -2.11 |  $J_{Fe-Cu}$ = 5.39
|  |  $J_{Cu-Mn}$ = -1.17 |  $J_{Cu-Mn}$ = -0.89 |  $J_{Cu-Mn}$ = 2.68
$J_{c}^{IP}$ |  $J_{Fe-Cu}$ = 1.29 |  $J_{Fe-Cu}$ = 1.18 |  $J_{Fe-Cu}$ = 1.16 |  $J_{Fe-Cu}$ = -8.31
|  |  $J_{Cu-Mn}$ = -2.42 |  $J_{Cu-Mn}$ = -2.91 |  $J_{Cu-Mn}$ = -0.77
Figure 4: Depicting the nature of interaction (by colored bonds) between the
nearest neighbour magnetic moments (on the basis of Table III) for the for (a)
x=0, (b) x=0.125, (c) x=0.25 (confg-I), (d) x=0.25 (confg-III). Yellow and red
colored bonds represent the ferromagnetic (FM) and anti-ferromagnetic (AFM)
interaction, respectively. Stacking of pyramids along $c$ is shown on the
right side for reference.

From Table 3, we observe that the $ab$-plane exchange interactions between Fe-
Fe, Cu-Cu and Fe-Mn remain antiferromagnetic in nature with Mn doping similar
to the parent compound. The largest interaction among all is Cu-Cu exchange
which is seen to decrease with Mn doping. This is in line with the
experimental observation that with Mn doping the magnetic transition
temperatures of both commensurate and incommensurate phases decreasezhang2021
. Looking at the exchange interactions along $c$ direction we find that for
parent compound there exists ferromagnetic interaction within the bi-pyramid
(JcBP, see Fig.4) whereas inter-bilayer (JcIP) it is antiferro which is
consistent with the experimentally observed commensurate magnetic
structuremorin2015 . When we dope 12.5$\%$ Mn at Fe sites (x=0.125), we
observe that inter-bilayer Cu-Mn exchange becomes ferromagnetic whereas the
corresponding Fe-Cu exchange remains antiferro. The magnitude of this inter-
bilayer Cu-Mn exchange is also larger than other Fe/Mn-Cu exchange along
$c$-direction. Therefore, this causes frustration in the commensurate magnetic
interactions forcing the moments to cant away. When we move to higher Mn
doping (i.e. x = 0.25, confg-I) we observe similar behaviour as in x = 0.125.
Interestingly, in confg-III of x = 0.25, we observe that along c direction the
exchange interactions within the bi-pyramid have become antiferromagnetic
whereas the inter-bilayer ones ferromagnetic; just opposite to the case in
parent compound. These observations have been depicted in Fig.4. In this case
the frustration is therefore released. This can be further understood from the
electronic structure which we present in the following section.

###  III.3 Electronic structure

We performed electronic structure calculations for undoped and Mn doped YBCFO
for their respective magnetic ground state within GGA+U+SO approximation as
shown in Fig.3. We present in Fig.5 the total density of states (TDOS) where
we observe insulating state for the undoped (x = 0.0) and 12.5$\%$ Mn doped (x
= 0.125) case (see Fig.5(a) and (b)) whereas for 25 $\%$ Mn doped case
(x=0.25) out of the four configurations we have studied, three (confg-I,II and
IV) are insulating and one (confg-III, Fig.3(e)) is metallic. In Fig.5(c) we
present the TDOS for the lowest energy configuration (i.e. confg-I, 3(c)) at
x=0.25 and in Fig.5(d) we present the same for highest energy configuration
(i.e. confg-III). We further observe that around the Fermi level Cu/Fe/Mn
states are present which strongly hybridized with oxygen states. To understand
the mechanism of different ground state spin orientations as discussed in the
previous section and also different transport behaviour (insulating/metallic)
we looked into the details of electronic structure around the Fermi level.

Figure 5: Atom projected total density of states (TDOS) calculated within
GGA+U+SO for (a) x=0, (b) x=0.125, (c) x=0.25 (confg-I) and (d) x=0.25 (confg-
III). Figure 6: Orbital projected partial density of states (PDOS) calculated
within GGA+U+SO for (a) x=0, (b) x=0.125, (c) x=0.25 (confg-I) and (d) x=0.25
(confg-II). Figure 7: 3D electron density plotted for the energy range : -0.25
eV to Fermi level (0 eV) showing the occupied orbitals at various sites for
(a) x=0, (b) x=0.125, (c) x=0.25 (confg-I) and (d) x=0.25 (confg-III).

We calculated the orbital projected partial density of states (PDOS) as well
as 3D electron densities which we present in Fig.6 and Fig.7, respectively. As
we see in Fig.6(a) the Cu d${}_{x^{2}-y^{2}}$ is the highest occupied orbital
in the parent compound (x = 0.0). Upon doping Mn at x = 0.125 we observe that
Mn d${}_{z^{2}}$ becomes the highest occupied orbital which hybridizes with Cu
d${}_{z^{2}}$ and oxygen p-orbital. Same situation is also observed at x =
0.25 with confg-I with even stronger hybridization between Mn d${}_{z^{2}}$,
Cu d${}_{z^{2}}$ and Oxygen pz orbital. This can be clearly seen in calculated
3D electron densities in the energy range -0.25eV to Fermi level (0 eV) which
are shown in Fig.7(a) to (c) for x = 0.0, 0.125 and 0.25 respectively.
Interestingly, in confg-III at x=0.25 (see Fig.6(d)) we observe that there is
a charge transfer from the occupied Mn d${}_{z^{2}}$ orbital to unoccupied Fe
d${}_{z^{2}}$ and Cu d${}_{x^{2}-y^{2}}$ orbitals making the system metallic.
From Fig.7(d) we can clearly see the occupied Cu d${}_{x^{2}-y^{2}}$ orbitals
similar to the parent compound. Thus we can conclude that with Mn doping Cu
d${}_{z^{2}}$ becomes the highest occupied orbital due to its hybridization
with Mn d${}_{z^{2}}$ in contrast to the parent compound where Cu
d${}_{x^{2}-y^{2}}$ orbital was the highest occupied orbital. The occupancy of
out of plane oriented Cu d${}_{z^{2}}$ orbital in the doped compounds in place
of the planar Cu d${}_{x^{2}-y^{2}}$ of parent compound facilitates the spin
moments canting towards $c$ direction upon doping Mn as observed in the
experiment. In confg-III of x = 0.25 as the Cu d${}_{x^{2}-y^{2}}$ gets
occupied through charge transfer from Mn, the spin moments prefer to lie in
the $ab$-plane as observed in our calculation (see Fig.3(e)). Therefore, the
orbitals play a major role in driving the spin moments out of plane in YBCFMO.

##  IV Conclusions

In summary, we carried out a thorough investigation to understand the
mechanism behind spin reorientation upon Mn doping at Fe sites in YBCFO using
density functional theory calculations. In this context, we studied the non-
collinear magnetic order, electronic structure and orbital state of doped and
undoped YBCFO. Through a detailed analysis of optimized crystal structures of
the undoped (x=0.0) and doped (x=0.125 and 0.25) compounds, we have found the
local Jahn-Teller distortions of CuO5, FeO5 and MnO5 pyramids which are
consistent with experimental results observed by Zhang et al.zhang2021 . From
our total energy calculations we have obtained the lowest energy
configurations for Fe/Cu distribution as well as Mn distribution in the three
compounds. Considering the respective lowest energy configurations we have
performed non-collinear magnetic calculations within GGA+U+SO approximation to
obtain the ground state magnetic order. Considering the commensurate magnetic
structure for all the three compounds, we observe that the spin moments indeed
prefer to align along crystallographic $c$ direction in the doped compounds
whereas in parent compound they are aligned along (111) direction. This is
consistent with the trend observed in the experimental measurements. Further,
our estimates of nearest neighbour magnetic exchange interactions in the
$ab$-plane and along $c$ direction reveal that the Cu-Cu exchange in
$ab$-plane which is AFM and the largest amongst all, decreases with Mn doping.
This is again in line with the experimental observation of lowering of
transition temperature upon doping. More interestingly, we find that Mn doping
gives rise to a frustrating inter-bilayer ferromagnetic exchange which
otherwise should have been antiferromagnetic in the commensurate magnetic
structure. This we believe may be one of the reasons behind spins canting away
from their corresponding orientation in the parent compound. The other
possible reason is the change in the nature of highest occupied Cu orbital. In
parent compound, the orbital is d${}_{x^{2}-y^{2}}$ which is a planar orbital
but in the doped compound due to hybridization with Mn d${}_{z^{2}}$, Cu
d${}_{z^{2}}$ becomes the highest occupied orbital. This change in orbital
state we believe drives the spin moments out of plane towards $c$ direction in
the doped compounds. We found one particular configuration of Mn distribution
(not lowest energy) in x=0.25 case which is metallic due to charge transfer
from Mn d${}_{z^{2}}$ to unoccupied Cu d${}_{x^{2}-y^{2}}$ in the conduction
band. Interestingly, we find that in this case also the spin moments prefer to
lie in the $ab$-plane. Therefore, the spin orientation in the doped compound
appears to be orbitally driven aided by the frustrating exchange interactions.

##  V Acknowledgement

MS would like to acknowledge the Ministry of Education, India for the research
fellowship. MS and TM acknowledge the National Supercomputing Mission (NSM)
for providing computing resources of ‘PARAM Ganga’ at the Indian Institute of
Technology Roorkee, which is implemented by C-DAC and supported by the
Ministry of Electronics and Information Technology (MeitY) and Department of
Science and Technology (DST), Government of India. TM acknowledeges Science
and Engineering Research Board (SERB), India for funding support through
MATRICS research grant (MTR/2020/000419).

## References

  * (1) Yoshinori Tokura, Science 312 1481-1482 (2006).
  * (2) W. Eerenstein, N. D. Mathur, and J. F. Scott, Nature 442 759-765 (2006).
  * (3) Buurma, A. J. C., G. R. Blake, T. T. M. Palstra, and Umut Adem. ”Multiferroic materials: physics and properties.” Reference Module in Materials Science and Materials Engineering (2016).
  * (4) D.I. Khomskii, Journal of Magnetism and Magnetic Materials 306 (2006).
  * (5) Daniel Khomskii, Physics 2, 20 (2009).
  * (6) T. Kimura, Y. Sekio, H. Nakmura, T. Siegrist, and A. P. Ramirez, Nature materials 7, 291 (2008).
  * (7) L. Er-Rakho, C. Michel, Ph. Lacorre, and B. Raveau, Journal of Solid State Chemistry 73, 531-535 (1988).
  * (8) B. Kundys, A. Maignan, and Ch. Simon, Applied Physics Letters 94, 072506 (2009).
  * (9) Yuji Kawamura, Takahiro Kai, Erika Satomi, Yukio Yasui, Yoshiaki Kobayashi, Masatoshi Sato, and Kazuhisa Kakurai, Journal of the Physical Society of Japan 79, 073705 (2010) .
  * (10) M. Morin, A. Scaramucci, M. Bartkowiak, E. Pomjakushina, G. Deng,1, D. Sheptyakov, L. Keller, J. Rodriguez-Carvajal, N. A. Spaldin, M. Kenzelmann, K. Conder, and M. Medarde, Physical Review B 91, 064408 (2015).
  * (11) Yen-Chung Lai, Chao-Hung Du, Chun-Hao Lai, Yu-Hui Liang, Chin-Wei Wang, Kirrily C Rule, Hung-Cheng Wu, Hung-Duen Yang, Wei-Tin Chen, G J Shu, and F-C Chou, Journal of Physics: Condensed Matter 29 145801 (2017).
  * (12) Dibyendu Dey, S. Nandy, T. Maitra, C. S. Yadav, and A. Taraphder Scientific Reports 8 2404 (2018).
  * (13) Andrea Scaramucci, Hiroshi Shinaoka, Maxim V. Mostovoy, Markus Müller, Christopher Mudry, Matthias Troyer, and Nicola A. Spaldin, Phys. Rev. X 8, 011005 (2018).
  * (14) Mickaël Morin, Emmanuel Canevt2, Adrien Raynad1, Marek Bartkowiak, Denis Sheptyakov, Voraksmy Ban, Michel Kenzelmann, Ekaterina Pomjakushina, Kazimierz Conder, and Marisa Medarde, Nature Communications 7 13758 (2016).
  * (15) Tian Shang, Emmanuel Canévet, Mickaël Morin1, Denis Sheptyakov, María Teresa Fernández-Díaz, Ekaterina Pomjakushina, and Marisa Medarde, Science Advances 4 eaau6386 (2018).
  * (16) S. Lal, S. K. Upadhyay, K. Mukherjee, and C. S. Yadav; Euro Phys. Letts. 117, 67006 (2017).
  * (17) M. J. Ruiz-Aragón, E. Morán, U. Amador, J. L. Martınez, N. H. Andersen , and H. Ehrenberg, Physical Review B 58 10 (1998).
  * (18) A W Mombru, K Prassides, C Christides, R Erwin, M Pissas, C Mitros, and D Niarchos, Journal of Physics: Condensed Matter 10 1247–1258 (1998).
  * (19) V.Caignaert, I. Mirebeau, F. Bourée, N. Nguyen, A. Ducouret, J-M. Greneche, and B.Raveau, Journal of Solid State Chemistry 114 24-35 (1995).
  * (20) Xiaodong Zhang, Arnau Romagueraa, Oscar Fabelob, Francois Fauthc, Javier Herrero-Martínc, José Luis García-Muñoza, Acta Materialia 206 116608 (2021).
  * (21) Hosho Katsura, Naoto Nagaosa, and Alexander V. Balatsky, Physical Review Letters 95, 057205 (2005).
  * (22) G. Kresse and J. Hafner, Physical Review B $\mathbf{47}$, 558 (1993).
  * (23) G. Kresse and J. Furthmüller, Physical Review B $\mathbf{54}$, 11169 (1996).
  * (24) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E.Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Physical Review Letters 100, 136406 (2008).
  * (25) D D Koelling and B N Harmon, Journal of Physics C: Solid State Physics $\mathbf{10}$, 3107 (1977).
  * (26) H. J. Xiang, E. J. Kan, S.-H. Wei, M.-H. Whangbo, and X. G. Gong, Phys. Rev. B 84, 224429 (2011).


# Computational Chemistry on Quantum Computers

Jerimiah B. Wright
Department of Chemistry
University of Texas
Austin, TX 78705
<EMAIL_ADDRESS>
All simulations studied used Qumquat, an experimental high-level quantum
programming language. https://github.com/patrickrall/Qumquat

###### Abstract

The purpose of this experiment was to use the known analytical techniques to
study the creation, simulation, and measurements of molecular Hamiltonians.
The techniques used consisted of the Linear Combination of Atomic Orbitals
(LCAO), the Linear Combination of Unitaries (LCU), and the Phase Estimation
Algorithm (PEA). The molecules studied were $H_{2}$ with and without spin, as
well as $He_{2}$ without spin. Hamiltonians were created under the LCAO basis,
and reconstructed using the Jordan-Winger transform in order to create a
linear combination of Pauli spin operators. The lengths of each molecular
Hamiltonian greatly increased from the $H_{2}$ without spin, to $He_{2}$. This
resulted in a reduced ability to simulate the Hamiltonians under ideal
conditions. Thus, only low orders of l = 1 and l = 2 were used when expanding
the Hamiltonian in accordance to the LCU method of simulation. The resulting
Hamiltonians were measured using PEA, and plotted against function of
$\frac{2\pi(K)}{N}$ and the probability distribution of each register. The
resolution of the graph was dependent on the amount of registers, N, being
used. However, the reduction of order hardly changed the image of the $H_{2}$
graphs.

##  1 Introduction

###  1.1 Computational Chemistry

Computational chemistry is the combination and application of chemical,
mathematical and computing skills, used to solve chemical problems. These
problems range from the study of molecular structure, electronic charge
density distributions, reactivity, and many more. These computational results
are usually used to complement experimental results, however as they can also
shed light on unobserved chemical phenomena, they are also used to design new
drugs and materials.

The Schroëdinger equation is one of the starting points for most quantum
chemical calculations. It can be used to predict or calculate the energy of
atoms as particles move around. However, for many body systems, the
Schroëdinger equation cannot be solved analytically. To this end, many chemist
use different mathematical methods to approximate solutions to the equation.
These methods range from simplified forms of the first-principles equations
that are easier or faster to solve, to methods that limit the size of the
system such as periodic boundary conditions. Because approximation methods are
required to achieve solutions to the Schroëdinger equation, the goal of
computational chemist has been to minimize the residual error of the
approximations while keeping the calculations tractable.

Performing any kind of chemical computation requires a method of representing
the system. These methods are known as basis sets, and are typically comprised
of a set of functions that are used to represent a wave function. One of the
more common basis sets is composed of the Linear Combination of Atomic
Orbitals. LCAO is a quantum superposition of atomic orbitals and a technique
for calculating molecular orbitals in quantum chemistry. Assuming the
molecular orbital $\phi_{i}$ is equal to the number of atomic orbitals
included in the linear expansion, the expression
$\phi_{i}=\sum_{j}c_{ij}\chi_{j}$ can be formed where where $\phi_{i}$ is a
molecular orbital represented as the sum of n atomic orbitals $\chi_{r}$, each
multiplied by a corresponding coefficient $\ c_{ij}$, and j represents which
atomic orbital is combined in the term. This method can be incredibly useful
when trying to describe molecular orbitals because of its simplicity.

###  1.2 The Molecular Hamiltonian

A Hamiltonian operator, Ĥ, is an operator corresponding to the sum of the
potential and kinetic energies of all particles in a system, $\ket{\psi}$.
Using this operator, the Schroëdinger equation can be expressed as
$\hat{H}\ket{\psi}=E_{a}\ket{\psi}$. In the case of molecules, this
Hamiltonian is known as the Molecular Hamiltonian. The most common
mathematical approximation is the Bohrn and Oppenheimer approximation in which
nuclear kinetic energy terms within the Hamiltonian are neglected. This
results in a Hamiltonian that only considers the kinetic energies of electrons
and the Coulomb interactions between them. This new Hamiltonian is known as
the Electronic Hamiltonian, and is shown in its second quantized form below.

$\hat{H}=\sum_{i,j}h_{i,j}a^{\dagger}_{i}a_{j}+\dfrac{1}{2}\sum_{i,j,k,l}h_{i,j,k,l}a^{\dagger}_{i}a^{\dagger}_{k}a_{j}a_{l}$
(1)

The Electronic Hamiltonian’s second quantized form tracks each orbital and
stores whether there is an electron present in each of them. In order to do
this annihilation and creation operators are used as they can act on
electronic states. The annihilation operator is denoted $a$, and it lowers the
number of particles in a given state by one. The creation operator is denoted
$a^{\dagger}$ and it increases the number of particles in a given state by
one, as well as being the complex conjugate of the annihilation operator. The
coefficients $h_{ij}$ and $h_{ijkl}$ are one-and two-electron overlap
integrals which can be computed classically.

###  1.3 Quantum Representation

Computational chemist often use computers to solve chemical problems like
molecular structure. Classical software packages are often limited as they
must store an electrons wave function of each particle as a probability
distribution. This becomes a problem as the memory required to store the wave
function in the form of bits scales exponentially with the number of electrons
in the molecule. With the introduction of quantum bits, or qubits, a particles
wave function can be more efficiently stored. This is because qubits, as the
name implies, are quantum allowing them to be described as wave functions
themselves.

In practical applications, the qubits used for computation must be
distinguishable. This becomes an issue when they are used to represent the
electrons of molecular systems as the electrons are indistinguishable
particles. In order to circumvent this issue, the Jordan-Winger transform is
used to express the annihilation operator $a_{j}$, and the creation operator
$a^{\dagger}_{j}$, in terms of the Pauli spin operators X,Y,Z and the identity
I, that correspond to the algebra of distinguishable $\nicefrac{{1}}{{2}}$
spin particles. The Jordan-Winger transform is given by the equation below.

$\begin{split}a_{j}\equiv I^{\otimes n-j-1}\otimes\hat{Q}^{-}\otimes
Z^{\otimes j}\\\ a^{\dagger}_{j}\equiv I^{\otimes
n-j-1}\otimes\hat{Q}^{+}\otimes Z^{\otimes j}\end{split}$ (2)

*Note that $\hat{Q}^{-}=\ket{0}\bra{1}=\frac{1}{2}(X+iY)$ and $\hat{Q}^{+}=\ket{1}\bra{0}=\frac{1}{2}(X-iY)$.

Once the Hamiltonian has been expressed as a sum of elementary Pauli spin
operators the dynamics can then be compiled into fundamental gate operations
using a host of well-known techniques. One such technique is the Linear
Combination of Unitaries, or LCU.

The LCU algorithm is used to solve the Hamiltonian Simulation problem, where
the goal is construct a unitary operation U that follows the equation shown
below.

$\norm{U-e^{-iHt}}\leq\epsilon$ (3)

*Note that $\epsilon$ represents the error parameter.

Suppose there is unitary that can can be written as a linear combination of
efficiently implementable unitary matrices Ui , i.e., $V=\sum_{i}a_{i}U_{i}$.
This unitary U then maps $\ket{i}\ket{\psi}$ to $\ket{i}U_{i}\ket{\psi}$. This
means that if the Hamiltonian itself is a linear combination of unitaries
where Ui is unitary for all i, then the map $e^{iHt}$, is also a linear
combination of unitaries.

$e^{-iHt}=\sum_{l=0}^{\infty}\frac{1}{l!}(-iHt)^{l}=\sum_{l=0}^{\infty}\frac{1}{l!}(-i(\sum_{i}a_{i}U_{i}t)^{l}$
(4)

When expanded out, each term is a scalar multiplied by a product of unitary
operators. Although this is an infinite sum, a good approximation can be
obtained by truncating the sum. Using this idea, a similar method can be
applied to simulate sparse Hamiltonians that are not restricted to being
Hermitian. The new operator V is defined as
$V:=e^{{-iH}/{m}}=e^{{-i\sum_{j}U_{j}}/{m}}$. It is important to note that the
coefficient $a_{i}$ can be absorbed int the unitary $U_{i}$ without any loss
in relative phase. Expanding the definition of the matrix exponential, and
truncating the sum after l = k results in the matrix $\tilde{V}$.

$\tilde{V}:=\sum_{l=0}^{k}\dfrac{1}{l!m^{l}}(-i\sum_{j}U_{j})^{l}$ (5)

##  2 Experimental

###  2.1 Creating the Hamiltonian

The molecules of interest in this experiment were diatomic hydrogen and
helium. A unique form of the electronic Hamiltonian (1) was used to simulate
the quantum chemistry of the molecules. This equation is shown below.

$\hat{H}=\sum_{i,j}h_{ij}b^{\dagger}_{i}b_{j}$ (6)

The Hamiltonian was created using the LCAO method to rewrite the molecular
orbitals ,$b^{\dagger}b$, in terms of separable atomic orbitals. By treating
each atomic orbital as a vector, ie. $\ket{1s_{A}}$, a matrix of inner
products, $\bra{1s_{A}}\ket{1s_{B}}$ could be produced. This matrix is known
as the overlap matrix.

Figure 1: The overlap between the 1s orbitals of Hydrogen

In order to write the the molecular orbitals as separable atomic orbitals, the
Gram-Schmidt process was used to create an orthonormal set of virtual orbitals
that do not overlap. These virtual orbitals were substituted into equation [6]
in order to produce a separable molecular Hamiltonian. The Jordan-Winger
transform was performed on the new molecular Hamiltonians, resulting in their
decomposition into Pauli spin operators.

###  2.2 Hamiltonian Simulation

Once the Hamiltonian has been expressed as a sum of elementary Pauli spin, the
LCU method could be used to simulate Hamiltonian dynamics. First, A unitary
matrix, A, was created to map $\ket{0}$ to
$\frac{1}{\sqrt{a}}\sum_{i}\sqrt{a_{i}}\ket{i}$. This was done using the
postselect feature in Qumquat, the Quantum Machine Learning and Quantum
Algorithms Toolkit. The Palui string Hamiltonians were then expanded in
accordance to equation [5]. The hydrogen molecule without spin was expanded
with orders of l = 1 and l = 2, while the hydrogen molecule with spin was
expanded only with an order of l = 2. The helium molecule without spin was
expanded with an order of l = 1. The resulting expand Pauli strings were
applied to the quantum circuit, conversting the eigenenergy into relative
phases. The relative phases can then be measred using the Phase Estimation
algorithm.

###  2.3 Measuring the Hamiltonian

The phase estimation algorithm, or PEA, was used to obtain the molecular
energies from the time evolution of the molecular Hamiltonian described in
section 2.2. Given a unitary matrix U and a quantum state $\ket{\psi}$ such
that U$\ket{\psi}=e^{2\pi i\theta}$, the algorithm estimates the value of
$\theta$ with high probability with additive error $\epsilon$. In the case of
calculating molecular energies, the unitary matrix U is the molecular
Hamiltonian produced after simulating its dynamics. The algorithm begins by
setting up two registers where the upper n qubits comprise the first register,
and the lower m qubits are the second register. A super position of the
initial state is created by applying n-bit Hadamard gate operations on the
first register, resulting in the state below.

$\frac{1}{2^{\nicefrac{{n}}{{2}}}}(\ket{0}+\ket{1})^{\otimes n}$ (7)

A controlled-U gate is applied to the second register only if its
corresponding control qubit is $\ket{1}$. Applying n many controlled-U gate
operations to the phase of an eigenstate of the Hamiltonian evolving dependent
on a register qubit, i.e.
$\ket{0}\ket{\psi_{n}}+e^{-iHt}\ket{1}\ket{\psi_{n}}=\ket{0}\ket{\psi_{n}}+e^{-iE_{n}t}\ket{1}\ket{\psi_{n}}$,
allows the state of the first register to be described by the following
equation.

$\frac{1}{2^{\nicefrac{{n}}{{2}}}}\sum_{k=0}^{N}e^{-iE_{n}t}$ (8)

* Note that N = $2^{n}-1$

When $E_{n}=2\pi(\theta-K)/t$ where $0\leq\theta\leq 1$ and K is an
unobservable integer, the unknown eigenvalue becomes encoded in the relative
phase of the register quantum state as $\ket{0}+e^{-2\pi i(\theta-K)}\ket{1}$.
The $\theta$ can then be estimated by applying an inverse quantum Fourier
transform, $\frac{1}{2^{\nicefrac{{n}}{{2}}}}\sum_{k=0}^{N}e^{\frac{-2\pi
in\theta}{N}}$, and performing a measurement in the computational basis on the
first register.

##  3 Results

$H_{2}(no-spin)=\blockarray{ccc}1S_{A}&1S_{B}\\\ \block{[cc]c}1S1S_{A}\\\
S11S_{B}\\\ $

$H_{2}(spin)=\blockarray{ccccc}1s_{A\uparrow}&1s_{A\downarrow}1s_{B\uparrow}1s_{B\downarrow}\\\
\block{[cccc]c}1S001s_{A\uparrow}\\\ S1001s_{A\downarrow}\\\
001S1s_{B\uparrow}\\\ 00S11s_{B\downarrow}\\\ $

$He_{2}(no-spin)=\blockarray{ccccc}1s_{A}&1s_{B}2s_{A}2s_{B}\\\
\block{[cccc]c}1S0S_{1}1s_{A}\\\ S1S_{1}01s_{B}\\\ 0S_{1}1S_{2}2s_{A}\\\
S_{1}0S_{2}12s_{B}\\\ $

Figure 2: Overlap Matrices for $H_{2}$ without spin, $H_{2}$ with spin, and
$He_{2}$ without spin

The LCAO approach was used in order to create a basis set in which the
molecular Hamiltonian was a linear combination of atomic orbitals. Writing the
Hamiltonian in this format aided in state preparation for Hamiltonian
Simulation.The molecules were written in terms of their overlap matrix (Figure
2), and their atomic orbitals were separated using the Gram-Schmit process.
The original orbitals were then written in terms of these new virtual
orbitals, in order to be used in equation [6]. Performing the Jordan-Winger
transform results in the Pauli decompositions of molecular Hamiltonians where
each varied in length increasing from $H_{2}$ to $He_{2}$(Table 1). The amount
of gate operations greatly increased the time it took to simulate and measure
the results of each Hamiltoniain.

Table 1: Table of Molecular Hamiltonians comprised of sums of Pauli operators $H_{2}$ (no-spin) |  |  $H_{2}$ (spin) |  |  $He_{2}$ (no-spin) |
---|---|---|---|---|---
1.5686986355290005 | II | 3.137397271058001 | IIII | 22.452642369057198 | IIII
-1.2843493177645002 | ZI | -1.2843493177645002 | ZIII | -1.511819721289001 | ZIII
-0.2843493177645001 | IZ | -0.2843493177645001 | IZII | -1.3117878918224433 | IZII
-0.4722596673392581 | XY | -1.2843493177645002 | IIZI | -5.029764413362022 | IIZI
0.4722596673392581 | YX | -0.2843493177645001 | IIIZ | -14.599270342583733 | IIIZ
|  | -0.4722596673392581 | XYII | -0.09446696524113857 | XYII
|  | 0.4722596673392581 | YXII | 0.09446696524113857 | YXII
|  | -0.4722596673392581 | IIXY | -0.9924233565701107 | XZYI
|  | 0.4722596673392581 | IIYX | 0.9924233565701107 | YZXI
|  |  |  | -1.1191380207487813 | IXXI
|  |  |  | -1.1191380207487813 | IYYI
|  |  |  | 1.8223341943756635 | XZZY
|  |  |  | -1.8223341943756635 | YZZX
|  |  |  | 3.026611588811937 | IXZX
|  |  |  | 3.026611588811937 | IYZY
|  |  |  | -7.9506036387717955 | IIXX
|  |  |  | -7.9506036387717955 | IIYY

Hamiltonian Simulation was done using the LCU method. Since the generated
Hamiltonians were already linear combinations of Pauli spin operators,
promoting them into equation [6] was simple. However, because of the length of
each Hamiltonian, only low orders of l could be studied, as the length would
increase exponentially, drastically increasing the time it took to simulate
the results.

The PEA algorithm was used to measure the eigenvalue of the of the Hamiltonian
operator by encoding the relative phase of the register quantum state as
$\ket{0}+e^{-i2\pi(\theta-K)/t}$ and applying an inverse quantum Fourier
transform, them measuring on the first register. Plotting these results as a
function of $\frac{2\pi(K)}{N}$ and probability distribution of each register,
produced the graphs shown in figures 3 and 4

(a) N = 100, l = 2

(b) N = 10, l = 2

(c) N = 10, l = 1

Figure 3: Plots of the $H_{2}$ molecule without spin under different
parameters.

(d) $H_{2}$ molecule with spin. N = 10, O = 2

(e) $He_{2}$ molecule with spin. N = 10, O = 1

Figure 4: Plots of the $H_{2}$ and $He_{2}$ molecules

##  4 Conclusions

The graphs shown in figures 3 and 4 demonstrate the results under different
values of N and l, where N represents the range of register used, and l
represents the order of the of the molecular Hamiltonian produced from
equation [5]. As no error mitigation was done, these results should be studied
qualitatively. As shown in figure 1, decreasing the range of N drastically
decreases the resolution of the graph, with the benefit of producing a graph
much quicker. Comparing 3(b) to 3(c) shows the result of changing the order of
the Hamiltonian expansion. For low values of l, the effect of changing the
order seems to be fairly negligible. This is explored further in the appendix.
Finally, the more complex Hamiltonians of $H_{2}$ with spin, and $He_{2}$
without spin are shown in figure 2a and 2b respectively. Despite Both graphs
having low values of N and l, they have distinct shapes. However, $H_{2}$
molecule with spin shares a nearly identical appearance to its no-spin counter
part. While these results must be interpreted qualitatively, they still help
prove the methodology of producing the Hamiltonians as legitimate through
reproducibility. If this method was flawless, it would be expected that the
spin variant would have a higher quality graph despite having the same N and l
values, due to a more complete description of the molecule. The $He_{2}$
molecule has a more unique shape indicating a clear difference between the two
molecules.

The purpose of this experiment was to use the known analytical techniques to
study the creation, simulation, and measurements of molecular Hamiltonians.
While this experiment was successful in implementing known techniques to
perform electronic structure calculations, it was extremely limited due to the
performance capabilities of classical computers. This greatly hindered the
ability for quantitative results due to using low resolution graphs to
identify the eigenvalues of each molecular Hamiltonian.

## References

  * [1] Whitfield, James D., Jacob Biamonte, and Alán Aspuru-Guzik.  simulation of electronic structure Hamiltonians using quantum computers.  2011\. S Molecular Physics 109(5): 735-750.
  * [2] Aarronson, Scott.  “Quantum Computing, Postselection, and Probabilistic Polynomial-Time.”  Arxiv.org, arxiv.org/pdf/quant-ph/0412187.pdf.
  * [3] Clark, Jim.  “ ELECTRONIC STRUCTURE AND ATOMIC ORBITALS.” Chemguide helping you understand Chemistry  www.chemguide.co.uk/basicorg/bonding/orbitals.html.
  * [4] Tuckerman, Mark E.  “CHEM-UA 127: Advanced General Chemistry I.” LINEAR COMBINATION OF ATOMIC ORBITALS,  www.nyu.edu/classes/tuckerman/adv.chem/lectures/lecture 14/lecture 14.pdf.
  * [5] Kothari, Robin  "Efficient algorithms in quantum query complexity" University of Waterloo’s institutional repository,  https://uwspace.uwaterloo.ca/bitstream/handle/10012/8625/Kothari Robin.pdf