# Theory for Planar Hall Effect in Organic Dirac Fermion System

Yuki Nakamura<EMAIL_ADDRESS>Takao Morinari
<EMAIL_ADDRESS>Course of Studies on Materials Science,
Graduate School of Human and Environmental Studies, Kyoto University, Kyoto
606-8501, Japan

###### Abstract

In a recent experiment on the interlayer magnetoresistance in the quasi-two-
dimensional organic salt, $\alpha$-(BEDT-TTF)2I3, it has been observed that at
low temperatures, interlayer tunneling attains phase coherence, leading to the
emergence of a three-dimensional electronic structure. Theoretically and
experimentally it has been suggested that the system exhibits characteristics
of a three-dimensional Dirac semimetal as a consequence of broken time-
reversal symmetry and inversion symmetry. Here, we perform a theoretical
calculation of the magnetoconductivity under an in-plane magnetic field and
demonstrate that the system displays a planar Hall effect. Our calculations
are based on a realistic model for $\alpha$-(BEDT-TTF)2I3 incorporating
interlayer tunneling and the tilt of the Dirac cone. Given that the planar
Hall effect is anticipated as a consequence of chiral anomaly, our findings
provide support for the classification of $\alpha$-(BEDT-TTF)2I3 as a three-
dimensional Dirac semimetal.

Massless Dirac and Weyl semimetals have been extensively studied recently
because of their unique and intriguing electrical properties.Murakami (2007);
Burkov and Balents (2011); Yang  _et al._ (2011); Xu  _et al._ (2011); Young
_et al._ (2012); Wang  _et al._ (2012, 2013); Neupane  _et al._ (2012);
Armitage  _et al._ (2018) The energy spectrum in these systems is
characterized by the touching of the valence band and conduction band at
discrete momentum points. The key distinction between the Dirac/Weyl semimetal
and the two-dimensional Dirac fermion system lies in the presence of broken
time-reversal symmetry and/or inversion symmetry. To realize the Dirac/Weyl
semimetal, it is necessary to break either time-reversal symmetry or inversion
symmetry, or both.

Organic charge-transfer salt, $\alpha$-(BEDT-TTF)2I3, has been studied as a
quasi-two-dimensional Dirac fermion system.Katayama  _et al._ (2006);
Kobayashi  _et al._ (2007); Kajita  _et al._ (2014) (Here, BEDT-TTF is
bis(ethylenedithio)tetrathiafulvalene.) One of the present authors
theoretically predictedMorinari (2020, ) that both time-reversal symmetry and
inversion symmetry are broken, and, as a result, the system becomes a three-
dimensional Dirac semimetal when the interlayer tunneling becomes phase
coherent at low temperatures. The phase coherence in the interlayer tunneling
is confirmed experimentallyTajima  _et al._ (2023) by the observation of the
peak structure in the interlayer magnetoresistance. Furthermore, the
observation of the negative magnetoresistance and the planar Hall effect (PHE)
has been reported recentlyTajima  _et al._ that is associated with chiral
anomalyNielsen and Ninomiya (1983); Son and Spivak (2013); Burkov (2015);
Nandy  _et al._ (2017); Burkov (2017); Xiong  _et al._ (2015); Huang  _et al._
(2015); Hirschberger  _et al._ (2016); Zhang  _et al._ (2016); Li  _et al._
(2016) in a Dirac semimetal.

In this Letter, we consider a model that includes inter-layer tunneling and
the tilt of the Dirac cone that exists in $\alpha$-(BEDT-TTF)2I3.Katayama  _et
al._ (2006); Kobayashi  _et al._ (2007) Based on the semiclassical Boltzmann
equation, we compute the magnetoconductivity under in-plane magnetic fields.
We show that the system exhibits a PHE using a set of realistic parameters for
$\alpha$-(BEDT-TTF)2I3.

In the absence of the interlayer tunneling, there are two Driac cones in the
$k_{x}$-$k_{y}$ plane.Katayama  _et al._ (2006); Kobayashi  _et al._ (2007) By
including the interlayer tunneling, each Dirac cone splits into two Dirac
cones. In contrast to systems where spin degeneracy is lifted due to the
breaking of time-reversal symmetry caused by magnetic correlations, the spin
remains degenerate in $\alpha$-(BEDT-TTF)2I3 because the time-reversal
symmetry breaking is not associated with magnetic correlations.Morinari (2020)
For the sake of simplicity, we neglect the spin degrees of freedom in the
follwoing analysis.

The Hamiltonian for two of the four Dirac cones is given by

$\displaystyle H(\bm{k})$ $\displaystyle=\hbar vk_{x}\tau_{x}+\hbar
vk_{y}\tau_{y}-2t_{2}\cos(ck_{z})\tau_{z}$
$\displaystyle\quad+\left[-2t_{1}\cos(ck_{z})+\hbar
uk_{x}\right]\tau_{0}+\varepsilon_{\rm D}.$ (1)

Here $k_{x}$ and $k_{y}$ are in-plane wave numbers measured from the Dirac
point and $k_{z}$ is the wave number perpendicular to the $k_{x}$-$k_{y}$
plane. We note that the position of the Dirac point in the plane is irrelevant
for the following calculation, though we need to include them to make clear
the presence of the symmetry breaking. The parameter $u$ describes the tilt of
the Dirac cone to the $k_{x}$ axis, and we neglect anisotropy in the Dirac
cone in the plane. $c$ is the lattice constant in the $c$-axis.
$\tau_{x},\tau_{y},\tau_{z}$ are the Pauli matrices and $\tau_{0}$ is the
$2\times 2$ identity matrix. $t_{1}$ and $t_{2}$ are the parameters for the
interlayer tunneling. $t_{1}$ is for the tunneling between the same molecules,
and $t_{2}$ is for the tunneling between the adjacent molecules along the
$a$-axis. Because of the parameter $t_{2}$ each of the Dirac cones at
$k_{z}=0$ splits into two Dirac cones residing at $k_{z}=\pm\pi/2$. The Dirac
cone is type-I in the $k_{x}$-$k_{y}$ plane,Tajima and Morinari (2018) so the
range of the parameter $u$ is $-v<u<v$. The other two Dirac cones are
described by Eq. (1) with $k_{x}\rightarrow-k_{x}$. We may assume
$t_{1}>t_{2}$ from the crystal structure of $\alpha$-(BEDT-TTF)2I3.Bender  _et
al._ (1984) In this case, the Dirac cone is type-IISoluyanov  _et al._ (2015)
in the $k_{z}$ direction. The parameter $\varepsilon_{\rm D}$ denotes the
energy of the Dirac point. We assign different values of $\varepsilon_{\rm D}$
to the two Dirac cones in the $k_{x}$-$k_{y}$ plane to incorporate the
symmetry breaking.

The energy dispersion is given by $E_{\bm{k}}^{\left(\pm\right)}=(\hbar
v/a)\tilde{E}_{\bm{k}}^{\left(\pm\right)}$ where

$\displaystyle\tilde{E}_{\bm{k}}^{\left(\pm\right)}=\pm\tilde{E}_{\bm{k}}-2\tilde{t_{1}}\cos(ck_{z})+\eta
ak_{x},$ (2)

with

$\tilde{E}_{\bm{k}}=\sqrt{a^{2}(k_{x}^{2}+k_{y}^{2})+4\tilde{t_{2}}^{2}\cos^{2}(ck_{z})}.$
(3)

Here, $a$ is the in-plane lattice constant. We take the same lattice constants
for $a$ and $b$ axes for simplicity. We defined the following dimensionless
parameters,

$\displaystyle\tilde{t_{1}}$ $\displaystyle=\frac{t_{1}}{\hbar v/a},$
$\displaystyle\tilde{t_{2}}$ $\displaystyle=\frac{t_{2}}{\hbar v/a},$
$\displaystyle\eta$ $\displaystyle=\frac{u}{v}.$ (4)

Taking $a=1.0\times 10^{-9}$ m and $v=5.0\times 10^{4}$ m/s, we find $\hbar
v/a=3.3\times 10^{-2}$ eV. Note that $\varepsilon_{\rm D}$ is note included in
$E_{\bm{k}}^{\left(\pm\right)}$.

Figure 1(a) shows the energy dispersion in the plane and Fig. 1(b) shows that
in the $k_{z}$ direction. We see that the Dirac cone is type-I in the
$k_{x}$-$k_{y}$ plane and type-II in the $k_{z}$ axis as stated above. Figure
1(c) shows the Fermi surface. If the Fermi energy is larger than $t_{1}$ and
$t_{2}$, the Fermi surface is a warped cylinder.Tajima  _et al._ (2023) For
$\alpha$-(BEDT-TTF)2I3, the Fermi energy is expected to be smaller than
$t_{1}$ and $t_{2}$.Tajima  _et al._ (2023) In this case, the Fermi surface
splits into a single electronic Fermi surface and two hole Fermi surfaces as
shown in Fig. 1(c). Because of the tilt parameter $\eta$, which is slightly
lower than one,Tajima and Morinari (2018) the Fermi surface is largely
deformed.

Figure 1:  (Color online) (a) Energy dispersion of the model described by the
Hamiltonian (1) in the $k_{x}$-$k_{y}$ plane and (b) along the $k_{z}$ axis.
The tilt parameter in the plane is $\eta=0.70$, and so the Dirac cone is
type-I. The interlayer hopping parameters are $\tilde{t_{1}}=0.10$ and
$\tilde{t_{2}}=0.05$, and so the Dirac cone is type-II in the $k_{z}$ axis.
(c) The Fermi surface around the two Dirac cones. The Fermi energy is set to
be zero and we set $\varepsilon_{\rm D}=0.03$. The Fermi surface consists of
three portions: the middle one is the electron Fermi surface and the other two
are the hole Fermi surfaces.

We calculate the magnetoconductivity using the semiclassical Boltzmann
equation employing the relaxation time approximation. The application of the
Boltzmann equation is justified when $\omega_{c}\tau<1$ with $\omega_{c}$
being the cyclotron frequency and $\tau$ being the scattering time. Therefore,
our result is limited to the regime of relatively weak magnetic fields. In the
presence of the electric field ${\bf E}$ and the magnetic field ${\bf B}$, the
quasiclassical equation of motion is given byDuval  _et al._ (2006); Xiao  _et
al._ (2006)

$\displaystyle\hbar\frac{{d{\bf{k}}}}{{dt}}$ $\displaystyle=$
$\displaystyle\frac{1}{{1+\frac{e}{\hbar}{\bf{B}}\cdot\bm{\Omega_{k}}}}$ (5)
$\displaystyle\times\left[{-e{{\bf{v}}_{\bf{k}}}\times{\bf{B}}-e{\bf{E}}-\frac{{{e^{2}}}}{\hbar}\left({{\bf{B}}\cdot{\bf{E}}}\right)\bm{\Omega_{k}}}\right],$
$\displaystyle\frac{{d{\bf{r}}}}{{dt}}$ $\displaystyle=$
$\displaystyle\frac{1}{{1+\frac{e}{\hbar}{\bf{B}}\cdot\bm{\Omega_{k}}}}$ (6)
$\displaystyle\times\left[{{{\bf{v}}_{\bf{k}}}+\frac{e}{\hbar}\left({\bm{\Omega_{k}}\cdot{{\bf{v}}_{\bf{k}}}}\right){\bf{B}}+\frac{e}{\hbar}{\bf{E}}\times\bm{\Omega_{k}}}\right],$

where $\Omega_{\bf{k}}$ is the Berry curvature. There are three contributions
arising from the Berry curvature: The second term in the angular bracket in
Eq. (6) leads to a chiral magnetic effect.Vilenkin (1980); Fukushima  _et al._
(2008) The third term leads to an anomalous Hall effect.Yang  _et al._ (2011)
We do not consider these two effects here and focus on the chiral anomaly
effect arising from the third term in the angular bracket in Eq. (5).

From the energy dispersion (2), the group velocity is given by

$\displaystyle\bm{v_{k}}^{\left(\pm\right)}=$ $\displaystyle
v\left(\pm\frac{ak_{x}}{\tilde{E}_{\bm{k}}}+\eta,\pm\frac{ak_{y}}{\tilde{E}_{\bm{k}}},\right.$
$\displaystyle\left.\mp\frac{4\frac{c}{a}\tilde{t_{2}}^{2}\sin(ck_{z})\cos(ck_{z})}{\tilde{E}_{\bm{k}}}+2\frac{c}{a}\tilde{t_{1}}\sin(ck_{z})\right)$
(7)

