# First-principles calculation of the electronic and optical properties of
Gd2FeCrO6 double perovskite: Effect of Hubbard U parameter

Subrata Das a,b, M. D. I. Bhuyan a and M. A. Basith a<EMAIL_ADDRESS>aDepartment of Electrical and Electronic Engineering, Bangladesh University of
Engineering and Technology, Dhaka, 1000, Bangladesh , b Nanotechnology
Research Laboratory, Department of Physics, Bangladesh University of
Engineering and Technology, Dhaka-1000, Bangladesh.

DOI: 10.1016/j.jmrt.2021.06.026

###### Abstract

We have synthesized Gd2FeCrO6 (GFCO) double perovskite which crystallized in
monoclinic structure with P21/n space group. The UV-visible and
photoluminescence spectroscopic analyses confirmed its direct band gap
semiconducting nature. Here, by employing experimentally obtained structural
parameters in first-principles calculation, we have reported the spin-
polarized electronic band structure, charge carrier effective masses, density
of states, electronic charge density distribution and optical absorption
property of this newly synthesized GFCO double perovskite. Moreover, the
effects of on-site d-d Coulomb interaction energy (Ueff) on the electronic and
optical properties were investigated by applying a range of Hubbard Ueff
parameter from 0 to 6 eV to the Fe-3d and Cr-3d orbitals within the
generalized gradient approximation (GGA) and GGA+U methods. Notably, when we
applied Ueff in the range of 1 to 5 eV, both the up-spin and down-spin band
structures were observed to be direct. The charge carrier effective masses
were also found to enhance gradually from Ueff = 1 eV to 5 eV, however, these
values were anomalous for Ueff = 0 and 6 eV. These results suggest that Ueff
should be limited within the range of 1 to 5 eV to calculate the structural,
electronic and optical properties of GFCO double perovskite. Finally we
observed that considering Ueff = 3 eV, the theoretically calculated optical
band gap $\sim$1.99 eV matched well with the experimentally obtained value
$\sim$2.0 eV. The outcomes of our finding imply that the Ueff value of 3 eV
most accurately localized the Fe-3d and Cr-3d orbitals of GFCO keeping the
effect of self-interaction error from the other orbitals almost negligible.
Therefore, we may recommend Ueff = 3 eV for first-principles calculation of
the electronic and optical properties of GFCO double perovskite that might
have potential in photocatalytic and related solar energy applications.

##  I Introduction

Over the past few decades, the double perovskite oxides A2BB′O6 (A = rare
earth metal ions, B = transition metal ions) have gained immense research
interest owing to their rich multifunctional characteristics and potential
applications in the next generation spintronic devices Vasala and Karppinen
(2015); Gray  _et al._ (2010); Gaikwad  _et al._ (2019); Das  _et al._ (2008);
Cerón  _et al._ (2019); Yin  _et al._ (2019). Recently, several double
perovskites have been reported to possess promising applicability in the field
of photocatalysis, photovoltaic devices and photo(electro)chemical energy
storage systems Lin  _et al._ (2021); Mohassel  _et al._ (2020); Kangsabanik
and Alam (2020); Yin  _et al._ (2019). Especially, A2FeCrO6 double perovskites
(A = Pr, Bi etc.) having two 3d transition elements at B and B′ sites have
demonstrated fascinating optoelectronic properties such as favorable band gap
energy, strong absorbance in the visible regime of the solar energy etc Wu
_et al._ (2020); Gaikwad  _et al._ (2019); Nechache  _et al._ (2015).
Therefore, it is intriguing to investigate other members of A2FeCrO6 double
perovskite family for enormous applications in electronics, photochemistry and
optical technologies.

Nevertheless, it is difficult to synthesize perfectly ordered structure of
A2FeCrO6 double perovskites because of the complementary ionic radii of Fe and
Cr ions Gray  _et al._ (2010); Nair  _et al._ (2014). Hence, the double
perovskites of this family is yet less explored as compared to analogous
double perovskite materials. Recently, we have successfully synthesized
nanoparticles of double perovskite Gd2FeCrO6 (GFCO) for the first time by
optimizing the synthesis condition of a citrate-based sol-gel technique and
extensively investigated their crystallographic and chemical structure as well
as magnetic and optical behaviors Bhuyan  _et al._ (2021). Interestingly, the
favorable surface morphology, optimal direct band gap of $\sim$2.0 eV and the
band edge positions of synthesized GFCO double perovskite have revealed its
promising potential for visible light driven photocatalysis and related
applications. Unfortunately, the complexity of experimental conditions and to
some extent, the unavailability of the required experimental set-up posed
difficulty to understand their optical and electronic characteristics at
atomic level. Such impediment can be overcome by performing density functional
theory (DFT) based first-principles calculation systematically Zhang  _et al._
(2016).

However, it should be noted that the standard DFT methods, for instance local
density approximation (LDA) and generalized gradient approximation (GGA) of
the exchange-correlation functional have some limitations in analyzing
correctly the electronic properties of strongly correlated systems like GFCO
Shenton  _et al._ (2017). To be specific, the LDA and GGA methods have
demonstrated systemic failures to explain the on-site Coulomb interactions of
highly localized electrons because of the erroneous electron self-interaction
Himmetoglu  _et al._ (2014). One of the notable deficiencies of standard
approximations is the underestimation of the optical band gap in
semiconducting double perovskites which might have potentiality in
photocatalytic and optoelectronic applications Terakura  _et al._ (1984);
Brown and Page (2020). These limitations of standard DFT can be reasonably
corrected by GGA+U method in which on-site Hubbard-like correction is applied
to the effective potential Lu and Liu (2014); Wang  _et al._ (2006); Zhang
_et al._ (2010). Notably, two free parameters, U and J are required to
effectively tune the on-site Coulomb and exchange interactions, respectively
Shenton  _et al._ (2017). In an approach proposed by Dudarev et. al. Dudarev
_et al._ (1998), these two parameters can be combined into a single Hubbard
Ueff correction parameter where Ueff = U-J Wang  _et al._ (2006).

Typically, within GGA+U calculation, the Ueff value is selected such that the
calculated band gap matches with the experimentally obtained one Shenton  _et
al._ (2017). However, a number of recent investigations Shenton  _et al._
(2017); Perdew (1985) have demonstrated that choosing Ueff parameter only to
match the band gap, may introduce spurious effects for strongly correlated
materials. For instance, Shenton et al. Shenton  _et al._ (2017) reported that
a Ueff value of 5 eV or larger is required to match the experimental band gap
of BiFeO3. However, the ordering of the Fe d orbitals at the conduction band
minimum of BiFeO3 inverts for Ueff $>$ 4 eV. Hence, careful consideration is
required to apply the most accurate Hubbard Ueff parameter in GGA so that the
theoretical band gap value closely matches with the experimental one without
introducing any significant error in the character of electronic band edges.
To the best of our knowledge, the influences of Hubbard Ueff parameter on the
optical and electronic properties of double perovskites possessing two 3d
transition elements like GFCO have not been extensively investigated yet.

Therefore, in the current work, we have extensively investigated the effect of
Hubbard Ueff parameter on the crystallographic parameters, spin-polarized
electronic band structure and optical properties of our recently synthesized
GFCO nanoparticles by performing first-principles calculation via both GGA and
GGA+U methods. To ensure reliability, we have employed our experimentally
obtained structural parameters for theoretical analysis. We observed that the
variation in Ueff had insignificant effect on the structural parameters of
GFCO. However, the character and curvature of electronic band edges could not
be determined accurately without applying Ueff i.e. for Ueff = 0 eV and also,
for Ueff $>$ 5 eV in GGA+U calculation. Finally, considering the theoretically
calculated optical band gap, we conjectured that a Ueff value of 3 eV is
reasonable to employ to the Fe-3d and Cr-3d orbitals of GFCO double perovskite
within this calculation.

##  II Experimental and computational details

Double perovskite GFCO nanoparticles were synthesized by adopting a standard
citrate-based sol-gel technique, as discussed in details in our previous work
Bhuyan  _et al._ (2021). The crystallographic phase, lattice parameters, bond
angles and bond lengths of as-synthesized GFCO were determined by Rietveld
refinement of the powder XRD data using FullProf computer program package
Rodriguez-Carvajal (1990). An ultraviolet-visible (UV-visible)
spectrophotometer (UV-2600, Shimadzu) was used to obtain absorbance spectrum
of the as-synthesized GFCO for wavelengths ranging from 200 to 800 nm. Steady-
state photoluminescence (PL) spectroscopy was conducted at room temperature by
Spectro Fluorophotometer (RF-6000, Shimadzu). In the present investigation,
based on our experimental findings, spin-polarized optical properties and
accurate electronic band structure of as-prepared GFCO nanoparticles were
determined theoretically by DFT based first-principles calculation.

Figure 1:  Plane-wave cutoff energy convergence for structural optimization.

The theoretical calculations were carried out using both generalized gradient
approximation (GGA) and GGA+U methods within the plane wave pseudopotential
(PWPP) framework as implemented in the Cambridge Serial Total Energy Package
(CASTEP) Segall  _et al._ (2002); Rozilah  _et al._ (2017). The
crystallographic structural parameters obtained from the Rietveld refined
powder XRD spectrum of GFCO Bhuyan  _et al._ (2021) were employed for DFT
calculation. Prior to calculation, the geometry was optimized via Brodyden-
Fletcher-Goldfarb-Shanno (BFGS) scheme applying energy of 10-5 eV/atom,
maximum force of 0.05 eV/Å and maximum stress of 0.1 GPa Fischer and Almlof
(1992). The Gd-4f85s25p66s2, Fe-3d64s2, Cr-3s23p63d54s1 and O-2s22p4 electrons
were treated as valence electrons. The plane-wave cutoff energy convergence
result for structural optimization is demonstrated in Fig. 1 where the dashed
line represents the default energy cutoff. As can be observed, 450 eV energy
cutoff was found sufficient to achieve the converged ground-state energy of
GFCO. Hence, the plane-wave basis set was employed with the optimized energy
cutoff of 450 eV. Moreover, Brillouin-zone integration were carried out with a
5$\times$5$\times$3 Monkhorst-Pack k-point mesh Monkhorst and Pack (1976).
Spin polarized mode was endorsed during self-consistent field (SCF)
calculations and a SCF tolerance of 2$\times$10-6 eV per atoms was used.

