# Piezo-optic effect of high-harmonic generation in semiconductors

Tomohiro Tamaya<EMAIL_ADDRESS>Takeo Kato  Institute for Solid
State Physics, University of Tokyo, Kashiwa, 277-8581, Japan

###### Abstract

We theoretically investigate the piezo-optic effect of high-harmonic
generation (HHG) in shear-strained semiconductors. By focusing on a typical
semiconductor, GaAs, we show that there is optical activity, meaning different
responses to right-handed and left-handed elliptically polarized electric
fields. We also show that this optical activity is more pronounced for higher
harmonics whose perturbative order exceeds the band-gap energy. These findings
point to a useful pathway for strain engineering of nonlinear optics to
control the reciprocity of HHG.

##  I Introduction

High-harmonic generation (HHG) is one of the most fundamental topics in
nonlinear opticsBoyd (2008); Yariv and Yeh (1984); Shen (1984). Experimental
progress from the perturbative to the nonperturbative regime in gaseous media
has paved the way for developing novel optical devices for, e.g., generating
short-wavelength attosecond pulses Corkum (1993); Protopapas  _et al._ (1997);
Brabec and Krausz (2000); Agostini and DiMauro (2004); Corkum and Krausz
(2007); Krausz and Ivanov (2009). Moreover, HHG in the nonperturbative regime
has been experimentally observed in solidsGhimire  _et al._ (2011); Schubert
_et al._ (2014); Luu  _et al._ (2015); Hohenleutner  _et al._ (2015); Vampa
_et al._ (2015); Ghimire and Reis (2019), and ensuing studies have opened up a
new field in condensed-matter science Ndabashimiye  _et al._ (2016); Lanin
_et al._ (2017); Liu  _et al._ (2017a); You  _et al._ (2017); Yoshikawa  _et
al._ (2017); Lanin  _et al._ (2017); Kim  _et al._ (2017); Jiang  _et al._
(2018); Langer  _et al._ (2018); G. Vampa and Y. S. You and H. Liu and S.
Ghimire and D. A. Reis (2018); Silva  _et al._ (2018); Hirori  _et al._
(2019); Cheng  _et al._ (2020). In contrast to gaseous media, HHG in solids
has various inherent properties that are rooted in the crystallinity of the
medium and may provide means of developing new optical devices that use HHG.
Thus, it is important to study the characteristics of HHG in various materials
and to devise a control method that can provide a possible route to novel
optical technology.

The most important aspects that determine the properties of HHG are the band
structures of the materials and the corresponding Bloch wavefunctions. The
Hamiltonian of the light-matter interaction is principally made up of these
elements and HHG is expected to yield unusual new features by appropriately
choosing those materials. A recent experimental study reported that HHG in
monolayer $\rm{MoS}_{2}$ was polarized perpendicular to the linearly polarized
pump fieldYou  _et al._ (2017), an effect that was mainly explained in terms
of the anomalous transverse intraband current arising from the material’s
Berry curvature. Thus, the properties of the Bloch wavefunctions in materials
directly affect the characteristics of HHG, and thus, exploring methods of
controlling these wavefunctions are crucial for applications of HHG.

One possible way to control the Bloch wavefunctions in materials is strain
engineeringSmith (1954); Llordes  _et al._ (2012); Lü  _et al._ (2016); Liu
_et al._ (2017b); Chen  _et al._ (2018); Zhao  _et al._ (2017); Zhu  _et al._
(2019); Steele  _et al._ (2019); Peng  _et al._ (2020); Conley  _et al._
(2013); He  _et al._ (2013); Frank  _et al._ (2011); Shi  _et al._ (2013);
Castellanos-Gomez  _et al._ (2015); Bertolazzi  _et al._ (2011); Castellanos-
Gomez  _et al._ (2012); Bao  _et al._ (2017); Chen  _et al._ (2020); Xiong
_et al._ (2020); Shi and Chen (2019); Wang  _et al._ (2015). Mechanical
deformation of a material modifies the Bloch wavefunctions by distorting the
crystal structure, and it can be used to control various physical properties
such as transport and optical response. We expect, for example, that shear
strain will rotate the direction of the generated current (see Fig. 1). This
rotation of the current direction indicates left-right symmetry breaking in
materials leading to emergent optical activity of HHG, i.e., different
responses to right-handed and left-handed elliptically polarized electric
fieldsBarron (2009); Mason (1982). In particular, the piezo-optic effect of
HHG referred to here could be used for applications such as spatially resolved
distortion measurement and mechanical control of HHG, which are considered
impossible for gaseous mediaWang  _et al._ (2020); Shao  _et al._ (2019).

Figure 1: (Color online) Schematic diagram of the current $\bm{J}(t)$
generated by a linearly polarized electric field. The left figure indicates
the generated current in the absence of strain, while the right figure
indicates the generated current when an external shear strain is applied. In
the absence of strain, current is generated only in the $z$-direction parallel
to the electric field of the incident light. When shear strain is applied to
the system, the direction of the current is distorted and $\bm{J}(t)$ has both
$z$\- and $x$-components.

In this paper, we theoretically investigate the characteristics of HHG in
shear-strained semiconductors. We construct a theoretical framework based on
the Luttinger-Kohn-Bir-Pikus HamiltonianLuttinger and Kohn (1955); Ahn  _et
al._ (1995); Dargys (2002); Pryor (1998); Tomić  _et al._ (2006); Chuang
(1991); Pfeffer and Zawadzki (1996); Sytnyk and Melnik (2018); Bir (1974);
Luque  _et al._ (2015), which provides us with a general platform to treat
various semiconductors under external shear strain. By performing a dynamical
simulation on GaAs, a typical III-V semiconductor, we find that external shear
strain indeed causes different responses to right-handed and left-handed
elliptically polarized electric fields. We also find that this non-reciprocity
is more pronounced for higher harmonics whose perturbative order exceeds the
band-gap energy.

The organization of this paper is as follows. In Sec. II, we introduce the
theoretical framework for HHG using the Luttinger-Kohn-Bir-Pikus Hamiltonian.
In Sec. III, we show numerical results for HHG emitted from shear-strained
GaAs. We also discuss the optical activity of HHG in detail by examining the
different optical responses to right- and left-handed elliptically polarized
light. Sec. IV summarizes the conclusions of this study.

##  II Formulation

The general formulation introduced here for HHG is applicable to various
semiconductors with a direct band gap at the $\Gamma$ point. First, we
describe the eight-band Luttinger-Kohn model in Sec. II.1; then, we extend it
to the strained case, called the Pikus-Bir Hamiltonian model, in Sec. II.2.
Next, we derive the light-matter interaction in terms of the Luttinger-Kohn-
Pikus-Bir model in Sec. II.3. Finally, we describe the time-dependent
Schrödinger equation for electrons and define the polarization currents in
Sec. II.4.

###  II.1 Luttinger-Kohn model

Let us consider a general microscopic Hamiltonian,

$\displaystyle
H=\frac{\bm{p}^{2}}{2m_{0}}+\sum_{i}\left[V_{i}(\bm{x})+\frac{\hbar}{4m_{0}^{2}c^{2}}(\nabla
V_{i}\times{\bm{p}})\cdot{\bm{\sigma}}\right],$ (1)

where $m_{0}$ is the electron mass, $\bm{p}$ is the momentum of the electron,
and $V_{i}({\bm{x}})=V(\bm{x}-\bm{R}_{i})$ is the periodic core potential of
atoms located at $\bm{R}_{i}$. The second term in brackets expresses the spin-
orbit coupling, where ${\bm{\sigma}}$ is the spin angular momentum. By
performing a band calculation, the Hamiltonian can be diagonalized as

$\displaystyle
H|\Psi_{n{\bm{k}}}\rangle=E_{n{\bm{k}}}|\Psi_{n{\bm{k}}}\rangle,$ (2)

where $n$ is the band index, ${\bm{k}}$ is the Bloch wavenumber,
$E_{n{\bm{k}}}$ is the energy dispersion of the $n$th band, and
$|\Psi_{n{\bm{k}}}\rangle$ is the Bloch wavefunction. Here, we will focus on
the bands near the band edge at the $\Gamma$ point (${\bm{k}}=0$) and restrict
them to eight bands composed of one conduction band ($n=1$) and three valence
bands, i.e., a heavy-hole band ($n=2$), a light-hole band ($n=3$), and a
split-off band ($n=4$), and their time-reversal counterparts ($n=5,6,7,8$).

We apply conventional ${\bm{k}}\cdot{\bm{p}}$ perturbation theoryHarrison
(2012); Kittel and Fong (1963); Peter and Cardona (2010); Voon and Willatzen
(2009) around the $\Gamma$ point using these eight bands. We rewrite the Bloch
wavefunction as
$\ket{\Psi_{n\bm{k}}}=e^{i{\bm{k}}\cdot{\bm{x}}}\ket{u_{n{\bm{k}}}}$. The
eigenvalue equation is rewritten as
$\tilde{H}\ket{u_{n{\bm{k}}}}=E_{n{\bm{k}}}\ket{u_{n{\bm{k}}}}$, where
$\tilde{H}$ is an effective Hamiltonian defined as

$\displaystyle\tilde{H}$ $\displaystyle\equiv
e^{-i{\bm{k}}\cdot{\bm{x}}}He^{i{\bm{k}}\cdot{\bm{x}}}\equiv\tilde{H}_{0}+\tilde{V},$
(3) $\displaystyle\tilde{H}_{0}$
$\displaystyle=\frac{{\bm{p}}^{2}}{2m_{0}}+\sum_{i}\left[V_{i}(\bm{x})+\frac{\hbar}{4m_{0}^{2}c^{2}}(\nabla
V_{i}\times{\bm{p}})\cdot{\bm{\sigma}}\right],$ (4) $\displaystyle\tilde{V}$
$\displaystyle=\frac{\hbar}{m_{0}}{\bm{k}}\cdot{\bm{p}}+\frac{\hbar^{2}\bm{k}^{2}}{2m_{0}}.$
(5)

The unperturbed Haimiltonian $\tilde{H}_{0}$ is diagonalized by the
wavefunction at the $\Gamma$ point, $\ket{u_{n}}\equiv\ket{u_{n{\bm{k}}=0}}$.
Following conventional $\bm{k}\cdot\bm{p}$ perturbation theory, we incorporate
the ${\bm{k}}$-dependence of the eigen wavefunctions by second-order
perturbation with respect to $\tilde{V}$, taking the effect of outside bands
other than the target bands into account. The resultant effective Hamiltonian
is Tomić  _et al._ (2006)