The Berry curvatureXiao  _et al._ (2010) is given by

$\displaystyle\bm{\Omega_{k}}^{\left(\pm\right)}=$
$\displaystyle\left(\mp\frac{2a^{2}c\tilde{\tau_{2}}k_{x}\sin(ck_{z})}{2\tilde{E}_{\bm{k}}^{3}},\mp\frac{2a^{2}c\tilde{\tau_{2}}k_{y}\sin(ck_{z})}{2\tilde{E}_{\bm{k}}^{3}},\right.$
$\displaystyle\quad\left.\pm\frac{2a^{2}\tilde{\tau_{2}}\cos(ck_{z})}{2\tilde{E}_{\bm{k}}^{3}}\right).$
(8)

Here, $\bm{v_{k}}^{\left(+\right)}$ and $\bm{\Omega_{k}}^{\left(+\right)}$ are
for the positive energy state, $\tilde{E}_{\bm{k}}^{\left(+\right)}$, and
$\bm{v_{k}}^{\left(-\right)}$ and $\bm{\Omega_{k}}^{\left(-\right)}$ are for
the negative energy state, $\tilde{E}_{\bm{k}}^{\left(-\right)}$.

From the Boltzmann equation, we obtain the formula for the
magnetoconductivityNandy  _et al._ (2017); Burkov (2017):

$\displaystyle\sigma_{xx}^{(\pm)}$ $\displaystyle=$
$\displaystyle\frac{2e^{2}\tau}{(2\pi)^{3}}\int
d^{3}\bm{k}\left(-\frac{\partial
f_{\rm{eq}}}{\partial\epsilon}\right)\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega_{k}}^{\left(\pm\right)}}$
(9)
$\displaystyle\quad\times\left[v_{x}^{\left(\pm\right)}+\frac{e}{\hbar}B_{x}(\bm{v_{k}^{\left(\pm\right)}}\cdot\bm{\Omega_{k}}^{\left(\pm\right)})\right]^{2},$
$\displaystyle\sigma_{xy}^{(\pm)}$ $\displaystyle=$
$\displaystyle\frac{2e^{2}\tau}{(2\pi)^{3}}\int
d^{3}\bm{k}\left(-\frac{\partial
f_{\rm{eq}}}{\partial\epsilon}\right)\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega_{k}}^{\left(\pm\right)}}$
(10)
$\displaystyle\quad\times\left[v_{x}^{\left(\pm\right)}+\frac{e}{\hbar}B_{x}(\bm{v_{k}}^{\left(\pm\right)}\cdot\bm{\Omega_{k}}^{\left(\pm\right)})\right]$
$\displaystyle\quad\times\left[v_{y}^{\left(\pm\right)}+\frac{e}{\hbar}B_{y}(\bm{v_{k}}^{\left(\pm\right)}\cdot\bm{\Omega_{k}}^{\left(\pm\right)})\right],$

where $f_{\rm{eq}}$ is the equilibrium Fermi-Dirac distribution function. We
compute the components of the positive energy state, denoted by superscript
$(+)$ and the negative energy state, denoted by superscript $(-)$, separately.
Here, the magnetic field is given by
$\bm{B}=(B_{x},B_{y},0)=(B\cos\phi,B\sin\phi,0)$. In order to obtain the total
magnetoconductivity, we take the sum of $\sigma_{xx}^{(+)}+\sigma_{xx}^{(-)}$
and $\sigma_{xy}^{(+)}+\sigma_{xy}^{(-)}$. We also calculate the contribution
from the other two Dirac cones. The splitting of each Dirac cone in the
$k_{z}$ direction results in a twofold multiplication factor.

The result is shown in Fig. 2. We subtract the constant value from
$\sigma_{xx}$, and only the oscillating component is shown in Fig. 2(a). Thus,
we show only the components of the planar Hall effect corresponding to
$\sigma_{xx}$ and $\sigma_{xy}$, denoted as $\sigma_{xx}^{\rm PHE}$ and
$\sigma_{xy}^{\rm PHE}$, respectively. They are plotted as the function of
$\phi$ for different values of $b=\left(a/\ell_{B}\right)^{2}$ with
$\ell_{B}=\sqrt{\hbar/eB}$ the magnetic length. $b$ is defined as the
dimensionless magnetic field parameter. At $B=1$ T, $b=1.5\times 10^{-3}$. The
unit of conductivity is $e^{2}\tau v/(2\pi^{3}\hbar ac)$. For the interlayer
tunneling parameters, $\tilde{t}_{1}$ and $\tilde{t}_{2}$, we take
$\tilde{t}_{1}=0.10$ and $\tilde{t}_{2}=0.05$. For the tilt parameter we take
$\eta=0.7$. This set of parameters is reasonable for $\alpha$-(BEDT-TTF)2I3.
We note that both $\sigma_{xx}^{\rm PHE}$ and $\sigma_{xy}^{\rm PHE}$ exhibit
the periodicity of $\pi$. This oscillating behavior can be associated with the
chiral anomaly.Nandy  _et al._ (2017); Burkov (2017) Qualitatively similar
behavior is observed in a recent experiment.Tajima  _et al._

Figure 2:  (Color online) (a) Longitudinal conductivity and (b) planar Hall
conductivity as the function of $\phi$ for different values of $b$. The unit
of conductivity is $e^{2}\tau v/(2\pi^{3}\hbar ac)$. The constant component is
subtracted from $\sigma_{xx}$. The parameters are $\tilde{t_{1}}=0.10$,
$\tilde{t_{2}}=0.05$, and $\eta=0.70$. For the Dirac point energies, we set
$\varepsilon_{\rm D}/(\hbar v/a)=0.5$ for two Dirac cones along the $k_{z}$
axis and $\varepsilon_{\rm D}/(\hbar v/a)=-0.4$ for the other two Dirac cones.

We also examined the magnetic field parameter $b$ dependence of the amplitude
of $\sigma_{xx}^{\rm PHE}$ and $\sigma_{xy}^{\rm PHE}$ as shown in Fig. 3(a).
We find that the amplitude varies quadratically with the magnetic field. If
there remains the effect associated with the tilt of the Dirac cone, we may
expect a linear dependence, but there is no such component. This is understood
by complete cancellation between the contribution from the Dirac cones with
opposite tilts and chiralities. We note that the amplitudes of
$\sigma_{xx}^{\rm PHE}$ is slightly larger than $\sigma_{xy}^{\rm PHE}$. This
behavior is qualitatively in agreement with experimental observations, where
the amplitude of $\sigma_{xx}^{\rm PHE}$ is ten times larger than that of
$\sigma_{xy}^{\rm PHE}$ at 3 T.Tajima  _et al._ The difference of the
amplitutdes is associated with the inteplay between the tilt parameter
dependence of the group velocity and the density of states. To make clear the
tilt parameter dependence, we calculate $\eta$ dependence of $\sigma_{xx}^{\rm
PHE}$ and $\sigma_{xy}^{\rm PHE}$ as shown in Fig. 3(b). When $\eta=0$, there
is no difference in the amplitudes of $\sigma_{xx}^{\rm PHE}$ and
$\sigma_{xy}^{\rm PHE}$. Their difference increases as we increase $\eta$.
However, the result depends on the choice of two values of $\varepsilon_{\rm
D}$. If we take a different set of values for $\varepsilon_{\rm D}$, we obtain
a different $\eta$ dependence. The energy dispersion exhibits particle-hole
symmetry; however, the integrands in Eqs. (9) and (10) do not. Consequently,
the $\eta$ dependence of $\sigma_{xx}^{\rm PHE}$ and $\sigma_{xy}^{\rm PHE}$
is non-trivial.

Figure 3:  (Color online) (a) The magnetic field dependence of the amplitudes
of $\sigma_{xx}^{\rm PHE}$ and $\sigma_{xx}^{\rm PHE}$ at $b=0.001$ using the
same parameters as in Fig. 2. The amplitude varies quadratically with the
magnetic field and there is no linear term. (b) The tilt parameter dependence
of the amplitudes of $\sigma_{xx}^{\rm PHE}$ and $\sigma_{xx}^{\rm PHE}$ at
$b=0.001$.

To conclude, we have shown that the magnetoconductivity exhibit PHE in a
realistic model for $\alpha$-(BEDT-TTF)2I3. Since $\alpha$-(BEDT-TTF)2I3 does
not show any indication of ferromagnetism,Konoike  _et al._ (2022) the
presence of PHE suggests the chiral anomaly effect that is associated with a
three-dimensional Dirac semimetal. While our analysis is confined to a small
magnetic field range due to the utilization of the semiclassical Boltzmann
equation, we anticipate the occurrence of the PHE at high magnetic fields,
provided there is no qualitative change between the low and high magnetic
field regimes. This seems to be consistent with the recent experiment.Tajima
_et al._ In conjunction with the experimental findings,Tajima  _et al._ (2023,
) our results provide strong support for the classification of $\alpha$-(BEDT-
TTF)2I3 as a three-dimensional Dirac semimetal under conditions of low
temperatures and high pressures.

###### Acknowledgements.

We thank N. Tajima for helpful discussions and sharing experimental data. The
research was supported by JSPS KAKENHI Grant Number 22K03533.

## References

  * Murakami (2007) S. Murakami, New J. Phys. 9, 356 (2007).
  * Burkov and Balents (2011) A. A. Burkov and L. Balents, Phys. Rev. Lett. 107 (2011).
  * Yang  _et al._ (2011) K.-Y. Yang, Y.-M. Lu, and Y. Ran, Phys. Rev. B 84 (2011).
  * Xu  _et al._ (2011) G. Xu, H. Weng, Z. Wang, X. Dai, and Z. Fang, Phys. Rev. Lett. 107, 186806 (2011).
  * Young  _et al._ (2012) S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane, E. J. Mele, and A. M. Rappe, Phys. Rev. Lett. 108, 140405 (2012).
  * Wang  _et al._ (2012) Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, and Z. Fang, Phys. Rev. B 85, 195320 (2012).
  * Wang  _et al._ (2013) Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Phys. Rev. B 88, 125427 (2013).
  * Neupane  _et al._ (2012) M. Neupane, S.-Y. Xu, L. A. Wray, A. Petersen, R. Shankar, N. Alidoust, C. Liu, A. Fedorov, H. Ji, J. M. Allred, Y. S. Hor, T.-R. Chang, H.-T. Jeng, H. Lin, A. Bansil, R. J. Cava, and M. Z. Hasan, Phys. Rev. B 85, 235406 (2012).
  * Armitage  _et al._ (2018) N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Katayama  _et al._ (2006) S. Katayama, A. Kobayashi, and Y. Suzumura, J. Phys. Soc. Jpn. 75, 054705 (2006).
  * Kobayashi  _et al._ (2007) A. Kobayashi, S. Katayama, Y. Suzumura, and H. Fukuyama, J. Phys. Soc. Jpn. 76, 034711 (2007).
  * Kajita  _et al._ (2014) K. Kajita, Y. Nishio, N. Tajima, Y. Suzumura, and A. Kobayashi, J. Phys. Soc. Jpn. 83, 072002 (2014).
  * Morinari (2020) T. Morinari, J. Phys. Soc. Jpn. 89, 073705 (2020).
  * (14) T. Morinari, arXiv:2308.11141 .
  * Tajima  _et al._ (2023) N. Tajima, Y. Kawasugi, T. Morinari, R. Oka, T. Naito, and R. Kato, J. Phys. Soc. Jpn. 92 (2023).
  * (16) N. Tajima, Y. Kawasugi, T. Morinari, R. Oka, T. Naito, and R. Kato, arXiv:2302.05616 .
  * Nielsen and Ninomiya (1983) H. Nielsen and M. Ninomiya, Phys. Lett. B 130, 389 (1983).
  * Son and Spivak (2013) D. T. Son and B. Z. Spivak, Phys. Rev. B 88, 104412 (2013).
  * Burkov (2015) A. A. Burkov, J. Phys.: Condens. Matter 27, 113201 (2015).
  * Nandy  _et al._ (2017) S. Nandy, G. Sharma, A. Taraphder, and S. Tewari, Phys. Rev. Lett. 119, 176804 (2017).
  * Burkov (2017) A. A. Burkov, Phys. Rev. B 96, 041110(R) (2017).
  * Xiong  _et al._ (2015) J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava, and N. P. Ong, Science 350, 413 (2015).
  * Huang  _et al._ (2015) X. Huang, L. Zhao, Y. Long, P. Wang, D. Chen, Z. Yang, H. Liang, M. Xue, H. Weng, Z. Fang, X. Dai, and G. Chen, Phys. Rev. X 5, 031023 (2015).
  * Hirschberger  _et al._ (2016) M. Hirschberger, S. Kushwaha, Z. Wang, Q. Gibson, S. Liang, C. A. Belvin, B. A. Bernevig, R. J. Cava, and N. P. Ong, Nat. Mater. 15, 1161 (2016).
  * Zhang  _et al._ (2016) C.-L. Zhang, S.-Y. Xu, I. Belopolski, Z. Yuan, Z. Lin, B. Tong, G. Bian, N. Alidoust, C.-C. Lee, S.-M. Huang, T.-R. Chang, G. Chang, C.-H. Hsu, H.-T. Jeng, M. Neupane, D. S. Sanchez, H. Zheng, J. Wang, H. Lin, C. Zhang, H.-Z. Lu, S.-Q. Shen, T. Neupert, M. Z. Hasan, and S. Jia, Nat. Commun. 7, 10735 (2016).
  * Li  _et al._ (2016) H. Li, H. He, H.-Z. Lu, H. Zhang, H. Liu, R. Ma, Z. Fan, S.-Q. Shen, and J. Wang, Nat. Commun. 7 (2016).
  * Tajima and Morinari (2018) N. Tajima and T. Morinari, J. Phys. Soc. Jpn. 87, 045002 (2018).
  * Bender  _et al._ (1984) K. Bender, I. Hennig, D. Schweitzer, K. Dietz, H. Endres, and H. J. Keller, Mol. Cryst. Liq. Cryst. 108, 359 (1984).
  * Soluyanov  _et al._ (2015) A. A. Soluyanov, D. Gresch, Z. Wang, Q. Wu, M. Troyer, X. Dai, and B. A. Bernevig, Nature 527, 495 (2015).
  * Duval  _et al._ (2006) C. Duval, Z. Horváth, P. A. Horváthy, L. Martina, and P. C. Stichel, Mod. Phys. Lett. B 20, 373 (2006).
  * Xiao  _et al._ (2006) D. Xiao, Y. Yao, Z. Fang, and Q. Niu, Phys. Rev. Lett. 97, 026603 (2006).
  * Vilenkin (1980) A. Vilenkin, Phys. Rev. D 22, 3080 (1980).
  * Fukushima  _et al._ (2008) K. Fukushima, D. E. Kharzeev, and H. J. Warringa, Phys. Rev. D 78, 074033 (2008).
  * Xiao  _et al._ (2010) D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Konoike  _et al._ (2022) T. Konoike, T. Terashima, S. Uji, Y. Hattori, and R. Kato, J. Phys. Soc. Japan 91 (2022).