Notably, to describe the exchange-correlation energy, at first we have used
the GGA (Ueff = 0 eV) method based on the Perdew-Burke-Ernzerhof functional
(PBE) within on-the-fly generated ultrasoft pseudopotentials (USP). Then a
spin-polarized calculation was carried out to confirm the exact ground state
of GFCO double perovskite Zhang  _et al._ (2010). Further, the GGA+U
calculation was performed with different values of Ueff to investigate the
effect of Hubbard Ueff parameter on the structural, electronic and optical
properties of the GFCO ground state Dudarev  _et al._ (1998). To be specific,
Ueff was varied from 1 to 6 eV for Fe-3d and Cr-3d orbitals whereas for Gd-4f
orbital, Ueff was kept fixed at 6 eV in accordance with previous
investigations Rozilah  _et al._ (2017); Chakraborty  _et al._ (2008).

The optical absorption coefficient was determined using the equation
$\leavevmode\nobreak\ \alpha=\sqrt[]{2}\omega\sqrt[]{\leavevmode\nobreak\
\sqrt[]{\varepsilon_{1}^{2}\left(\omega\right)+\varepsilon_{2}^{2}\left(\omega\right)}-\varepsilon_{1}\left(\omega\right)}$,
where $\varepsilon_{1}\left(\omega\right)$ and
$\varepsilon_{2}\left(\omega\right)$ denote frequency dependent real and
imaginary parts of dielectric function, $\omega$ represents the photon
frequency Saha  _et al._ (2000). $\varepsilon_{1}\left(\omega\right)$ was
calculated from the $\varepsilon_{2}\left(\omega\right)$ by the Kramers-Kronig
relationship Zhang  _et al._ (2008).

##  III Results and discussion

###  III.1 Crystal structure

Table 1: Lattice parameters, monoclinic angle and unit cell volume of
Gd2FeCrO6 for different values of Ueff obtained via first-principles
calculation along with the corresponding experimental values.

| Experimental value | Ueff= 0 eV | Ueff = 1 eV | Ueff = 2 eV | Ueff = 3 eV | Ueff = 4 eV | Ueff = 5 eV | Ueff = 6 eV
---|---|---|---|---|---|---|---|---
a ($\AA$) | 5.359 | 5.397 | 5.442 | 5.448 | 5.457 | 5.466 | 5.470 | 5.502
b ($\AA$) | 5.590 | 5.601 | 5.676 | 5.682 | 5.689 | 5.698 | 5.701 | 5.750
c ($\AA$) | 7.675 | 7.688 | 7.781 | 7.798 | 7.808 | 7.820 | 7.821 | 7.886
$\beta$ (∘) | 89.958 | 89.691 | 90.004 | 90.003 | 89.998 | 89.991 | 89.990 | 89.989
Volume($\AA$3) | 229.92 | 232.43 | 240.35 | 241.42 | 242.46 | 243.73 | 243.85 | 249.54

In our previous investigation Bhuyan  _et al._ (2021), we have extensively
investigated the crystallographic structure of our as-synthesized GFCO
nanoparticles at room temperature by performing Rietveld refinement analysis
of their powder XRD pattern. Notably, it was observed that GFCO crystallizes
in monoclinic structure with P21/n space group. The lattice parameters of GFCO
unit cell were found to be $a$ = 5.359(1) $\AA$, $b$ = 5.590(2) $\AA$, $c$ =
7.675(3) $\AA$, monoclinic angle $\beta$ = 89.958(1)∘ with cell volume 229.920
$\AA^{3}$. In the present work, we have performed first-principles calculation
to determine the structural parameters of GFCO by GGA (Ueff = 0 eV) and GGA+U
(Ueff = 1 to 6 eV) approaches. The calculated lattice constants $a$, $b$ and
$c$, monoclinic angle ($\beta$) along with unit cell volume are tabulated in
Table 1. For comparison, we have also included the experimentally obtained
structural parameters in the table. Noticeably, the lattice constants and
monoclinic angles obtained via first-principles calculation were found to be
slightly larger than the experimental values. Nevertheless, this mismatch
remained nominal i.e. within 3% of the experimental results which is
consistent with a number of previous investigations Shenton  _et al._ (2017);
Mann  _et al._ (2016). For instance, the lattice parameters of CO2-metal
organic framework calculated with Hubbard U corrections were reported to
remain within 3% of experimental values Mann  _et al._ (2016). Moreover, it
can be observed that the calculated lattice parameters and unit cell volume
enhanced with the increment of Ueff from 0 to 6 eV Shenton  _et al._ (2017).
In contrast, the monoclinic angle obtained by the GGA+U method decreased with
increasing Ueff in the range of 1 to 6 eV. Such observations are in well
agreement with previous investigations of related materials Shenton  _et al._
(2017); Neaton  _et al._ (2005).

Figure 2:  Experimentally obtained (a) Tauc plot for direct optical band gap
estimation and (b) steady-state photoluminescence spectrum of Gd2FeCrO6
nanoparticles Bhuyan  _et al._ (2021).

###  III.2 Experimentally obtained optical properties

The optical characteristics of as-synthesized GFCO nanoparticles were
extensively investigated by obtaining their UV-visible absorbance spectrum and
was reported previously Bhuyan  _et al._ (2021). The absorbance data was
employed to calculate the optical band gap of synthesized GFCO double
perovskite using Tauc relation Tauc  _et al._ (1966). The generated Tauc plot
for estimating the direct optical band gap of GFCO nanoparticles is shown in
Fig. 2(a). The abscissa intercept of the tangent to the linear region of the
curve demonstrated that the optical band gap value is $\sim$2.0 eV.

To further ensure the estimated direct optical band gap of GFCO perovskite,
the steady-state PL spectrum of the synthesized material was recorded for an
excitation wavelength of 230 nm Bhuyan  _et al._ (2021). The position of the
PL peak in Fig. 2(b) confirmed that the band gap value of GFCO is $\sim$1.98
eV which closely matches with the direct band gap value obtained from the Tauc
plot (Fig. 2(a)). Therefore, both the UV-visible and PL spectroscopic analyses
revealed that our as-synthesized nanostructured GFCO is a direct band-gap
double perovskite with a band gap value of $\sim$2.0 eV. The experimentally
observed optical properties of GFCO perovskite suggest the semiconducting
nature of GFCO and most importantly, demonstrates its ability to absorb light
of the visible spectrum of the solar illumination.

Figure 3:  Electronic band structure of Gd2FeCrO6 for (a) Ueff = 0 eV, (b)
Ueff = 1 eV, (c) Ueff = 3 eV, and (d) Ueff = 6 eV. Black and blue curves
represent up-spin and down-spin orientations, respectively. The energy ranges
from -3 to 3 eV and the zero is set to the Fermi energy EF.

###  III.3 Theoretical investigation of electronic properties

####  III.3.1 Electronic band structure

The theoretical investigation was initiated by calculating the spin-polarized
electronic band structure of GFCO within the GGA method (Ueff = 0 eV) as shown
in Fig. 3(a). The dotted horizontal line between the valence and conduction
bands represents the Fermi level Hou  _et al._ (2014). As can be observed, for
up-spin orientation, we obtained a direct electronic band gap of 0.84 eV
whereas for down-spin, the band gap is found to be indirect having a value of
0.27 eV. Notably, both of these values are much smaller than the direct
optical band gap $\sim$2.0 eV of as-synthesized GFCO as confirmed by UV-
visible and PL spectroscopic analyses. It is well known that electronic band
gap of any material would be larger than its optical band gap Bredas (2014),
hence we may infer that the band structure formed for Ueff = 0 eV is
incorrect. Therefore, we have determined the electronic band structure of GFCO
via GGA+U method with Ueff = 1 to 6 eV. The electronic band structures
obtained for Ueff = 1, 3, 6 eV are presented in Fig. 3(b), (c) and (d),
respectively and the band structures calculated for Ueff = 2, 4, 5 eV are
provided in Fig. S1 of Electronic Supplementary Information (ESI). As can be
seen in these two figures, the gap in the up-spin band enlarges with
increasing Ueff which can be attributed to the enhanced localization of the
Fe-3d and Cr-3d orbitals due to increased Ueff Shenton  _et al._ (2017). It is
worth noting that for Ueff = 1 to 5 eV, both the valence band maximum (VBM)
and conduction band minimum (CBM) were within the A symmetry point indicating
direct band structure.

However, at Ueff = 6 eV, we obtained an indirect up-spin band structure which
is inconsistent with all the previous cases. It should be noted that so far
both experimental and theoretical calculations provided strong evidence in
support of direct band structure of GFCO. Therefore, the indirect band
structure obtained for Ueff = 6 eV calls into question the applicability of
employing large Hubbard parameter (i.e. Ueff $>$ 5 eV) in first-principles
calculations of GFCO double perovskite Shenton  _et al._ (2017). Moreover, it
is interesting to note that the theoretically calculated up-spin band gap
values were within the range of 2.46 to 2.84 eV which ensures the
semiconducting nature of GFCO double perovskite as was also evident from the
UV-visible and PL spectroscopic analyses (Fig. 2).