$\displaystyle\braket{u_{n}}{\tilde{H}^{\rm{eff}}_{0}}{u_{n^{\prime}}}=\left(\begin{array}[]{cc}H^{\bm{k}}_{uu}&H^{\bm{k}}_{ul}\\\
H^{\bm{k}}_{lu}&H^{\bm{k}}_{ll}\end{array}\right),$ (8)

where $H^{\bm{k}}_{uu}$, $H^{\bm{k}}_{ul}$, $H^{\bm{k}}_{lu}$, and
$H^{\bm{k}}_{ll}$ are 4 $\times$ 4 submatrices. The submatrix
$H^{\bm{k}}_{uu}$ has the form,

$\displaystyle
H^{\bm{k}}_{uu}=\left(\begin{array}[]{cccc}E_{CB}&-\sqrt{3}T&\sqrt{2}U&-U\\\
-\sqrt{3}T^{*}&E_{HH}&\sqrt{2}S&-S\\\
\sqrt{2}U&\sqrt{2}S^{*}&E_{LH}&-\sqrt{2}Q\\\ -U&-S^{*}&-\sqrt{2}Q&E_{SO}\\\
\end{array}\right),$ (13)

while the submatrix $H_{ll}$ is defined as $H_{ll}=H_{uu}^{*}$. The
submatrices, $H_{ul}$ and $H_{lu}$, are expressed as

$\displaystyle
H^{\bm{k}}_{ul}=\left(\begin{array}[]{cccc}0&0&-T^{*}&-\sqrt{2}T^{*}\\\
0&0&-R&-\sqrt{2}R\\\ T^{*}&R&0&\sqrt{3}S\\\
\sqrt{2}T^{*}&\sqrt{2}R&-\sqrt{3}S&0\\\ \end{array}\right),$ (18)

and $H_{lu}=H_{ul}^{\dagger}$. The diagonal elements of $H_{uu}$ and $H_{ll}$
are defined as

$\displaystyle E_{CB}$ $\displaystyle=E_{g}+O,$ (19) $\displaystyle E_{HH}$
$\displaystyle=-(P+Q),$ (20) $\displaystyle E_{LH}$ $\displaystyle=-(P-Q),$
(21) $\displaystyle E_{SO}$ $\displaystyle=-(P+\Delta_{SO}).$ (22)

The subscripts CB, HH, LH, and SO stand for conduction, heavy-hole, light-hole
and split-off bands, respectively, and $E_{g}$ and $\Delta_{SO}$ are the band-
gap energy and the split-off energy due to the spin-orbit interaction, where
$E_{g}=1.42$ eV and $\Delta_{SO}=0.34$ eV in GaAs.

In the absence of external strain, the matrix elements are

$\displaystyle O$
$\displaystyle=\frac{\hbar^{2}}{2m_{0}}\gamma_{C}\left(k_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}\right),$
(23) $\displaystyle P$
$\displaystyle=\frac{\hbar^{2}}{2m_{0}}\gamma_{1}\left(k_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}\right),$
(24) $\displaystyle Q$
$\displaystyle=\frac{\hbar^{2}}{2m_{0}}\gamma_{2}\left(k_{x}^{2}+k_{y}^{2}-2k_{z}^{2}\right),$
(25) $\displaystyle R$
$\displaystyle=\frac{\hbar^{2}}{2m_{0}}\sqrt{3}\left[\gamma_{2}(k_{x}^{2}-k_{y}^{2})-2i\gamma_{3}k_{x}k_{y}\right],$
(26) $\displaystyle S$
$\displaystyle=\frac{\hbar^{2}}{2m_{0}}\sqrt{6}\gamma_{3}\left(k_{x}-ik_{y}\right)k_{z},$
(27) $\displaystyle T$
$\displaystyle=\frac{1}{\sqrt{6}}P_{0}\left(k_{x}+ik_{y}\right),$ (28)
$\displaystyle U$ $\displaystyle=\frac{1}{\sqrt{3}}P_{0}k_{z}.$ (29)

Here, $k_{x}$, $k_{y}$, and $k_{z}$ denote components of the Bloch wavevector
along the [100], [010], and [001] crystallographic directions, respectively,
and $\gamma_{0}$, $\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$ are the
Luttinger parameters. We set the Luttinger parameters of GaAs to be
$\gamma_{C}=0.5$, $\gamma_{1}=2.7$, $\gamma_{2}=-0.1$, and $\gamma_{3}=0.7$
following Ref. Sytnyk and Melnik, 2018. The dipole matrix element (the Kane
matrix element) is defined as

$P_{0}=-i\left(\frac{\hbar}{m_{0}}\right)\braket{s;\sigma}{p_{\lambda}}{\lambda;\sigma},$
(30)

where $\lambda=x,y,z$. The value of $P_{0}$ can be absorbed into the
definition of the Rabi frequency introduced later.

###  II.2 Strain-induced effect in semiconductors

Strain in the crystal is expressed by displacement of the lattice vectors from
those of the unstrained crystal, ${\bm{x}}_{i}$ ($i=x,y,z$):

$\displaystyle\delta x^{j}_{i}=\sum_{i}\delta_{ij}x^{j}_{i},$ (31)

where $x^{j}_{i}$ is the $j$-th component of ${\bm{x}}_{i}$ ($j=x,y,z$), and
$\delta_{ij}$ ($i,j=x,y,z$) denote components of the strain tensor. The effect
of the strain can be incorporated into the ${\bm{k}}\cdot{\bm{p}}$ band
structure calculations by adding an extra perturbation term to the unstrained
potentialBir (1974). Thus, the change in the matrix elements in the presence
of the strain is obtained as $O\rightarrow O+\delta O$, $P\rightarrow P+\delta
P$, and so on, whereLuttinger and Kohn (1955); Ahn  _et al._ (1995); Dargys
(2002); Pryor (1998); Tomić  _et al._ (2006); Chuang (1991); Pfeffer and
Zawadzki (1996); Sytnyk and Melnik (2018); Bir (1974); Luque  _et al._ (2015)

$\displaystyle\delta O$
$\displaystyle=+a_{c}\left(\delta_{xx}+\delta_{yy}+\delta_{zz}\right),$ (32)
$\displaystyle\delta P$
$\displaystyle=-a_{v}(\delta_{xx}+\delta_{yy}+\delta_{zz}),$ (33)
$\displaystyle\delta Q$
$\displaystyle=-\frac{b_{v}}{2}\left(\delta_{xx}+\delta_{yy}+\delta_{zz}\right),$
(34) $\displaystyle\delta R$
$\displaystyle=-\frac{\sqrt{3}}{2}b_{v}\left(\delta_{xx}-\delta_{yy}\right)+id_{v}\epsilon_{xy},$
(35) $\displaystyle\delta S$
$\displaystyle=-\frac{d_{v}}{\sqrt{2}}\left(\delta_{zx}-i\delta_{yz}\right),$
(36) $\displaystyle\delta T$
$\displaystyle=-\frac{1}{\sqrt{6}}P_{0}\sum_{j}\left(\delta_{xj}+i\delta_{yj}\right)k_{j},$
(37) $\displaystyle\delta U$
$\displaystyle=-\frac{1}{\sqrt{3}}P_{0}\sum_{j}\delta_{zj}k_{j}{.}$ (38)

Here, $a_{c}$ and $a_{v}$ are the conduction- and valence-band hydrostatic
deformation potentials of the host material, and $b_{v}$ and $d_{v}$ are the
shear deformation potentials along the [001] and [111] directions of the host
material, respectively. Here, we will set the deformation potentials of GaAs
to be $a_{c}=-9.3$ eV, $a_{v}=-0.7$ eV, $b_{v}=2.0$ eV, and $d_{v}=5.4$ eV,
following Ref. Pryor, 1998.

###  II.3 Light-matter interaction

Next, let us consider a bulk crystal of GaAs that is driven by elliptically
polarized electric fields. Here, we take the $z$-axis ([001] direction) to be
the major axis and the $x$-axis ([100] directions) to be the minor axis. Then,
the vector potential of the elliptically polarized electric field
${\bm{A}}(t)$ can be defined as

$\displaystyle{\bm{A}}(t)$ $\displaystyle=(A_{x}(t),0,A_{z}(t))$
$\displaystyle=A_{0}f(t)(\eta\sin{\omega t},0,\cos{\omega t}),$ (39)

where $\eta$ and $A_{0}$ are the ellipticity and the amplitude of the electric
field, respectively, and $f(t)$ is the envelop function defined as

$\displaystyle{f(t)=\exp\left(-\frac{(t-t_{0})^{2}}{\tau^{2}}\right).}$ (40)

Here, we set $t_{0}=24\pi/\omega$ and $\tau=4\pi/\omega$.

We introduce the light-matter interaction through the vector potential:

$\displaystyle H_{0}+H_{\rm ex}$
$\displaystyle=\frac{1}{2m_{0}}\left(\bm{p}-\frac{e}{c}\bm{A}(t)\right)^{2}$
$\displaystyle+\sum_{i}\left[V_{i}(\bm{x})+\frac{\hbar}{4m_{0}^{2}c^{2}}(\nabla
V_{i}\times{\bm{p}})\cdot{\bm{\sigma}}\right].$ (41)

Here, we have assumed that the term caused by the replacement
${\bm{p}}\rightarrow{\bm{p}}-e/c{\bm{A}}(t)$ in the spin-orbit interaction is
small enough to be neglected. Thus, the Hamiltonian for the light-matter
interaction is

$\displaystyle H_{\rm
ex}=-\frac{e}{m_{0}c}{\bm{A}}(t)\cdot{\bm{p}}+\frac{e^{2}}{2m_{0}c^{2}}{\bm{A}}(t)^{2}.$
(42)

The second term in $H_{\rm ex}$ can be eliminated by performing a unitary
transformation $H_{\rm ex}\rightarrow U^{-1}_{1}H_{\rm ex}U_{1}$, where

$\displaystyle
U_{1}=\exp\left[\frac{ie^{2}}{2m_{0}c^{2}\hbar^{2}}\int^{t}_{0}dt^{\prime}\,{\bm{A}}(t^{\prime})^{2}\right].$
(43)

$H_{\rm ex}$ is then rewritten as $\tilde{H}_{\rm ex}$, which operates on the
eigenstate $\ket{u_{n{\bm{k}}}}$:

$\displaystyle\tilde{H}_{\rm ex}$ $\displaystyle\equiv
e^{-i{\bm{k}}\cdot{\bm{x}}}H_{\rm ex}e^{i{\bm{k}}\cdot{\bm{x}}}$
$\displaystyle=-\frac{e}{m_{0}c}\left[{\bm{A}}(t)\cdot\hbar\bm{k}+\bm{A}(t)\cdot{\bm{p}}\right].$
(44)

Here, we can also eliminate the first term in brackets through the unitary
transformation $\tilde{H}_{\rm ex}\rightarrow U_{2}^{-1}\tilde{H}_{\rm
ex}U_{2}$, whose matrix elements are expressed as

$\displaystyle\bra{u_{n{\bm{k}}}}U_{2}\ket{u_{n^{\prime}{\bm{k}}^{\prime}}}$
$\displaystyle\hskip
14.22636pt=\exp\left(-i\frac{e}{m_{0}c}\int_{0}^{t}dt^{\prime}\,{\bm{k}}\cdot{\bm{A}}(t^{\prime})\right)\delta_{{\bm{k}},{\bm{k}}^{\prime}}\delta_{n,n^{\prime}}.$
(45)

Thus, the Hamiltonian of the light-matter interaction becomes

$\displaystyle\tilde{H}_{\rm ex}$
$\displaystyle=-\frac{e}{m_{0}c}\bm{A}(t)\cdot{\bm{p}}.$ (46)

Figure 2: (Color online) (a) Schematic diagram of HHG in shear-strained GaAs
caused by a linearly polarized electric field ($\eta=0$). The polarization
axis of the light rotates as it propagates through the shear-strained
material. As a result, the electric field of the emitted light includes
components along the major axis (the $z$-direction) and the minor axis (the
$x$-direction). (b) Numerical result of HHG spectra in shear-strained GaAs
caused by linearly polarized electric field ($\eta=0$) in the case of
$\delta_{zx}=0.01$. The blue and red lines show the HHG spectra emitted along
the major axis ($z$-axis) and minor axis ($x$-axis), respectively.

To express the light-matter interaction in a simple matrix form, we will
return to the elemental basis set of $\ket{s\\!\uparrow}$,
$\ket{p_{x}\\!\uparrow}$, $\ket{p_{y}\\!\uparrow}$, $\ket{p_{z}\\!\uparrow}$,
$\ket{s\\!\downarrow}$, $\ket{p_{x}\\!\downarrow}$,
$\ket{p_{y}\\!\downarrow}$, and $\ket{p_{z}\\!\downarrow}$, and we will
redefine them as $\ket{v_{m}}$ $(m=1,2,\cdots,8)$ for simplicity of notation.
Note that the eigen wavefunctions at the $\Gamma$ point, $\ket{u_{n}}$’s
($n=1,2,\cdots,8)$, are expressed by a linear combination of $\ket{v_{m}}$
(see Appendix A for the explicit forms). In this basis set, almost all of the
matrix elements in $\tilde{H}_{\rm{ex}}$ are zero because of parity symmetry.
The nonzero matrix elements are given as

$\displaystyle\braket{s\sigma}{\tilde{H}_{\rm{ex}}}{p_{x}\sigma}$
$\displaystyle=(\braket{p_{x}\sigma}{\tilde{H}_{\rm{ex}}}{s\sigma})^{*}$
$\displaystyle=-i\hbar\Omega_{R0}f(t)\cos\omega t,$ (47)
$\displaystyle\braket{s\sigma}{\tilde{H}_{\rm{ex}}}{p_{z}\sigma}$
$\displaystyle=(\braket{p_{z}\sigma}{\tilde{H}_{\rm{ex}}}{s\sigma})^{*}$
$\displaystyle=-i\hbar\Omega_{R0}\eta f(t)\sin\omega t,$ (48)

for $\sigma=\uparrow,\downarrow$, where we have defined the Rabi frequency as
$\Omega_{R0}=(e/c\hbar^{2})P_{0}A_{0}$. Thus, we have derived the light-matter
interaction Hamiltonian in a simple matrix form using the elemental basis set
$\ket{v_{m}}$. Here, we set $E_{g}=7\hbar\omega$ and $\Omega_{R0}=4\omega$,
respectively.

###  II.4 Dynamical simulation

We solve the time-dependent Schrödinger equation,

$\displaystyle i\hbar\frac{\partial}{\partial t}\ket{u_{{\bm{k}}}{(t)}}$
$\displaystyle=\tilde{H}_{{\rm tot}}^{{\bm{k}}}\ket{u_{\bm{k}}{(t)}},$ (49)

where $\tilde{H}_{{\rm tot}}^{{\bm{k}}}=\tilde{H}_{{\rm
eff}}^{{\bm{k}}}+\tilde{H}_{{\rm ex}}$ is the total Hamiltonian. In the
simulation, we employed the atomic basis $\ket{v_{m}}$ and expanded the
wavefunction $\ket{u_{\bm{k}}{(t)}}$ as

$\displaystyle\ket{u_{\bm{k}}{(t)}}=\sum_{m=1}^{8}a_{m{\bm{k}}}(t)\ket{v_{m}}.$
(50)

Using this basis set, the matrix elements of the light-matter interaction
$\tilde{H}_{{\rm ex}}$ are given by Eqs. (47)-(48), while those of the system
Hamiltonian are given as

$\displaystyle(\tilde{H}_{\rm
eff}^{\bm{k}})_{mm^{\prime}}\equiv\bra{v_{m}}\tilde{H}_{\rm
eff}^{\bm{k}}\ket{v_{m^{\prime}}}$ $\displaystyle\hskip
14.22636pt=\sum_{n,n^{\prime}=1}^{8}\braket{v_{m}}{u_{n}}\braket{u_{n}}{\tilde{H}_{\rm
eff}^{\bm{k}}}{u_{n^{\prime}}}\braket{u_{n^{\prime}}}{v_{m^{\prime}}}.$ (51)

Here, $\braket{u_{n}}{\tilde{H}_{\rm eff}^{\bm{k}}}{u_{n^{\prime}}}$ is as in
Eq. (8), and $(U)_{nm}=\braket{u_{n}}{v_{m}}$ is the unitary matrix for the
basis transformation, whose explicit forms are in Appendix A. Thus, the time-
dependent Schrödinger equation finally becomes

$\displaystyle
i\hbar\frac{da_{m{\bm{k}}}}{dt}=\sum_{m^{\prime}=1}^{8}(\tilde{H}_{\rm
eff}^{\bm{k}}+\tilde{H}_{\rm ex})_{mm^{\prime}}a_{m^{\prime}{\bm{k}}}(t).$
(52)

The generated currents along the [001] and [100] directions are calculated as

$\displaystyle J_{x}(t)$ $\displaystyle=-\frac{1}{c}\Braket{\frac{\partial
H_{I}}{\partial A_{x}}}$
$\displaystyle\propto-i\sum_{\bm{k}\sigma}\left[a_{s\sigma{\bm{k}}}(t)^{*}a_{p_{x}\sigma{\bm{k}}}(t)-{\rm
c.c.}\right]$ (53) $\displaystyle J_{z}(t)$
$\displaystyle=-\frac{1}{c}\Braket{\frac{\partial H_{I}}{\partial A_{z}}}$
$\displaystyle\propto-i\sum_{{\bm{k}}\sigma}\left[a_{s\sigma{\bm{k}}}(t)^{*}a_{p_{z}\sigma{\bm{k}}}(t)-{\rm
c.c.}\right]$ (54)

The HHG spectra in GaAs along the [001] and [100] directions are calculated as
$I_{z}=\left|\omega{\cal J}_{z}(\omega)\right|^{2}$ and
$I_{x}=\left|\omega{\cal J}_{x}(\omega)\right|^{2}$, where ${\cal
J}_{z}(\omega)$ and ${\cal J}_{x}(\omega)$ are the Fourier transforms of the
generated currents $J_{z}(t)$ and $J_{x}(t)$. We numerically solve the time-
dependent differential equation (52) under the initial conditions where
$\ket{u_{\bm{k}}(t=0)}=\ket{u_{n}}$ for occupied valence bands
($n=2,3,4,6,7,8$) and sum up the currents with respect to these six initial
conditions and the Bloch wavenumber ${\bm{k}}$.

##  III RESULTS AND DISCUSSION

In this section, we discuss the characteristics of HHG originating from shear-
strained GaAs. First, let us examine the numerical results of HHG caused by a
linearly polarized electric field in Sec. III.1. Here, we identify a rotation
of the polarization axis of the emitted light. Next, let us examine the
numerical results of HHG caused by elliptically polarized light in Sec. III.2.
These results indicate a breakdown in reciprocity of HHG in the shear-strained
material. Sec. III.3 discusses the physical interpretation of the numerical
results.

###  III.1 Linearly polarized electric fields

Now, let us consider the case of linearly polarized electric fields ($\eta=0$)
in shear-strained GaAs:

$\displaystyle\delta_{ij}=\left\\{\begin{array}[]{ll}\delta_{zx}\neq
0,&(i,j)=(z,x),\\\ 0,&({\rm otherwise}).\end{array}\right.$ (57)

Fig. 2 (a) shows a schematic diagram of the effect of strain on HHG
($\delta_{zx}\neq 0$). When the electric field of the incident light is
polarized in the $z$-direction, current is generated along the major ($z$) and
minor ($x$) axes (see also Fig. 1). As a result, the electric fields of the
emitted light also include an $x$-component, resulting in a rotation of the
polarized light.

Figure 3: (Color online) (a) Schematic diagram of HHG for shear-strained GaAs
caused by elliptically polarized electric fields ($\eta\neq 0$). (b) Numerical
result for the calculated ellipticity dependence of HHG in the absence of
strain. The blue and red lines show the ellipticity dependences of HHG emitted
along to the major axis ($z$-axis) and minor axis ($x$-axis), respectively.
The finite shear strain changes these dependences and is expected to cause
optical activity, i.e., different responses to the right-handed ($\eta>0$) and
left-handed ($\eta<0$) elliptically polarized electric fields.  Figure 4:
(Color online) Calculated ellipticity dependences of high-order harmonic
intensities focusing on the third ((a1) and (b1)), fifth ((a2) and (b2)),
seventh ((a3) and (b3)), and ninth harmonics ((a4) and (b4)). Figures
(a1)-(a4) show the harmonic intensities emitted along the major axis
($z$-axis), while figures (b1)-(b4) show those along the minor axis
($x$-axis). The red, green, and blue lines indicate the ellipticity
dependences for $\delta_{zx}{=0}$, 0.005, and 0.01, respectively. These
figures indicate that the right-handed ($\eta>0$) and left-handed ($\eta<0$)
elliptically polarized electric fields yield the different HHG spectra. This
breakdown of reciprocity becomes significant for the high-order harmonics
above the band-gap energy (n $\geq$ $E_{g}/\hbar\omega_{0}=7$).

Fig. 2 (b) shows the calculated HHG spectra for shear-strained GaAs
($\delta_{zx}=0.01$). The blue and red curves denote the intensity of the
emitted light parallel to the major ($z$) and minor ($x$) axes, respectively.
We note that a few megapascal pressure yields one percent displacement of the
lattice vectors in GaAs ($\delta_{zx}=0.01$), that can be realized in an
experiment by using the boat technique or the liquid encapsulated Czochralski
techniqueYonenaga and Sumino (1992). The spectra for the $z$-direction (the
major axis) exhibit several features characteristics of HHG; the peaks
correspond to $n\omega$ for odd $n$, and their heights first decay
exponentially with respect to $n$ in the perturbative regime ($n<7$),
eventually reaching a plateau in the intermediate regime ($7\leq n\leq 13$),
and finally collapsing exponentially again for large $n$ ($>13$). The spectra
for the $x$-direction (the minor axis) are similar to those in the
$z$-direction except that their intensity is much smaller. Here, it is
remarkable that the ratio of the intensity in the $x$-direction relative to
that of the $z$-direction is enhanced around the band-gap energy, that is,
$n\omega\approx E_{\rm g}$, where $E_{\rm g}$ is the band gap energy. To
discuss this strain-induced effect in detail, the next subsection examines the
case of elliptically polarized electric fields.

###  III.2 Elliptically polarized electric fields

Next, let us discuss the case of elliptically polarized electric fields
($\eta\neq 0$). Figure 3 (a) shows a schematic diagram of the strain-induced
effect on GaAs driven by elliptically polarized electric fields
($\delta_{zx}\neq 0$). The $z$\- and $x$-components of the electric fields of
the emitted light are modified from those of the incident light. The
intensities of the emitted light in the two directions are plotted as a
function of ellipticity $\eta$. As a reference, Fig. 3 (b) plots the
ellipticity dependences of the seventh harmonics emitted along to the $z$-axis
(blue line) and $x$-axis (red line) for unstrained GaAs ($\delta_{zx}=0$). In
this figure, we can identify that the HHG intensity in the $z$-direction (the
major axis) has a single Gaussian-like peak at $\eta=0$, while the intensity
in the $x$-direction (the minor axis) has two peaks at finite values of
$\eta$. These features have already been identified in the previous work
Tamaya  _et al._ (2016a); the double peaks in the $x$-direction grow with
increasing field intensity and become especially pronounced in the semimetal
regimeTamaya  _et al._ (2016b). This phenomenon was observed in an experiment
on HHG using graphene and $\rm{MoS_{2}}$Yoshikawa  _et al._ (2017).

Fig. 4 shows high-order harmonic intensities of the third ((a1) and (b1)),
fifth ((a2) and (b2)), seventh ((a3) and (b3)), and ninth harmonics ((a4) and
(b4)) as a function of ellipticity $\eta$ for three values of $\delta_{zx}$.
The red, green, and blue lines indicate the HHG spectra for
$\delta_{zx}=0.01$, $0.005$, and $0$, respectively. For unstrained GaAs
($\delta_{zx}=0$), the HHG intensity is symmetric with respect to an inversion
of ellipticity ($\eta\rightarrow-\eta$). This inversion symmetry, however, is
broken for strained GaAs ($\delta_{zx}=0.005,0.01$). The peak position in the
$z$-direction (the major axis) shifts toward positive $\eta$ (see figures
(a1)-(a4)). As the order of HHG increases, the peak shift becomes more
significant, and its height gradually decreases in comparison with the
unstrained case. For the $x$-direction (the minor axis), the heights of the
two peaks become different in the strained case (see figures (b1)-(b4)). From
the numerical results of the seventh and ninth HHG (see figures (b3) and
(b4)), we find that one of the two peaks disappears for sufficiently large
strain.

###  III.3 Breakdown of reciprocity

The breakdown of the reciprocity relation between the right-handed ($\eta>0$)
and left-handed ($\eta<0$) elliptically polarized electric fields obtained in
this work can be understood in terms of symmetry reduction of the crystal
structure. The external shear strain changes the crystal structure of GaAs
from cubic (zinc blende structure) into monoclinic, inducing off-diagonal
elements in the dielectric tensor $\epsilon_{zx}$Boyd (2008); Powell (2010).
The emergence of $\epsilon_{zx}$ directly means optical activity, i.e.,
different responses to left-handed and right-handed elliptically polarized
electric fields. It should be noted that the breakdown in reciprocity becomes
rather pronounced when the emitted photon energy exceeds the band-gap energy,
which is $E_{g}=7\hbar\omega_{0}$ in our calculation. The origin of this
feature is conjectured to be as follows. When the emitted photon energy is
larger than the band-gap energy, the number of excitation channels relevant to
HHG largely increases in comparison with those for the low-order harmonics
below the band-gap energy. The increase in the available channels contributes
to emergence of the plateau structure in the HHG spectraTamaya and Kato
(2019), as well as oscillatory behavior as a function of field strengthXia
_et al._ (2020). Thus, we suppose that the sensitivity of the high-order
harmonics above the band-gap energy to the external strain effect is caused by
an increase in the excitation channels. This breakdown of reciprocity in HHG
that is sensitive to shear strain may be used for applications such as
spatially-resolved measurement of lattice deformation. The features revealed
here may also be useful for developing a mechanical control of HHG
ellipticity.

Note that the strain effect (piezo-optic effect) should appear in other
semiconductors, such as AlAs and InAs, in accordance with the same formulation
based on the Luttinger-Kohn-Bir-Pikus model. By performing numerical
calculations with different Luttinger parameters, we can easily identify
similar properties of HHG in these materials. We also suppose that the piezo-
optic effect obtained here is not specific to materials covered by the
Luttinger-Kohn-Bir-Pikus model and that it appears in various systems, because
the symmetry reduction of the lattice structure can be induced by shear
strain. For a general discussion, we need to extend our theory to take the
details of the band structure of the materials as well as the corresponding
Bloch wavefunctions into account. This will be left as a future problem.

Finally, we comment on how the non-perturbative effect in HHG is important for
obtaining a large piezo-optical effect. According to the previous workTamaya
_et al._ (2016a); Yoshikawa  _et al._ (2017), the ellipticity dependence of
HHG parallel to the minor axis ($x$-axis) is not so pronounced in the
multiphoton absorption (perturbative) regime. Therefore, to obtain a
sufficiently large piezo-optical effect, we have to apply a strong field so
that non-perturbative regimes, which were identified as the AC Zener or the
semimetal regimes in Ref. Tamaya  _et al._ , 2016b, are realized. Because the
amplitude of the perturbation is roughly estimated by the ratio between the
field intensity (the Rabi frequency) and the band-gap energy, it is expected
that the piezo-optical effect discussed here can be observed more clearly by
using narrow-gap and zero-gap semiconductors, such as InSb and graphene.

##  IV Conclusion

We theoretically investigated the shear-strain effect of HHG in GaAs. By
constructing a theoretical framework based on the Luttinger-Kohn-Bir-Pikus
model, we calculated the spectra of HHG emissions parallel to the major and
minor axes. Our numerical results for linearly polarized incident light
implied that shear-stained materials have optical activity (non-reciprocity),
i.e., different responses to right-handed and left-handed elliptically
polarized electric fields. To verify this conjecture, we calculated the
ellipticity dependence of HHG emitted from shear-strained GaAs. Consequently,
we found a breakdown in reciprocity with respect to inversion of the
ellipticity. We also found that this breakdown is much more pronounced for
higher-order harmonics that exceed the band-gap energy. These features can be
understood in terms of a reduction in the symmetry of the crystal structure of
GaAs inducing off-diagonal elements in the dielectric tensor $\epsilon_{zx}$.
The conclusions presented in this paper are generally applicable to various
semiconductors because the symmetry reduction can be induced by the shear-
strain effect. Our study thus provides a foundation for strain engineering of
nonlinear optics, including HHG, that is impossible in gaseous media.

T. T. and T. K. gratefully acknowledge support from the Japan Society for the
Promotion of Science (JSPS KAKENHI Grants No. JP19K14624 and No. JP20K03831).

##  Appendix A Basis transformation

The eigen wavefunctions at the $\Gamma$ point,
$\ket{\Psi_{n{\bm{k}}=0}}=\ket{u_{n}}$ ($n=1,2,\cdots,8$) are given as a
superposition of atomic orbitals. In the presence of the spin-orbit
interaction, they are categorized by $\ket{J,{J}_{z}}$, where $J$ is the total
orbital angular momentum and $J_{z}$ is the $z$-component of the angular
momentum:

$\displaystyle\ket{u_{1}}$
$\displaystyle\equiv\Ket{\frac{1}{2},+\frac{1}{2}}=\ket{s\\!\uparrow},$
$\displaystyle\ket{u_{2}}$
$\displaystyle\equiv\Ket{\frac{3}{2},+\frac{3}{2}}=\frac{i}{\sqrt{2}}(\ket{p_{x}\\!\uparrow}+i\ket{p_{y}\\!\uparrow}),$
$\displaystyle\ket{u_{3}}$
$\displaystyle\equiv\Ket{\frac{3}{2},+\frac{1}{2}}=\frac{i}{\sqrt{6}}(\ket{p_{x}\\!\downarrow}+i\ket{p_{y}\\!\downarrow}-2\ket{p_{z}\\!\uparrow}),$
$\displaystyle\ket{u_{4}}$
$\displaystyle\equiv\Ket{\frac{1}{2},+\frac{1}{2}}=\frac{i}{\sqrt{3}}(\ket{p_{x}\\!\downarrow}+i\ket{p_{y}\\!\downarrow}+\ket{p_{z}\\!\uparrow}).$

Here, we define $\ket{s\sigma}$, $\ket{p_{x}\sigma}$, $\ket{p_{y}\sigma}$, and
$\ket{p_{z}\sigma}$ to be the $s$-, $p_{x}$-, $p_{y}$-, and $p_{z}$-like
wavefunctions, respectively, for spin components $\sigma=\uparrow$ or
$\downarrow$. The remaining set of Bloch basis states are expressed as

$\displaystyle\ket{u_{5}}$
$\displaystyle\equiv\Ket{\frac{1}{2},-\frac{1}{2}}=-\ket{s\\!\downarrow},$
(58) $\displaystyle\ket{u_{6}}$
$\displaystyle\equiv\Ket{\frac{3}{2},-\frac{3}{2}}=-\frac{i}{\sqrt{2}}(\ket{p_{x}\\!\downarrow}-i\ket{p_{y}\\!\downarrow}),$
(59) $\displaystyle\ket{u_{7}}$
$\displaystyle\equiv\Ket{\frac{3}{2},-\frac{1}{2}}=\frac{i}{\sqrt{6}}(\ket{p_{x}\\!\uparrow}-i\ket{p_{y}\\!\uparrow}+2\ket{p_{z}\\!\downarrow}),$
(60) $\displaystyle\ket{u_{8}}$
$\displaystyle\equiv\Ket{\frac{1}{2},-\frac{1}{2}}=\frac{i}{\sqrt{3}}(\ket{p_{x}\\!\uparrow}-i\ket{p_{y}\\!\uparrow}-\ket{p_{z}\\!\downarrow}).$
(61)

From these expressions, one can easily obtain the matrix element of the
unitary operation, $(U)_{nm}=\braket{u_{n}}{v_{m}}$.

## References

  * Boyd (2008) R. W. Boyd,  _Nonlinear Optics_ , 3rd ed. (2008).
  * Yariv and Yeh (1984) A. Yariv and P. Yeh,  _Optical Waves in Crystals_ (1984).
  * Shen (1984) Y. R. Shen,  _The Principles of Nonlinear Optics_ (1984).
  * Corkum (1993) P. B. Corkum, Phys. Rev. Lett. 71, 1994 (1993).
  * Protopapas  _et al._ (1997) M. Protopapas, C. H. Keitel, and P. L. Knight, Rep. Prog. Phys. 60, 389 (1997).
  * Brabec and Krausz (2000) T. Brabec and F. Krausz, Rev. Mod. Phys. 72, 545 (2000).
  * Agostini and DiMauro (2004) P. Agostini and L. F. DiMauro, Rep. Prog. Phys. 67, 813 (2004).
  * Corkum and Krausz (2007) P. B. Corkum and F. Krausz, Nature Phys. 3, 381 (2007).
  * Krausz and Ivanov (2009) F. Krausz and M. Ivanov, Rev. Mod. Phys. 81, 163 (2009).
  * Ghimire  _et al._ (2011) S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and D. A. Reis, Nature Phys. 7, 138 (2011).
  * Schubert  _et al._ (2014) O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, Nature Photon. 8, 119 (2014).
  * Luu  _et al._ (2015) T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T. Hassan, and E. Goulielmakis, Nature 521, 498 (2015).
  * Hohenleutner  _et al._ (2015) M. Hohenleutner, F. Langer, O. Schubert, M. Knorr, U. Huttner, S. Koch, M. Kira, and R. Huber, Nature 523, 572 (2015).
  * Vampa  _et al._ (2015) G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, and P. B. Corkum, Nature 522, 462 (2015).
  * Ghimire and Reis (2019) S. Ghimire and D. A. Reis, Nature Phys. 15, 10 (2019).
  * Ndabashimiye  _et al._ (2016) G. Ndabashimiye, S. Ghimire, M. Wu, D. A. Browne, K. J. Schafer, M. B. Gaarde, and D. A. Reis, Nature 534, 520 (2016).
  * Lanin  _et al._ (2017) A. A. Lanin, E. A. Stepanov, A. B. Fedotov, and A. M. Zheltikov, Optica 4, 516 (2017).
  * Liu  _et al._ (2017a) H. Liu, Y. Li, Y. S. You, S. Ghimire, T. F. Heinz, and D. A. Reis, Nature Phys. 13, 262 (2017a).
  * You  _et al._ (2017) Y. S. You, D. A. Reis, and S. Ghimire, Nature Phys. 13, 345 (2017).
  * Yoshikawa  _et al._ (2017) N. Yoshikawa, T. Tamaya, and K. Tanaka, Science 356, 736 (2017).
  * Kim  _et al._ (2017) H. Kim, S. Han, Y. W. Kim, S. Kim, and S.-W. Kim, ACS Photonics 4, 1627 (2017).
  * Jiang  _et al._ (2018) S. Jiang, J. Chen, H. Wei, C. Yu, R. Lu, and C. D. Lin, Phys. Rev. Lett. 120, 253201 (2018).
  * Langer  _et al._ (2018) F. Langer, C. Schmid, S. Schlauderer, M. Gmitra, J. Fabian, P. Nagler, C. Schüller, T. Korn, P. Hawkins, J. Steiner, U. Huttner, S. Koch, M. Kira, and R. Huber, Nature 557, 76 (2018).
  * G. Vampa and Y. S. You and H. Liu and S. Ghimire and D. A. Reis (2018) G. Vampa and Y. S. You and H. Liu and S. Ghimire and D. A. Reis, Opt. Express 26, 12210 (2018).
  * Silva  _et al._ (2018) R. E. F. Silva, I. V. Blinov, A. N. Rubtsov, O. Smirnova, and M. Ivanov, Nature Photon. 12, 266 (2018).
  * Hirori  _et al._ (2019) H. Hirori, P. Xia, Y. Shinohara, T. Otobe, Y. Sanari, H. Tahara, N. Ishii, J. Itatani, K. L. Ishikawa, T. Aharen, M. Ozaki, A. Wakamiya, and Y. Kanemitsu, APL Mater. 7, 041107 (2019).
  * Cheng  _et al._ (2020) B. Cheng, N. Kanda, T. N. Ikeda, T. Matsuda, P. Xia, T. Schumann, S. Stemmer, J. Itatani, N. P. Armitage, and R. Matsunaga, Phys. Rev. Lett. 124, 117402 (2020).
  * Smith (1954) C. S. Smith, Physical review 94, 42 (1954).
  * Llordes  _et al._ (2012) A. Llordes, A. Palau, J. Gázquez, M. Coll, R. Vlad, A. Pomar, J. Arbiol, R. Guzman, S. Ye, V. Rouco, _et al._ , Nature materials 11, 329 (2012).
  * Lü  _et al._ (2016) X. Lü, Y. Wang, C. C. Stoumpos, Q. Hu, X. Guo, H. Chen, L. Yang, J. S. Smith, W. Yang, Y. Zhao, _et al._ , Advanced Materials 28, 8663 (2016).
  * Liu  _et al._ (2017b) G. Liu, L. Kong, J. Gong, W. Yang, H.-k. Mao, Q. Hu, Z. Liu, R. D. Schaller, D. Zhang, and T. Xu, Advanced Functional Materials 27, 1604208 (2017b).
  * Chen  _et al._ (2018) B. Chen, T. Li, Q. Dong, E. Mosconi, J. Song, Z. Chen, Y. Deng, Y. Liu, S. Ducharme, A. Gruverman, _et al._ , Nature materials 17, 1020 (2018).
  * Zhao  _et al._ (2017) J. Zhao, Y. Deng, H. Wei, X. Zheng, Z. Yu, Y. Shao, J. E. Shield, and J. Huang, Science advances 3, eaao5616 (2017).
  * Zhu  _et al._ (2019) C. Zhu, X. Niu, Y. Fu, N. Li, C. Hu, Y. Chen, X. He, G. Na, P. Liu, H. Zai, _et al._ , Nature communications 10, 1 (2019).
  * Steele  _et al._ (2019) J. A. Steele, H. Jin, I. Dovgaliuk, R. F. Berger, T. Braeckevelt, H. Yuan, C. Martin, E. Solano, K. Lejaeghere, S. M. Rogge, _et al._ , Science 365, 679 (2019).
  * Peng  _et al._ (2020) Z. Peng, X. Chen, Y. Fan, D. J. Srolovitz, and D. Lei, Light: Science & Applications 9, 1 (2020).
  * Conley  _et al._ (2013) H. J. Conley, B. Wang, J. I. Ziegler, R. F. Haglund Jr, S. T. Pantelides, and K. I. Bolotin, Nano letters 13, 3626 (2013).
  * He  _et al._ (2013) K. He, C. Poole, K. F. Mak, and J. Shan, Nano letters 13, 2931 (2013).
  * Frank  _et al._ (2011) O. Frank, M. Mohr, J. Maultzsch, C. Thomsen, I. Riaz, R. Jalil, K. S. Novoselov, G. Tsoukleri, J. Parthenios, K. Papagelis, _et al._ , Acs Nano 5, 2231 (2011).
  * Shi  _et al._ (2013) H. Shi, H. Pan, Y.-W. Zhang, and B. I. Yakobson, Physical Review B 87, 155304 (2013).
  * Castellanos-Gomez  _et al._ (2015) A. Castellanos-Gomez, V. Singh, H. S. van der Zant, and G. A. Steele, Annalen der Physik 527, 27 (2015).
  * Bertolazzi  _et al._ (2011) S. Bertolazzi, J. Brivio, and A. Kis, ACS nano 5, 9703 (2011).
  * Castellanos-Gomez  _et al._ (2012) A. Castellanos-Gomez, M. Poot, G. A. Steele, H. S. Van Der Zant, N. Agraït, and G. Rubio-Bollinger, Advanced materials 24, 772 (2012).
  * Bao  _et al._ (2017) S. Bao, D. Kim, C. Onwukaeme, S. Gupta, K. Saraswat, K. H. Lee, Y. Kim, D. Min, Y. Jung, H. Qiu, _et al._ , Nature communications 8, 1 (2017).
  * Chen  _et al._ (2020) Y. Chen, Y. Lei, Y. Li, Y. Yu, J. Cai, M.-H. Chiu, R. Rao, Y. Gu, C. Wang, W. Choi, _et al._ , Nature 577, 209 (2020).
  * Xiong  _et al._ (2020) P. Xiong, F. Zhang, X. Zhang, S. Wang, H. Liu, B. Sun, J. Zhang, Y. Sun, R. Ma, Y. Bando, _et al._ , Nature communications 11, 1 (2020).
  * Shi and Chen (2019) J. Shi and L.-Q. Chen, “Strain engineering in functional materials,” (2019).
  * Wang  _et al._ (2015) Y. Wang, X. Lu, W. Yang, T. Wen, L. Yang, X. Ren, L. Wang, Z. Lin, and Y. Zhao, Journal of the American Chemical Society 137, 11144 (2015).
  * Barron (2009) L. D. Barron,  _Molecular light scattering and optical activity_ (Cambridge University Press, 2009).
  * Mason (1982) S. F. Mason,  _Molecular optical activity and the chiral discriminations_ (Cambridge University Press, 1982).
  * Wang  _et al._ (2020) Z. Wang, S. Jiang, G. Yuan, T. Wu, C. Li, C. Qian, C. Jin, C. Yu, W. Hua, and R. Lu, SCIENCE CHINA Physics, Mechanics & Astronomy 63, 257311 (2020).
  * Shao  _et al._ (2019) T.-J. Shao, Y. Xu, X.-H. Huang, and X.-B. Bian, Physical Review A 99, 013432 (2019).
  * Luttinger and Kohn (1955) J. M. Luttinger and W. Kohn, Phys. Rev. 97, 869 (1955).
  * Ahn  _et al._ (1995) D. Ahn, S. J. Yoon, S. L. Chuang, and C. Chang, Journal of Applied Physics 78, 2489 (1995), https://doi.org/10.1063/1.360103 .
  * Dargys (2002) A. Dargys, Phys. Rev. B 66, 165216 (2002).
  * Pryor (1998) C. Pryor, Phys. Rev. B 57, 7190 (1998).
  * Tomić  _et al._ (2006) S. Tomić, A. G. Sunderland, and I. J. Bush, J. Mater. Chem. 16, 1963 (2006).
  * Chuang (1991) S. L. Chuang, Phys. Rev. B 43, 9649 (1991).
  * Pfeffer and Zawadzki (1996) P. Pfeffer and W. Zawadzki, Phys. Rev. B 53, 12813 (1996).
  * Sytnyk and Melnik (2018) D. Sytnyk and R. Melnik, “The luttinger-kohn theory for multiband hamiltonians: A revision of ellipticity requirements,” (2018), arXiv:1808.06988 [cond-mat.mtrl-sci] .
  * Bir (1974) G. L. i. Bir,  _Symmetry and strain-induced effects in semiconductors_ (1974).
  * Luque  _et al._ (2015) A. Luque, A. Panchak, A. Mellor, A. Vlasov, A. Martí, and V. Andreev, Solar Energy Materials and Solar Cells 141, 39 (2015).
  * Harrison (2012) W. A. Harrison,  _Electronic structure and the properties of solids: the physics of the chemical bond_ (Courier Corporation, 2012).
  * Kittel and Fong (1963) C. Kittel and C.-y. Fong,  _Quantum theory of solids_ , Vol. 5 (Wiley New York, 1963).
  * Peter and Cardona (2010) Y. Peter and M. Cardona,  _Fundamentals of semiconductors: physics and materials properties_ (Springer Science & Business Media, 2010).
  * Voon and Willatzen (2009) L. C. L. Y. Voon and M. Willatzen,  _The kp method: electronic properties of semiconductors_ (Springer Science & Business Media, 2009).
  * Yonenaga and Sumino (1992) I. Yonenaga and K. Sumino, Journal of applied physics 71, 4249 (1992).
  * Tamaya  _et al._ (2016a) T. Tamaya, A. Ishikawa, T. Ogawa, and K. Tanaka, Phys. Rev. B 94, 241107(R) (2016a).
  * Tamaya  _et al._ (2016b) T. Tamaya, A. Ishikawa, T. Ogawa, and K. Tanaka, Phys. Rev. Lett. 116, 016601 (2016b).
  * Powell (2010) R. C. Powell,  _Symmetry, group theory, and the physical properties of crystals_ , Vol. 172 (Springer, 2010).
  * Tamaya and Kato (2019) T. Tamaya and T. Kato, Phys. Rev. B 100, 081203(R) (2019).
  * Xia  _et al._ (2020) P. Xia, T. Tamaya, C. Kim, F. Lu, T. Kanai, N. Ishii, J. Itatani, H. Akiyama, and T. Kato, “High-harmonic generation in gaas beyond the perturbative regime,” (2020), arXiv:2004.04492 [cond-mat.mes-hall] .


# Offline Recommender System Evaluation under Unobserved Confounding

Olivier Jeunen  ShareChatEdinburghUK and  Ben London  Amazon MusicSeattleWAUSA

(2023)

###### Abstract.

Off-Policy Estimation (OPE) methods allow us to learn and evaluate decision-
making policies from logged data. This makes them an attractive choice for the
offline evaluation of recommender systems, and several recent works have
reported successful adoption of OPE methods to this end. An important
assumption that makes this work is the absence of unobserved confounders:
random variables that influence both actions and rewards at data collection
time. Because the data collection policy is typically under the practitioner’s
control, the unconfoundedness assumption is often left implicit, and its
violations are rarely dealt with in the existing literature.

This work aims to highlight the problems that arise when performing off-policy
estimation in the presence of unobserved confounders, specifically focusing on
a recommendation use-case. We focus on policy-based estimators, where the
logging propensities are learned from logged data. We characterise the
statistical bias that arises due to confounding, and show how existing
diagnostics are unable to uncover such cases. Because the bias depends
directly on the _true_ and unobserved logging propensities, it is non-
identifiable. As the unconfoundedness assumption is famously untestable, this
becomes especially problematic. This paper emphasises this common, yet often
overlooked issue. Through synthetic data, we empirically show how naïve
propensity estimation under confounding can lead to severely biased metric
estimates that are allowed to fly under the radar. We aim to cultivate an
awareness among researchers and practitioners of this important problem, and
touch upon potential research directions towards mitigating its effects.

††journalyear: 2023††copyright: none††conference: ACM RecSys; September 19,
2023; Singapore††booktitle: CONSEQUENCES Workshop at ACM RecSys
(CONSEQUENCES@RecSys’23), September 19, 2023, Singapore

##  1\. Introduction & Motivation

Inferring cause and effect from observational data is not a straightforward
task, due to the presence of “ _unobserved confounders_ ” (Pearl, 2009).
Indeed, if an unobserved variable exists that influences both the treatment
and its outcome, this can easily lead to biased estimates of the treatment
effect. In the extreme case where the estimate changes sign, this is known as
“Simpson’s Paradox” (Julious and Mullee, 1994). Bottou et al. describe how to
avoid such unpleasantries in counterfactual learning scenarios, by carefully
modelling the data-generating process (Bottou et al., 2013). Jadidinejad et
al. describe how an instance of Simpson’s Paradox is prevalent to occur in the
offline evaluation of recommender systems—as the test data is often influenced
by an unknown recommendation policy that takes the role of the confounder
(Jadidinejad et al., 2021). Their proposed solution to this problem involves
_Inverse Propensity Score_ (IPS) weighting (Owen, 2013, Ch. 9), with
propensities that are estimated from logged data. Estimating propensities is
common practice in cases where the true propensities are unknown (Yang et al.,
2018; Ai et al., 2018), and empirical results indicate that IPS with estimated
propensities can even lead to favourable variance (Hanna et al., 2019).
Whether to simplify experimental setup or to deal with missing information in
publicly available datasets, these works often make the limiting assumption
that the propensities are independent of user or context (see (Jadidinejad et
al., 2021, §6.3) or (Yang et al., 2018, §3.3)). In the very likely case that
these assumptions are violated, this implies the presence of unobserved
confounders.111An analogous problem occurs in ranking applications when
contextual independence is incorrectly assumed for position bias estimates
(Fang et al., 2019; Jeunen, 2023). Existing diagnostics for validating logged
bandit feedback cannot detect these issues (Li et al., 2015; London and
Joachims, 2022), as we will show formally.

Unobserved confounding is a well-known issue in general off-policy
reinforcement learning, and several methods have been proposed to deal with it
in the recent literature. They typically leverage additional data (such as
interventions (Wang et al., 2021; Gasse et al., 2021) or instrumental
variables (Xu et al., 2023)) or make assumptions about the nature of
confounding variables (Namkoong et al., 2020; Bennett et al., 2021; Kausik et
al., 2023). Analogously, instrumental variables have been leveraged to test
for the unconfoundedness assumption (de Luna and Johansson, 2014), and other
statistical methods have been proposed to assess the sensitivity of results to
potential confounders (Liu et al., 2013). Nevertheless, in the absence of
additional tools, unconfoundedness is a famously _untestable_ assumption,
rendering its effects especially troublesome. Focusing on the off-policy
bandit setting with a guiding example in recommendation, we aim to answer:

“ _Can we reliably select the optimal policy from a set of competing policies,
under unobserved confounding?_ ”

##  2\. Off-Policy Estimation in the Presence of Unobserved Confounders

Figure 1. Probabilistic Graphical Model (PGM) for our setup.

[Figure described in text.]Figure described in text.

Throughout this work, we denote random variables as $X$, with specific
instances as $x\in\mathcal{X}$. A contextual bandit problem consists of
_contexts_ $X$ (e.g. user and context features), _actions_ $A$ (e.g. item
recommendations), and _rewards_ $R$ (e.g. clicks, streams, revenue). Rewards
are causally influenced by both contexts and actions, as illustrated by the
edges in the causal graph shown in Figure 1. A contextual _policy_ $\pi$
determines which actions are selected (or sampled), thereby inducing a
probability distribution over $A$, which is often denoted with the shorthand
$\pi(a|x)\coloneqq\mathsf{P}(A=a|X=x;\Pi=\pi)$. A policy’s effectiveness is
measured by the expected reward obtained when selecting actions according to
that policy: $\mathbb{E}_{x}\mathbb{E}_{a\sim\pi(\cdot|x)}[R]$. In a
recommendation application, this value can be estimated by deploying the
policy in an online experiment. However, since such experiments are typically
costly, and we may have many policies to evaluate, we would rather obtain
reward estimates by other means.

Suppose there is an existing deployed policy, called the _logging policy_
$\pi_{0}$, with which we collect a dataset
$\mathcal{D}_{0}\coloneqq\\{(a_{i},r_{i})_{i=1}^{N}\\}$. We will assume, as is
often the case, that the logging policy, and the contextual covariates it
uses, are unobservable, as indicated by the dashed nodes and edges in Figure
1. Our goal is to leverage data logged under $\pi_{0}$ to estimate the
expected reward under $\pi$—a problem often referred to as _off-policy_
estimation (Vasile et al., 2020; Saito and Joachims, 2021).

One simple off-policy estimator is the _Direct Method_ (DM). DM computes the
expected reward under $\pi$ (Eq. 1) using a model $\widehat{R}^{A}_{\rm
DM}(a)$ that estimates the reward for every available action. Since we assume
that contextual covariates are unavailable, the best we can do for
$\widehat{R}^{A}_{\rm DM}(a)$ is to naïvely count the observed rewards for
every action (Eq. 2).  (1) $\widehat{R}_{\rm
DM}(\pi)=\sum_{a\in\mathcal{A}}\widehat{R}^{A}_{\rm DM}(a)\pi(a)$ (2)
$\widehat{R}^{A}_{\rm
DM}(a)=\frac{\sum_{(a_{i},r_{i})\in\mathcal{D}_{0}}\mathbbm{1}\\{a_{i}=a\\}\cdot
r_{i}}{\sum_{(a_{i},r_{i})\in\mathcal{D}_{0}}\mathbbm{1}\\{a_{i}=a\\}}$

Unfortunately, this estimator is biased, for two reasons:  (a) it does not
take into account the covariates $X$ (i.e., the model is mis-specified), and
(b) it ignores the selection bias from the logging policy $\pi_{0}$,
influencing the estimates in Eq. 2.

In theory, we can bypass both the model mis-specification and selection bias
problems by leveraging the _ideal_ IPS estimator (Eq. 3), which is provably
unbiased. Importantly, ideal IPS requires access to both the contextual
covariates and the exact action probabilities (_propensities_) under
$\pi_{0}$—which we assume are unavailable. Accordingly, we will adopt the
common practice of using estimated logging propensities $\widehat{\pi}_{0}$
for IPS (Eq. 4). As the estimated propensities cannot properly consider all
covariates, this leads to unobserved confounding.

(3) $\widehat{R}_{\rm ideal-
IPS}(\pi)=\frac{1}{\left|\mathcal{D}_{0}\right|}\sum_{({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}x_{i}},a_{i},r_{i})\in\mathcal{D}_{0}}r_{i}\frac{\pi(a_{i})}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pi_{0}(a_{i}|x_{i})}}$