# Learning to Bid in Forward Electricity Markets
Using a No-Regret Algorithm

Arega Getaneh Abatea, Dorsa Majdib, Jalal Kazempoura, and Maryam Kamgarpourc a
Department of Wind and Energy Systems, Technical University of Denmark, Kgs.
Lyngby, Denmark  b Sharif University of Technology, Tehran, Iran  c Sycamore
lab, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland  {ageab,
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

It is a common practice in the current literature of electricity markets to
use game-theoretic approaches for strategic price bidding. However, they
generally rely on the assumption that the strategic bidders have prior
knowledge of rival bids, either perfectly or with some uncertainty. This is
not necessarily a realistic assumption. This paper takes a different approach
by relaxing such an assumption and exploits a no-regret learning algorithm for
repeated games. In particular, by using the _a posteriori_ information about
rivals’ bids, a learner can implement a no-regret algorithm to optimize
her/his decision making. Given this information, we utilize a multiplicative
weight-update algorithm, adapting bidding strategies over multiple rounds of
an auction to minimize her/his regret. Our numerical results show that when
the proposed learning approach is used the social cost and the market-clearing
prices can be higher than those corresponding to the classical game-theoretic
approaches. The takeaway for market regulators is that electricity markets
might be exposed to greater market power of suppliers than what classical
analysis shows.

###### Index Terms:

Online learning, Strategic bidding, No-regret algorithm, Market power,
Diagonalization.

##  I Introduction

###  I-A Background and Motivation

Strategic bidding is a key task in the daily operation of electricity
suppliers, so-called bidders. There is a vast literature, mostly built upon
game-theoretic approaches such as bi-level programming, allowing bidders to
derive their strategic bids in terms of price or quantity or both [1, 2, 3].
One critical assumption in such works is that the bidder has knowledge of the
decision making process of rival bidders, market price formation, and thereby
the market outcomes for every bidder. While some works model the uncertainty
on rival bids [4], they still assume that the bidder perfectly knows how the
market is being cleared.

These assumptions do not necessarily hold in real-world electricity auctions
for several reasons, such as (i) bidding decisions are private to bidders,
(ii) bidders typically only have access to public post-auction information
such as aggregate demand and supply curves and market-clearing prices, and
(iii) electricity auctions often involve multiple rounds and different set of
bidders in each round. Therefore, formulating bidding problems with such
unrealistic assumptions may eventually result in suboptimal decisions [5, 6]
or outcomes that do not match the real world.

This paper adopts an approach based on learning in repeated games [7], and
leverages no-regret learning from the online optimization literature [8, 9] to
address this problem. There has been a growing body of work that leverages
online and reinforcement learning for bidding in electricity markets. In the
following section, we highlight the related works.

###  I-B Literature Review

In [6], the authors investigate the sequential optimization of bidders to
maximize payoffs using an adversarial multiarmed bandit-based real-time
bidding scheme in an electricity market. For strategic bidding in electricity
markets, a deep reinforcement learning method is developed in [10], where
physical non-convex operating characteristics of generators are accounted for.
Considering unobservable information, [11] develops a multi-agent policy
gradient algorithm aimed at approximating the Nash equilibrium among multiple
strategic bidders in the day-ahead electricity market. The authors conclude
that their algorithm enables all the market participants to find superior
bidding strategies with increased profit gains. Reference [12] considers
online learning in the general class of repeated electricity market auctions,
and presents an approach to incorporate extra information available when a bid
is not accepted, to estimate the utility function of a bidder. A learning-
based approach is developed in [13] using a multi-armed bandit algorithm for
real-time pricing in a demand response program. Finally, [14] utilizes online
learning approaches to address the problem of optimal bidding for virtual
trading in two-settlement electricity markets.

To the best of our knowledge, none of the above work, and more broadly
existing studies on learning in electricity markets, have considered the day-
ahead electricity market auction in light of public post-market clearing data.
In particular, in several current day-ahead markets, after each round of the
auction, each bidder not only observes the outcomes of her/his submitted bids
but also gets access to (i) hourly market-clearing price, (ii) hourly buy and
sell volumes, and (iii) hourly aggregate supply and demand curves. Such
feedback information can be used by the bidder to assess potential outcomes of
her/his alternative price bids111Please see the day-ahead aggregate supply and
demand curves published by European market operators, such as EPEX [15] and
Nord Pool [16].. We refer to such an information setting as full information
feedback.

While advanced learning algorithms can be designed and tailored to a specific
electricity market, our goal in this paper is to use an off-the-shelf online
learning algorithm under full information feedback to understand the effect of
online learning on social welfare and the learner’s utility.

###  I-C Our Contributions

Our learning approach aims to serve as a benchmark to analyze how much market
power bidders could potentially exert in a day-ahead market. To this end, we
employ full information feedback and a simple market setting for which we can
compute the Nash equilibrium. This enables us to compare social welfare and
bidders’ payoffs under the learning strategy with that corresponding to the
Nash equilibrium strategy. The learning algorithm employed is the Hedge
algorithm [17] and enables the bidder over time, that is, over auction rounds,
to effectively learn how to place a strategic price bid in a day-ahead
electricity market without making any assumption on knowing the rivals’ true
costs or their next actions a priori.

To compare the learning results with existing approaches in the literature, we
consider each bidder solving a bi-level program, perfectly knowing the best
response of rival bidders, leading to an equilibrium problem with equilibrium
constraints (EPEC)[18]. This problem can then be solved by a diagonalization
method [19, 1]. In our case study, we verify that the approach converges to a
Nash equilibrium strategy and hence, serves as a benchmark for our learning
approach.

Our study illustrates the extent of social welfare suboptimality when an
online learning approach is employed. Furthermore, it illustrates the cases
for which a given bidder can improve her/his utilities by employing such a
learning approach. Last, it also highlights the limitations of the no-regret
criteria when used to infer bidders’ utilities in an electricity market
auction.

The rest of the paper is organized as follows. Section II provides
preliminaries on repeated electricity market auctions. Section III lays out
the proposed no-regret algorithm. Section IV formulates the benchmark model
and its solution method. Section V presents illustrative case studies. Section
VI concludes the paper.

##  II Repeated Electricity Market Auction

We study a day-ahead electricity market and consider the electricity market
auction as our stage game. The set of participants consists of the power
generators referred to as Bidders
$\ell\in\mathcal{N}=\\{1,\ldots,|\mathcal{N}|\\}$ and an auctioneer or Market
Operator whose goal is to procure the total market demand, modeled by a price-
inelastic demand $Q\in\mathbb{R}_{+}$, from the bidders at the minimum cost.
We assume each bidder has a private true cost function
$C_{\ell}:\mathbb{X}_{\ell}\rightarrow\mathbb{R}_{+}$,
$\mathbb{X}_{\ell}\subseteq\mathbb{R}_{+}$. We further assume that
$0\in\mathbb{X}_{\ell}$ and $C_{\ell}(0)=0$ and each bidder’s cost function is
quadratic, therefore is of the form
$C_{\ell}(x)=\frac{1}{2}c_{\ell}x^{2}+d_{\ell}x$ [20]. Each bidder $\ell$ has
a finite strategy set
$\mathcal{K}_{\ell}=\left\\{1,\ldots,\left|\mathcal{K}_{\ell}\right|\right\\}$
that consists of her/his true cost function and bid functions of the form
$b_{\ell}^{k},$ where
$b_{\ell}^{k}(x)=\frac{1}{2}c^{k}_{\ell}x^{2}+d^{k}_{\ell}x$. Without loss of
generality, $1\in\mathcal{K}_{\ell}$ corresponds to the true cost function.

Let $T\\!\in\\!\mathbb{N}$ be the time horizon length, i.e., the number of
rounds the underlying auction is repeated. To illustrate the performance of
the algorithm in a benchmark setting, we assume that for all rounds, the total
demand $Q$, the true costs, and the strategy set $\mathcal{A}_{\ell}$ remain
unchanged. Let $k_{\ell,t}\in\mathcal{K}_{\ell}$ denote the strategy of bidder
$\ell$ at time $t\leq T$ which is only revealed to the market operator and ex-
post published by the market operator anonymously as a part of the aggregate
supply curve. Given the strategy profile
$\mathcal{B}_{t}=\left\\{b_{\ell}^{k_{\ell,t}}\right\\}_{\ell\in\mathcal{N}}$,
a mechanism defines an allocation rule
$x_{\ell}^{*}(\mathcal{B}_{t})\in\mathbb{X}_{\ell}$, and a payment rule
$p_{\ell}(\mathcal{B}_{t})\in\mathbb{R}$ for each bidder $\ell$. In many
auctions, the allocation rule is determined by an optimization problem

$\displaystyle J(\mathcal{B}_{t})=$
$\displaystyle\min_{x\in\mathbb{X}_{t}}\sum_{\ell\in\mathcal{N}}b_{\ell}^{k_{\ell,t}}\left(x_{\ell}\right)$
(1a) $\displaystyle\text{ s.t.
}\sum_{\ell\in\mathcal{N}}x_{\ell}\in\mathbb{S},$ (1b)

where
$\mathbb{X}=\prod_{\ell\in\mathcal{N}}\mathbb{X}_{\ell}\subset\mathbb{R}_{+}^{|\mathcal{N}|}$
is the decision set for the optimization problem, and the set
$\mathbb{S}\subset\mathbb{R}_{+}$ corresponds to the market constraints.

In an electricity market auction, these constraints in their simplest form may
correspond to the capacity constraints of bidders and the total demand of the
market. Hence, the market-clearing optimization problem solved by the market
operator can be expressed as

$\displaystyle\min_{x\in\mathbb{X}}\quad$
$\displaystyle\sum_{\ell\in\mathcal{N}}\frac{1}{2}c_{\ell}^{k_{\ell,t}}x_{\ell}^{2}+d_{\ell}^{k_{\ell,t}}x_{\ell}$
(2a) s.t. $\displaystyle\sum_{\ell}x_{\ell}=Q,$ (2b) $\displaystyle 0\leq
x_{\ell}\leq\bar{x}_{\ell}\ \forall{\ell\in\mathcal{N}}.$ (2c)

Let the optimal solution of (2) be denoted by $x^{*}(\mathcal{B}_{t})$. Let
$\lambda^{*}(\mathcal{B}_{t})\\!\in\\!\mathbb{R}$ denote the Lagrange
multiplier associated with (2b). In the context of the electricity market
auction, $\lambda^{*}(\mathcal{B}_{t})$ is called the marginal price and is
announced after each round $t$ to all bidders. The payment for each bidder
$\ell$ is
$p_{\ell}(\mathcal{B}_{t})=\lambda^{*}(\mathcal{B}_{t})x_{\ell}^{*}$. The
utility of bidder $\ell$ is linear in the payment received;
$u_{\ell,t}=u_{\ell}(\mathcal{B}_{t})=p_{\ell}(\mathcal{B}_{t})-C_{\ell}\left(x_{\ell}^{*}\right)$.
A bidder whose bid is not accepted i.e., $x_{\ell}^{*}=0$, is not paid, and
$u_{\ell,t}=0$.

While Nash equilibrium offers a theoretical solution concept, real-world
bidders are profit-maximizing entities operating within constraints of privacy
and limited information. Consequently, it may not be practical to assume that
they can compute or would be willing to select their Nash equilibrium
strategy. Instead, a more pragmatic assumption is that bidders choose their
strategies using some adaptive algorithm based on observed auction data. In
the current practice of European day-ahead electricity markets, market
operators publish ex-post hourly aggregate supply and demand curves (see,
e.g., [15] and [16]). This information could be leveraged to enable each
bidder to retrospectively calculate her/his payoff for any chosen bid. The
possibility of a posteriori evaluating the payoff function is referred to as
_full information feedback_ in online learning, in contrast to the so-called
bandit feedback where the bidder evaluates her/his payoff only for the
submitted bid.

Remark: The considered setup abstracts several realistic considerations,
including time-varying demands, different sets of possible bids and bidders in
each round of the auction, different types of bid curves such as piece-wise
constant, and merit order dispatch approaches. With our simplifications, we
are able to benchmark the learning approach and provide a deeper understanding
of the potential of exploiting posterior information feedback in the market
under consideration.

##  III No-Regret Learning Approach

###  III-A Regret and connections to game

From the perspective of a learning bidder $\ell$, the problem of learning in a
repeated auction can be cast as an instance of online learning. In such a
setting, the learner is faced with a sequence of time-varying utility
functions and she/he aims to maximize her/his cumulative utility. Note that in
our auction setting, the time variations on the utility of bidder $\ell$ are
induced due to the bidders changing strategies at each auction round.

Recall $\mathcal{B}_{t}$ is the strategy profile chosen at auction round $t$
and let $\mathcal{B}_{-\ell,t}=$
$\left\\{b_{j}^{k_{j,t}}\right\\}_{j\in\mathcal{N}\backslash\\{\ell\\}}$. The
regret of bidder $\ell$ at time $T$ is defined as follows:

######  Definition 1 (Regret).

The regret of bidder $\ell$ at time $T$ is

$\displaystyle R_{\ell,T}=\max_{k\in\mathcal{K}_{\ell}}\
\sum_{t=1}^{T}u_{\ell}\left(\left\\{\mathcal{B}_{-\ell,t},b_{\ell}^{k}\right\\}\right)-\sum_{t=1}^{T}u_{\ell,t}.$
(3)

Observe that after $T$ rounds, $R_{\ell,T}$ measures the difference between
two quantities. The first term quantifies the maximum payoff bidder $\ell$
could have made had the bidder known the sequence of rival bids ahead of time,
and had the bidder chosen the best fixed strategy in $\mathcal{K}_{\ell}$. The
second term represents the cumulative utility derived from the strategy
employed by bidder $\ell$ throughout the rounds. An algorithm for bidder
$\ell$ is no-regret if $R_{\ell,T}/T\rightarrow 0$ as $T\rightarrow\infty$.

Note that while no-regret algorithms can be defined for the case in which the
sequence of received payoffs is arbitrary, they have additional connections to
equilibria in a game setting. In particular, if every participant bids
according to a no-regret algorithm, the empirical distribution of
participants’ bids converges to a coarse-correlated equilibrium, a relaxation
of the Nash equilibria, of the one-shot game [21].

###  III-B No-regret algorithm

We consider the case in which, after each round of the auction, the bidder can
accurately calculate her/his payoff had the bidder submitted alternative bids,
namely, the full information feedback. In practice, this can be done based on
the information published by the market operator after the day-ahead market
clearing. This feedback allows the bidder to observe a vector of rewards. To
exploit this information, we implement the Hedge algorithm, modeling the
bidder’s reward function based on historical rounds of the auction. This
algorithm enables the bidder to sequentially and adaptively learn how to
maximize her/his _cumulative_ payoff by strategically selecting her/his bids
across multiple auction rounds.

Algorithm 1 Hedge algorithm for bidder $\ell$

1:Inputs: Strategy set $\mathcal{K}_{\ell}$ with $|\mathcal{K}_{\ell}|=K$,
parameter $\eta$

2:Initialize weights:
$\mathbf{w}_{\ell,1}=\frac{1}{K}(1,\ldots,1)\in\mathbb{R}^{K}$

3:for auction rounds $t=1$ to $T$, do

4:Sample bidding action $k=(c_{\ell}^{k},d_{\ell}^{k})\in\mathcal{K}_{\ell}$
randomly, such that $k\sim\mathbf{w}_{\ell,t}$.

5:Observe auction outcomes $\\{x_{\ell,t},\lambda_{t},u_{\ell,t}\\}$ by
solving (2).

6:Compute the reward for every action using the aggregate supply curve
$\forall{i}\in\mathcal{K}_{\ell}$.

7:Update weights: For $\forall i\in\mathcal{K}_{\ell}$,
$[\mathbf{w}_{\ell,t+1}]_{i}=\frac{[\mathbf{w}_{\ell,t}]_{i}\exp(-\eta(1-[\mathbf{u}_{\ell,t}]_{i}))}{\sum\limits_{j=1}^{K}[\mathbf{w}_{\ell,t}]_{j}\exp(-\eta(1-[\mathbf{u}_{\ell,t}]_{j}))}.$

8:end for

This no-regret Hedge algorithm is standard in online learning with full
information feedback [22]. The steps are outlined in Algorithm 1. In Step $1$,
we provide the input data for bidder $\ell$, including her/his set of actions
(i.e., price bid options)
$\mathcal{K}_{\ell}=\\{(c_{\ell}^{1},d_{\ell}^{1}),...,(c_{\ell}^{K},d_{\ell}^{K})\\}$,
where $K=\left|\mathcal{K}_{\ell}\right|$ is the number of bidding options of
bidder $\ell$, and $\eta$ is the learning rate.

Step $2$ initializes the uniform weight vector for each bidding option
$\mathbf{w}_{\ell,1}\in\mathbb{R}^{K}$. The main iteration procedure of the
algorithm is shown in Steps $3$–$7$, which are repeated at each auction round
$t$. In Step $4$, given the weights $\mathbf{w}_{\ell,t}$, the algorithm
chooses a bidding option $k_{\ell}(t)=(c_{\ell}^{k},d_{\ell}^{k})$ sampled
from a distribution defined by the weights. This sampled action is the one
that the bidder submits to the auction in round $t$. The distribution is
initialized uniformly and is updated iteratively as will be described, based
on the received utilities. The rationale behind the algorithm is to update
weights $\mathbf{w}_{\ell,t}$ over rounds, such that bidding options with
potentially higher cumulative payoffs get greater weights.

Based on the bids played $k=(c_{\ell}^{k},d_{\ell}^{k}\in\mathcal{K_{\ell}})$
in Step $4$, in Step $5$, bidder $\ell$ observes the outcomes of the auction
for round $t$. These observations include the market-clearing price
$\lambda_{t}$ (public information), the bidder’s allocation $x_{\ell,t}$, and
the payoff $u_{\ell,t}$ corresponding to the chosen action
$k\in\mathcal{K}_{\ell}$ at round $t$. Assuming a posteriori information on
the rival bids, in Step $6$, bidder $\ell$ can calculate the payoffs that
could have been obtained if alternate bids had been submitted in round $t$. In
Step $7$, the algorithm updates the weights of each bidding strategy for
bidder $\ell$ based on their performance in round $t$. The update is based on
the exponential weight rule, so as to decrease the weight of poor performing
actions severely. This expression is then normalized by dividing it by the sum
of the updated weights of all actions, ensuring the weight vectors define a
probability distribution on the set of available actions of the bidder. The
idea is to incrementally shift the weights towards more successful actions
(with higher utilities), thereby nudging the bidder towards improved bidding
choice over time. With $\eta=\sqrt{\frac{8\log(K)}{T}}$ the Hedge algorithm is
known to obtain the best regret rates in online learning [22].

##  IV Benchmark

In the current literature of day-ahead electricity markets, sequential bidding
decisions are mostly addressed through bi-level programming, representing
single-leader single-follower or single-leader multi-follower Stackelberg
games. As our benchmark, we consider a multi-leader single-follower game. That
is, a set of strategic bidders aims to maximize their payoffs, while the
market operator seeks to maximize social welfare using an economic dispatch
model. In such hierarchical games, achieving a Nash equilibrium between the
leader(s) and the follower(s) is generally challenging [23, 24, 25]. As bi-
level programming is widely used in the literature for the bidding strategy
problem, we adopt it as a benchmark for the proposed no-regret learning
algorithm. In this setting, the assumption is that bidders know the market
operator’s response before their actions [4].

The general solution method for bi-level programming is to transform the
problem into an equivalent single-level optimization problem, which can also
be interpreted as a mathematical program with equilibrium constraints (MPEC)
[26], in a single-leader case. To achieve this, as it is common in the
literature [27], we replace the follower’s optimization problem with its
Karush–Kuhn–Tucker (KKT) optimality conditions. To accommodate multiple
leaders, the MPEC model can then be extended to a multi-leader single (multi)
followers problem, known as an EPEC [18]. Solving EPECs is challenging due to
their inherent non-convexity, coupled constraints, and the interplay of
multiple equilibria within a single problem. In the literature, iterative ways
exist to solve this class of problems such as the diagonalization method [19,
1]. This iterative approach is what we use, and will be outlined briefly.

Algorithm 2 Diagonalization algorithm to solve the EPEC

1:Initialize:
$\mathcal{B}_{0}=\\{c_{\ell}^{k},d_{\ell}^{k}\quad\forall\ell\\};$ maximum
number of iterations $\mathcal{I}$; convergence criterion $\varepsilon$.

2:for iterations $i=1$ to $\mathcal{I}$ do

3: for bidders $\ell=1$ to $|\mathcal{N}|$ do

4: Solve bidder $\ell$’s MIQP given rival strategies.

5: Update $\mathcal{B}_{\ell,i}$;

6: end for

7: if
$\|\mathcal{B}_{\ell,i}-\mathcal{B}_{\ell,i-1}\|\leq\epsilon\quad\forall\ell\in|\mathcal{N}|$
then

8: Convergence achieved, stop the algorithm.

9: end if

10: if $i=\mathcal{I}$ then

11: Convergence failed, stop the algorithm.

12: end if

13:end for

We consider $|\mathcal{N}|$ strategic bidders who can be leaders and take
their decisions first by anticipating the market operator’s response. We also
consider the market operator who clears the market as a common follower and
reacts to the decisions of the upper-level strategic bidders. The objective of
each bidder $\ell$ is to maximize her/his own payoff by anticipating how the
market operator reacts to her/his decision. On the other hand, the market
operator in the lower-level problem minimizes the social cost, i.e., the total
cost of meeting demand, accounting for the bid received by $\ell$ and her/his
rivals $-\ell$, subject to market constraints. The bi-level problem for each
bidder $\ell$ is expressed below

$\displaystyle\small\max_{\Xi^{\text{Upper}}}\quad$ $\displaystyle\lambda
x_{\ell}-\big{(}\frac{1}{2}c_{\ell}x_{\ell}^{2}+d_{\ell}x_{\ell}\big{)}$ (4a)
s.t. $\displaystyle
c^{k}_{\ell}\in\mathcal{C}_{\ell},d_{\ell}^{k}\in\mathcal{D}_{\ell}$ (4b)
$\displaystyle\min_{\Xi^{\text{Lower}}}\quad$
$\displaystyle\frac{1}{2}c_{\ell}^{k}x_{\ell}^{2}+d_{\ell}^{k}x_{\ell}+\sum_{-\ell}\big{(}\frac{1}{2}c_{-\ell}^{k}x_{-\ell}^{2}+d_{-\ell}^{k}x_{-\ell}\big{)}$
(4c) s.t. $\displaystyle x_{\ell}+\sum_{-\ell}x_{-\ell}=Q:\quad\lambda\quad$
(4d) $\displaystyle 0\leq
x_{\ell}\leq\bar{x}_{\ell}:\quad\underline{\mu}_{\ell},\bar{\mu}_{\ell}\quad$
(4e) $\displaystyle 0\leq
x_{-\ell}\leq\bar{x}_{-\ell}:\quad\underline{\mu}_{-\ell},\bar{\mu}_{-\ell}\quad\forall-\ell,$
(4f)

where $\Xi^{\text{Upper}}\\!=\\!\\{c_{\ell}^{k},d_{\ell}^{k}\\}$ is the set of
decision variables of the upper-level bidder,
$\Xi^{\text{Lower}}\\!=\\!\\{x_{\ell},x_{-\ell}\\}$ is the set of primal
variables for the lower-level problem, and
$\Xi^{\text{dual}}\\!=\\!\\{\lambda,\bar{\mu}_{\ell},\underline{\mu}_{\ell},\bar{\mu}_{-\ell},\underline{\mu}_{-\ell}\\}$
represents the corresponding dual variables of the lower-level problem. The
first term in the objective function (4a) of bidder $\ell$ is her/his revenue,
i.e., the market-clearing price $\lambda$ multiplied by the production
quantity $x_{\ell}$. The second term is the true quadratic cost function for
producing $x_{\ell}$. The upper-level constraint (4b) defines the potential
strategic actions of the bidder, such that the bidder can bid her/his cost
function with coefficients $c^{k}_{\ell}$ and $d^{k}_{\ell}$, which are not
necessarily identical to true coefficients $c_{\ell}$ and $d_{\ell}$. The
lower-level objective function (4c) minimizes the social cost from the
perspective of the market operator. By this, we implicitly assume that bidder
$\ell$, who solves bi-level program (4), perfectly knows the bidding decisions
of her/his rivals $-\ell$.

Notice that from the perspective of a single bidder, an optimizer of the above
optimization problem would correspond to her/his _best response_ action.
Computing the best response action in general is intractable. Thus, this bi-
level problem can be reformulated as an MPEC by replacing (4c)-(4f) by its KKT
conditions. After reformulations proposed by [4], the problem is cast as a
mixed-integer quadratic program (MIQP). The resulting EPEC can be solved using
a diagonalization algorithm, an iterative process whose details are given in
Algorithm 2. Note that the iterations imply that bidders one by one solve (or
approximate) their best-response actions. In the case in which they solve the
MIQP problem exactly in each iteration, and that the algorithm converges, it
follows that the resulting bidding strategy is a Nash equilibrium.

Table I: Bidding options. While {$c_{1},d_{1}$} represents the true cost
function of each bidder, {$c_{o},d_{o}$}, $\forall{o}=\\{2,...,10\\}$ are
her/his available nine options for placing a strategic bid.

Bidder | $c_{1}$ | $c_{2}$ | $c_{3}$ | $c_{4}$ | $c_{5}$ | $c_{6}$ | $c_{7}$ | $c_{8}$ | $c_{9}$ | $c_{10}$ | $d_{1}$ | $d_{2}$ | $d_{3}$ | $d_{4}$ | $d_{5}$ | $d_{6}$ | $d_{7}$ | $d_{8}$ | $d_{9}$ | $d_{10}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$1$ | 0.070 | 0.080 | 0.090 | 0.100 | 0.120 | 0.075 | 0.085 | 0.095 | 0.150 | 0.170 | 9 | 10 | 11.5 | 14 | 13 | 10 | 12 | 13 | 15 | 11
$2$ | 0.020 | 0.050 | 0.060 | 0.150 | 0.250 | 0.025 | 0.150 | 0.900 | 0.130 | 0.310 | 10 | 15 | 12 | 17 | 11 | 12 | 14 | 13 | 16 | 14
$3$ | 0.030 | 0.040 | 0.060 | 0.120 | 0.140 | 0.095 | 0.080 | 0.090 | 0.210 | 0.270 | 12 | 14 | 13 | 16 | 15 | 13 | 12 | 15 | 14 | 12
$4$ | 0.008 | 0.010 | 0.075 | 0.240 | 0.310 | 0.080 | 0.090 | 0.050 | 0.110 | 0.140 | 12 | 14 | 17 | 15 | 17 | 14 | 11 | 17 | 15 | 12
$5$ | 0.010 | 0.090 | 0.100 | 0.210 | 0.970 | 0.020 | 0.130 | 0.075 | 0.190 | 0.095 | 11 | 13 | 11 | 17 | 20 | 12 | 11 | 15 | 17 | 20

Remark: In cases where the number of bidding options for each bidder is small,
an alternative approach to solving (4) would be an enumeration process. For
instance, consider a scenario with three bidding options for every bidder. In
every iteration, with fixed rival bids, each bidder $\ell$ solves three
quadratic programs (QPs) (2) (which could be solved in parallel)—one for each
bidding option, calculates her/his payoff and chooses the best response
accordingly. This approach requires each bidder to solve one QP per bid option
in every iteration, thus bypassing the need to solve a single MIQP in every
iteration of Algorithm 2.

##  V Numerical Study

We consider five bidders, namely Bidders $1$ to $5$, participating in an hour-
ahead electricity market. Since no inter-temporal constraint is enforced in
the market-clearing optimization problem (2), the auction for different hours
can be conducted separately and independently. This enables us to focus on the
bidding strategy for a specific hour and base our analysis on the hour-ahead
market [28]. We fix the total load $Q$ to be $1148.4$ MW, which is the total
demand in the Danish bidding zone DK$1$ reported by Nord Pool for hour
$9$:$00$-$10$:$00$ on August $27$, $2023$. Table I presents the action sets.
For each bidder, we consider ten options for coefficients {$c$, $d$}, where
{$c_{1},d_{1}$} refers to her/his true quadratic cost function. The remaining
nine bidding options {$c_{2},d_{2}$} to {$c_{10},d_{10}$} let the bidder place
a strategic price bid. These nine options are generated by gradually
increasing the quadratic and linear terms of the cost function. The capacity
of bidders is set identically to $700$ MW. All source codes are provided in
[29].

###  V-A Cases Based on Bidding Approaches

We first start by defining two benchmark cases and then introduce six
learning-oriented cases, resulting in eight cases. The two benchmark cases are
as follows:

  * •

Case $a$: Best Response: All Bidders $1$ to $5$ have perfect information of
rival bids and use this to compute their best response. For this case, we use
the diagonalization approach in Algorithm 2.

  * •

Case $b$: Trustful: All Bidders $1$ to $5$ submit their true generation cost
{$c_{1},d_{1}$}. This results in a perfectly competitive market, for which we
solve the market-clearing optimization problem (2) given true costs.

We now define six learning-oriented cases. In all these six cases, Bidder $5$
employs either the proposed Hedge algorithm 1 (cases $c$, $e$, and $g$) to
learn her/his strategic price bids over auction rounds or bids randomly (cases
$d$, $f$ and $h$). Conversely, rival bidders (Bidders $1$-$4$) utilize their
true cost bids (cases $c$ and $d$), learn using the Hedge algorithm (cases $e$
and $f$), and bid randomly (cases $g$ and $h$):

  * •

Case $c$: Trustful vs Hedge: We assume Bidders $1$-$4$ are naive and provide
their true cost bids in every round, while Bidder $5$ learns by employing the
Hedge algorithm. That means, only Bidder 5 is allowed to choose her/his bids
from the ten options provided in Table I and learns how to bid strategically
based on the Hedge algorithm.

  * •

Case $d$: Trustful vs Random: In this case, instead of using the Hedge
algorithm, Bidder 5 bids randomly, while rival bidders consistently submit
their true costs.

  * •

$\rm{Case}$ $e$: Hedge vs Hedge: Both Bidder 5 and her/his rival bidders
(Bidders $1$-$4$) utilize the Hedge algorithm to simultaneously learn their
strategic bids.

  * •

Case $f$: Hedge vs Random: Bidders $1$-$4$ employ the Hedge algorithm, while
Bidder $5$ bids randomly.

  * •

Case $g$: Random vs Hedge: Bidders $1$-$4$ choose their bids uniformly at
random in each round of the auction while Bidder $5$ exploits the proposed
Hedge algorithm.

  * •

Case $h$: Random vs Random: All bidders sample bids uniformly at random in
each round of the auction from bidding options in Table I.

Table II: Social cost [€] for various cases in increasing order.

Cases | Trustful | Trustful vs Hedge | Trustful vs Random | Best Response | Hedge vs Hedge | Hedge vs Random | Random vs Hedge | Random vs Random
---|---|---|---|---|---|---|---|---
Social cost [€] | $19,419$ | $21,085$ | $21,332$ | $24,408$ | $26,162$ | $28,799$ | $29,292$ | $35,763$

For all the six learning cases, we use the optimization problem (2) to
simulate the auction environment and determine the bidders’ bidding strategies
over $T=200$ rounds of the auction. Consequently, we run these $200$ auction
rounds $15$ times to calculate the average outcomes across different runs.
Hereafter, solid lines in our figures show the average values over $15$ runs,
whereas shaded areas cover the mean value plus and minus the standard
deviation. Unlike the six learning-oriented cases ($c,d,e,f,g,h$), the
benchmarks (cases $a,b$) are not built upon data gathered from learning
algorithms, since bidders either bid their true costs or submit their a priori
computed best response bids.

###  V-B Market Outcomes: System-level Results

For all eight cases, we report the social cost, denoting the social welfare,
and the market-clearing price. Recall that the social cost is the optimal
value of objective function (2a), and the market-clearing price is the optimal
value of the dual variable associated with (2b). We note that the EPEC problem
solved by the diagonalization Algorithm 2 converged after four iterations with
$\epsilon=0.0004$. Furthermore, through enumeration, we verified that the
converged bid profiles were a Nash equilibrium of the game.

Table II provides the social cost of the eight cases, arranged in ascending
order. We discuss the results of this table together with those plotted in
Figure 1, illustrating the evolution of social cost (left panel) and that of
the market-clearing price (right panel) over $200$ auction rounds. The social
cost results in Table II correspond to those in the last round of the upper
plot of Figure 1. As expected, the social cost is lower in the $\rm{Trustful}$
than the cases in which Bidder $5$ is acting strategically through the Hedge
algorithm or randomly. Compared to the $\rm{Best}$ $\rm{Response}$ where
bidders are using their Nash equilibrium strategy, the social cost increases
when one or all bidders adopt a learning approach. Furthermore, with the
increase of randomness (zero, one, and all bidders applying $\rm{Random}$
strategy), the social cost increases.

Given that implementing the Hedge algorithm and computing the best response
strategies require progressively more information than playing a random
strategy, the results above quantify the value of information to the bidders,
and suggest that having more information increases market efficiency.

Figure 1: Evolution of social cost (left panel) and market-clearing price
(right panel) over $200$ auction rounds.

Figure 2: Evolution of average regret (left panel) and average payoff (right
panel) for Bidder $5$ over $200$ auction rounds.

###  V-C Market Outcomes: Individual Results

Figure 2, left panel, illustrates the average regret of Bidder $5$ in the six
learning-oriented cases. As expected, when rivals stay with a given approach
(Hedge algorithm, or randomization), Bidder $5$’s regret decreases when
employing the Hedge algorithm in comparison to playing randomly.

Figure 2, right panel, illustrates the payoff of Bidder $5$. From this figure,
we make three observations below.

First, if rivals bid randomly, Bidder $5$ earns a higher payoff, regardless of
her/his bidding strategy, compared to other cases where rivals do not submit
randomly. When Bidder $5$ employs the Hedge algorithm while rival bidders act
randomly (the red curve), her/his payoffs are higher than if she/he were to
also use a random strategy (yellow curve).

Second, if rivals employ the Hedge algorithm, it is slightly better for Bidder
$5$ to use the Hedge algorithm too (blue curve) compared to bidding randomly
(gray curve). While this is expected, it is interesting to observe that in
these two cases, the average regret for Bidder $5$ is decreasing despite the
fact that her/his payoffs are decreasing. While the decreasing trend of regret
by the learning algorithm is expected, in this auction setting higher regret
does not mean a lower payoff (see orange line) and lower regret does not mean
a higher payoff (see blue). Hence, this suggests that regret may not be the
best benchmark for the problem at hand and it is more important to directly
measure the bidders’ utilities and social costs.

Third, if rival bidders submit bids based on their $\rm{Trustful}$ bids,
Bidder $5$ earns a higher payoff by bidding her/his true costs (red dashed
line) compared to when Bidder $5$ employs the Hedge algorithm (green) or bids
randomly (brown). The results show that Bidder $5$’s payoff, when using the
Hedge algorithm (green), while her/his rivals consistently use their
$\rm{Trustful}$ bidding, slowly converges to all $\rm{Trustful}$ bidding (red
dashed line). This particular observation is due to the fact that by
enumeration, we could verify that Bidder $5$’s best response to others’
playing their true costs is to play her/his true cost.

##  VI Conclusion

This paper presented a no-regret learning approach to analyze the extent of
market power that bidders can exert by learning their strategic bids in the
day-ahead electricity market. To this end, the Hedge algorithm, which is known
to obtain the best regret rates in online learning, is utilized. The proposed
learning approach shows promising results by relying solely on the posterior
information published by the market operator, without making any assumptions
about rival bids a priori in the bidding decision making process.

From the system perspective, we observed that electricity markets might be
exposed to greater market power compared to the case in which bidders provide
their best responses, as classical game theoretic approaches would result in,
but use a learning approach. We further observed the less information bidders
have to choose their actions, the higher the social welfare will be.

From a bidder’s perspective, we illustrate the decreasing regret of a single
bidder over auction rounds. However, bidders’ payoffs by using the Hedge
algorithm highly depend on rivals’ strategies. In particular, unless other
bidders are playing a fixed strategy, a given bidder has a higher payoff if
she/he uses the proposed online learning approach. However, if other bidders
play a fixed strategy over auction rounds (submitting either their true costs
or their best response strategies), then there is clearly a single fixed
strategy that outperforms no-regret learning and in this case, the no-regret
algorithm converged to the corresponding strategy. Furthermore, from our case
study, we observed that a bidder earns a higher payoff when she/he is the only
one who employs a learning approach, and rival bidders submit their bids
randomly, compared to the cases where everyone is learning simultaneously.
Lastly, we observed that a decrease in regret in one player may not imply an
increase in the payoffs of the player and vice versa. This motivates use of
other metrics such as policy regret to address learning in repeated
electricity market auctions.

Our formulations and case studies serve as a preliminary evaluation of the
potential of the online learning approach in electricity markets. There are
several directions for future research. First, a large number of simulations
with different bidders and market constraints are needed to generalize the
observations made in our case study to realistic market settings. Second, in
order to simulate the full information feedback considered in our case study,
the aggregate supply curves published a posteriori need to be utilized. Third,
in a realistic market setting, contextual information such as weather
forecasts could complement a bidder’s information and thus, it is important to
leverage this information in online learning to improve a bidder’s
performance. Fourth, the auction environment can have multiple trading stages
such as spot and balancing markets.

## Acknowledgments

Arega Getaneh Abate has received funding from the European Union’s Horizon
$2020$ research and innovation program under the Marie Skłodowska-Curie grant
agreement No $899987$. Maryam Kamgarpour gratefully acknowledges the project
funding under “UrbanTwin” with the financial support of the ETH-Domain Joint
Initiative program in the Strategic Area Energy, Climate and Sustainable
Environment.

## References

  * [1] X. Hu and D. Ralph, “Using EPECs to model bilevel games in restructured electricity markets with locational prices,” _Operations Research_ , vol. 55, no. 5, pp. 809–827, 2007.
  * [2] B. F. Hobbs, C. B. Metzler, and J. S. Pang, “Strategic gaming analysis for electric power systems: An MPEC approach,” _IEEE Transactions on Power Systems_ , vol. 15, no. 2, pp. 638–645, 2000.
  * [3] S. Martín, Y. Smeers, and J. A. Aguado, “A stochastic two settlement equilibrium model for electricity markets with wind generation,” _IEEE Transactions on Power Systems_ , vol. 30, no. 1, pp. 233–245, 2015.
  * [4] C. Ruiz and A. J. Conejo, “Pool strategy of a producer with endogenous formation of locational marginal prices,” _IEEE Transactions on Power Systems_ , vol. 24, no. 4, pp. 1855–1866, 2009.
  * [5] P. Auer, N. Cesa-Bianchi, and P. Fischer, “Finite-time analysis of the multiarmed bandit problem,” _Machine Learning_ , vol. 47, pp. 235–256, 2002\.
  * [6] Y. Wang, B. Zhang, J. Ma, and Q. Jin, “Earning while learning: An adversarial multi-armed bandit based real-time bidding scheme in deregulated electricity market,” _IEEE Transactions on Network Science and Engineering_ , vol. 9, no. 6, pp. 3991–4000, 2022.
  * [7] J. H. Nachbar, “Prediction, optimization, and learning in repeated games,” _Econometrica_ , vol. 65, no. 2, pp. 275–309, 1997.
  * [8] F. Orabona, “A modern introduction to online learning,” _arXiv preprint arXiv:1912.13213_ , 2022.
  * [9] S. Shalev-Shwartz _et al._ , “Online learning and online convex optimization,” _Foundations and Trends® in Machine Learning_ , vol. 4, no. 2, pp. 107–194, 2012.
  * [10] Y. Ye, D. Qiu, M. Sun, D. Papadaskalopoulos, and G. Strbac, “Deep reinforcement learning for strategic bidding in electricity markets,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 2, pp. 1343–1355, 2020.
  * [11] Y. Du, F. Li, H. Zandi, and Y. Xue, “Approximating Nash equilibrium in day-ahead electricity market bidding with multi-agent deep reinforcement learning,” _Journal of Modern Power Systems and Clean Energy_ , vol. 9, no. 3, pp. 534–544, 2021.
  * [12] O. Karaca, P. G. Sessa, A. Leidi, and M. Kamgarpour, “No-regret learning from partially observed data in repeated auctions,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 14–19, 2020.
  * [13] Z. Zhao, A. L. Liu, and Y. Chen, “Electricity demand response under real-time pricing: A multi-armed bandit game,” in _Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC)_ , 2018, pp. 748–756.
  * [14] S. Baltaoglu, L. Tong, and Q. Zhao, “Algorithmic bidding for virtual trading in electricity markets,” _IEEE Transactions on Power Systems_ , vol. 34, no. 1, pp. 535–543, 2018.
  * [15] European Energy Exchange. (2023) EPEX Spot aggregate curves. [Online]. Available: https://www.epexspot.com/en/market-data
  * [16] Nord Pool. (2024) Aggregate market data: Nord Pool day-ahead aggregate curves. [Online]. Available: https://data.nordpoolgroup.com/auction/day-ahead/aggregated-bid-curves?deliveryDate=latest&bidCurveArea=BALTIC
  * [17] Y. Freundand and R. E. Schapire, “Adaptive game playing using multiplicative weights,” _Games and Economic Behavior_ , vol. 29, pp. 79–103, 1999.
  * [18] A. Ehrenmann, _Equilibrium problems with equilibrium constraints and their application to electricity markets_. Cambridge University, UK, 2004, Ph.D. Dissertation.
  * [19] B. H. Ahn and W. W. Hogan, “On convergence of the PIES algorithm for computing equilibria,” _Operations Research_ , vol. 30, pp. 281–300, 1982\.
  * [20] P. G. Sessa, I. Bogunovic, M. Kamgarpour, and A. Krause, “No-regret learning in unknown games with correlated payoffs,” in _Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 32, 2019.
  * [21] N. Cesa-Bianchi and G. Lugosi, _Prediction, learning, and games_. Cambridge University Press, 2006.
  * [22] Y. Freund and R. E. Schapire, “A decision-theoretic generalization of on-line learning and an application to boosting,” _Journal of Computer and System Sciences_ , vol. 55, no. 1, pp. 119–139, 1997.
  * [23] S. Leyffer and T. Munson, “Solving multi-leader–common-follower games,” _Optimisation Methods & Software_, vol. 25, no. 4, pp. 601–623, 2010.
  * [24] V. Conitzer and D. Korzhyk, “Commitment to correlated strategies,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 25, no. 1, 2011, pp. 632–637.
  * [25] N. Basilico, S. Coniglio, N. Gatti, A. Marchesi _et al._ , “Bilevel programming approaches to the computation of optimistic and pessimistic single-leader-multi-follower equilibria,” _Leibniz International Proceedings in Informatics_ , vol. 75, pp. 1–14, 2017.
  * [26] Z. Q. Luo, J. S. Pang, and D. Ralph, _Mathematical programs with equilibrium constraints_. Cambridge University Press, 1996.
  * [27] D. Pozo, E. Sauma, and J. Contreras, “Basic theoretical foundations and insights on bilevel models and their applications to power systems,” _Annals of Operations Research_ , vol. 254, pp. 303–334, 2017.
  * [28] R. Chen, I. C. Paschalidis, M. C. Caramanis, and P. Andrianesis, “Learning from past bids to participate strategically in day-ahead electricity markets,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 5, pp. 5794–5806, 2019.
  * [29] A. Getaneh Abate, D. Majdi, J. Kazempour, and M. Kamgarpour. (2023) Electronic companion: Learning to bid in forward electricity markets using a no-regret algorithm. [Online]. Available: https://github.com/AregaGetaneh/Day_ahead-bidding-strategy