In the case of down-spin orientation, for all values of Ueff (1 to 6 eV), both
the CBM and VBM were obtained within the G symmetry point indicating again
direct band structure. Interestingly, the band gap value increased
monotonically from Ueff = 1 to 5 eV but an anomalous decrease can be observed
at Ueff = 6 eV (Fig. 3(d)). Therefore, it might be conjectured that the
optimized value of Ueff for first-principles calculation of GFCO would be less
than 6 eV Shenton  _et al._ (2017).

Figure 4:  Variation in absolute charge carrier effective mass as a function
of Ueff. me∗ and mh∗ are the electron and hole effective masses, respectively
in units of the electron rest mass, mo

For understanding the carrier transport in the material, we have quantified
the curvature at band extrema by calculating the technologically important
charge carrier effective masses using following expression Reunchan  _et al._
(2016).

$m^{*}=\hbar^{2}\left(\frac{\mathrm{d}^{2}E}{\mathrm{d}k^{2}}\right)^{-1}$ (1)

Here, E is the band-edge energy as a function of wave-vector k. The electron
effective mass (me∗) was calculated by parabolic fitting of the E-k curve
within the small region of wave-vector near the CBM Reunchan  _et al._ (2016);
Das  _et al._ (2019). The hole effective mass (mh∗) was estimated by analyzing
the region near the VBM using similar approach Reunchan  _et al._ (2016). The
variations of me∗ and mh∗ of GFCO double perovskite as a function of Ueff are
illustrated in Fig. 4 both for up-spin and down-spin orientations. For up-spin
band structure, we can observe a reduction in me∗ from 15.5mo to 10.6mo for
Ueff = 0 eV to Ueff = 1 eV suggesting an increase in curvature at the CBM.
Between Ueff of 1 to 5 eV, the values of me∗ changed nominally ($\sim$1.3mo).
Further, for Ueff = 6 eV, an enhancement of 3.5mo can be noticed which
corresponds to the reduction in curvature at CBM. Clearly, in Fig. 4, we can
observe a similar dependence of mh∗ on Ueff both for up and down-spin
orientations. Moreover, for down-spin band structure, a notable increase can
be noticed in me∗ from 26.1mo to 81mo for Ueff = 5 eV to Ueff = 6 eV. However,
the variation was comparatively smaller ($\sim$ 11.2mo) in the range of Ueff =
0 to 5 eV. It is worth mentioning that while varying the values of Ueff for
electronic band structure calculation, we had kept the structural parameters
fixed at experimentally obtained values. Hence, we can confirm that the
variations in me∗ and mh∗ with Ueff are solely due to the electronic effects.
Such variation trend again justifies that it would be worthwhile to consider
Ueff within the range 1 to 5 eV for GFCO double perovskite.

####  III.3.2 Density of states

Figure 5:  Calculated total density of states (TDOSs) and partial density of
states (PDOSs) of Gd-4f, Fe-3d, Cr-3d and O-2p orbitals for both up-spin and
down-spin channels. The panels (a)–(d) show the DOSs for Ueff = 0, 1, 3, 6 eV,
respectively. The zero is set to the Fermi energy.

In order to resolve the contribution of each individual orbital to the
electronic bands of GFCO double perovskite, we have calculated the total
density of states (TDOS) and partial density of states (PDOS) for Gd-4f,
Fe-3d, Cr-3d and O-2p orbitals via GGA (Ueff = 0 eV) and GGA+U methods (Ueff =
1 to 6 eV). The DOSs obtained for Ueff = 0, 1, 3, 6 eV are presented in Fig.
5(a), (b), (c) and (d), respectively and the DOSs calculated for Ueff = 2, 4,
5 eV are provided in Fig. S2 of ESI.

Beginning with the character of the conduction band (E-EF$>$ 0 eV) obtained
for up-spin orientation, we can make the following observations. At Ueff = 0
eV, Fig. 5(a), the conduction band lying at around 1 eV was made up of a
hybridization of Fe-3d and O-2p orbitals. This band disappeared after
considering the effect of Hubbard Ueff which corresponds to the decrement in
me∗ at Ueff = 1 eV (Fig. 5(b)). Notably, the conduction band (at around 2.5
eV) obtained for Ueff = 1 eV had primarily the characteristics of Cr-3d with a
minor contribution from O-2p. Further, as shown in Fig. 5(c) and ESI Fig. S2,
with the increase of Ueff up to 5 eV, this conduction band shifted to higher
energy resulting in the enlargement of band gap. Also, the contribution of
Cr-3d to the up-spin conduction band enhanced with increasing Ueff.
Interestingly, at Ueff = 6 eV (Fig. 5(d)), we observed that the computed PDOS
for Gd-4f orbital radically altered and the conduction band arose mainly due
to the hybridization of Gd-4f with a minor contribution from Cr-3d orbital.
Such anomaly is another indication for the limitation of Ueff $>$ 5 eV to
explain the electronic band structure of GFCO double perovskite accurately.

Now, if we analyze the valence bands (E-EF$<$ 0 eV) of up-spin electronic band
structure shown in Fig. 5(a), it can be observed that the valence band at Ueff
= 0 eV was composed of Cr-3d, O-2p as well as Fe-3d orbitals. However, when we
employed Ueff, the contribution of Fe-3d got diminished. To be specific, the
up-spin valence band for Ueff = 1 eV was made up of the hybridization of Cr-3d
and O-2p states (Fig. 5(b)). From Fig. 5(c), (d) and ESI Fig. S2, it can be
noticed that with increasing Ueff, the contribution of Cr-3d state gradually
decreased leaving only O-2p to dominate the valence band at Ueff = 6 eV which
can be associated with the enhancement of mh∗.

Further Fig. 5(a) demonstrates that at Ueff = 0 eV, the conduction band
(E-EF$>$ 0 eV) of down-spin structure was dominated mostly by Fe-3d orbital
with a negligible contribution from O-2p. As can be noticed in Fig. 5(b), (c)
and ESI Fig. S2, the characteristics of the conduction band did not change
much due to the effect of Hubbard Ueff up to 5 eV. Only a gradual shifting of
conduction band to a higher energy can be noticed with increasing Ueff as
expected Shenton  _et al._ (2017). However, when we employed a Ueff of 6 eV
(Fig. 5(d)), a new flat conduction band appeared at around 1.25 eV owing to
the sole contribution of Gd-4f state. Furthermore, in the case of the down-
spin valence band (E-EF$<$ 0 eV), a significant influence of Ueff can be
observed. Without considering the Coulomb repulsion effect (i.e. for Ueff = 0
eV in Fig. 5(a)), we obtained the valence band near the Fermi level which is
attributed to the hybridization of Fe-3d and O-2p orbitals. After applying
Ueff = 1 eV (Fig. 5(b)), this band disappeared followed by the emergence of a
new band at around -1.5 eV which shifted gradually to higher energy with
increasing Ueff of up to 6 eV (Fig. 5(b)-(d) and ESI Fig. S2). Notably, the
formation of this band has the contribution from only O-2p orbital.

####  III.3.3 Mulliken population analysis

Table 2:  Mulliken effective charges of individual atoms, bond populations and bond lengths of Gd2FeCrO6 for different values of Ueff obtained via Mulliken population analysis. | Ueff = 0 eV | Ueff =1 eV | Ueff = 2 eV | Ueff = 3 eV | Ueff = 4 eV | Ueff = 5 eV | Ueff = 6 eV
---|---|---|---|---|---|---|---
Atom | Mulliken effective charge (e) |
Gd | 1.52 | 1.49 | 1.49 | 1.49 | 1.49 | 1.50 | 1.50
Fe | 0.59 | 0.78 | 0.79 | 0.80 | 0.82 | 0.82 | 0.84
Cr | 0.53 | 0.58 | 0.59 | 0.61 | 0.62 | 0.63 | 0.65
O | -0.69 | -0.72 | -0.73 | -0.73 | -0.74 | -0.74 | -0.75
Bond | Bond population |
Gd-O | 0.1450 | 0.1587 | 0.1600 | 0.1625 | 0.1637 | 0.1662 | 0.1700
Fe-O | 0.3200 | 0.3100 | 0.3166 | 0.3166 | 0.3166 | 0.3166 | 0.3066
Cr-O | 0.3533 | 0.3200 | 0.3200 | 0.3200 | 0.3200 | 0.3133 | 0.3100
O-O | -0.0350 | -0.0333 | -0.0333 | -0.0333 | -0.0333 | -0.0333 | -0.0333
Bond | Bond length (Å) |
Gd-O | 2.475 | 2.501 | 2.504 | 2.505 | 2.507 | 2.507 | 2.524
Fe-O | 2.023 | 2.051 | 2.061 | 2.065 | 2.065 | 2.069 | 2.083
Cr-O | 2.026 | 2.036 | 2.042 | 2.049 | 2.050 | 2.057 | 2.079

The effective atomic charge, bond population and bond length in a crystalline
solid can be obtained from Mulliken population analysis which provides insight
into the distribution of electrons among different parts of the atomic bonds,
covalency of bonding as well as bond strength Mulliken (1955); Segall  _et
al._ (1996). Mulliken effective charge, Q($\alpha$) of a particular atom
$\alpha$ can be calculated using the following expression Mulliken (1955).

$Q(\alpha)=\sum_{k}\omega_{k}\sum\sum_{\mu}^{on\;\alpha}\sum_{v}P_{\mu\nu}(k)S_{\mu\nu}(k)$
(2)