(4) $\widehat{R}_{\rm estim-
IPS}(\pi)=\frac{1}{\left|\mathcal{D}_{0}\right|}\sum_{(a_{i},r_{i})\in\mathcal{D}_{0}}r_{i}\frac{\pi(a_{i})}{\widehat{\pi}_{0}(a_{i})}$

Using the fact that the ideal IPS estimator is unbiased, we can quantify the
bias of the estimated IPS estimator as:

(5) $\mathbb{E}[\widehat{R}_{\rm estim-
IPS}(\pi)]-\mathbb{E}_{a\sim\pi}[R]=\mathbb{E}[\widehat{R}_{\rm estim-
IPS}(\pi)]-\mathbb{E}[\widehat{R}_{\rm ideal-
IPS}(\pi)]=\mathbb{E}\left[R\pi(A|X)\left(\frac{1}{\widehat{\pi}_{0}(A)}-\frac{1}{\pi_{0}(A|X)}\right)\right].$

To further illustrate our point, we resort to Pearl’s do-calculus framework
(Pearl, 2009). What OPE methods wish to estimate is the expected value of the
reward given that a new policy _intervenes_ on the action distribution. When
unobserved confounders are present, this _interventional_ quantity is _not_
equal to the _observational_ quantity we can estimate from logged data:
$\mathop{\mathbb{E}}\left[R|A=a\right]\neq\mathop{\mathbb{E}}\left[R|{\rm
do}(A=a)\right]$. Instead, we would require the “backdoor adjustment” to
obtain:

(6) $\mathop{\mathbb{E}}\left[R|{\rm
do}(A=a)\right]=\sum_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}x\in\mathcal{X}}}\mathop{\mathbb{E}}\left[R|A=a,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}X=x}\right].$

It should be clear that without access to $X$, this estimand is non-
identifiable, and this problem is not easily solved.

##  3\. Existing Diagnostics for Logging Propensities do not Uncover
Confounding Bias

Several diagnostics have been proposed in the literature to detect data
quality issues with logged bandit feedback. In particular, they try to uncover
cases where the two classical assumptions of the IPS estimator do not hold (Li
et al., 2015; London and Joachims, 2022):  (1) either the empirical action
frequencies in the data do not match those implied by the logged propensities,
or  (2) the logging policy does not have full support over the action space.
Note that the presence of unobserved confounders does _not_ automatically
violate these assumptions. As a result, the diagnostics that were proposed
will _not_ detect confounding bias. Logging propensities can be estimated by
empirically counting logged actions, as shown in Eq. 7. In doing so, we obtain
unbiased estimates of the true marginal action probabilities. Indeed,
$\lim_{N\to\infty}\widehat{\pi_{0}}(a)=\mathsf{P}(A=a|\Pi=\pi_{0})$.

Li et al. propose the use of _arithmetic_ and _harmonic_ mean tests to compare
empirical action frequencies with the logging propensities (Li et al., 2015).
As we _define_ the logging propensities to be equal to the empirical action
frequencies, it should be clear that this test will trivially pass.
Alternatively, London and Joachims propose to use the average importance
weight as a control variate, whose expectation should equal 1 for any target
policy $\pi$ (London and Joachims, 2022). Here as well, because the marginal
propensities are unbiased (Eq. 7), we can show that the control variate
remains unbiased as well (Eq. 8).