# Transport and optical properties of the chiral semiconductor Ag3AuSe2

Juyeon Won  Department of Materials Science and Engineering and Materials
Research Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL
61801, USA  Soyeun Kim  Department of Physics, University of Illinois at
Urbana-Champaign, Urbana, 61801 IL, USA  Materials Research Laboratory,
University of Illinois at Urbana-Champaign, Urbana, 61801 IL, USA  Martin
Gutierrez-Amigo  Department of Physics, University of the Basque Country
(UPV/EHU), Apartado 644, 48080 Bilbao, Spain  Centro de F’isica de Materiales
(CSIC-UPV/EHU), Manuel de Lardizabal Pasealekua 5, 20018 Donostia/San
Sebastián, Spain  Simon Bettler  Department of Physics, University of Illinois
at Urbana-Champaign, Urbana, 61801 IL, USA  Materials Research Laboratory,
University of Illinois at Urbana-Champaign, Urbana, 61801 IL, USA  Bumjoo Lee
Center for Correlated Electron Systems, Institute for Basic Science, Seoul
08826, Republic of Korea  Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Republic of Korea  Jaeseok Son  Center for Correlated
Electron Systems, Institute for Basic Science, Seoul 08826, Republic of Korea
Department of Physics and Astronomy, Seoul National University, Seoul 08826,
Republic of Korea  Tae Won Noh  Center for Correlated Electron Systems,
Institute for Basic Science, Seoul 08826, Republic of Korea  Department of
Physics and Astronomy, Seoul National University, Seoul 08826, Republic of
Korea  Ion Errea  Donostia International Physics Center, P. Manuel de
Lardizabal 4, 20018 Donostia-San Sebastian, Spain  Centro de F’isica de
Materiales (CSIC-UPV/EHU), Manuel de Lardizabal Pasealekua 5, 20018
Donostia/San Sebastián, Spain  Fisika Aplikatua 1 Saila, Gipuzkoako
Ingeniaritza Eskola, University of the Basque Country (UPV/EHU), Europa Plaza
1, 20018 Donostia/San Sebastián, Spain  Maia G. Vergniory  Donostia
International Physics Center, P. Manuel de Lardizabal 4, 20018 Donostia-San
Sebastian, Spain  Max Planck Institute for Chemical Physics of Solids , 01187
Dresden, Germany  Peter Abbamonte  Department of Physics, University of
Illinois at Urbana-Champaign, Urbana, 61801 IL, USA  Materials Research
Laboratory, University of Illinois at Urbana-Champaign, Urbana, 61801 IL, USA
Fahad Mahmood  Department of Physics, University of Illinois at Urbana-
Champaign, Urbana, 61801 IL, USA  Materials Research Laboratory, University of
Illinois at Urbana-Champaign, Urbana, 61801 IL, USA  Daniel P. Shoemaker
<EMAIL_ADDRESS>Department of Materials Science and Engineering and
Materials Research Laboratory, University of Illinois at Urbana-Champaign,
Urbana, IL 61801, USA