where Pμν denotes an element of the density matrix and Sμν refers to the
overlap matrix. And the bond population, P($\alpha\beta$) between two atoms
$\alpha$ and $\beta$ can be expressed as Mulliken (1955)-

$P(\alpha\beta)=\sum_{k}\omega_{k}\sum\sum_{\mu}^{on\;\alpha}\sum_{v}^{on\;\beta}2P_{\mu\nu}(k)S_{\mu\nu}(k)$
(3)

Table II provides the calculated Mulliken effective charges of individual
atoms, bond populations and bond lengths between different atoms in GFCO
crystal structure. Noticeably, for all values of Ueff, the Mulliken effective
charges of the individual Gd, Fe, Cr and O atoms are found to be reasonably
smaller than their formal ionic charges which are +3, +3, +3, and -2,
respectively. Such difference between Mulliken effective and formal ionic
charges is an indication of the existence of mixed ionic and covalent bonding
in GFCO Segall  _et al._ (1996). It should be noted that small Mulliken
effective charge of an atom is associated with its high level of covalency and
vice versa Mulliken (1955); Rozilah  _et al._ (2017). Therefore it can be
inferred that GFCO double perovskite includes chemical bonding with prominent
covalency. Further, in Table II, an enhancement can be observed in the
effective charges of Fe, Cr and O atoms after employing Ueff in the first-
principles calculation (Ueff =1 to 6 eV). This outcome indicates that the
degree of bond covalency in GFCO reduced to an extent due to the effect of on-
site Coulomb interaction.

Table II also presents the bond populations and bond lengths of Gd–O, Fe–O,
Cr–O and O-O bonds in GFCO double perovskite as obtained for Ueff = 0 to 6 eV.
It is noteworthy that a large positive value of bond population is associated
with high degree of covalency whereas a small bond population indicates high
degree of ionicity in the chemical bond Ching and Xu (1999). In the present
investigation, the bond populations of the Gd–O, Fe–O and Cr–O were determined
to be positive whereas the bond population of O-O was negative for all values
of Ueff. This result suggests that no bonds had formed between the O atoms in
GFCO double perovskite Chakraborty  _et al._ (2008); Yaakob  _et al._ (2015).
Moreover, the calculated bond populations of Fe-O and Cr-O are found to be
considerably larger than that of Gd-O. Such observation implies that the Fe-O
and Cr-O bonds possess higher degree of covalency as compared to Gd-O bonds.
It is also worth noticing that the bond populations of Fe-O and Cr-O
calculated by GGA+U method (Ueff =1 to 6 eV) are lower than the values
obtained by GGA method (Ueff = 0 eV). This result provides further evidence
for the influence of Coulomb repulsion to reduce the covalency of Fe-O and
Cr-O bonding in GFCO as was observed before by analyzing the calculated
Mulliken effective charges. Furthermore, Table II demonstrates that the bond
lengths of Fe-O and Cr-O are reasonably smaller than that of Gd-O which can be
attributed to their larger bond population and consequently, higher degree of
covalency as compared to Gd-O bonding.

####  III.3.4 Electron charge density

Figure 6:  Electronic charge density along z-axis of Gd2FeCrO6 for (a) Ueff =
0 eV, (b) Ueff = 1 eV, (c) Ueff = 3 eV, and (d) Ueff = 6 eV.

For further understanding of the chemical bonding nature, we have determined
the electron charge density distribution of GFCO double perovskite along the
z-axis by GGA (Ueff = 0 eV) and GGA+U (Ueff = 1, 3 and 6 eV) methods as
illustrated in Fig. 6. It is worth noting that a typical covalent bond between
two atoms involves overlapping of electron clouds from both of them and the
electrons remain concentrated in the overlapping region Phillips (1968). In
Fig. 6, for all values of Ueff, a larger overlap of electron cloud can be
observed between Fe/Cr and O atoms as compared to Gd and O atoms. Such
observation implies that the covalent bonds between Fe/Cr and O in GFCO are
considerably stronger than the bond between Gd and O as was also revealed by
Mulliken population analysis.

Moreover, it can be clearly seen that no electron clouds are concentrated in
the area between one of the two Gd atoms and O which implies the formation of
an ionic bond between these two atoms Craig  _et al._ (1954); Rozilah  _et
al._ (2017). Notably, the charge sharing between Fe/Cr-O can be attributed to
the hybridization among Fe/Cr-3d and O-2p orbitals and the Gd-O bond formation
is associated with the hybridization between Gd-4f and O-2p orbitals which was
demonstrated by the DOS curves (Fig. 5).

Furthermore, if we meticulously observe Fig. 6(a) and (b), it can be noticed
that the overlapping of electron clouds between Fe/Cr and O atoms got
reasonably narrower after considering the effect of on-site Coulomb
interaction (Ueff) in first-principles calculation. Afterward, the electron
charge density in the overlapped region between Fe/Cr and O atoms enhanced
with increasing Ueff (see Fig. 6(c) and (d)) which is consistent with the
outcome of Mulliken population analysis.

Figure 7:  (a) Variation of theoretically obtained absorption coefficient of
GFCO perovskite as a function of wavelength for different Ueff. (b) The
absorption coefficient vs. Ueff for some fixed values of the wavelength.

###  III.4 Light absorption property

The absorption coefficient provides valuable information about a material’s
light-harvesting ability. The optical absorption coefficient of GFCO has been
evaluated by first-principles calculation via GGA (Ueff = 0 eV) and GGA+U
(Ueff = 1 to 6 eV) approaches using a polarization technique which includes
the electric field vector as an isotropic average in all directions Brik  _et
al._ (2014). To gain additional distinguishable absorption peaks, a small
smearing value of 0.5 eV was used. Fig. 7(a) illustrates the calculated
absorption coefficients of GFCO double perovskite as a function of wavelength
to demonstrate the effect of Hubbard Ueff parameter on its light absorption
property. In Fig. 7(a), for all values of Ueff, two absorption peaks can be
clearly observed in the UV region which indicates the strong UV light
absorption capacity of GFCO. Noticeably, the stronger absorption coefficient
peak lies at around 120 nm and it undergoes a slight red-shift with increasing
Ueff. On contrary, the weaker absorption peak at around 320 nm is slightly
blue-shifted for higher values of Ueff. For comparison, we have also provided
the absorbance spectrum of GFCO obtained via UV-visible spectroscopy in ESI
Fig. S3 Bhuyan  _et al._ (2021). It is worth noticing that the experimentally
obtained spectrum has two additional bands in the visible region along with
the two bands we obtained theoretically in the UV regime. This might be
related to the fact that the DFT calculations were performed for 0 K whereas
the experiment was conducted at room temperature. Such difference in
conditions can be attributed to the discrepancy between the optical and
theoretical absorbance spectra of GFCO Honglin  _et al._ (2014); Harun  _et
al._ (2020).

Fig. 7(b) shows the absorption coefficient vs. Ueff curves of GFCO for some
fixed values of wavelength ranging from 150 to 600 nm. We can observe that for
the wavelength of 150 nm, the absorption coefficient attains its minimum and
maximum values for Ueff = 4 eV and Ueff = 5 eV, respectively suggesting the
underestimation and overestimation of the coefficient at these Ueff values.
Similarly, it can be found that except at Ueff = 3 eV, the absorption
coefficient is either overestimated or underestimated at all other values of
Ueff for any of the specific wavelengths such as 200, 250, 300 nm etc. This
intriguing observation is highlighted by the rectangle in Fig. 7(b).

###  III.5 Comparison of experimental and theoretical optical band gaps

Further, we have employed the calculated absorption coefficients to
theoretically calculate the optical band gap values of GFCO double perovskite
using Tauc relation Tauc  _et al._ (1966). Fig. 8 shows the variation in
theoretically calculated direct band gap values as a function of Ueff.
Noticeably, a direct band gap value of 0.5 eV was obtained by GGA method (Ueff
= 0 eV) which is significantly smaller than the experimental value $\sim$2.0
eV. Further, with the increase of Ueff up to 5 eV, an almost linear increase
can be observed in the direct optical band gap values of GFCO double
perovskite. However, an anomalous decrease can be noticed for a further
increase of Ueff to 6 eV. It is intriguing to note that the direct optical
band gap ($\sim$1.99 eV) obtained for Ueff = 3 eV matches well with the
experimentally obtained one ($\sim$2.0 eV) which is marked by a circle in Fig.
8. Also, for Ueff = 6 eV, we found the calculated band gap value ($\sim$2.08
eV) to be quite close to the experimental one. But as we had demonstrated
before, the character and curvature of the conduction and valence bands
inaccurately changed for Ueff $>$ 5 eV. Hence, fitting Ueff to the band gap
alone may provide erroneous results for such double perovskite materials
especially in the cases where the role of the band edge character is crucial.

Finally, it is fascinating to note that our predictions for all the properties
of GFCO (i.e. structural, electronic and optical) considering Ueff = 3 eV
closely matched with the experimental results without any over or
underestimation of band gap values. This implies that the Ueff value of 3 eV
most accurately localized the Fe-3d and Cr-3d orbitals of GFCO. Moreover, the
almost accurate estimation of band gap values suggests that the effect of
self-interaction error from the other orbitals of GFCO was almost negligible
Brown and Page (2020); Perdew (1985). Therefore we recommend Ueff = 3 eV for
the GGA+U calculation of electronic and optical properties of GFCO double
perovskite.

Figure 8:  Variation in theoretically calculated direct optical band gap as a
function of Ueff. The red circle represents the experimentally obtained
optical band gap value.

##  IV Conclusions