(7)
$\widehat{\pi}_{0}(a)=\frac{1}{|\mathcal{D}_{0}|}\sum_{(a_{i},r_{i})\in\mathcal{D}_{0}}\\!\\!\\!\\!\mathbbm{1}\\{a_{i}=a\\}\underset{\lim\limits_{N\to\infty}}{=}\mathsf{P}(A=a|\Pi=\pi_{0})=\sum_{x\in\mathcal{X}}\mathsf{P}(A=a|X=x,\Pi=\pi_{0})\mathsf{P}(X=x).$

###### Theorem 3.1.

When unobserved confounders are present and logging propensities are estimated
from empirical data (thus ignoring the confounders), the expected value of the
importance weights equals 1 for any target policy:

$\mathop{\mathbb{E}}_{\begin{subarray}{c}x\sim\mathsf{P}(X)\\\
a\sim\mathsf{P}(A|X=x,\Pi=\pi_{0})\end{subarray}}\left[\frac{\pi(a)}{\widehat{\pi}_{0}(a)}\right]=1.$

###### Proof.

(8)
$\begin{split}\mathop{\mathbb{E}}_{\begin{subarray}{c}x\sim\mathsf{P}(X)\\\
a\sim\mathsf{P}(A|X=x,\Pi=\pi_{0})\end{subarray}}\left[\frac{\pi(a)}{\widehat{\pi}_{0}(a)}\right]&=\sum_{a\in\mathcal{A}}\sum_{x\in\mathcal{X}}\frac{\pi(a)}{\widehat{\pi}_{0}(a)}\mathsf{P}(A=a|X=x,\Pi=\pi_{0})\mathsf{P}(X=x)\\\
&=\sum_{a\in\mathcal{A}}\frac{\pi(a)}{\widehat{\pi}_{0}(a)}\sum_{x\in\mathcal{X}}\mathsf{P}(A=a|X=x,\Pi=\pi_{0})\mathsf{P}(X=x)\\\
&\underset{\lim\limits_{N\to\infty}}{=}\sum_{a\in\mathcal{A}}\pi(a)\frac{\sum_{x\in\mathcal{X}}\mathsf{P}(A=a|X=x,\Pi=\pi_{0})\mathsf{P}(X=x)}{\sum_{x\in\mathcal{X}}\mathsf{P}(A=a|X=x,\Pi=\pi_{0})\mathsf{P}(X=x)}=\sum_{a\in\mathcal{A}}\pi(a)=1\qed\\\
\end{split}$

∎

As such, existing diagnostics are unable to detect issues of unobserved
confounding. This implies that the self-normalised IPS (SNIPS) estimator and
its extensions that adopt the above control variate to reduce the variance of
the IPS estimator, would exhibit the same bias as estimated IPS when
unobserved confounders are present (Swaminathan and Joachims, 2015; Joachims
et al., 2018).

##  4\. Empirical Validation of the Effects of Unobserved Confounding on
Synthetic Data

We now describe a guiding example, and provide a notebook that implements the
methods described earlier at github.com/olivierjeunen/confounding-
consequences-2023/. Consider a setting with two possible actions and a binary
covariate $X=\\{x_{0},x_{1}\\}$, following the distribution in Table 1(a)
(parameterised with $\alpha\in\left[\frac{1}{2},1\right]$). Rewards are
Bernoulli-distributed (Table 1(b)). The logging policy is contextual, taking a
suboptimal action with probability $\epsilon\in[0,1]$ (Table 1(c)). We can map
this to an intuitive setting: action $a_{1}$ is of general appeal to the
entire population (i.e. $R\perp\\!\\!\\!\perp X|A=a_{1}$); whereas action
$a_{0}$, on the other hand, is specifically appealing to a more niche user-
base (i.e. $\mathbb{E}[R|X=x_{0},A=a_{0}]>\mathbb{E}[R|X=x_{0},A=a_{1}]$, but
$\mathsf{P}(X=x_{0})<\mathsf{P}(X=x_{1})$). Estimates for logging propensities
can be obtained by empirical counting, as in Eq. 7. The expected value for
these estimated context-independent propensities is shown in Table 1(d).

Naïve propensity estimation methods suffer from confounding bias. We simulate
an off-policy estimation setup where we wish to evaluate deterministic
policies $\pi_{a}(a)\equiv 1$. We obtain $N=2\cdot 10^{6}$ samples from the
synthetic distribution described in Table 1, and compute the confounded
estimate $\widehat{R}_{\rm estim-IPS}$, as well as the unobservable ideal IPS
estimate $\widehat{R}_{\rm ideal-IPS}$. We vary both the level of selection
bias $\epsilon$ (over the x-axis), and the confounding distribution $\alpha$
(over columns) in Fig. 2, where the y-axis shows the estimated difference in
rewards from policies $\pi_{a_{1}}$ and $\pi_{a_{0}}$. We shade the positive
and negative regions in the plot to clearly visualise when an off-policy
estimator allows us to _correctly_ identify the optimal policy, or when it
does not. We observe that the IPS estimator with estimated propensities fails
considerably, in that it will incorrectly identify $\pi_{a_{0}}$ as the
reward-maximising policy. Only when $\epsilon$ is sufficiently high (i.e.
approaching a uniform logging policy for $\epsilon=0.5$, and hence no
confounding is present), $\widehat{R}_{\rm estim-IPS}$ is able to correctly
identify $\pi_{a_{1}}$. This shows that, even in simplified settings, the
estimates we obtain from IPS with confounded estimates lead to misleading
conclusions. Furthermore, existing diagnostics cannot detect these problems
when they occur.

$\mathsf{P}(X=x_{0})$ | $\mathsf{P}(X=x_{1})$
---|---
$1-\alpha$ | $\alpha$
(a) Covariate distribution