###### Abstract

Previous band structure calculations predicted Ag3AuSe2 to be a semiconductor
with a band gap of approximately 1 eV. Here, we report single crystal growth
of Ag3AuSe2 and its transport and optical properties. Single crystals of
Ag3AuSe2 were synthesized by slow-cooling from the melt, and grain sizes were
confirmed to be greater than 2 mm using electron backscatter diffraction.
Optical and transport measurements reveal that Ag3AuSe2 is a highly resistive
semiconductor with a band gap of and activation energy around 0.3 eV. Our
first-principles calculations show that the experimentally-determined band gap
lies between the predicted band gaps from GGA and hybrid functionals. We
predict band inversion to be possible by applying tensile strain. The
sensitivity of the gap to Ag/Au ordering, chemical substitution, and heat
treatment merit further investigation.

##  I Introduction

Ag3AuSe2 was first reported reported in 1971 as the mineral fischesserite from
Předbořice, Czech Republic, by Johan, et al.Johan  _et al._ (1971) It was the
first gold-containing selenide reported, and its crystal structure was
identified as isostructural to the chiral mineral petzite Ag3AuTe2 in space
group $I4_{1}32$. Moderate reflectivity in the visible range and a Moh’s
hardness of two were corroborated on another natural sample in 2004.Bindi and
Cipriani (2004) Initial work on the silver-gold chalcogenides used
differential thermal analysis (DTA) to show that they undergo order-disorder
transitions to body-centered cubic Ag2S phases upon heating, with Ag3AuSe2
transforming at 270∘C.Tavernier  _et al._ (1967) The tellurite Ag3AuTe2
undergoes two transitions at 220∘C and 320∘C, with the intermediate crystal
structure unknown, even though high-temperature X-ray diffraction (XRD) data
were collected.Smit  _et al._ (1970) The sulfide uytenbogaardtite, on the
other hand, has a more complex distorted tetragonal $P4_{1}22$ structure at
room temperatureBarton  _et al._ (1978) and a transition at 185∘C.Tavernier
_et al._ (1967)