We have employed the first-principles based GGA and GGA+U methods to calculate
the spin-polarized electronic band structure, Mulliken bond population,
electron charge density distribution and optical characteristics of our newly
synthesized direct band-gap semiconductor Gd2FeCrO6 (GFCO) double perovskite
for a range of Ueff between 0 and 6 eV, applied to the Fe-3d and Cr-3d
orbitals. The structural parameters of the monoclinic GFCO crystal varied
nominally with Ueff. To the contrary, the variation of Ueff demonstrated
significant effect on the electronic band structure which indicates the
importance of employing reasonable value of Ueff to correct the over-
delocalization of the Fe/Cr-3d states. For Ueff $>$ 5 eV, the computed partial
density of states for Gd-4f orbital radically altered which had significantly
changed the band structure. In particular, we observed that the character and
curvature of the conduction and valence bands largely varied for Ueff $>$ 5 eV
leading to enormous changes in calculated charge carrier effective masses.
Notably, in the case of Ueff = 3 eV, the theoretically calculated direct
optical band gap $\sim$1.99 eV matched well with the experimental value
$\sim$2.0 eV. These findings justify that it might be worthwhile to employ
Ueff = 3 eV to accurately calculate the structural, electronic and optical
properties of GFCO double perovskite. The outcome of this investigation might
be useful for a keen understating of the electronic and optical properties of
this newly synthesized double perovskite and related material systems for
photocatalytic applications via band gap engineering. This study also reveals
the importance of conducting systematic analysis of the influence of on-site
Coulomb interaction on band gaps as well as on the electronic structure as a
whole for related other double perovskite materials for which experimental
data are still not available.

## Acknowledgements

The computational facility provided by the Institute of Information and
Communication Technology (IICT), Bangladesh University of Engineering and
Technology (BUET) is sincerely acknowledged. The authors would also like to
acknowledge the Committee for Advanced Studies and Research (CASR), BUET.

## Data availability

The raw and processed data required to reproduce these findings cannot be
shared at this time due to technical or time limitations.

## Supplementary Information

Additional electronic supplementary information (ESI) is available. See DOI:
10.1016/j.jmrt.2021.06.026

## References

  * Vasala and Karppinen (2015) S. Vasala and M. Karppinen, Progress in solid state chemistry 43, 1 (2015).
  * Gray  _et al._ (2010) B. Gray, H. N. Lee, J. Liu, J. Chakhalian, and J. Freeland, Applied Physics Letters 97, 013105 (2010).
  * Gaikwad  _et al._ (2019) V. M. Gaikwad, M. Brahma, R. Borah, and S. Ravi, Journal of Solid State Chemistry 278, 120903 (2019).
  * Das  _et al._ (2008) H. Das, U. V. Waghmare, T. Saha-Dasgupta, and D. Sarma, Physical review letters 100, 186402 (2008).
  * Cerón  _et al._ (2019) J. G. Cerón, J. A. Rodríguez, A. Rosales-Rivera, J. C. Farfán, J. C. Vasquez, D. L. Téllez, J. Roa-Rojas, _et al._ , Journal of Materials Research and Technology 8, 3978 (2019).
  * Yin  _et al._ (2019) W.-J. Yin, B. Weng, J. Ge, Q. Sun, Z. Li, and Y. Yan, Energy & Environmental Science 12, 442 (2019).
  * Lin  _et al._ (2021) C. Lin, Y. Zhao, Y. Liu, W. Zhang, C. Shao, and Z. Yang, Journal of Materials Research and Technology 11, 1645 (2021).
  * Mohassel  _et al._ (2020) R. Mohassel, M. Amiri, A. K. Abbas, A. Sobhani, M. Ashrafi, H. Moayedi, and M. Salavati-Niasari, Journal of Materials Research and Technology 9, 1720 (2020).
  * Kangsabanik and Alam (2020) J. Kangsabanik and A. Alam, The Journal of Physical Chemistry Letters 11, 5148 (2020).
  * Wu  _et al._ (2020) H. Wu, Z. Pei, W. Xia, Y. Lu, K. Leng, and X. Zhu, Journal of Alloys and Compounds 819, 153007 (2020).
  * Nechache  _et al._ (2015) R. Nechache, C. Harnagea, S. Li, L. Cardenas, W. Huang, J. Chakrabartty, and F. Rosei, Nature Photonics 9, 61 (2015).
  * Nair  _et al._ (2014) H. S. Nair, R. Pradheesh, Y. Xiao, D. Cherian, S. Elizabeth, T. Hansen, T. Chatterji, and T. Brückel, Journal of applied physics 116, 123907 (2014).
  * Bhuyan  _et al._ (2021) M. Bhuyan, S. Das, and M. Basith, Journal of Alloys and Compounds 878, 160389 (2021).
  * Zhang  _et al._ (2016) R. Zhang, J. Zhao, Y. Yang, Z. Lu, and W. Shi, Computational Condensed Matter 6, 5 (2016).
  * Shenton  _et al._ (2017) J. K. Shenton, D. R. Bowler, and W. L. Cheah, Journal of Physics: Condensed Matter 29, 445501 (2017).
  * Himmetoglu  _et al._ (2014) B. Himmetoglu, A. Floris, S. De Gironcoli, and M. Cococcioni, International Journal of Quantum Chemistry 114, 14 (2014).
  * Terakura  _et al._ (1984) K. Terakura, T. Oguchi, A. Williams, and J. Kübler, Physical Review B 30, 4734 (1984).
  * Brown and Page (2020) J. J. Brown and A. J. Page, The Journal of Chemical Physics 153, 224116 (2020).
  * Lu and Liu (2014) D. Lu and P. Liu, The Journal of chemical physics 140, 084101 (2014).
  * Wang  _et al._ (2006) L. Wang, T. Maxisch, and G. Ceder, Physical Review B 73, 195107 (2006).
  * Zhang  _et al._ (2010) S. Zhang, O. Y. Kontsevoi, A. J. Freeman, and G. B. Olson, Physical Review B 82, 224107 (2010).
  * Dudarev  _et al._ (1998) S. Dudarev, G. Botton, S. Savrasov, C. Humphreys, and A. Sutton, Physical Review B 57, 1505 (1998).
  * Perdew (1985) J. P. Perdew, International Journal of Quantum Chemistry 28, 497 (1985).
  * Rodriguez-Carvajal (1990) J. Rodriguez-Carvajal, in  _FULLPROF: a program for Rietveld refinement and pattern matching analysis, Abstracts of the Satellite meeting on powder diffraction of the XV congress of the IUCr, Toulouse, France_ , Vol. 127 (1990).
  * Segall  _et al._ (2002) M. Segall, P. J. Lindan, M. a. Probert, C. J. Pickard, P. J. Hasnip, S. Clark, and M. Payne, Journal of physics: condensed matter 14, 2717 (2002).
  * Rozilah  _et al._ (2017) R. Rozilah, M. Yaakob, Z. Mohamed, and A. Yahya, Materials Research Express 4, 066103 (2017).
  * Fischer and Almlof (1992) T. H. Fischer and J. Almlof, The Journal of Physical Chemistry 96, 9768 (1992).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Physical review B 13, 5188 (1976).
  * Chakraborty  _et al._ (2008) M. Chakraborty, P. Pal, and B. R. Sekhar, Solid state communications 145, 197 (2008).
  * Saha  _et al._ (2000) S. Saha, T. Sinha, and A. Mookerjee, Physical Review B 62, 8828 (2000).
  * Zhang  _et al._ (2008) X. Zhang, M. Guo, W. Li, and C. Liu, Journal of Applied Physics 103, 063721 (2008).
  * Mann  _et al._ (2016) G. W. Mann, K. Lee, M. Cococcioni, B. Smit, and J. B. Neaton, The Journal of chemical physics 144, 174104 (2016).
  * Neaton  _et al._ (2005) J. Neaton, C. Ederer, U. Waghmare, N. Spaldin, and K. Rabe, Physical Review B 71, 014113 (2005).
  * Tauc  _et al._ (1966) J. Tauc, R. Grigorovici, and A. Vancu, physica status solidi (b) 15, 627 (1966).
  * Hou  _et al._ (2014) Y. Hou, H. Xiang, and X. Gong, Physical Review B 89, 064415 (2014).
  * Bredas (2014) J.-L. Bredas, Materials Horizons 1, 17 (2014).
  * Reunchan  _et al._ (2016) P. Reunchan, A. Boonchun, and N. Umezawa, Physical Chemistry Chemical Physics 18, 23407 (2016).
  * Das  _et al._ (2019) S. Das, M. Bhowal, and S. Dhar, Journal of Applied Physics 125, 075705 (2019).
  * Mulliken (1955) R. Mulliken, The Journal of Chemical Physics 23, 1841 (1955).
  * Segall  _et al._ (1996) M. Segall, R. Shah, C. J. Pickard, and M. Payne, Physical Review B 54, 16317 (1996).
  * Ching and Xu (1999) W. Ching and Y.-N. Xu, Physical Review B 59, 12815 (1999).
  * Yaakob  _et al._ (2015) M. Yaakob, M. Taib, L. Lu, O. Hassan, and M. Yahya, Materials Research Express 2, 116101 (2015).
  * Phillips (1968) J. Phillips, Physical Review 166, 832 (1968).
  * Craig  _et al._ (1954) D. P. Craig, A. Maccoll, R. Nyholm, L. Orgel, and L. Sutton, Journal of the Chemical Society (Resumed) , 332 (1954).
  * Brik  _et al._ (2014) M. Brik, O. Parasyuk, G. Myronchuk, and I. Kityk, Materials Chemistry and Physics 147, 155 (2014).
  * Honglin  _et al._ (2014) L. Honglin, L. Yingbo, L. Jinzhu, and Y. Ke, Journal of alloys and compounds 617, 102 (2014).
  * Harun  _et al._ (2020) K. Harun, N. A. Salleh, B. Deghfel, M. K. Yaakob, and A. A. Mohamad, Results in Physics 16, 102829 (2020).


]Presented XXXXX; received XXXXX; accepted XXXXX; published online XXXXX