$\mathbf{\mathbb{E}[R|X,A]}$ | $a_{0}$ | $a_{1}$
---|---|---
$x_{0}$ | $1.0$ | $0.7$
$x_{1}$ | $0.0$ | $0.7$
(b) Reward distribution

$\mathbf{\pi_{0}(a|x)}$ | $a_{0}$ | $a_{1}$
---|---|---
$x_{0}$ | $1.0-\epsilon$ | $\epsilon$
$x_{1}$ | $\epsilon$ | $1.0-\epsilon$
(c) Logging policy distribution

$\mathbf{\widehat{\pi}_{0}(a)}$
---
$a_{0}$ | $(1-\epsilon)\alpha+\epsilon(1-\alpha)$
$a_{1}$ | $(1-\epsilon)(1-\alpha)+\epsilon\alpha$
(d) Logging propensity estimates

Table 1. Data distributions for a guiding example that highlights issues with
off-policy estimation under unobserved confounding. Figure 2. Estimated
differences in reward for various estimators. IPS with estimated propensities
suffers from unobserved confounding.

##  5\. Conclusions & Outlook

Unobserved confounders lead to biased estimates, both for DM- and IPS-based
methods. This problem has received considerable attention in the research
literature for general offline reinforcement learning use-cases, but the
literature dealing with these issues in recommendation settings remains
scarce. Our work highlights that this is problematic—especially in cases where
propensities are estimated under simplifying independence assumptions. In
doing so, we add to the literature identifying problematic practices that
might hamper progress in the field (Jeunen, 2019; Krichene and Rendle, 2020;
Sun, 2023; Jeunen et al., 2023).

## References

  * (1)
  * Ai et al. (2018) Qingyao Ai, Keping Bi, Cheng Luo, Jiafeng Guo, and W. Bruce Croft. 2018.  Unbiased Learning to Rank with Unbiased Propensity Estimation. In _The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval_ _(SIGIR ’18)_. ACM, 385–394.  https://doi.org/10.1145/3209978.3209986
  * Bennett et al. (2021) Andrew Bennett, Nathan Kallus, Lihong Li, and Ali Mousavi. 2021\.  Off-policy Evaluation in Infinite-Horizon Reinforcement Learning with Latent Confounders. In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ _(Proceedings of Machine Learning Research, Vol. 130)_ , Arindam Banerjee and Kenji Fukumizu (Eds.). PMLR, 1999–2007.  https://proceedings.mlr.press/v130/bennett21a.html
  * Bottou et al. (2013) Léon Bottou, Jonas Peters, Joaquin Quiñonero-Candela, Denis X. Charles, D. Max Chickering, Elon Portugaly, Dipankar Ray, Patrice Simard, and Ed Snelson. 2013\.  Counterfactual Reasoning and Learning Systems: The Example of Computational Advertising.  _Journal of Machine Learning Research_ 14, 101 (2013), 3207–3260.  http://jmlr.org/papers/v14/bottou13a.html
  * de Luna and Johansson (2014) Xavier de Luna and Per Johansson. 2014.  Testing for the Unconfoundedness Assumption Using an Instrumental Assumption.  _Journal of Causal Inference_ 2, 2 (2014), 187–199.  https://doi.org/doi:10.1515/jci-2013-0011
  * Fang et al. (2019) Zhichong Fang, Aman Agarwal, and Thorsten Joachims. 2019\.  Intervention Harvesting for Context-Dependent Examination-Bias Estimation. In _Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval_ _(SIGIR’19)_. ACM, 825–834.  https://doi.org/10.1145/3331184.3331238
  * Gasse et al. (2021) Maxime Gasse, Damien Grasset, Guillaume Gaudron, and Pierre-Yves Oudeyer. 2021.  Causal Reinforcement Learning using Observational and Interventional Data.  arXiv:2106.14421 [cs.LG]
  * Hanna et al. (2019) Josiah Hanna, Scott Niekum, and Peter Stone. 2019\.  Importance Sampling Policy Evaluation with an Estimated Behavior Policy. In _Proceedings of the 36th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 97)_ , Kamalika Chaudhuri and Ruslan Salakhutdinov (Eds.). PMLR, 2605–2613.  https://proceedings.mlr.press/v97/hanna19a.html
  * Jadidinejad et al. (2021) Amir H. Jadidinejad, Craig Macdonald, and Iadh Ounis. 2021\.  The Simpson’s Paradox in the Offline Evaluation of Recommendation Systems.  _ACM Trans. Inf. Syst._ 40, 1, Article 4 (sep 2021), 22 pages.  https://doi.org/10.1145/3458509
  * Jeunen (2019) Olivier Jeunen. 2019\.  Revisiting Offline Evaluation for Implicit-Feedback Recommender Systems. In _Proceedings of the 13th ACM Conference on Recommender Systems_ _(RecSys ’19)_. ACM, 596–600.  https://doi.org/10.1145/3298689.3347069
  * Jeunen (2023) Olivier Jeunen. 2023\.  A Probabilistic Position Bias Model for Short-Video Recommendation Feeds. In _Proceedings of the 17th ACM Conference on Recommender Systems_ _(RecSys ’23)_. ACM.
  * Jeunen et al. (2023) Olivier Jeunen, Ivan Potapov, and Aleksei Ustimenko. 2023\.  On (Normalised) Discounted Cumulative Gain as an Offline Evaluation Metric for Top-$n$ Recommendation.  arXiv:2307.15053 [cs.IR]
  * Joachims et al. (2018) Thorsten Joachims, Adith Swaminathan, and Maarten de Rijke. 2018\.  Deep Learning with Logged Bandit Feedback. In _International Conference on Learning Representations_.  https://openreview.net/forum?id=SJaP_-xAb
  * Julious and Mullee (1994) Steven A Julious and Mark A Mullee. 1994.  Confounding and Simpson’s paradox.  _Bmj_ 309, 6967 (1994), 1480–1481.
  * Kausik et al. (2023) Chinmaya Kausik, Yangyi Lu, Kevin Tan, Yixin Wang, and Ambuj Tewari. 2023.  Offline Policy Evaluation and Optimization under Confounding.  arXiv:2211.16583 [stat.ML]
  * Krichene and Rendle (2020) Walid Krichene and Steffen Rendle. 2020.  On Sampled Metrics for Item Recommendation. In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_ _(KDD ’20)_. ACM, 1748–1757.  https://doi.org/10.1145/3394486.3403226
  * Li et al. (2015) Lihong Li, Shunbao Chen, Jim Kleban, and Ankur Gupta. 2015\.  Counterfactual Estimation and Optimization of Click Metrics in Search Engines: A Case Study. In _Proceedings of the 24th International Conference on World Wide Web_ _(WWW ’15 Companion)_. ACM, 929–934.  https://doi.org/10.1145/2740908.2742562
  * Liu et al. (2013) Weiwei Liu, S Janet Kuramoto, and Elizabeth A Stuart. 2013\.  An introduction to sensitivity analysis for unobserved confounding in nonexperimental prevention research.  _Prev Sci_ 14, 6 (Dec. 2013), 570–580.
  * London and Joachims (2022) Ben London and Thorsten Joachims. 2022.  Control Variate Diagnostics for Detecting Problems in Logged Bandit Feedback. In _CONSEQUENCES+REVEAL Workshop at ACM RecSys ’22_ _(CONSEQUENCES+REVEAL ’22)_.
  * Namkoong et al. (2020) Hongseok Namkoong, Ramtin Keramati, Steve Yadlowsky, and Emma Brunskill. 2020.  Off-policy Policy Evaluation For Sequential Decisions Under Unobserved Confounding. In _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (Eds.), Vol. 33. Curran Associates, Inc., 18819–18831.  https://proceedings.neurips.cc/paper_files/paper/2020/file/da21bae82c02d1e2b8168d57cd3fbab7-Paper.pdf
  * Owen (2013) Art B. Owen. 2013\.  _Monte Carlo theory, methods and examples_.
  * Pearl (2009) Judea Pearl. 2009\.  _Causality_.  Cambridge university press.
  * Saito and Joachims (2021) Yuta Saito and Thorsten Joachims. 2021.  Counterfactual Learning and Evaluation for Recommender Systems: Foundations, Implementations, and Recent Advances. In _Proc. of the 15th ACM Conference on Recommender Systems_ _(RecSys ’21)_. ACM, 828–830.  https://doi.org/10.1145/3460231.3473320
  * Sun (2023) Aixin Sun. 2023\.  Take a Fresh Look at Recommender Systems from an Evaluation Standpoint. In _Proceedings of the 46th International ACM SIGIR Conference on Research and Development in Information Retrieval_ _(SIGIR ’23)_. ACM, 2629–2638.  https://doi.org/10.1145/3539618.3591931
  * Swaminathan and Joachims (2015) Adith Swaminathan and Thorsten Joachims. 2015.  The Self-Normalized Estimator for Counterfactual Learning. In _Advances in Neural Information Processing Systems_ , C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett (Eds.), Vol. 28. Curran Associates, Inc.  https://proceedings.neurips.cc/paper_files/paper/2015/file/39027dfad5138c9ca0c474d71db915c3-Paper.pdf
  * Vasile et al. (2020) Flavian Vasile, David Rohde, Olivier Jeunen, and Amine Benhalloum. 2020.  A Gentle Introduction to Recommendation as Counterfactual Policy Learning. In _Proceedings of the 28th ACM Conference on User Modeling, Adaptation and Personalization_ _(UMAP ’20)_. ACM, 392–393.  https://doi.org/10.1145/3340631.3398666
  * Wang et al. (2021) Lingxiao Wang, Zhuoran Yang, and Zhaoran Wang. 2021\.  Provably Efficient Causal Reinforcement Learning with Confounded Observational Data. In _Advances in Neural Information Processing Systems_ , M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (Eds.), Vol. 34. Curran Associates, Inc., 21164–21175.
  * Xu et al. (2023) Yang Xu, Jin Zhu, Chengchun Shi, Shikai Luo, and Rui Song. 2023.  An Instrumental Variable Approach to Confounded Off-Policy Evaluation. In _Proceedings of the 40th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 202)_ , Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (Eds.). PMLR, 38848–38880.  https://proceedings.mlr.press/v202/xu23x.html
  * Yang et al. (2018) Longqi Yang, Yin Cui, Yuan Xuan, Chenyang Wang, Serge Belongie, and Deborah Estrin. 2018\.  Unbiased Offline Recommender Evaluation for Missing-Not-at-Random Implicit Feedback. In _Proceedings of the 12th ACM Conference on Recommender Systems_ _(RecSys ’18)_. ACM, 279–287.  https://doi.org/10.1145/3240323.3240355