The first studies on synthetic fischesserite by Smit, et al.Smit  _et al._
(1970) focused on confirming the previously-reported phase transitions by DTA
and XRD.

WiegersWiegers (1976) investigated the electronic and ionic conductivity of
Ag3AuSe2 with a claimed band gap from diffuse reflectance spectroscopy of 0.9
eV, but the samples were processed at 450∘C and no evidence of their purity
(XRD or otherwise) was provided. The remaining work on synthetic Ag3AuSe2
include a Mössbauer study showing that the ordered silver-gold chalcogenides
contain monovalent Au+ and stronger bonding than the binary
chalcogenides,Wagner  _et al._ (1992) and electromotive force measurements of
thermodynamic equations of state that confirmed the stability of these
compounds but did not precisely agree with the order-disorder transition
temperatures. Osadchii and Echmaeva ; Feng and Taskinen (2014)

Ag3AuSe2 and Ag3AuTe2 are chiral and isostructural at room temperature with
the space group of $I4_{1}32$. The chirality and possible narrow band gap (or
metallicity) of these materials gives rise to possible non-reciprocal and
topological phenomena. In this report, we investigate the room-temperature
phase of Ag3AuSe2.

More recently, the chirality and possible narrow band gap (or possibility of
band inversion and topological phenomena) has led renewed computational focus
on Ag3AuSe2. Faizan used density functional theory (DFT) to evaluate the
promise of Ag3AuSe2 for thermoelectric applications, and found the band gap to
be about 1.0 eV using the generalized gradient approximation (GGA).Faizan  _et
al._ (2016, 2017) However, published resistivity measurements (with the
aforementioned question of sample purity)Wiegers (1976) gave a lower
activation energy of around 0.5 eV. Calculations using the local density
approximation (LDA), by Fang, et al. are expected to underestimate the band
gap and found a gap of about 0.2 eV. Fang  _et al._ (2002) Bradlyn and
collaboratorsBradlyn  _et al._ (2016) predicted it to host multifold
topological crossings. Since 2016 multifold crossing in chiral topological
semimetals have attracted great attention owing to their large Chern numbers
Bradlyn  _et al._ (2016); Chang  _et al._ (2017), long Fermi arcs Schröter
_et al._ (2019, 2020); Sanchez  _et al._ (2019) and optical quantized
responses in terms of Chern numbers.Flicker  _et al._ (2018) Better
understanding the band separation in Ag3AuSe2 is required to determine whether
closing this narrow gap could lead to topological insulating behavior was
examined by Sanchez-MartinezSánchez-Martínez  _et al._ (2019). DFT-GGA
calculations showed that decrease in lattice parameters result in band gap
widening. In addition, they found that compression of 2% yields a meV-range
band gap, which may be suitable for dark matter detector applications. All
previous computational studies have pointed toward a direct band gap, but no
single crystal measurements have been made, nor have diffraction data
confirmed the phase purity of synthetic samples. Here, we report single
crystal synthesis, transport and optical properties of Ag3AuSe2 and a compare
computed electronic structures of the experimental and computationally-relaxed
phases.