# Magnetic Reconnection on a Klein Bottle

Luke Xia<EMAIL_ADDRESS>University of California, Irvine 92697, USA  M.
Swisdak<EMAIL_ADDRESS>Institute for Research in Electronics and Applied
Physics, University of Maryland, College Park, Maryland 20742, USA

([)

###### Abstract

We present a new boundary condition for simulations of magnetic reconnection
based on the topology of a Klein bottle. When applicable, the new condition is
computationally cheaper than fully periodic boundary conditions, reconnects
more flux than systems with conducting boundaries, and does not require
assumptions about regions external to the simulation as is necessary for open
boundaries. The new condition reproduces the expected features of
reconnection, but cannot be straightforwardly applied in systems with
asymmetric upstream plasmas.

##  I Introduction

During reconnection oppositely directed components of a magnetic field come
together at an X-point and undergo a topological reordering that enables the
transfer of magnetic energy to the surrounding plasma. In the collisionless
plasmas of space and astrophysics, reconnection has particular importance due
to its role in catalyzing energy release over a wide range of scales in
flares, jets, and other phenomena.

Unsurprisingly then, numerical simulations of reconnection have been a subject
of significant interest. Any simulation entails certain fundamental choices,
among them the boundary conditions imposed on the computational domain.
Idealized simulations – i.e., those not modeling a particular event – often
begin in a configuration based on the Harris sheet [1], a steady-state
solution of the Vlasov equation in which the magnetic field has the simple
profile

$\mathbf{B}=B_{0}\tanh(y/w_{0})\,\mathbf{\hat{x}},$ (1)

with $w_{0}$ a constant. (For definiteness, we assume the simulation domain
occupies the space $-L_{x}/2\leq x\leq L_{x}/2$, $-L_{y}/2\leq y\leq L_{y}/2$.
A fully three-dimensional domain is not necessary for reconnection to occur
and so the $z$ coordinate will be suppressed.) The reversal of the field
across $y=0$ implies the existence of a current sheet, with the current
pointing in the $-z$ direction.

The form of the magnetic field in the Harris sheet imposes constraints on the
boundary conditions. Periodic boundaries are typically used in the horizontal
($x$) direction, but the anti-symmetric nature of $\mathbf{B}$,
$B_{x}(y=-L_{y}/2)=-B_{x}(y=L_{y}/2)$, precludes the use of periodic
boundaries in the vertical ($y$) direction. Two options are usually
considered. The first is the inclusion of hard conducting walls at the top and
bottom boundaries of the domain. Appropriate boundary conditions are imposed
on the electromagnetic fields (e.g., Dirichlet for $E_{x}$, $B_{y}$, and
$E_{z}$ and von Neumann for the other components) while particles specularly
reflect from the walls. The influential GEM Challenge [2] used these
conditions and so hereinafter we will use its name to denote them. For the GEM
Challenge boundary conditions the computational domain is topologically
equivalent (homeomorphic) to a cylinder.

The second common option is to add a second current sheet of opposite
orientation so that

$\mathbf{B}=B_{0}\left[\tanh\left(\frac{y-L_{y}/4}{w_{0}}\right)-\tanh\left(\frac{y+L_{y}/4}{w_{0}}\right)+1\right]\,\mathbf{\hat{x}}$
(2)

This choice allows periodic boundary conditions in the vertical direction
because $B_{x}(y=-L_{y}/2)=B_{x}(y=L_{y}/2)$. (For both choices, one usually
also requires $w_{0}\ll L_{y}$ so that $\partial B_{x}/\partial y$, and hence
the current, is negligible for $y=\pm L_{y}/2$.) The implementation of fully
periodic boundaries is typically straightforward since no special conditions
need to be placed on the electromagnetic fields or particles at any point in
the grid. In this case, the computational domain is homeomorphic to a torus.

Either choice comes with limitations. Simulations using the GEM Challenge
boundary conditions will stagnate when the reconnected flux accumulating in
the magnetic islands presses against the hard boundaries (see Figure 3 below).
Fully periodic boundary conditions employing the field described in equation 2
can mitigate this stagnation by offsetting the X-points on the two current
layers by $L_{x}/2$, thus allowing the flux accumulating in the island of one
layer to drive inflow towards the X-line on the other. However, keeping the
amount of flux available to reconnect the same for both types of boundary
conditions (equivalently, keeping the aspect ratios of the current sheets the
same) requires a doubling of the size of the computational domain, and hence
the computational cost, in the fully periodic case. Less-common alternatives,
e.g., open boundaries or driven reconnection [3, 4, 5, 6], increase the
computational complexity and can only partially mitigate these issues (often
while raising new ones).

Here we present a novel type of periodic boundary condition in which the
computational domain is homeomorphic to a Klein bottle. Making this choice
allows the simulation to enjoy the benefits of the two major options discussed
above: the computational expense of the GEM Challenge conditions coupled with
the self-driving of the fully periodic domain. We show that these boundary
conditions preserve the important properties of magnetic reconnection. Section
II discusses our computational methods, including the implementation of the
new boundary condition. Results from the simulations are presented in Section
III, while Section IV offers our conclusions.

##  II Computational Methods

###  II.1 Code

The simulations are performed with p3d, a massively parallel particle-in-cell
(PIC) code with a long history of simulating magnetic reconnection [7]. It
employs units in which a reference density $n_{0}$ and magnetic field strength
$B_{0}$ define the units of length, the proton inertial length
$d_{p0}=c/\omega_{p0}$ (where $\omega_{p0}=\sqrt{4\pi n_{0}e^{2}/m_{p}}$ is
the proton plasma frequency), and time, the proton cyclotron time
$\Omega_{p0}^{-1}=m_{p}c/eB_{0}$. Velocities are measured in units of a
reference Alfvén speed $v_{A0}=\sqrt{B_{0}^{2}/4\pi m_{p}n_{0}}$ while
electric fields and temperatures are normalized to $v_{A0}B_{0}/c$ and
$m_{p}v_{A0}^{2}$, respectively. The simulations presented here follow
particles in a 2x3v (also known as 2.5D) phase space in which no variations
are allowed in the third spatial dimension (i.e., $\partial/\partial z=0$).

With the exception of a slight alteration in the size of the computational
domain, the parameters are chosen to match those of the GEM Challenge [2]. In
order to reduce the separation between the spatial scales of the protons and
the electrons, and hence the computational expense, the electron mass is taken
to be $m_{e}/m_{i}=0.04$, which implies the electron inertial length
$d_{e}=0.2d_{p}$. For similar computational reasons, the ratio of the speed of
light to the Alfvén speed (equivalently, the ratio of the proton plasma and
cyclotron frequencies) is taken to be 20. The initial electron and proton
temperatures are $1/12$ and $5/12$, respectively. The timestep is chosen to
satisfy the CFL condition and, as a consequence, particles can travel at most
one grid cell per timestep.

For later reference, it is useful to discuss the layout of the computational
domain in more detail than usual. To parallelize a computation, p3d breaks the
domain into subdomains in configuration space, each of which is assigned to a
single processor. The processors are organized into a rectangular grid of size
$P=P_{x}\times P_{y}\times P_{z}$ and each processor is further subdivided
into $N_{x}\times N_{y}\times N_{z}$ grid points. For notational simplicity we
will assume $N_{z}=P_{z}=1$ in what follows, but the extension to multiple
processors and gridpoints in the third dimension is straightforward. Each
processor has an overall label that can range from $0$ to $P-1$ as well as a
label in each dimension: $p_{x}\in\\{0,1,2,...,P_{x}-1\\}$,
$p_{y}\in\\{0,1,2,...,P_{y}-1\\}$, and $p_{z}\in\\{0,1,2,...,P_{z}-1\\}$. The
overall label and the triplet $(p_{x},p_{y},p_{z})$ are unique to each
processor and it is straightforward, given the dimensions of the computational
domain, to convert between the two. For example, the processor at the bottom-
left of the domain, $P=0$, has triplet $(0,0,0)$. The processor to its
immediate right, $P=1$, has triplet $(1,0,0)$ and the processor immediately
above is $P=P_{x}$ with triplet $(0,1,0)$. For communication purposes, every
processor needs to know the identities of the other processors bordering it,
and while this calculation is simple for those on the interior of the
computational domain, it can be more complicated for those on the boundaries.

###  II.2 Klein Bottles and Processor Mapping

Figure 1: Immersion of a two-dimensional Klein bottle in three-dimensional
space.

The novel boundary condition presented here is based on the topology of a
Klein bottle, a two-dimensional non-orientable manifold. (The perhaps-more-
familiar Möbius strip is another example of a non-orientable manifold.) A
Klein bottle is a surface with no distinct front or back side that can be
obtained by identifying and gluing together specific edges of a rectangle in a
non-trivial way, resulting in a structure that cannot be embedded smoothly in
three-dimensional space without self-intersections or singularities (see
Figure 1 for a representation). Higher-dimensional analogues of Klein bottles
can be similarly constructed. For example, mapping the edges of a rectangular
prism in a particular way produces a system homeomorphic to a Klein geometry
in $\mathbb{R}^{4}$.

Figure 2: Topological mappings of (i) a torus (ii) a Klein bottle (iii) a
phase-shifted Klein bottle. The horizontal boundaries are periodic. Squares
represent processors and those of the same color abut in the vertical
direction. The arrow demonstrates how a vector transforms across the boundary.

In numerical simulations, the identification and gluing-together process is
equivalent to describing how processors abut one another when crossing a
boundary. This process is given a pictorial representation in Figure 2. In
this visualization, small squares represent individual processors; squares
painted the same color neighbor one another in the vertical direction. (Only
the top and bottom boundaries are of interest, as the horizontal direction is
assumed to be periodic in each panel.) In the configuration shown, with
$P=P_{x}\times P_{y}\times P_{z}=16\times 8\times 1=128$, the processors are
numbered from $0$ to $127$, with those on the bottom row having
$P\in\\{0,1,...15\\}$. The fact that the doubly periodic boundary conditions
described in Section I are topologically equivalent to a torus is demonstrated
in the first panel. Travelling downward from the bottom row brings one to the
top row, where the processors have $P\in\\{112,113,...127\\}$. In other words,
processors next to the boundary (on either side) with triplets of general form
$(i,j,0)$ border processors with triplets of the form $(i,P_{y}-(j+1),0)$. The
overplotted arrow illustrates how a vector transforms when moving between the
top and bottom rows of processors. Note that interpreting the arrow as a
magnetic field vector makes it clear why a system with a single current sheet
– across which the field must switch direction – cannot use these boundary
conditions.

To map a two-dimensional sheet to a Klein bottle, one connects the left and
right sides together as usual. However, processors on the top and bottom
boundaries now adjoin counterparts at mirrored positions on the $x$-axis, as
shown in Figure 2(ii). A processor on the top or bottom boundary with triplet
of the form $(i,j,0)$ borders a processor with triplet
$(P_{x}-(i+1),P_{y}-(j+1),0)$. (Note the set of neighboring gridpoints within
processors are also mirrored, so that the leftmost gridpoint on, for example,
the bottom row of processor 1 is considered to be directly above the rightmost
gridpoint in the top row of processor 126.) As a consequence, vectors flip
orientation across the boundary, as shown by the overplotted arrow.
Implementation of this boundary condition requires several algorithmic
changes. Processors on the top and bottom rows have a new set of neighboring
processors and the $x$ and $z$ components of any vector quantity, e.g.,
particle velocities and electromagnetic fields, must be flipped across the
boundary. (In practice, processors in p3d interact with neighbors through
guard cells. Any changes to the signs of gridded quantities required by the
boundary conditions occurs when these cells are updated.) This means that only
one current sheet needs to be included in the domain, as the reversal in
$B_{x}$ is naturally accommodated.

However, the basic Klein boundary condition shown in Figure 2(ii) is not ideal
for reconnection simulations. The reason becomes clear after consideration of
the magnetic configuration at late times. Reconnection proceeding at the
X-line (assumed to be at the center of the domain), produces an island of
magnetic flux centered on the left/right (they are equivalent) boundary. This
island grows, eventually becoming large enough to approach the top and bottom
boundaries. In runs employing the GEM Challenge boundary condition, this
island presses against the conducting walls, producing a back-pressure that
resists further reconnection and eventually stifles the simulation. The same
issue bedevils the Klein boundary condition, although in that case the back-
pressure is not provided by conducting walls but by the island itself – growth
of the island in, for instance, the bottom-left of the domain, presses against
the island in the upper-right.

Fortunately, a straightforward adjustment that mimics the behavior of doubly
periodic boundary conditions can alleviate this issue by, in essence, creating
self-driving reconnection. Arranging for the island growing on the left and
right sides of the domain to emerge from the top and bottom boundaries in the
middle of the domain extends the time that reconnection will occur. The phase
shift is produced by displacing the set of neighboring processors for either
the top or bottom boundary by $L_{x}/2$, forming the configuration displayed
in Figure 2(iii). Returning to the example of a $16\times 8$-processor domain,
each processor with $P\in\\{0,1,\dots 7\\}$ now neighbors the corresponding
processor $P\in\\{119,118,\dots 112\\}$, and vice versa, while each processor
$P\in\\{8,9,\dots 15\\}$ abuts the corresponding processor
$P\in\\{127,126,\dots,120\\}$. More generally, a processor on the border with
triplet $(i,j,0)$ directly borders a processor with triplet
$([P_{x}/2-(i+1)]\,\bmod P_{x},P_{y}-(j+1),0)$. As with the unshifted Klein
boundaries, the set of neighboring gridpoints are again mirrored across the
boundary so, for example, the leftmost gridpoint in the bottom row of
processor 1 is considered to be directly above the rightmost gridpoint of
processor 118.

##  III Simulation Results

We now present results from simulations using the various types of boundary
conditions discussed above. The computational domain measures $L_{x}\times
L_{y}=51.2d_{i}\times 12.8d_{i}$ for the GEM Challenge and shifted Klein cases
and is doubled in $L_{y}$ for the doubly periodic case. (In other words, each
reconnecting current sheet has the same aspect ratio in all runs.) The
reconnecting component of the magnetic field has the initial form shown in
either equation 1 or 2 while the density has the asymptotic value of
$0.2n_{0}$ and varies as necessary to ensure force balance in the $y$
direction. Section III.1 considers the case where the asymptotic fields are
anti-parallel, as in the original GEM Challenge, while Section III.2 adds a
$B_{z}$ component that reduces the shear angle between the asymptotic fields.

###  III.1 Anti-Parallel Reconnection

Figure 3 shows a comparison between runs with the GEM Challenge and the
shifted Klein boundary conditions with the three panels in each column
depicting the simulations at $t\Omega_{cp}=20$, $40$, and $60$, respectively.
In each, the out-of-plane current density is shown in color and is overplotted
with magnetic field lines. For the GEM Challenge boundaries (left column), the
islands of reconnected magnetic flux grow freely at first (top panel) until
they begin to interact with the walls at the top and bottom of the domain
(middle panel). The increasing pressure in the island exerts a growing force
on plasma ejected from the X-point and suppresses further reconnection (bottom
panel).

Figure 3: Out-of-plane electron current density and magnetic field lines at
$t\Omega_{cp}=20$, $40$, and $60$ for the GEM Challenge boundary conditions
(left) and the Klein boundary conditions (right).

In contrast, the panels in the right column of Figure 3 show results from an
otherwise-identical simulation that employs the shifted Klein boundary
conditions described above. By the time of the middle panel the effects of the
boundary conditions are already noticeable in the magnetic field lines since,
unlike in the GEM Challenge simulation, they are free to pass through the top
and bottom boundaries. By late time (bottom panel), the islands of reconnected
flux have grown substantially and, due to the shift at the boundaries,
continue to drive reconnection at the X-line.

The lack of hard walls at the top and bottom boundaries makes the Klein
boundary conditions similar to the doubly-periodic case. Figure 4 presents a
comparison of several quantities at $t\Omega_{cp}=60$, a time of robust
reconnection, from runs with both boundary conditions. The panels from the
Klein boundaries (left column) show the entire domain, while those from the
doubly periodic case (right column) show only one of the two current sheets.
The panels depict, from top to bottom, $J_{ez}$, $v_{ex}$, $v_{ix}$, $B_{z}$,
and $E_{y}$. Each horizontal pair of images uses the same color scale and so
may be directly compared. The two columns are quite similar, suggesting that
the Klein boundary conditions do not significantly affect the evolution of the
system.

Figure 4: Comparison of results with the Klein (left column) and doubly
periodic (right column) boundary conditions at $t\Omega_{cp}=60$. (A)-(B):
Out-of-plane electron current density ($J_{ez}$); (C)-(D) Horizontal electron
velocity ($v_{ex}$); (E)-(F): Horizontal ion velocity ($v_{ix}$); (G)-(H):
$B_{z}$; (I)-(J): $E_{y}$

We next compare all three runs in Figure 5, which plots the reconnected flux,
$\Delta\psi$, and the reconnection rate, $d(\Delta\psi)/dt$, versus time for
the different cases: GEM Challenge in blue dashes, doubly periodic in red dot-
dashes, and shifted Klein in black. Contours of $\psi$ are the magnetic field
lines shown in Figures 3 and 4 and in the first and third cases $\Delta\psi$
is merely the difference in $\psi$ between the X- and O-points. In the doubly
periodic case the curve plots the average of the differences observed in the
two current sheets. Each of the simulations exhibit an initial overshoot where
the reconnection rate spikes to $\approx 0.25$, but in the GEM Challenge case
that spike is followed by a steady decline, with reconnection essentially
ending after $t\Omega_{cp}\approx 50$ when $\Delta\psi$ plateaus. In contrast,
the doubly periodic and shifted Klein runs both settle into a state with
reasonably constant reconnection between $t\Omega_{cp}=40$ and
$t\Omega_{cp}=80$ with $d(\Delta\psi)/dt\approx 0.15$. The large excursions
near $t\Omega_{cp}=60$ correspond to the time when the once-reconnected field
lines surrounding the magnetic islands interact with an X-point a second time,
either that on the opposite current sheet for the doubly periodic case or the
original one for the shifted Klein (see Figure 3). It is interesting to note
that the reconnection rate after these fluctuations roughly remains unchanged.

The elapsed time (measured on a clock) for the three computations differed by
less than $3\%$. However, the GEM Challenge and shifted Klein cases were run
on the same number of processors, while the doubly periodic case, for which
the computational domain was doubled in size, used twice as many. Hence, while
the doubly periodic and shifted Klein cases found very similar physical
results, the former incurred twice the computational expense of the latter.
(For the relatively small simulations considered here, p3d exhibits near-
perfect weak scaling and so any extra overhead incurred by doubling the number
of processors is minimal.)

Figure 5: Reconnected flux (top), $\Delta\psi$, and reconnection rate
(bottom), $d(\Delta\psi)/dt$ versus time for the GEM Challenge (blue dash),
doubly periodic (red dot-dash), and shifted Klein (black) boundary conditions.

###  III.2 Guide Field Reconnection

Adding a spatially constant guide field, $B_{g}$ to the initial conditions of
the GEM Challenge and doubly periodic cases is straightforward. The Klein case
is more subtle. Recall that the $x$ and $z$ components of vectors are flipped
when crossing the Klein boundary. As a consequence, a spatially constant guide
field would exhibit a discontinuous jump – a gridpoint at the top or bottom
boundary with $B_{z}=B_{g}$ would be adjacent to one with $B_{z}=-B_{g}$. Such
a configuration would generate a large current and lead to instability. This
behavior arises because Klein bottles are non-orientable, which means there is
ambiguity in the definition of "out of the page" and "into the page". (This is
perhaps most easily seen with the Möbius strip. If one takes an arrow pointing
normal to the surface and translates it once around the strip to the starting
location, the head will point in the opposite direction.)

To sidestep this obstacle, we can slightly modify the boundary condition.
Write $B_{z}$ as the sum of a constant background $B_{z0}$ and a perturbation:
$B_{z}=B_{z0}+\delta B_{z}$. Then, rather than apply the Klein boundary to the
full magnitude, only change the sign of the perturbation $\delta B_{z}$ when
it crosses the boundary. Note that in the run discussed in Section III.1,
$B_{z0}=0$ so this procedure is identical to modifying the full $B_{z}$. It
might be argued that, for consistency, such a modification should be applied
to any $z$ component (e.g., $E_{z}$ or $J_{z}$), but in the initial states
considered here all other quantities have no uniform component and so making
such an adjustment would have no effect.

We now compare results from three simulations with $B_{g0}=1$. Figure 6 shows
the evolution of the GEM Challenge boundary conditions (left) and the Klein
boundaries with the modification discussed above. As is typical, the addition
of a guide field breaks a symmetry of the system and produces differences
between the separatrices, while the reconnection layer becomes more likely to
form plasmoids, as can be seen in the middle panels. As with the anti-parallel
case shown in Figure 3, the GEM Challenge system is restricted by the walls
and stagnates by the last panel. The Klein boundaries again allow flux to pass
through them as can be seen by the field lines intersecting the top and bottom
of the domain.

A comparison between the Klein and doubly periodic runs shows more significant
differences than in the anti-parallel case. Figure 7 compares the two in the
same format as Figure 4, but at an earlier time, $t\Omega_{cp}=45$. The reason
for this choice will be discussed below. The plasmoids in the Klein case
distort the reconnecting layer, but the overall structure is generally
similar. (Note that the color bar is the same for each pair of images so
direct comparisons can be made.) The largest differences come in panels G and
H, which show $B_{z}$ in the two runs. While the fields in the islands of
reconnected flux are roughly comparable, bands of low-strength (blue) field
are apparent just upstream of the current sheet in the doubly periodic case.
These features represent the last vestiges of unreconnected plasma being
pushed towards the X-line by an island expanding on the other (unseen) sheet.
These features do not appear in the Klein boundary case, although hints are
visible at the the top and bottom of the domain. Still, the boundary
conditions are clearly beginning to retard the inflow.

The effects of the various boundary conditions can also be seen in Figure 8,
which shows the reconnected fluxes and reconnection rates for the three runs
in the same format as Figure 5. After the peak at $t\Omega_{cp}\approx 25$,
the GEM Challenge run essentially stagnates, with the reconnection rate
falling to zero. The other two cases exhibit a similar earlier peak and
decline, but both have a later period of increasingly fast reconnection with
the Klein case peaking at $t\Omega_{cp}\approx 45$, the time shown in Figure
7. After this point, however, the Klein system stagnates while the doubly
periodic system continues to enjoy robust reconnection. It is worth noting
that simulations are usually halted once current sheets begin interacting with
each other, so the lack of agreement after $t\Omega_{cp}\approx 50-55$ is not
significant from a practical standpoint.

The elapsed times for the three computations were again similar, differing by
less than $4\%$. But, as previously, the doubly periodic case incurred twice
the computational expense due to its larger domain.

Figure 6: Out-of-plane electron current density and magnetic field lines at
$t\Omega_{cp}=20$, $40$, and $60$ for the GEM Challenge boundary conditions
(left) and the Klein boundary conditions (right) in runs with $B_{g}=1$.
Figure 7: Comparison of results with the Klein (left column) and doubly
periodic (right column) boundary conditions at $t\Omega_{cp}=45$ in runs with
$B_{g}=1$. (A)-(B): Out-of-plane electron current density ($J_{ez}$); (C)-(D)
Horizontal electron velocity ($v_{ex}$); (E)-(F): Horizontal ion velocity
($v_{ix}$); (G)-(H): $B_{z}$; (I)-(J): $E_{y}$ Figure 8: Reconnected flux
(top), $\Delta\psi$, and reconnection rate (bottom), $d(\Delta\psi)/dt$ versus
time for the GEM Challenge (blue dash), doubly periodic (red dot-dash), and
shifted Klein (black) boundary conditions in runs with a guide field
$B_{g}=1$.

##  IV Discussion

We have shown that a novel boundary condition based on the topology of a Klein
bottle can be used in numerical simulations of magnetic reconnection. The new
boundary condition combines the advantages of other common options by being
both self-reinforcing (as in systems with doubly periodic boundary conditions)
and computationally compact (as seen in the single current sheet of the GEM
Challenge boundaries).

The extension of Klein boundary conditions to full three-dimensional (3x3v)
simulations should be straightforward. The additional dimension can be
physically pictured as adding a "thickness" to the Klein bottle in Figure 1.
From a computational perspective, the most significant change is to the
identification of neighboring processors across the Klein boundary.
Neighboring processors in $z$ must change in the manner shown in panel (ii) of
Figure 2 – recall that the $z$ component of vectors flips across the Klein
boundary – but do not have to undergo the shift of panel (iii).

However, regardless of the dimensionality, Klein boundary conditions do have
limitations. Symmetric reconnection. in which the asymptotic plasmas on the
two sides of the current sheet share the same characteristics, is an
idealization, albeit a common one, that sidesteps many of the complexities of
actual systems. In reality, asymmetries can exist in any or all of the plasma
density, plasma temperature, strength of the reconnecting component of the
field, strength of the guide field, etc., and at some locations (e.g., the
terrestrial magnetopause), reconnection is nearly always asymmetric. Such
complications are straightforward to include with GEM Challenge or doubly
periodic boundary conditions because the two sides of a reconnecting current
sheet represent distinct plasmas. For Klein boundaries there is no distinction
between the two sides, and so simulating asymmetric reconnection while
maintaining the associated benefits (low computational cost and non-
stagnation) is not straightforward.

Finally, the topological equivalence of the doubly periodic system to a torus
implies the existence of certain conservation laws. Specifically, consider the
surface integral of Faraday’s Law, $\partial\mathbf{B}/\partial
t=K\bm{\nabla\times}\mathbf{E}$, where the proportionality constant $K$
depends on the choice of units. Performing a surface integral and invoking
Stokes’ theorem transforms the right-hand side into a boundary integral, which
vanishes because a torus has no boundary. Hence, the magnetic flux through the
plane of the simulation must be constant, which can be verified numerically. A
Klein bottle also has no boundary, but it is non-orientable and so Stokes’
theorem does not apply. As a consequence, the magnetic flux through the plane
of the simulation can, and does, change when Klein boundary conditions are
employed. While this difference does not appear to affect the dynamics of
reconnection, further research will be necessary to explore whether it has
physical significance.

###### Acknowledgements.

The authors acknowledge the support of NASA grants 80NSSC22K0352 and
18-DRIVE18_2-0029, "Our Heliospheric Shield", 80NSSC22M0164). The authors also
acknowledge UMD’s TREND program, funded through NSF PHY-2150399. The authors
would like to thank P. R. Colarco and E. A. Colarco for providing visual
inspiration. The simulations were carried out at the National Energy Research
Scientific Computing Center (NERSC). The data used to perform the analysis and
construct the figures for this paper are preserved at
https://doi.org/10.5281/zenodo.12785961. Further data from the associated runs
are stored at the NERSC High Performance Storage System and are available upon
request.