##  II Methods

Single crystals of Ag3AuSe2 with grains larger than 2mm were synthesized. Ag
(99.9% metals basis), Au (99.98% metals basis), and Se (99.9999% metals basis)
powders were mixed in a stoichiometric ratio inside an Ar filled glove box and
sealed in an evacuated silica tube with 7 mm inner diameter. The ampule was
placed in a muffle furnace and heated to 750∘C at 1.5∘C/min, annealed for 24
h, cooled to 700∘C at 0.03∘/min, then cooled to room temperature over 24 h.
The resulting product was a solid ingot with metallic luster. A portion of the
ingot was crushed and analyzed using powder X-ray diffraction (PXRD) on a
Bruker D8 diffractometer with a Cu X-ray source in reflection geometry. The
crystal structure was refined using the Rietveld method with GSAS-IIToby and
Von Dreele (2013).

Ingots of Ag3AuSe2were manually polished with SiC grinding papers and alumina
slurries down to 0.3 $\mu$m. Electron backscatter diffraction (EBSD) was
performed using a Thermo Scios 2 DualBeam scanning electron microscope (SEM).
We found that Ag3AuSe2 surface oxidizes within minutes of exposure in air,
which made it difficult to obtain clear Kikuchi patterns. To avoid this issue,
mechanically polished samples were milled with 7keV Ar+ ions (Gatan PECS II)
immediately before transferring into the SEM chamber. Resistivity measurements
were carried out using the 2-point contact probe method in a Quantum Design
PPMS DynaCool. A rod of Ag3AuSe2 (1.5 $\times$ 2.8 $\times$ 0.7 mm) was
mounted with Kapton tape and 2 gold wire contacts were made by silver epoxy.

Spectroscopic ellipsometry measurements were performed using M-2000 and IR-
VASE ellipsometers (J. A. Woollam Co.). Optical conductivity was measured in
the UV-visible (0.74-6 eV) and mid-near IR (70-560 meV) energy ranges. Both
ellipsometeers were calibrated with a SiO2 (25 nm)/Si wafer prior to scanning.
The experimental $\Psi$ and $\Delta$ parameters were obtained at $70^{\circ}$
incident angle, then converted into the real and imaginary parts of the
optical constants to satisfy the Kramers-Kronig relation. The optical
conductivity determined at $60^{\circ}$ and $80^{\circ}$ did not change
significantly, as expected for a cubic material without optical anisotropy.

First-principles density functional theory (DFT) simulations were performed
using Quantum Espresso (QE) Giannozzi  _et al._ (2017) and the Vienna Ab
Initio Simulation Package (VASP) Kresse and Furthmüller (1996, 1996) with
projector-augmented wave pseudopotentials. We used a variety of
approximations: (1) The local density approximation (LDA), (2) the generalized
gradient approximation with the Perdew-Burke-Ernzerhof parameterization (PBE)
Perdew  _et al._ (1996), (3) the modified Becke-Johnson (mBJ) method Tran and
Blaha (2009) and (4) the Heyd–Scuseria–Ernzerhof (HSE) approximation with the
HSE06 parametrization Heyd  _et al._ (2003). Structural relaxations were
performed with QE using PBE and a plane-wave basis with a kinetic energy
cutoff of 100 Ry and a 7$\times$7$\times$7 $k$-mesh, which is stopped when
pressures are below 0.1 kbar. Band calculations were performed using VASP in
the presence of spin-orbit coupling, an energy cutoff of 500 eV and a
10$\times$10$\times$10 (5$\times$5$\times$5) $k$-mesh for LDA, PBE and mBJ
(HSE06).

##  III Results and Discussion

###  III.1 Growth and structure

Figure 1:  The high-temperature $\alpha$-Ag2Se structure with partial Ag
occupancy is shown in (a). Ag3AuSe2 is formed upon cooling by ordering Ag
displacements and Au substitution on the $\alpha$-Ag2Se to form a unit cell
that is 8 times the volume. The full unit cell of Ag3AuSe2 is shown in (b).
(c) shows the most apparent chiral $4_{1}$ chain viewed down $<$100$>$.

The structure of Ag3AuSe2 seems complex, but can be understood as a set of
modification from the high-temperature $\alpha$-Ag2Se structure. Fig.1 (a)
shows the body-centered-cubic $Im\overline{3}m$ $\alpha$-Ag2Se structure in
which Se2- ions sit at lattice points, while Ag+ ions are statistically
distributed over 4 positions on and surrounding the face centers. Upon
cooling, the metal cations in Ag3AuSe2 order into one of the available
displaced sites, while also distributing the Au+ evenly throughout the cell.
This new distorted cell must be doubled in each direction to accommodate the
Ag/Au ordering. The resulting structure of Ag3AuSe2 (I4132 ($a=9.96$ Å $Z=8$)
is shown in Fig. 1(b), viewed down the $a$ axis. The chiral $4_{1}$ ordering
of the cell is best viewed down the $<$100$>$ direction, forming a spiral of
Ag and Au in Fig. 1(c).

Figure 2:  Rietveld refinement to X-ray powder diffraction pattern of a
crushed ingot of Ag3AuSe2 measured in reflection geometry.

Crushed ingots of Ag3AuSe2 gave PXRD data shown in Fig. 2, which clearly
indicate phase purity and the 110 peak at $2\theta=12.6^{\circ}$ that
indicates the doubled $I4_{1}32$ ordering described in Fig. 1. As expected for
a cubic material and first noted by Johan,Johan  _et al._ (1971) Ag3AuSe2 does
not cleave leaving large surfaces, with a cross section of the fractured ingot
shown in Fig. 3(a). The material forms large grains after cooling from the
melt, so determining the grain structure from polished surfaces can be a
challenge because there may be only one or few grains in cross section, such
as the polished face in Fig. 3(b). Therefore, EBSD is necessary to confirm the
orientation of various grains in the face. A typical cross section under EBSD
mapping is shown in Fig. 3(c).

Figure 3:  Fractured and polished (1200 grit SiC) surfaces of Ag3AuSe2 are
shown in (a) and (b), respectively. EBSD grain mapping in (c) indicates grain
sizes larger than 2 mm. The green and purple regions have approximate
$<$101$>$ and $<$112$>$ directions normal to the sample surface.

###  III.2 Electronic properties

Figure 4:  Two-point electrical resistivity (a) of Ag3AuSe2 shows
semiconducting behavior, while a linear fit to the Arrhenius plot in (b) gives
an activation energy $E_{a}=0.32$ eV.

The electrical resistivity and corresponding Arrhenius plot of Ag3AuSe2 are
shown in Fig.4. We find that Ag3AuSe2 is highly resistive and exhibits a
semiconductor-like temperature dependence, with an Arrhenius activation energy
$E_{a}=0.32$ eV, which is close to the previously reported value of 0.5 eV
from Wiegers.Wiegers (1976) Given the different microstructures (crystal
versus powder) and annealing temperatures (750 versus 450∘C) between our study
and that work, there may be significant differences between the defect
populations that scatter carriers in the two resistivity measurements, but the
activation energies are the same order of magnitude. The resistivity of
Ag3AuSe2 is orders of magnitude greater than that of Ag3AuTe2 in the same
structure type.Young  _et al._ (2000)

Figure 5:  Optical conductivity $\sigma_{1}$ (a) and the square of optical
absorption coefficient $\alpha$ (b) for Ag3AuSe2 are determined by
spectroscopic ellipsometry. In (b), the direct band gap edge $E_{g}$ of 0.27
($\pm$0.05) eV shown, which is estimated using the linear relation of
$\alpha^{2}$ and ($E-E_{g}$) of semiconductors Fox and Fox (2001).

In order to investigate the electronic band structure of Ag3AuSe2, we
performed spectroscopic ellipsometry measurements. Fig. 5(a) shows the
experimental optical conductivity $\sigma_{1}$ of Ag3AuSe2 measured at room
temperature. No sign of free carrier response (intraband transition) was
captured in the $\sigma_{1}$ down to the lowest frequency measured, consistent
with the measured high resistivity and finite gap predicted from previous
calculations. Fang  _et al._ (2002); Faizan  _et al._ (2016, 2017) The band
gap edge of Ag3AuSe2 was extracted from the squared optical absorption
coefficient $\alpha^{2}$ as in Fig. 5(b). From the Fermi’s golden rule and the
joint density of states for parabolic bands, a frequency dependence of
absorption near a direct band gap ($E_{\text{G}}$) is represented as Fox and
Fox (2001):

$\alpha\propto\sqrt{\omega-E_{\text{G}}},$ (1)

for $\omega>E_{\text{G}}$. Previous band structure studies of Ag3AuSe2, and
our own, support parabolic curvature of the band edges of the direct gap at
the Brillouin zone center $\Gamma$.Fang  _et al._ (2002); Wiegers (1976);
Faizan  _et al._ (2016) Using the above relation, the direct band gap of
Ag3AuSe2 is estimated to be 0.27 eV, as illustrated in Fig. 5(b).
Uncertainties in the exact onset of the rise in $alpha^{2}$ and the imperfect
geometry of the resistivity sample measured in Fig. 4 lead to a rough
agreement of both measurements to $\sim 0.3$ eV.

Above an excitation energy of 0.27 eV the optical conductivity gradually
increases and shows two clear peaks near 1.8 and 3.2 eV. These peaks arise
from higher energy interband transitions, i.e. optical excitation between the
valence and conduction bands that are further away from the Fermi energy. Note
that there may be additional interband transitions due to the observed kink-
like features in $\sigma_{1}$ at 1.3 eV, 4 eV and 5.5 eV.

Figure 6:  DFT-calculated electronic band structures of Ag3AuSe2 as obtained
from PBE and mBJ functionals. The general features of the band structures are
equivalent for both functionals, but the direct band gap is predicted to be 71
meV for PBE and 907 meV for the hybrid mBJ functional.

The DFT-calculated band structures of Ag3AuSe2 from PBE and mBJ functionals
are shown in Fig. 6. Although the gaps differ by about 1 eV, the shape and
features are consistent for both functionals. A direct gap is observed at
$\Gamma$ and bands are symmetric around $P_{1}$. Gaps of Ag3AuSe2 were
calculated using multiple approximations for both the experimentally-refined
and DFT-PBE relaxed structures. The obtained gaps are listed in Table 1. The
gaps span a wide range, from fully overlapped metallic bands to 931 meV. The
differences between experimentally-refined and DFT-relaxed structures for a
given functional are typically no larger than 30 meV, whereas the differences
between LDA/PBE and the hybrid functionals mBJ/HSE06 for a given structure is
much larger, over 800 meV. It is important to note that although the 30 meV
differences are comparatively small, subtle differences in the atomic
positions will have a large effect on transport if the band gap is of the same
order of magnitude (tens of meV). These small effects are also evident in
shifts of the bands due to increased overlap with compressive strain. The PBE
band gaps of the Ag3AuSe2 DFT-PBE-relaxed structures with varying lattice
strain are plotted in Fig. 7, showing an approach to metallicity as the
tensile strain reaches 3%, which is a large value of strain to be reached by
mechanical compression, but may be attainable by chemical substitution of the
Se2- anion in particular.

Figure 7:  The DFT-PBE band gap of the relaxed Ag3AuSe2 structure decreases with tensile strain, leading to band overlap around 3% strain.  | Experimental structure | DFT-PBE structure
---|---|---
LDA | No gap | No gap
PBE | 71 meV | 50 meV
mBJ | 931 meV | 928 meV
HSE06 | 907 meV | 872 meV
Table 1:  DFT-calculated band gaps of Ag3AuSe2 for experimentally-refined
structures and DFT-PBE-relaxed structures. Four pairs of ab-initio
calculations using the LDA, PBE, mBJ and HSE06 approximations for the exchange
correlation energy are used.

The noticeably larger band gap found by the hybrid mBJ and HSE functionals
points to additional disorder that may lead to defect states that dominate the
electronic properties of Ag3AuSe2. Percent-level antisite defects in the Ag/Au
ordering that occurs at the 270∘C order-disorder transition would be
challenging to quantify, and would likely remain coherent with the large-
grained microstructure. Additional single-crystal diffraction may shed light
on this possibility. Additionally, more detailed investigation of the optical
properties, including angle-resolved photoemission spectroscopy, should
correlate larger features of the conduction band with first-principles
results, rather than states exclusively near band edges. Understanding these
details for the entire class of fischesserite-type phases will reveal their
potential for nonlinear optics and topological responses.

##  IV Conclusions

We have established that large single crystals of Ag3AuSe2 can be grown by
slow-cooling from the melt. Ellipsometry and resistivity measurements point to
an activation energy and optical band gap around 0.3 eV, which is between the
PBE-estimated band gap of 0.07 eV (which should be an underestimation) and the
mBJ and HSE06 band gaps of about 0.9 eV, which we expect to be accurately
predicted. Given the sensitivity of the band gap to strain, and the necessary
order-disorder transition that Ag3AuSe2 experiences on cooling, it follows
that subtle tuning of the defect population, surface treatment, and
stoichiometry of Ag3AuSe2 and the related fischesserite-petzite family of
materials can lead to intentional tuning of the band gap. It remains to be
seen whether a band inversion can be obtained upon application of pressure of
chemical substitution with the telluride.

##  V Acknowledgments

This study was supported by the Center for Quantum Sensing and Quantum
Materials, an Energy Frontier Research Center funded by the U. S. Department
of Energy, Office of Science, Basic Energy Sciences under Award DE-SC0021238.
The authors acknowledge the use of facilities and instrumentation at the
Materials Research Laboratory Central Research Facilities, University of
Illinois, partially supported by NSF through the University of Illinois
Materials Research Science and Engineering Center DMR-1720633. MGV, IE and MGA
acknowledge the Spanish Ministerio de Ciencia e Innovación (grant number
PID2019109905GB-C21) and Programa Red Guipuzcoana de Ciencia, Tecnología e
Innovación 2021 No. 2021-CIEN-000070-01 Gipuzkoa Next. MGV thanks support from
the DeutscheForschungsgemeinschaft (DFG, German Research Foundation) GA
3314/1-1 – FOR5249 (QUAST). SB acknowledges support through the Early Postdoc
Mobility Fellowship from the Swiss National Science Foundation (Grant number
P2EZP2 191885). BL, JSS, and TWN were supported by the Institute for Basic
Science (IBS) in Korea (Grant No. IBS-R009-D1).

## References

  * Johan  _et al._ (1971) Z. Johan, P. Picot, R. Pierrot, and M. Kvacek, Bull. Soc. Fr. Mineral. Crist. 94, 381 (1971).
  * Bindi and Cipriani (2004) L. Bindi and C. Cipriani, Canad. Mineral. 42, 1733 (2004).
  * Tavernier  _et al._ (1967) B. H. Tavernier, J. Vervecken, P. Messien, and M. Baiwir, Z. Anorg. Allg. Chem. 356, 77 (1967).
  * Smit  _et al._ (1970) T. J. M. Smit, E. Venema, J. Wiersma, and G. A. Wiegers, J. Solid State Chem. 2, 309 (1970).
  * Barton  _et al._ (1978) M. D. Barton, C. Kieft, E. A. J. Burke, and I. S. Oen, Canad. Mineral. 16, 651 (1978).
  * Wiegers (1976) G. A. Wiegers, J. Less Common Met. 48, 269 (1976).
  * Wagner  _et al._ (1992) F. E. Wagner, J. A. Sawicki, J. Friedl, J. A. Mandarino, and D. C. Harris, Canad. Mineral. 30, 327 (1992).
  * (8) E. G. Osadchii and E. A. Echmaeva,  92, 640, publisher: De Gruyter.
  * Feng and Taskinen (2014) D. Feng and P. Taskinen, J. Alloys Compd. 583, 176 (2014).
  * Faizan  _et al._ (2016) M. Faizan, G. Murtaza, S. Azam, S. A. Khan, A. Mahmood, and A. Yar, Mater. Sci. Semicon. Proc. 52, 8 (2016).
  * Faizan  _et al._ (2017) M. Faizan, H. Ullah, S. H. Khan, S. M. Ramay, S. A. S. Qaid, A. Mahmood, and M. Ali, Int. J. Mod. Phys. B 31, 1750253 (2017).
  * Fang  _et al._ (2002) C. Fang, R. Groot, and G. Wiegers, J. Phys. Chem. Solids 63, 457 (2002).
  * Bradlyn  _et al._ (2016) B. Bradlyn, J. Cano, Z. Wang, M. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Science 353 (2016).
  * Chang  _et al._ (2017) G. Chang, S.-Y. Xu, B. J. Wieder, D. S. Sanchez, S.-M. Huang, I. Belopolski, T.-R. Chang, S. Zhang, A. Bansil, H. Lin, and M. Z. Hasan, Phys. Rev. Lett. 119, 206401 (2017).
  * Schröter  _et al._ (2019) N. B. Schröter, D. Pei, M. G. Vergniory, Y. Sun, K. Manna, F. De Juan, J. A. Krieger, V. Süss, M. Schmidt, P. Dudin, _et al._ , Nature Physics 15, 759 (2019).
  * Schröter  _et al._ (2020) N. B. Schröter, S. Stolz, K. Manna, F. De Juan, M. G. Vergniory, J. A. Krieger, D. Pei, T. Schmitt, P. Dudin, T. K. Kim, _et al._ , Science 369, 179 (2020).
  * Sanchez  _et al._ (2019) D. S. Sanchez, I. Belopolski, T. A. Cochran, X. Xu, J.-X. Yin, G. Chang, W. Xie, K. Manna, V. Süß, C.-Y. Huang, N. Alidoust, D. Multer, S. S. Zhang, N. Shumiya, X. Wang, G.-Q. Wang, T.-R. Chang, C. Felser, S.-Y. Xu, S. Jia, H. Lin, and M. Z. Hasan, Nature 567, 500 (2019).
  * Flicker  _et al._ (2018) F. Flicker, F. De Juan, B. Bradlyn, T. Morimoto, M. G. Vergniory, and A. G. Grushin, Physical Review B 98, 155145 (2018).
  * Sánchez-Martínez  _et al._ (2019) I. Sánchez-Martínez, M.-Áand Robredo, A. Bidaurrazaga, A. Bergara, F. d. Juan, A. G. Grushin, and M. G. Vergniory, J. Phys.: Mater. 3, 014001 (2019).
  * Toby and Von Dreele (2013) B. H. Toby and R. B. Von Dreele, J. Appl. Cryst. 46, 544 (2013).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, J. Phys. Cond. Mat. 29, 465901 (2017).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Tran and Blaha (2009) F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
  * Heyd  _et al._ (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003).
  * Young  _et al._ (2000) D. P. Young, C. L. Brown, P. Khalifah, R. J. Cava, and A. P. Ramirez, J. Applied Phys. 88, 5221 (2000).
  * Fox and Fox (2001) A. Fox and D. Fox,  _Optical Properties of Solids_ (Oxford University Press, 2001).