## References

  * Harris [1962] E. G. Harris, Nuovo Cim. 23, 115 (1962).
  * Birn  _et al._ [2001] J. Birn, J. F. Drake, M. A. Shay, B. N. Rogers, R. E. Denton, M. Hesse, M. Kuznetsova, Z. W. Ma, A. Bhattacharjee, A. Otto, and P. L. Pritchett, J. Geophys. Res. 106, 3715 (2001).
  * Horiuchi, Pei, and Sato [2001] R. Horiuchi, W. Pei, and T. Sato, Earth Planets Space 53, 439 (2001).
  * Daughton, Scudder, and Karimabadi [2006] W. Daughton, J. Scudder, and H. Karimabadi, Phys. Plasmas 13, 072101 (2006), 10.1063/1.2218817.
  * Divin  _et al._ [2007] A. V. Divin, M. I. Sitnov, M. Swisdak, and J. F. Drake, Geophys. Res. Lett. 34, L09109 (2007), 10.1029/2007GL029292.
  * Roytershteyn  _et al._ [2010] V. Roytershteyn, W. Daughton, S. Dorfman, Y. Ren, H. Ji, M. Yamada, H. Karimabadi, L. Yin, B. J. Albright, and K. J. Bowers, Phys. Plasmas 17 (2010), 10.1063/1.3399787.
  * Zeiler  _et al._ [2002] A. Zeiler, D. Biskamp, J. F. Drake, B. N. Rogers, M. A. Shay, and M. Scholer, J. Geophys. Res. 107, 1230 (2002).