# Liouvillian exceptional points of an open driven two-level system

Nikhil Seshadri  Harvard University, Cambridge, MA 02138, USA  Anqi Li
Department of Chemistry & Biochemistry, University of California San Diego, La
Jolla, CA 92093, USA  Michael Galperin<EMAIL_ADDRESS>Department of
Chemistry & Biochemistry, University of California San Diego, La Jolla, CA
92093, USA

###### Abstract

We study the applicability of the Liouvillian exceptional points (LEPs)
approach to nanoscale open quantum systems. A generic model of the driven two-
level system in a thermal environment is analyzed within the nonequilibrium
Green’s function (NEGF) and Bloch quantum master equation (QME) formulations.
We derive the latter starting from the exact NEGF Dyson equations and
highlight the qualitative limitations of the LEP treatment by examining the
approximations employed in its derivation. We find that non-Markov character
of evolution in open quantum systems does not allow for the introduction of
the concept of exceptional points for a description of their dynamics.
Theoretical analysis is illustrated with numerical simulations.

##  I Introduction

Non-Hermitian quantum mechanics Moiseyev (2011) is an accepted way of treating
open quantum systems which is employed in many fields of theoretical research
from optics, opto-mechanics, and polaritonics, to quantum field theory,
molecular physics, and quantum transport. Complex values of operator spectra
in these considerations reflect the non-stationary character of system states
with the balance between gain and loss accounted for by the imaginary parts of
eigenvalues. The most non-trivial physics (such as unidirectional transport,
anomalous lasing and absorption, and chiral modes) takes place at and in the
vicinity of the degeneracies of the complex eigenvalues – exceptional points
(EPs).

Experimentally, EP behavior has been observed mostly in optics Miri and Alù
(2019), in the setting of a chaotic optical microcavity Lee  _et al._ (2009),
optical coupled systems with a complex index potential Rüter  _et al._ (2010),
and photonic lattices Hahn  _et al._ (2016). EP systems were suggested as a
platform for development of topological optoelectronics Doppler  _et al._
(2016); Ergoktas  _et al._ (2022). Recently, observations of EPs in single-
spin systems (nitrogen-vacancy centers in diamonds) were also reported Wu  _et
al._ (2019). The sensitivity of EP system responses to parameter changes led
to suggestions of employing EP systems as optical Wiersig (2020); Yang  _et
al._ (2023) and quantum Liang  _et al._ (2023) sensors. Decoherence
enhancement observed in the vicinity of EPs Zhang  _et al._ (2018); Naghiloo
_et al._ (2019) opens a way for the exploration of EPs for quantum information
processing. EP physics was also observed in polaritonic systems (exciton-
polaritons in semiconductor microcavities) Gao  _et al._ (2015) and in thermal
transport (chiral heat transport) Xu  _et al._ (2023).

The majority of theoretical considerations use effective non-Hermitian
Hamiltonians as operators describing EP physics Günther  _et al._ (2007);
Rotter (2009); Toroker and Peskin (2009); Uzdin  _et al._ (2011); Heiss
(2012); Garmon  _et al._ (2012); Delga  _et al._ (2014); Rotter and Bird
(2015); Yang  _et al._ (2020); Engelhardt and Cao (2022); Ferrier  _et al._
(2022); Engelhardt and Cao (2023); Li  _et al._ (2023). These operators are
formed by adding complex absorbing potentials (retarded and/or advanced
projections of self-energies) to Hermitian system Hamiltonians. Their
degeneracies, the Hamiltonian EPs (HEPs), are the focus of these studies.

Another non-Hermitian operator describing the evolution of open quantum
systems is the Liouvillian. Its degeneracies, Liouvillian EPs (LEPs), were
also discussed recently Pick  _et al._ (2019); Minganti  _et al._ (2019,
2020); Arkhipov  _et al._ (2020); Chimzak  _et al._ (2023) . Analytical
studies comparing HEPs and LEPs conclude that the two types of EPs have
essentially different properties and that they become equivalent only in the
semiclassical limit. Similar to HEPs, Liouvillian based analysis predicts non-
trivial behavior at or in vicinity of LEPs. For example, LEPs were shown to
represent a threshold between diffusive and ballistic motion in a 1d quantum
Lorentz gas Hashimoto  _et al._ (2016a, b). Enhancement of decoherence rate
Chen  _et al._ (2021); Khandelwal  _et al._ (2021); Larson and Qvarfort
(2023), possibility of chiral state transfer Chen  _et al._ (2022), and
optimization of steering towards a predesigned target state Kumar  _et al._
(2022) are predicted in the presence of LEPs. Finally, recent experiment
demonstrated enhanced performance of the single-ion quantum heat engine from
the LEPs Bu  _et al._ (2023).

Recently, we studied applicability of the concept of HEPs in nanoscale open
quantum systems Mukamel  _et al._ (2023). Utilizing a model of two vibrational
modes in a cavity we compared standard nonequilibrium Green’s function (NEGF)
with HEP based predictions. We derived the latter from the former and
discussed approximations required to reduce exact NEGF to approximate HEP
description. In particular, we showed that HEP disregards lesser and greater
projections of self-energy due to intra-system interactions while keeping its
retarded projection which makes the HEP treatment inconsistent and may lead to
qualitative failures. Another limiting factor of the HEP approach is its
Markov character.

Here, we present analysis of LEP based considerations starting from exact NEGF
treatment and exploring approximations necessary to reduce the latter to the
approximate LEP description. The two most basic and widely employed models for
LEP analysis are the driven two-level system (TLS) Am-Shallem  _et al._
(2015); Hatano (2019); Perina Jr  _et al._ (2022) and oscillator Tay (2023) in
a generic environment. We use the TLS as a model for comparison between NEGF
and LEP methods. Similar to our findings in Ref. Mukamel  _et al._ , 2023, LEP
is also limited by its Markov character. Nevertheless, contrary to the HEP,
the Liouvillian based treatment disregards the retarded projection of the
self-energies while keeping their lesser and greater projections. Some
limitations in the applicability of LEP methods to nanoscale open quantum
systems are illustrated with simulations comparing the NEGF and Bloch quantum
master equation (QME) results for driven TLS in a thermal environment.

In Section II we introduce the model and present its NEGF treatment. We then
utilize NEGF as a starting point for derivation of the Bloch QME and its
generalization which accounts for dissipation and discuss approximations
necessary to reduce exact NEGF treatment to approximate Redfield/Lindblad QME.
Section III compares results of numerical simulations performed within the
NEGF formulation and within the two types of the Bloch QME formulations.
Conclusions are drawn in Section IV.

##  II Driven TLS in a thermal environment

###  II.1 Model

We consider a two-level system that is driven by external classical field
$E(t)$ and dissipated by a thermal bath. The latter is continuum of Bose modes
$\\{\alpha\\}$. The Hamiltonian of this model is

$\hat{H}(t)=\hat{H}^{S}(t)+\hat{H}^{B}+\hat{V}^{SB}$ (1)

where $\hat{H}^{S}(t)$ and $\hat{H}^{B}$ describe decoupled system and bath,
respectively. $\hat{V}^{SB}$ is the system-bath coupling. Explicit expressions
for each of the terms are given by

$\begin{split}\hat{H}^{S}(t)&=\sum_{i=1,2}\varepsilon_{i}\hat{d}_{i}^{\dagger}\hat{d}_{i}-\mu
E(t)\left(\hat{d}_{1}^{\dagger}\hat{d}_{2}+\hat{d}_{2}^{\dagger}\hat{d}_{1}\right)\\\
\hat{H}^{B}&=\sum_{\alpha}\omega_{\alpha}\hat{b}_{\alpha}^{\dagger}\hat{b}_{\alpha}\\\
\hat{V}^{SB}&=\sum_{i,j=1,2}\sum_{\alpha}\left(V_{ij,\alpha}\left[\hat{d}_{i}^{\dagger}\hat{d}_{j}\right]^{\dagger}\hat{b}_{\alpha}+V_{\alpha,ij}\,\hat{b}^{\dagger}_{\alpha}\left[\hat{d}_{i}^{\dagger}\hat{d}_{j}\right]\right)\end{split}$
(2)

Here, $\hat{d}_{i}^{\dagger}$ ($\hat{d}_{i}$) and $\hat{b}_{\alpha}^{\dagger}$
($\hat{b}_{\alpha}$) creates (annihilates) an electron in level $i$ and an
excitation in mode $\alpha$, respectively. $\mu$ is the transition dipole
moment. The driving function is taken to be harmonic

$E(t)=E_{0}\cos(\omega_{0}t)$ (3)

In the following analysis, we assume $\varepsilon_{1}<\varepsilon_{2}$ and
consider coupling to the thermal bath in the rotating-wave approximation
(RWA); that is, $V_{21,\alpha}=V_{\alpha,21}=0$. We note that the RWA is
central for the derivation of the Bloch QME.

Figure 1:  Sketch of a model for optically driven two-level system (S) in
thermal environment (B).

###  II.2 NEGF formulation

Within the NEGF formulation, the central quantity of interest is the single-
particle Green’s function of the system defined on the Keldysh contour

$G_{ij}(\tau_{1},\tau_{2})\equiv-i\langle
T_{c}\,\hat{d}_{i}(\tau_{1})\,\hat{d}_{j}^{\dagger}(\tau_{2})\rangle.$ (4)

Here, $T_{c}$ is the contour ordering operator, $\tau_{1,2}$ are the contour
variables, and the creation (annihilation) operator
$\hat{d}_{j}^{\dagger}(\tau_{2})$ ($\hat{d}_{i}(\tau_{1})$) is in the
Heisenberg picture. Knowledge of $G_{ij}(\tau_{1},\tau_{2})$ allows for the
calculation of characteristics of the system and its responses to external
perturbations. In particular, in the single-electron subspace of the problem,
the system density matrix is given by the lesser projection of the Green’s
function (4) taken at equal times

$\rho_{ij}(t)=-i\,G^{<}_{ij}(t,t)$ (5)

This relation is central for comparison between the NEGF and Bloch quantum
master equation (QME) results.

The dynamics of the system is described by the Dyson equation for the Green’s
function (4)

$\displaystyle
i\frac{\partial}{\partial\tau_{1}}G_{ij}(\tau_{1},\tau_{2})=\delta_{i,j}\,\delta(\tau_{1},\tau_{2})$
(6)
$\displaystyle+\sum_{n=1,2}H_{in}^{S}(t_{1})\,G_{nj}(\tau_{1},\tau_{2})+\int_{c}d\tau\,\Sigma_{in}(\tau_{1},\tau)\,G_{nj}(\tau,\tau_{2})$

where $t_{1}$ is physical time corresponding to contour variable $\tau_{1}$
and $\Sigma(\tau_{1},\tau)$ is the self-energy due to the coupling of the
system to the bath. While the exact expression for the latter is not
accessible due to the many-body character of the system-bath coupling
$\hat{V}^{SB}$, an appropriate level of theory for future comparison with the
Bloch QME can be achieved by a second order diagrammatic expansion. Within
this (Hartree-Fock) approximation, the expression for the self-energy is (see
Appendix A for derivation)

$\displaystyle\Sigma_{ij}(\tau_{1},\tau_{2})=\delta(\tau_{1},\tau_{2})\,\sum_{n_{1},n_{2}}\int_{c}d\tau\,\rho_{n_{1}n_{2}}(t)$
$\displaystyle\qquad\qquad\,\,\,\times\bigg{[}\Pi_{ji,n_{2}n_{1}}(\tau_{1},\tau)+\Pi_{n_{1}n_{2},ij}(\tau,\tau_{1})\bigg{]}$
(7)
$\displaystyle+i\sum_{n_{1},n_{2}}\bigg{[}\Pi_{n_{1}i,n_{2}j}(\tau_{1},\tau_{2})+\Pi_{jn_{2},in_{1}}(\tau_{2},\tau_{1})\bigg{]}G_{n_{1}n_{2}}(\tau_{1},\tau_{2})$

Here,

$\Pi_{n_{1}n_{2},n_{3}n_{4}}(\tau_{1},\tau_{2})\equiv\sum_{\alpha}V_{n_{1}n_{2},\alpha}\,F^{(0)}_{\alpha}(\tau_{1},\tau_{2})\,V_{\alpha,n_{3}n_{4}}$
(8)

is the thermal bath-induced effective interaction between transitions
$n_{1}n_{2}$ and $n_{3}n_{4}$, and

$F^{(0)}_{\alpha}(\tau_{1},\tau_{2})\equiv-i\left\langle
T_{c}\,\hat{b}_{\alpha}(\tau_{1})\,\hat{b}_{\alpha}^{\dagger}(\tau_{2})\right\rangle_{0}$
(9)

is the Green’s function of free phonon mode $\alpha$ in the bath.

###  II.3 Bloch QME

The derivation of an approximate Redfield/Lindblad QME starts from the exact
equation-of-motion (EOM) for the density matrix given by (5), which is derived
within the NEGF formulation. The EOM is (see Appendix B for derivation)

$\displaystyle\frac{d}{dt}\rho_{ij}(t)$
$\displaystyle=i\,\omega_{ji}\,\rho_{ij}(t)-i\,\mu\,E(t)\left[\rho_{i\bar{j}}(t)-\rho_{\bar{i}j}(t)\right]$
(10)
$\displaystyle+\sum_{n,n_{1},n_{2}}\int_{-\infty}^{t}dt^{\prime}\,\bigg{[}\Pi^{>}_{ni,n_{1}n_{2}}(t-t^{\prime})\,G^{(2)\,<}_{n_{1}n_{2},nj}(t^{\prime},t)-\Pi^{<}_{ni,n_{1}n_{2}}(t-t^{\prime})\,G^{(2)\,>}_{n_{1}n_{2},nj}(t^{\prime},t)$
$\displaystyle\qquad\qquad\qquad\,\,\,\,\,-\Pi^{>}_{jn,n_{1}n_{2}}(t-t^{\prime})\,G^{(2)\,<}_{n_{1}n_{2},in}(t^{\prime},t)+\Pi^{<}_{jn,n_{1}n_{2}}(t-t^{\prime})\,G^{(2)\,>}_{n_{1}n_{2},in}(t^{\prime},t)$
$\displaystyle\qquad\qquad\qquad\,\,\,\,\,+G^{(2)\,<}_{ni,n_{1}n_{2}}(t,t^{\prime})\,\Pi^{>}_{n_{1}n_{2},nj}(t^{\prime}-t)-G^{(2)\,>}_{ni,n_{1}n_{2}}(t,t^{\prime})\,\Pi^{<}_{n_{1}n_{2},nj}(t^{\prime}-t)$
$\displaystyle\qquad\qquad\qquad\,\,\,\,\,-G^{(2)\,<}_{jn,n_{1}n_{2}}(t,t^{\prime})\,\Pi^{>}_{n_{1}n_{2},in}(t^{\prime}-t)+G^{(2)\,>}_{jn,n_{1}n_{2}}(t,t^{\prime})\,\Pi^{<}_{n_{1}n_{2},in}(t^{\prime}-t)\bigg{]}$

where $\bar{i},\bar{j}=2(1)$ for $i,j=1(2)$,
$\omega_{ji}\equiv\varepsilon_{j}-\varepsilon_{i}$, and

$G^{(2)}_{n_{1}n_{2},n_{3}n_{4}}(\tau_{1},\tau_{2})\equiv-i\langle
T_{c}\,\left[\hat{d}_{n_{1}}^{\dagger}\hat{d}_{n_{2}}\right](\tau_{1})\,\left[\hat{d}_{n_{3}}^{\dagger}\hat{d}_{n_{4}}\right]^{\dagger}(\tau_{2})\rangle$
(11)

is the two-particle Green’s function. Reducing the exact EOM (10) to the
Redfield/Lindblad QME requires approximating its right side with a Markov
dynamics. The Redfield/Lindblad QME can be obtained from the Green’s function
Dyson equation by employing a Kadanoff-Baym-like ansatz Haug and Jauho (2008)

$G^{(2)\,\gtrless}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{2})\approx\theta(t_{1}-t_{2})\,e^{-i\omega_{21}(t_{1}-t_{2})}\,G^{(2)\,\gtrless}_{n_{1}n_{2},n_{3}n_{4}}(t_{2},t_{2})+\theta(t_{2}-t_{1})\,e^{-i\omega_{43}(t_{1}-t_{2})}\,G^{(2)\,\gtrless}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{1})$
(12)

where $\theta(\ldots)$ is the Heaviside step function. Employing this ansatz
leads to the Bloch equations (see Appendix C for derivation)

$\frac{d}{dt}\begin{pmatrix}\rho_{11}(t)\\\ \rho_{22}(t)\\\ \rho_{12}(t)\\\
\rho_{21}(t)\end{pmatrix}=-i\begin{bmatrix}-iW_{2\leftarrow 1}&iW_{1\leftarrow
2}&\mu E(t)&-\mu E(t)\\\ iW_{2\leftarrow 1}&-iW_{1\leftarrow 2}&-\mu E(t)&\mu
E(t)\\\ \mu E(t)&-\mu E(t)&-\omega_{21}-iW_{d}&0\\\ -\mu E(t)&\mu
E(t)&0&\omega_{21}-iW_{d}\end{bmatrix}\begin{pmatrix}\rho_{11}(t)\\\
\rho_{22}(t)\\\ \rho_{12}(t)\\\ \rho_{21}(t)\end{pmatrix}$ (13)

Here,

$\begin{split}W_{2\leftarrow
1}&\equiv\Gamma_{12,12}(\omega_{21})\,N(\omega_{21}),\\\ W_{1\leftarrow
2}&\equiv\Gamma_{12,12}(\omega_{21})\,\left[1+N(\omega_{21})\right],\end{split}$
(14)

are the population transfer rates, and

$\begin{split}&W_{d}\equiv\frac{W_{2\leftarrow 1}+W_{1\leftarrow 2}}{2}\\\
&\quad\,\,\,+\frac{\Gamma_{11,11}(0)+\Gamma_{22,22}(0)}{2}\left[1+2\,N(0)\right]\end{split}$
(15)

is the dephasing rate. In Eqs. (14)-(15)

$\Gamma_{n_{1}n_{2},n_{3}n_{4}}(\omega)\equiv
2\pi\sum_{\alpha}V_{n_{1}n_{2},\alpha}\,V_{\alpha,n_{3}n_{4}}\,\delta(\omega-\omega_{\alpha})$
(16)

is the dissipation matrix.

Finally, by using (3) with $\omega_{0}-\omega_{21}\equiv\Delta\ll\omega_{0}$,
going into rotating frame of the field

$\tilde{\rho}_{12}(t)\equiv e^{-i\omega_{0}t}\,\rho_{12}(t),$ (17)

and introducing the spin operators

$\begin{split}\tilde{S}_{x}(t)&\equiv\tilde{\rho}_{21}(t)+\tilde{\rho}_{12}(t)\\\
\tilde{S}_{y}(t)&\equiv
i\left[\tilde{\rho}_{21}(t)-\tilde{\rho}_{12}(t)\right]\\\
S_{z}(t)&\equiv\rho_{22}(t)-\rho_{11}(t)\end{split}$ (18)

one can employ the rotating wave approximation (RWA) to express the Bloch QME
(13) as an EOM for the spin operator

$\begin{split}\frac{d}{dt}\begin{pmatrix}\tilde{S}_{x}(t)\\\
\tilde{S}_{y}(t)\\\
S_{z}(t)\end{pmatrix}&=\begin{bmatrix}-\frac{1}{T_{2}}&\Delta&0\\\
-\Delta&-\frac{1}{T_{2}}&\mu E_{0}\\\ 0&-\mu
E_{0}&-\frac{1}{T_{1}}\end{bmatrix}\begin{pmatrix}\tilde{S}_{x}(t)\\\
\tilde{S}_{y}(t)\\\ S_{z}(t)\end{pmatrix}+\begin{pmatrix}0\\\ 0\\\
\frac{S_{z}^{0}}{T_{1}}\end{pmatrix}\end{split}$ (19)

Here,

$\begin{split}&\frac{1}{T_{1}}\equiv W_{2\leftarrow 1}+W_{1\leftarrow
2},\qquad\frac{1}{T_{2}}\equiv W_{d},\\\ &S_{z}^{0}\equiv\frac{W_{2\leftarrow
1}-W_{1\leftarrow 2}}{W_{2\leftarrow 1}+W_{1\leftarrow 2}}\end{split}$ (20)

###  II.4 Generalized Bloch QMEs

While deriving the Bloch QME (13) from the exact EOM (10), one looses proper
non-Markov evolution and disregards dissipation. Note that while the former is
common for Hamiltonian and Liouvillian EP formulations Mukamel  _et al._
(2023), the latter is specific to Liouvillian EPs. Indeed, the Hamiltonian EP
formulation disregards the lesser/greater projections of the self-energy,
while the ansatz (12) misses the retarded projection, however, the generalized
Kadanoff-Baym ansatz (GKBA) in the NEGF literature Lipavský  _et al._ (1986);
Haug and Jauho (2008) does preserve information about dissipation. To
construct the Liouville space analog we follow procedure originally introduced
in Ref. Esposito and Galperin, 2009. This leads to (see Appendix D for
derivation)

$\begin{split}&G^{(2)\,\gtrless}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{2})\approx\\\
&\qquad
i\sum_{e,f}\bigg{[}\mathcal{G}^{r}_{n_{1}n_{2},ef}(t_{1}-t_{2})\,G^{(2)\,\gtrless}_{ef,n_{3}n_{4}}(t_{2},t_{2})\\\
&\qquad\quad\,\,-G^{(2)\,\gtrless}_{n_{1}n_{2},ef}(t_{1},t_{1})\,\mathcal{G}^{a}_{ef,n_{3}n_{4}}(t_{1}-t_{2})\bigg{]}\end{split}$
(21)

where

$\begin{split}\mathcal{G}^{r}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{2})&\equiv-i\theta(t_{1}-t_{2})\llangle
n_{2}n_{1}\rvert\,\mathcal{U}_{eff}(t_{1},t_{2})\,\lvert n_{4}n_{3}\rrangle\\\
\mathcal{G}^{a}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{2})&\equiv+i\theta(t_{2}-t_{1})\llangle
n_{2}n_{1}\rvert\,\mathcal{U}^{\dagger}_{eff}(t_{2},t_{1})\,\lvert
n_{4}n_{3}\rrangle\end{split}$ (22)

are the retarded and advanced Green’s functions in Liouville space,
respectively, and

$\mathcal{U}_{eff}(t_{1},t_{2})\equiv
T\,\exp\left[-i\int_{t_{2}}^{t_{1}}dt\,\mathcal{L}_{eff}(t)\right]$ (23)

is the Liouville space effective evolution operator, where
$\mathcal{L}_{eff}(t)$ defines time evolution in the system subspace of the
problem. We note that Eq.(21) is an approximation. The approximation is
introduced by employing projection operator (41) in exact expressions (39)
which makes (21) a second order contribution in infinite diagrammatic
expansion of the coupled system-bath evolution in strength of the system-bath
coupling.

We employ parts of the Redfield/Lindblad Liouvillian matrix in the right side
of Eq.(13), as the system evolution generator. In particular, retaining only
free evolution (i.e disregarding driving $\mu E(t)$ and dissipation
$W_{2\leftarrow 1}$, $W_{1\leftarrow 2}$, and $W_{d}$) reduces (21) to (12).

Keeping the dissipation, using (21) in (10), and assuming the Born-Markov
approximation leads to a generalized version of the Bloch QME, which retains
the same form (13), although with renormalized ($\bar{W}$) dissipation rates

$\displaystyle\bar{W}_{2\leftarrow
1}\equiv-i\int\frac{d\omega}{2\pi}\,\Pi^{<}_{12,12}(\omega)\,\mbox{Im}\left[\frac{1}{\omega-\omega_{21}+iW_{d}}\right]$
$\displaystyle\bar{W}_{1\leftarrow
2}\equiv-i\int\frac{d\omega}{2\pi}\,\Pi^{>}_{12,12}(\omega)\,\mbox{Im}\left[\frac{1}{\omega-\omega_{21}+iW_{d}}\right]$
$\displaystyle\bar{W}_{d}\equiv-\int\frac{d\omega}{2\pi}\,\bigg{[}\frac{\Pi^{<}_{11,11}(\omega)+\Pi^{>}_{22,22}(\omega)}{\omega+iW_{d}}$
(24)
$\displaystyle-\frac{\Pi^{>}_{11,11}(\omega)+\Pi^{<}_{22,22}(\omega)}{\omega-
iW_{d}}$
$\displaystyle+\frac{V^{R}_{11,1}\,\left[V^{R}\right]^{-1}_{1,11}\,\Pi^{>}_{12,12}(\omega)+V^{R}_{22,1}\,\left[V^{R}\right]^{-1}_{1,22}\,\Pi^{<}_{12,12}(\omega)}{\omega-\omega_{21}+i\delta}$
$\displaystyle+\frac{V^{R}_{11,2}\,\left[V^{R}\right]^{-1}_{2,11}\,\Pi^{>}_{12,12}(\omega)+V^{R}_{22,2}\,\left[V^{R}\right]^{-1}_{2,22}\,\Pi^{<}_{12,12}(\omega)}{\omega-\omega_{21}+i\left(W_{2\leftarrow
1}+W_{1\leftarrow 2}\right)}\bigg{]}$

Here, $\mathbf{V}^{R}$ is the right eigenvector of the Liouvillian matrix.
Note that while also keeping driving terms in the effective evolution is
possible, we will not pursue this direction because the accepted approach
regarding the derivation of the standard Bloch QME requires one to disregard
the driving term when deriving dissipators of the Liouvillian. Note also, that
using the Liouville space generalized Kadanoff-Baym ansatz on the Keldysh
anti-contour Esposito and Galperin (2010) would lead to the same form of the
generalized Bloch equation.

Finally, one can choose to solve the time-nonlocal (non-Markov) version of the
QME. Using (21) in (10) without the Born-Markov assumption leads to

$\begin{split}&\frac{d}{dt}\rho_{11}(t)=\mu\,E(t)\,2\,\mbox{Im}\left[\rho_{12}(t)\right]\\\
&+2\,\mbox{Im}\int_{-\infty}^{t}dt^{\prime}\,\bigg{(}\Pi^{<}_{12,12}(t-t^{\prime})\,\rho_{11}(t^{\prime})\\\
&\qquad-\Pi^{>}_{12,12}(t-t^{\prime})\left[1-\rho_{11}(t^{\prime})\right]e^{(i\omega_{21}-W_{d})(t-t^{\prime})}\bigg{)}\\\
&\frac{d}{dt}\rho_{12}(t)=i\omega_{12}\rho_{12}(t)-i\mu\,E(t)\left[2\,\rho_{11}(t)-1\right]\\\
&-i\int_{-\infty}^{t}dt^{\prime}\bigg{(}\left[\Pi_{11,11}^{>}(t-t^{\prime})+\Pi_{22,22}^{<}(t-t^{\prime})\right.\\\
&\quad\left.+\Pi_{11,11}^{<}(t^{\prime}-t)+\Pi_{22,22}^{>}(t^{\prime}-t)\right]e^{(i\omega_{21}-W_{d})(t-t^{\prime})}\\\
&+\sum_{i=1,2}\left[V^{R}_{11,i}\,e^{-i\lambda_{i}(t-t^{\prime})}\,\left[V^{R}\right]^{-1}_{i,11}\,\Pi^{>}_{12,12}(t^{\prime}-t)\right.\\\
&\quad\left.+V^{R}_{22,i}\,e^{-i\lambda_{i}(t-t^{\prime})}\,\left[V^{R}\right]^{-1}_{i,22}\,\Pi^{<}_{12,12}(t^{\prime}-t)\right]\bigg{)}\rho_{12}(t^{\prime})\end{split}$
(25)

Here, $\lambda_{i}$ are the eigenvalues of the Liouvillian matrix. We note
that non-Markov version of the Bloch QME, Eq.(25), accounts for broadening of
system states induced by their hybridization with the bath which is completely
missed by the standard Bloch QME, Eq.(13). At the same time, this result is
still an approximation (it is only second order in infinite hybridization
expansion). That is, while for moderate coupling strengths Eq.(25) can produce
relatively accurate results, for significant system-bath coupling strengths
the approximation may fail.

Below we use the Bloch equation (13) and its generalizations (II.4) and (25)
to discuss the concept of exceptional points for a Liouville operator.
Following Ref. Am-Shallem  _et al._ , 2015, we evaluate the time dependence of
the $z$-projection of the spin operator $S_{z}(t)$ and use it in eigeinmode
analysis

$S_{z}(t)=\sum_{k}d_{k}\,e^{-i\omega_{k}t}$ (26)

Degeneracies of the complex eigenmodes $\omega_{k}$ represent LEPs. As
discussed in Ref. Fuchs  _et al._ , 2014, the latter can be approximately
found from the points of divergence of the absolute values of the coefficients
$\lvert d_{k}\rvert$, although extended analysis is needed for further
characterization. We will use the parameters found for LEPs in Ref. Am-Shallem
_et al._ , 2015 as a starting point for our consideration.

Figure 2:  Time dependence of the $z$-component of spin operator, $S_{z}(t)$.
Panel (a) presents the results of simulations employing the standard Bloch
QME, Eq.(19) (dotted line, red); Bloch QME with dissipation included, Eq.
(II.4) (dashed line, black); and non-Markov Bloch QME with dissipation
included, Eq.(25) (dash-dotted line, blue). Panel (b) shows the results of the
NEGF Hartree-Fock simulations, Eqs. (5)-(II.2).

##  III Numerical results

We now evaluate $S_{z}(t)$ within different methodologies and use the results
of simulations to obtain exceptional points for the Liouville operator.

Unless stated otherwise, the parameters of the simulations are the following.
The energy levels of the system are $\varepsilon_{1}=0$ and
$\varepsilon_{2}=1$, the laser detuning
$\Delta\equiv\omega_{0}-\omega_{21}=0.00102$, and the coupling to driving
field $\mu\,E_{0}=0.001$. For simplicity we take
$V_{\alpha,11}=V_{\alpha,22}=0$, so that the dephasing rates are
$1/T_{1}=2/T_{2}=0.1$. The temperature of the bath is assumed to be zero.
Simulations were performed on a time grid of $200$ points with step $1$. We
confirmed that simulations on a grid of $2000$ points with step $0.1$ yield
similar results.

For non-Markov simulations we employ the bath spectral function

$J(\omega)\equiv\left(\frac{\omega}{\omega_{21}}\right)^{2}\,\exp\left[2\left(1-\frac{\omega}{\omega_{21}}\right)\right]$
(27)

and the bath dephasing rate is defined as

$\Gamma_{12,12}(\omega)\equiv\Gamma_{12,12}(\omega)\,J(\omega)=J(\omega)/T_{1}$
(28)

Fast Fourier transform is performed on a grid of 10001 points and utilizes the
FFTW library Frigo and Johnson (2005). The NEGF simulations are performed by
employing the procedure first introduced in Ref. Stan  _et al._ (2009).

Figure 2 shows time dependence of the $z$-projection of the spin operator
after employing the various approaches. Note that the differences in shapes of
the curves $S_{z}(t)$ reflect differences in underlying eigenmode
compositions. Note that the differences in long time value of the projection
are due to renormalization of the dissipation parameters (II.4) and thus are
of secondary importance. Comparing QME and NEGF results (panels (a) and (b),
respectively) we note oscillating behavior of the NEGF stationary state and
difference in magnitude of the signal. The reason for the discrepancy are
assumptions made when deriving the Bloch QME and its generalizations: 1. the
rotating wave approximation in external driving and 2\. neglecting effect of
the driving term on dissipator super-operator (i.e. dissipator is derived as
if there is no driving). Within the NEGF, the driving term is taken into
account exactly.

We now turn to the exceptional points analysis of the resulting time series.
Because we use the parameters of Ref. Am-Shallem  _et al._ , 2015 we know that
our standard Bloch QME simulations are performed in vicinity of LEP of second
order. Thus, the divergence of coefficient $\lvert d_{k}\rvert$ indicates the
presence of an exceptional point. Instead of the harmonic inversion analysis
employed in Ref. Am-Shallem  _et al._ (2015) for eigenmode decomposition, we
use its filter diagonalization variant Mandelshtam and Taylor (1997). For the
parameters of the simulations, the latter method appears to be more stable.

Figure 3 presents the eigenmode analysis for the time series $S_{z}(t)$
obtained within different Bloch QME schemes. The divergence of the expansion
coefficient $\lvert d_{1}\rvert$ and the disappearance of the eigenmode
difference in the analysis of the standard Bloch QME results presented in
panel (a) indicate the presence of a second order LEP at $\mu\,E_{0}\sim
0.025$. Panel (b) shows similar analysis for generalized Bloch QME with
included dissipation. Similar to the standard Bloch QME three eigenmodes are
present in the region away from the LEP. QME rates renormalization, Eq.(II.4),
leads to a shift of the position of the LEP, which now occurs at
$\mu\,E_{0}\sim 0.015$. Result of analysis for the non-Markov QME is shown in
panel (c). We find four different eigenmodes in this case. Careful comparison
with Markov consideration of panel (b) shows absence of exceptional points:
one can see that difference in eigenmodes does not disappear.

Figure 3:  Eigenmode analysis for the $z$-component of spin operator
$S_{z}(t)$, Eq.(26). Shown are absolute values of the coefficient of expansion
$\lvert d_{1}\rvert$ (blue line, left axis) and difference between eigenmodes
$\lvert\omega_{1}-\omega_{2}\rvert$ (red line, right axis) for results of
simulations performed within (a) standard Bloch QME, (b) generalized Bloch QME
with dissipation included, and (c) non-Markov generalized Bloch QME with
dissipation included.

Figure 4 shows two eigenmodes which become degenerate at exceptional point.
One sees that for the standard (panel a) and the generalized (panel b) Markov
QME weak coupling to the driving field (below the LEP) corresponds to
situation where real parts of the eigenmodes coincide while imaginary parts
are different. Stronger couplings (above the LEP) correspond to zero
difference in imaginary parts and different real parts. Note that similar
behavior at LEP yields transition between diffusive and ballistic motion
Hashimoto  _et al._ (2016a, b) and enhancement of decoherence rate Chen  _et
al._ (2021); Khandelwal  _et al._ (2021); Larson and Qvarfort (2023). Behavior
of eigenmodes for results obtained within non-Markov QME (panel c) is more
complicated. No degeneracy is observed between the modes. Similarly, eigenmode
analysis for the NEGF results yields a large number of modes ($\sim 50$) with
no LEPs present.

The absence of exceptional points in the results of non-Markov evolution is
expected because the EOM for $S_{z}(t)$ is not generated by the time-
independent Liouvillian anymore. One can understand the absence of the LEPs in
this case from a purely mathematical perspective. Indeed, even if one starts
from a time-dependent characteristic for a LEP (for example, for LEP2 one
expects to have $\rho(t)\sim(d_{1}+d_{2}t)e^{-i\lambda t}$) the first step of
time evolution will annihilate the LEP time dependence due to convolution of
the density operator with the time-dependent function, Eq.(25). Indeed, taking
the integral in Eq.(25) with memory kernel which depends on time in a
complicated way does not preserve original form of $\rho(t)$. This can be
easily seen by expanding the kernel in Fourier series and performing time
integration.

Failure of the concept of the Liouvillian exceptional point for non-Markov
evolution is even more obvious when analyzing the more rigorous NEGF
formulation, Eqs. (6)-(8). Indeed, in the right side of the Dyson equation one
has product of two Green’s functions: one from Eq.(6), the other from
Eq.(II.2). In principle, one could start from Eqs. (6)-(8) and apply the
generalized Kadanoff-Baym ansatz to these expressions. This would yield an
analog of QME which differs from (but is more accurate than) the Bloch QME.
Such equation would contain $\rho^{2}$ in its right side which obviously
indicates that the form $\rho(t)=(d_{1}+d_{2}t)e^{-i\lambda t}$ does not
survive in non-Markov formulation.

We note that the central parameters for the accuracy of Markov approximation
are the characteristic times of bath $t_{c}^{B}$ and system $t_{c}^{S}$
dynamics: Markov approximation is accurate when $t_{c}^{B}\ll t_{c}^{S}$. For
the model, $t_{c}^{B}$ is defined by the bandwidth $W_{B}$, temperature $T$,
and structure of the bath spectral function $J(\omega)$:
$t_{c}^{B}\sim\hbar/W,\hbar/k_{B}T$. The characteristic time of system
dynamics is defined by the intra-system energy parameters (inter-level
separation $\omega_{21}$, driving frequency $\omega_{0}$, detuning $\Delta$)
and by dissipation rate due to coupling to the bath (e.g., $\Gamma_{12,12}$):
$t_{c}^{S}\sim
2\pi/\omega_{21},2\pi/\omega_{0},2\pi/\Delta,\hbar/\Gamma_{12,12}$.

Figure 4:  Eigenmodes $\omega_{1}$ and $\omega_{2}$ vs. coupling to driving
field $\mu\,E_{0}$. Shown are real (blue solid and blue dashed lines) and
imaginary (red dashed and red dash-dotted lines) parts of the eigenmodes for
$S_{z}(t)$ simulations performed within (a) standard Bloch QME, (b)
generalized Bloch QME with dissipation included, and (c) non-Markov
generalized Bloch QME with dissipation included.

##  IV Conclusion

We discuss the concept of Liouvillian exceptional points (LEPs) used in the
description of dynamics of open quantum systems. The discussion is focused on
a model of driven two-level system coupled to a thermal bath. Starting with
exact NEGF formulation of the problem and implementing set of approximations
we derive standard Bloch QME and its generalizations. The latter include
dissipation (retarded self-energy contribution). One of the generalizations is
non-Markov.

We compare this approach with our recent publication Mukamel  _et al._ (2023)
where similar analysis for the Hamiltonian exceptional points (HEPs) was
carried out. We note that both HEP and LEP approximations rely on Markov
description of the system evolution. In terms of neglected self-energies,
standard HEP and LEP considerations are complementary: while HEP disregards
lesser and greater projections of self-energy, standard LEP misses its
retarded projection (dissipation).

By performing simulations for parameters previously shown to provide
exceptional points Am-Shallem  _et al._ (2015) we find that generalized Bloch
QME which includes information about dissipation and treats evolution as
Markov process is capable to provide LEPs although for adjusted parameters.
The non-Markov character of evolution does not permit introduction of the
concept of LEPs. In particular, neither the non-Markov Bloch QME formulation
nor the NEGF formulation is capable of producing the LEPs. This inability of
using LEPs for description of non-Markov evolution is quite general. The
concept of the Liouvillian exceptional points can be introduced only for
Markovian dynamics.

We note that while the RWA should be used in derivation of the Bloch QME,
within the NEGF treatment the approximation may be relaxed. Such more general
consideration will not affect the conclusions. Indeed, the inability to
introduce LEPs directly follows from the fact that the time-dependent
characteristic expected for a LEP does not survive non-Markov time evolution
of the system. Absence of the RWA will only change a form of a time-dependent
function convolution with which will destroy the expected LEP time dependence.
Similarly, as long as system evolution is non-Markov the conclusions hold for
any driving frequency or in absence of external driving.

Finally, we stress that our work does not challenge existing experimental
observations, some of which are mentioned in introduction. We discuss
theoretical treatments used for explanation of those experiments, and indicate
possible pitfalls of the theory. For example, in many cases, theoretical
treatments utilizing Markov description and employing exceptional points
analysis will predict an abrupt ‘phase transition’ when crossing the
exceptional point. In reality (i.e. within a more accurate theoretical
analysis), the transition between two different regimes will be smooth. The
importance of the difference between the two (approximate and more accurate)
theoretical descriptions and whether the approximate (Markov) treatment may
lead to qualitative failures depends on the observable of interest.

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation
under Grant No. 2154323.

##  Appendix A Derivation of Eq.(II.2)

We start from the definition of the single-particle Green’s function, Eq.(4).
Taking derivative in the first contour variable yields

$\begin{split}&i\frac{\partial}{\partial\tau_{1}}G_{ij}(\tau_{1},\tau_{2})=\delta_{i,j}\,\delta(\tau_{1},\tau_{2})\\\
&+\sum_{n=1,2}H^{S}_{in}(t_{1})\,G_{nj}(\tau_{1},\tau_{2})\\\
&-i\sum_{n,\alpha}\bigg{[}V_{ni,\alpha}\,\langle
T_{c}\,\hat{d}_{n}(\tau_{1})\hat{b}_{\alpha}(\tau_{1})\,\hat{d}_{j}^{\dagger}(\tau_{2})\rangle\\\
&\qquad\,\,\,+\langle
T_{c}\,\hat{b}_{\alpha}^{\dagger}(\tau_{1})\hat{d}_{n}(\tau_{1})\,\hat{d}_{j}^{\dagger}(\tau_{2})\rangle\,V_{\alpha,in}\bigg{]}\end{split}$
(29)

First order expansion of the scattering operator in the rightmost term of the
expression leads to

$\begin{split}&i\frac{\partial}{\partial\tau_{1}}G_{ij}(\tau_{1},\tau_{2})=\delta_{i,j}\,\delta(\tau_{1},\tau_{2})\\\
&+\sum_{n=1,2}H^{S}_{in}(t_{1})\,G_{nj}(\tau_{1},\tau_{2})-i\sum_{n,n_{1},n_{2}}\int_{c}d\tau\\\
&\quad\bigg{[}\,\,\Pi_{ni,n_{1}n_{2}}(\tau_{1},\tau)\,\langle
T_{c}\,\hat{d}_{n}(\tau_{1})\,\hat{d}^{\dagger}_{n_{1}}(\tau)\hat{d}_{n_{2}}(\tau)\,\hat{d}_{j}^{\dagger}(\tau_{2})\rangle_{0}\\\
&\quad+\Pi_{n_{1}n_{2},in}(\tau,\tau_{1})\,\langle
T_{c}\,\hat{d}_{n}(\tau_{1})\,\hat{d}^{\dagger}_{n_{2}}(\tau)\hat{d}_{n_{1}}(\tau)\,\hat{d}_{j}^{\dagger}(\tau_{2})\rangle_{0}\bigg{]}\end{split}$
(30)

Here, $\Pi$ is defined in (8) and subscript $0$ indicates evolution driven by
system Hamiltonian. Employing Wick’s theorem to decouple multi-time
correlation functions in the last term on the right side and dressing the
result yields the Hartree-Fock approximation, Eq.(II.2).

##  Appendix B Derivation of Eq.(10)

Here we derive exact EOM for density matrix, Eq. (10), starting from EOM for
the Green’s function (4).

We start with writing the left and right EOMs for the lesser projection of the
Green’s function (4)

$\displaystyle i\frac{\partial}{\partial t_{1}}G^{<}_{ij}(t_{1},t_{2})$
$\displaystyle=\varepsilon_{i}G^{<}_{ij}(t_{1},t_{2})-\mu
E(t_{1})G_{\bar{i}j}^{<}(t_{1},t_{2})$ (31)
$\displaystyle+i\sum_{n,\alpha}\left[V_{ni,\alpha}\langle\hat{d}_{j}^{\dagger}(t_{2})\hat{d}_{n}(t_{1})\hat{b}_{\alpha}(t_{1})\rangle+V_{\alpha,in}\langle\hat{d}_{j}^{\dagger}(t_{2})\hat{d}_{n}(t_{1})\hat{b}^{\dagger}_{\alpha}(t_{1})\rangle\right]$
$\displaystyle-i\frac{\partial}{\partial t_{2}}G^{<}_{ij}(t_{1},t_{2})$
$\displaystyle=\varepsilon_{j}G^{<}_{ij}(t_{1},t_{2})-\mu
E(t_{2})G^{<}_{i\bar{j}}(t_{1},t_{2})$ (32)
$\displaystyle+i\sum_{n,\alpha}\left[V_{jn,\alpha}\langle\hat{b}_{\alpha}(t_{2})\hat{d}_{n}^{\dagger}(t_{2})\hat{d}_{i}(t_{1})\rangle+V_{\alpha,nj}\langle\hat{b}_{\alpha}^{\dagger}(t_{2})\hat{d}_{n}^{\dagger}(t_{2})\hat{d}_{i}(t_{1})\rangle\right]$

Taking $t_{1}=t_{2}\equiv t$ and subtracting (31) from (32) yields

$\displaystyle-i\frac{d}{dt}G_{ij}^{<}(t,t)=\omega_{ji}G^{<}_{ij}(t,t)-\mu
E(t)\left[G^{<}_{i\bar{j}}(t,t)-G^{<}_{\bar{i}j}(t,t)\right]$
$\displaystyle\qquad+\sum_{n,\alpha}\left[V_{ni,\alpha}\,{G}^{<}_{\alpha,nj}(t,t)+{G}^{>}_{jn,\alpha}(t,t)\,V_{\alpha,in}\right.$
(33)
$\displaystyle\left.\qquad\qquad\,\,-V_{jn,\alpha}\,{G}^{>}_{\alpha,in}(t,t)-{G}_{ni,\alpha}^{<}(t,t)\,V_{\alpha,nj}\right]$

Here,

$\begin{split}{G}_{\alpha,m_{1}m_{2}}(\tau_{1},\tau_{2})&\equiv-i\bigg{\langle}T_{c}\,\hat{b}_{\alpha}(\tau_{1})\left[\hat{d}_{m_{1}}^{\dagger}\hat{d}_{m_{2}}\right]^{\dagger}(\tau_{2})\bigg{\rangle}\\\
{G}_{m_{1}m_{2},\alpha}(\tau_{1},\tau_{2})&\equiv-i\bigg{\langle}T_{c}\,\left[\hat{d}_{m_{1}}^{\dagger}\hat{d}_{m_{2}}\right](\tau_{1})\,\,\hat{b}_{\alpha}^{\dagger}(\tau_{2})\bigg{\rangle}\end{split}$
(34)

are the mixed system-bath Green’s function which satisfy the Dyson equations

$\begin{split}&{G}_{\alpha,m_{1}m_{2}}(\tau_{1},\tau_{2})=\\\
&\,\,\,\,\,\sum_{n_{1},n_{2}}\int_{c}d\tau^{\prime}\,F^{(0)}_{\alpha}(\tau_{1},\tau^{\prime})\,V_{\alpha,n_{1}n_{2}}\,G^{(2)}_{n_{1}n_{2},m_{1}m_{2}}(\tau^{\prime},\tau_{2})\\\
&{G}_{m_{1}m_{2},\alpha}(\tau_{1},\tau_{2})=\\\
&\,\,\,\,\,\sum_{n_{1},n_{2}}\int_{c}d\tau^{\prime}\,G^{(2)}_{m_{1}m_{2},n_{1}n_{2}}(\tau_{1},\tau^{\prime})\,V_{n_{1}n_{2},\alpha}\,F^{(0)}_{\alpha}(\tau^{\prime},\tau_{2})\end{split}$
(35)

Green’s functions $F^{(0)}_{\alpha}$ and $G^{(2)}$ are defined in Eqs. (9) and
(11), respectively.

Taking lesser and greater projections of (35), setting $t_{1}=t_{2}\equiv t$,
and substituting resulting expressions into (B) leads to (10).

##  Appendix C Derivation of Eq.(13)

Here we derive Bloch equations, Eq. (13), starting from exact EOM for density
matrix, Eq. (10).

Substituting the Kadanoff-Baym ansatz Eq. (12) into (10) gives

$\displaystyle\frac{d}{dt}\rho_{ij}(t)\approx i\omega_{ji}\rho_{ij}(t)-i\mu
E(t)\left[\rho_{i\bar{j}}(t)-\rho_{\bar{i}j}(t)\right]$ (36)
$\displaystyle-i\sum_{n,n_{1},n_{2}}\int_{-\infty}^{t}dt^{\prime}\,\bigg{[}\Pi_{ni,n_{1}n_{2}}^{>}(t-t^{\prime})\,e^{i\omega_{jn}(t-t^{\prime})}\,\left\langle\left[\hat{d}_{j}^{\dagger}\hat{d}_{n}\hat{d}_{n_{1}}^{\dagger}\hat{d}_{n_{2}}\right](t^{\prime})\right\rangle-\Pi_{ni,n_{1}n_{2}}^{<}(t-t^{\prime})\,e^{i\omega_{jn}(t-t^{\prime})}\,\left\langle\left[\hat{d}_{n_{1}}^{\dagger}\hat{d}_{n_{2}}\hat{d}_{j}^{\dagger}\hat{d}_{n}\right](t^{\prime})\right\rangle$
$\displaystyle\qquad\qquad\qquad\quad\,\,\,-\Pi_{jn,n_{1}n_{2}}^{>}(t-t^{\prime})\,e^{i\omega_{ni}(t-t^{\prime})}\left\langle\left[\hat{d}_{n}^{\dagger}\hat{d}_{i}\hat{d}_{n_{1}}^{\dagger}\hat{d}_{n_{2}}\right](t^{\prime})\right\rangle+\Pi_{jn,n_{1}n_{2}}^{<}(t-t^{\prime})\,e^{i\omega_{ni}(t-t^{\prime})}\left\langle\left[\hat{d}_{n_{1}}^{\dagger}\hat{d}_{n_{2}}\hat{d}_{n}^{\dagger}\hat{d}_{i}\right](t^{\prime})\right\rangle$
$\displaystyle\qquad\qquad\qquad\quad\,\,\,+\Pi_{n_{1}n_{2},nj}^{>}(t^{\prime}-t)\,e^{i\omega_{in}(t^{\prime}-t)}\,\left\langle\left[\hat{d}_{n_{2}}^{\dagger}\hat{d}_{n_{1}}\hat{d}_{n}^{\dagger}\hat{d}_{i}\right](t^{\prime})\right\rangle-\Pi_{n_{1}n_{2},nj}^{<}(t^{\prime}-t)\,e^{i\omega_{in}(t^{\prime}-t)}\,\left\langle\left[\hat{d}_{n}^{\dagger}\hat{d}_{i}\hat{d}_{n_{2}}^{\dagger}\hat{d}_{n_{1}}\right](t^{\prime})\right\rangle$
$\displaystyle\qquad\qquad\qquad\quad\,\,\,-\Pi_{n_{1}n_{2},in}^{>}(t^{\prime}-t)\,e^{i\omega_{nj}(t^{\prime}-t)}\,\left\langle\left[\hat{d}_{n_{2}}^{\dagger}\hat{d}_{n_{1}}\hat{d}_{j}^{\dagger}\hat{d}_{n}\right](t^{\prime})\right\rangle+\Pi_{n_{1}n_{2},in}^{<}(t^{\prime}-t)\,e^{i\omega_{nj}(t^{\prime}-t)}\,\left\langle\left[\hat{d}_{j}^{\dagger}\hat{d}_{n}\hat{d}_{n_{2}}^{\dagger}\hat{d}_{n_{1}}\right](t^{\prime})\right\rangle\bigg{]}$

where self-energy $\Pi$ is defined in Eq. (8).

In the single electron subspace of the problem

$\begin{split}&\langle\hat{d}_{1}^{\dagger}\hat{d}_{1}\hat{d}_{1}^{\dagger}\hat{d}_{1}\rangle=\langle\hat{d}_{1}^{\dagger}\hat{d}_{2}\hat{d}_{2}^{\dagger}\hat{d}_{1}\rangle=\langle\hat{d}_{1}^{\dagger}\hat{d}_{1}\rangle\equiv\rho_{11},\\\
&\langle\hat{d}_{2}^{\dagger}\hat{d}_{1}\hat{d}_{1}^{\dagger}\hat{d}_{1}\rangle=\langle\hat{d}_{2}^{\dagger}\hat{d}_{2}\hat{d}_{2}^{\dagger}\hat{d}_{1}\rangle=\langle\hat{d}_{2}^{\dagger}\hat{d}_{1}\rangle\equiv\rho_{12}\\\
&\langle\hat{d}_{1}^{\dagger}\hat{d}_{1}\hat{d}_{1}^{\dagger}\hat{d}_{2}\rangle=\langle\hat{d}_{1}^{\dagger}\hat{d}_{2}\hat{d}_{2}^{\dagger}\hat{d}_{2}\rangle=\langle\hat{d}_{1}^{\dagger}\hat{d}_{2}\rangle\equiv\rho_{21}\\\
&\langle\hat{d}_{2}^{\dagger}\hat{d}_{2}\hat{d}_{2}^{\dagger}\hat{d}_{2}\rangle=\langle\hat{d}_{2}^{\dagger}\hat{d}_{1}\hat{d}_{1}^{\dagger}\hat{d}_{2}\rangle=\langle\hat{d}_{2}^{\dagger}\hat{d}_{2}\rangle\equiv\rho_{22}\end{split}$
(37)

with all other averages zero.

Substituting (37) in (36), employing fast bath approximation

$\begin{split}&\rho_{11}(t^{\prime})\approx\rho_{11}(t),\quad\rho_{12}(t^{\prime})\approx\rho_{12}(t)\,e^{+i\omega_{21}(t^{\prime}-t)},\\\
&\rho_{22}(t^{\prime})\approx\rho_{22}(t),\quad\rho_{21}(t^{\prime})\approx\rho_{21}(t)\,e^{-i\omega_{21}(t^{\prime}-t)},\end{split}$
(38)

and neglecting bath-induced couplings between populations and coherences leads
to (13).

##  Appendix D Derivation of Eq.(21)

Within the single-electron subspace of the problem there is a simple one-to-
one correspondence between the single-particle and many-body states of the
system. This correspondence allows one to express the lesser and greater
projections of the two-particle GF (11) as

$\begin{split}&G^{(2)\,<}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{2})\equiv-i\langle\left[\hat{d}_{n_{3}}^{\dagger}\hat{d}_{n_{4}}\right]^{\dagger}(t_{2})\,\left[\hat{d}_{n_{1}}^{\dagger}\hat{d}_{n_{2}}\right](t_{1})\rangle\\\
&\quad=-i\theta(t_{1}-t_{2})\llangle
n_{2}n_{1}\rvert\,\mathcal{U}(t_{1},t_{2})\,\lvert\,\rho_{SB}(t_{2})\,n_{4}n_{3}\rrangle\\\
&\quad\,\,\,\,-i\theta(t_{2}-t_{1})\llangle
n_{3}n_{4}\rvert\,\mathcal{U}(t_{2},t_{1})\,\lvert\,n_{1}n_{2}\,\rho_{SB}(t_{1})\rrangle\\\
&G^{(2)\,>}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{2})\equiv-i\langle\left[\hat{d}_{n_{1}}^{\dagger}\hat{d}_{n_{2}}\right](t_{1})\,\left[\hat{d}_{n_{3}}^{\dagger}\hat{d}_{n_{4}}\right]^{\dagger}(t_{2})\rangle\\\
&\quad=-i\theta(t_{1}-t_{2})\llangle
n_{2}n_{1}\rvert\,\mathcal{U}(t_{1},t_{2})\,\lvert\,n_{4}n_{3}\,\rho_{SB}(t_{2})\rrangle\\\
&\quad\,\,\,\,-i\theta(t_{2}-t_{1})\llangle
n_{3}n_{4}\rvert\,\mathcal{U}(t_{2},t_{1})\,\lvert\,\rho_{SB}(t_{1})\,n_{1}n_{2}\rrangle\end{split}$
(39)

where rightmost sides of the expressions are written in the Liouville space
notation $\lvert n_{1}n_{2}\rrangle\equiv\lvert n_{1}\rangle\langle
n_{2}\rvert$, $\hat{\rho}_{SB}$ is the total (system and bath) density
operator, and

$\mathcal{U}(t_{1},t_{2})\equiv
T\,\exp\left[-i\int_{t_{2}}^{t_{1}}dt\,\mathcal{L}(t)\right]$ (40)

is the Liouville space evolution operator.

By decoupling the system and the bath dynamics with projection super-operator

$\mathcal{P}\equiv\sum_{e,f}\lvert\,ef\,\rho_{B}^{eq}\,\rrangle\,\llangle\,ef\,I_{B}\rvert$
(41)

and introducing retarded and advanced Green’s functions in Liouville space

$\begin{split}&\mathcal{G}^{r}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{2})\equiv\\\
&\quad-i\theta(t_{1}-t_{2})\llangle
n_{2}n_{1}\,I_{B}\rvert\,\mathcal{U}(t_{1},t_{2})\,\lvert
n_{4}n_{3}\,\rho_{B}^{eq}\rrangle\\\
&\mathcal{G}^{a}_{n_{1}n_{2},n_{3}n_{4}}(t_{1},t_{2})\equiv\\\
&\quad+i\theta(t_{2}-t_{1})\llangle
n_{3}n_{4}\,I_{B}\rvert\,\mathcal{U}(t_{1},t_{2})\,\lvert
n_{1}n_{2}\,\rho_{B}^{eq}\rrangle\end{split}$ (42)

one can rewrite the exact expressions (39) in the form of generalized
Kadanoff-Baym ansatz in the Liouville space, Eq.(21).

Note that expressions for the Green’s functions in the Liouville space given
by Eq.(22), are equivalent to the their definition (42) via

$\mathcal{U}_{eff}(t_{1},t_{2})\equiv\llangle
I_{B}\rvert\,\mathcal{U}(t_{1},t_{2})\,\lvert\rho_{B}^{eq}\rrangle$ (43)

## References

  * Moiseyev (2011) N. Moiseyev,  _Non-Hermitian Quantum Mechanics_ (Cambridge University Press, Cambridge, 2011).
  * Miri and Alù (2019) M.-A. Miri and A. Alù, Science 363, eaar7709 (2019).
  * Lee  _et al._ (2009) S.-B. Lee, J. Yang, S. Moon, S.-Y. Lee, J.-B. Shim, S. W. Kim, J.-H. Lee, and K. An, Phys. Rev. Lett. 103, 134101 (2009).
  * Rüter  _et al._ (2010) C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Nature Phys. 6, 192 (2010).
  * Hahn  _et al._ (2016) C. Hahn, Y. Choi, J. Woong Yoon, S. Ho Song, C. Hwan Oh, and P. Berini, Nat. Commun. 7, 12201 (2016).
  * Doppler  _et al._ (2016) J. Doppler, A. A. Mailybaev, J. Böhm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, Nature 537, 76 (2016).
  * Ergoktas  _et al._ (2022) M. S. Ergoktas, S. Soleymani, N. Kakenov, K. Wang, T. B. Smith, G. Bakan, S. Balci, A. Principi, K. S. Novoselov, S. K. Ozdemir, and C. Kocabas, Science 376, 184 (2022).
  * Wu  _et al._ (2019) Y. Wu, W. Liu, J. Geng, X. Song, X. Ye, C.-K. Duan, X. Rong, and J. Du, Science 364, 878 (2019).
  * Wiersig (2020) J. Wiersig, Photon. Res. 8, 1457 (2020).
  * Yang  _et al._ (2023) M. Yang, H.-Q. Zhang, Y.-W. Liao, Z.-H. Liu, Z.-W. Zhou, X.-X. Zhou, J.-S. Xu, Y.-J. Han, C.-F. Li, and G.-C. Guo, Science Advances 9, eabp8943 (2023).
  * Liang  _et al._ (2023) C. Liang, Y. Tang, A.-N. Xu, and Y.-C. Liu, Phys. Rev. Lett. 130, 263601 (2023).
  * Zhang  _et al._ (2018) J. Zhang, B. Peng, Ş. K. Özdemir, K. Pichler, D. O. Krimer, G. Zhao, F. Nori, Y.-x. Liu, S. Rotter, and L. Yang, Nat. Photonics 12, 479 (2018).
  * Naghiloo  _et al._ (2019) M. Naghiloo, M. Abbasi, Y. N. Joglekar, and K. W. Murch, Nature Phys. 15, 1232 (2019).
  * Gao  _et al._ (2015) T. Gao, E. Estrecho, K. Y. Bliokh, T. C. H. Liew, M. D. Fraser, M. Brodbeck, M. Kamp, C. Schneider, S. Höfling, Y. Yamamoto, F. Nori, Y. S. Kivshar, A. G. Truscott, D. R. G., and N. A. Ostrovskaya, Nature 526, 554 (2015).
  * Xu  _et al._ (2023) G. Xu, X. Zhou, Y. Li, Q. Cao, W. Chen, Y. Xiao, L. Yang, and C.-W. Qiu, Phys. Rev. Lett. 130, 266303 (2023).
  * Günther  _et al._ (2007) U. Günther, I. Rotter, and B. F. Samsonov, J. Phys. A 40, 8815 (2007).
  * Rotter (2009) I. Rotter, J. Phys. A 42, 153001 (2009).
  * Toroker and Peskin (2009) M. C. Toroker and U. Peskin, J. Phys. B 42, 044013 (2009).
  * Uzdin  _et al._ (2011) R. Uzdin, A. Mailybaev, and N. Moiseyev, J. Phys. A: Math. and Theor. 44, 435302 (2011).
  * Heiss (2012) W. D. Heiss, J. Phys. A 45, 444016 (2012).
  * Garmon  _et al._ (2012) S. Garmon, I. Rotter, N. Hatano, and D. Segal, Int. J. Theor. Phys. 51, 3536 (2012).
  * Delga  _et al._ (2014) A. Delga, J. Feist, J. Bravo-Abad, and F. J. Garcia-Vidal, Journal of Optics 16, 114018 (2014).
  * Rotter and Bird (2015) I. Rotter and J. P. Bird, Rep. Prog. Phys. 78, 114001 (2015).
  * Yang  _et al._ (2020) C. Yang, X. Wei, J. Sheng, and H. Wu, Nat. Commun. 11, 4656 (2020).
  * Engelhardt and Cao (2022) G. Engelhardt and J. Cao, Phys. Rev. B 105, 064205 (2022).
  * Ferrier  _et al._ (2022) L. Ferrier, P. Bouteyre, A. Pick, S. Cueff, N. Dang, C. Diederichs, A. Belarouci, T. Benyattou, J. Zhao, R. Su, J. Xing, Q. Xiong, and H. Nguyen, Phys. Rev. Lett. 129, 083602 (2022).
  * Engelhardt and Cao (2023) G. Engelhardt and J. Cao, Phys. Rev. Lett. 130, 213602 (2023).
  * Li  _et al._ (2023) Z.-Z. Li, W. Chen, M. Abbasi, K. W. Murch, and K. B. Whaley, Phys. Rev. Lett. 131, 100202 (2023).
  * Pick  _et al._ (2019) A. Pick, S. Silberstein, N. Moiseyev, and N. Bar-Gill, Phys. Rev. Res. 1, 013015 (2019).
  * Minganti  _et al._ (2019) F. Minganti, A. Miranowicz, R. W. Chhajlany, and F. Nori, Phys. Rev. A 100, 062131 (2019).
  * Minganti  _et al._ (2020) F. Minganti, A. Miranowicz, R. W. Chhajlany, I. I. Arkhipov, and F. Nori, Phys. Rev. A 101, 062112 (2020).
  * Arkhipov  _et al._ (2020) I. I. Arkhipov, A. Miranowicz, F. Minganti, and F. Nori, Phys. Rev. A 102, 033715 (2020).
  * Chimzak  _et al._ (2023) G. Chimzak, A. Kowalewska-Kudlaszyk, E. Lange, K. Bartkiewicz, and J. Peiřina Jr., Sci. Rep. 13, 5859 (2023).
  * Hashimoto  _et al._ (2016a) K. Hashimoto, K. Kanki, S. Garmon, S. Tanaka, and T. Petrosky, Progress of Theoretical and Experimental Physics 2016, 053A02 (2016a).
  * Hashimoto  _et al._ (2016b) K. Hashimoto, K. Kanki, S. Tanaka, and T. Petrosky, in  _Non-Hermitian Hamiltonians in Quantum Physics_ , edited by F. Bagarello, R. Passante, and C. Trapani (Springer International Publishing, Cham, 2016) pp. 263–279.
  * Chen  _et al._ (2021) W. Chen, M. Abbasi, Y. N. Joglekar, and K. W. Murch, Phys. Rev. Lett. 127, 140504 (2021).
  * Khandelwal  _et al._ (2021) S. Khandelwal, N. Brunner, and G. Haack, PRX Quantum 2, 040346 (2021).
  * Larson and Qvarfort (2023) J. Larson and S. Qvarfort, Open Systems & Information Dynamics 30, 2350008 (2023).
  * Chen  _et al._ (2022) W. Chen, M. Abbasi, B. Ha, S. Erdamar, Y. N. Joglekar, and K. W. Murch, Phys. Rev. Lett. 128, 110402 (2022).
  * Kumar  _et al._ (2022) P. Kumar, K. Snizhko, Y. Gefen, and B. Rosenow, Phys. Rev. A 105, L010203 (2022).
  * Bu  _et al._ (2023) J.-T. Bu, J.-Q. Zhang, G.-Y. Ding, J.-C. Li, J.-W. Zhang, B. Wang, W.-Q. Ding, W.-F. Yuan, L. Chen, i. m. c. K. Özdemir, F. Zhou, H. Jing, and M. Feng, Phys. Rev. Lett. 130, 110402 (2023).
  * Mukamel  _et al._ (2023) S. Mukamel, A. Li, and M. Galperin, J. Chem. Phys. 158, 154106 (2023), publisher: American Institute of Physics.
  * Am-Shallem  _et al._ (2015) M. Am-Shallem, R. Kosloff, and N. Moiseyev, New J. Phys. 17, 113036 (2015).
  * Hatano (2019) N. Hatano, Mol. Phys. 117, 2121 (2019).
  * Perina Jr  _et al._ (2022) J. Perina Jr, A. Miranowicz, G. Chimczak, and A. Kowalewska-Kudlaszyk, Quantum 6, 883 (2022).
  * Tay (2023) B. A. Tay, Physica A 620, 128736 (2023).
  * Haug and Jauho (2008) H. Haug and A.-P. Jauho,  _Quantum Kinetics in Transport and Optics of Semiconductors_ , second, substantially revised edition ed. (Springer, Berlin Heidelberg, 2008).
  * Lipavský  _et al._ (1986) P. Lipavský, V. Špička, and B. Velický, Phys. Rev. B 34, 6933 (1986).
  * Esposito and Galperin (2009) M. Esposito and M. Galperin, Phys. Rev. B 79, 205303 (2009).
  * Esposito and Galperin (2010) M. Esposito and M. Galperin, J. Phys. Chem. C 114, 20362 (2010).
  * Fuchs  _et al._ (2014) J. Fuchs, J. Main, H. Cartarius, and G. Wunner, J. Phys. A: Math. and Theor. 47, 125304 (2014).
  * Frigo and Johnson (2005) M. Frigo and S. G. Johnson, Proc. IEEE 93, 216 (2005).
  * Stan  _et al._ (2009) A. Stan, N. E. Dahlen, and R. van Leeuwen, J. Chem. Phys. 130, 224101 (2009).
  * Mandelshtam and Taylor (1997) V. A. Mandelshtam and H. S. Taylor, J. Chem. Phys. 107, 6756 (1997).


#  Metrically differentiable set-valued functions
and their local linear approximants

Nira Dyn111Tel-Aviv University, Israel , Elza Farkhi11footnotemark: 1 , Alona
Mokhov222Afeka Tel-Aviv Academic College of Engineering, Israel

###### Abstract

A new notion of metric differentiability of set-valued functions at a point is
introduced in terms of right and left limits of special set-valued metric
divided differences of first order. A local metric linear approximant of a
metrically differentiable set-valued function at a point is defined and
studied. This local approximant may be regarded as a special realization of
the set-valued Euler approximants of M. S. Nikolskii and the directives of Z.
Artstein. Error estimates for the local metric linear approximant are
obtained. In particular, second order approximation is derived for a class of
“strongly” metrically differentiable set-valued maps.

Dedicated to Professor Dany Leviatan on the occasion of his 80th birthday

Key words: compact sets in Euclidean spaces, set-valued functions, metric
divided differences of first order, metric differentiability, local metric
linear approximants

##  1 Introduction

A differential calculus of set-valued functions (SVFs, set-valued maps,
multifunctions) is playing a significant role in different areas of applied
mathematics. Various notions of set-valued derivatives, motivated by different
needs and applications, have been introduced and studied in the last five
decades. The earlier papers focus on differentiabilty of multifunctions with
convex images (see e.g. [4, 5, 7]). Different concepts of set-valued
differentiability (such as Hukuhara-derivatives, contingent derivatives,
graphical derivatives, co-derivatives, generalized differentiability) are
introduced, discussed, and studied, for instance, in [3, 13, 16, 18, 19]. In
[12, 15, 21] the notions of affine, semi-affine, eclipsing, quasi-affine
mappings are suggested for the role of set-valued ”linear” approximants.
Derivatives of set-valued functions also arise in probability and statistics,
where another approach still might be required (see e.g. [14]).

In a series of works (e.g. [6, 8, 9, 10, 11]) we have developed a metric
approach to the approximation of SVFs, mapping a compact interval to general
(not necessarily convex) compact sets in ${\mathbb{R}}^{n}$. This approach,
inspired by the ideas from [1], includes metric tools such as metric chains,
metric linear combinations, metric integral. The metric tools are applied to
adapt classical sample-based approximation operators to multifunctions. In [6,
8, 9] we obtain approximation order at most $1$ for SVFs which are continuous
or of bounded variation. In [11] we define the notion of metric divided
differences of any order $r\in{\mathbb{N}}$ and obtain high order
approximation by metric piecewise-polynomial interpolation, without
introducing a notion of derivative for set-valued functions.

In this paper, we introduce the concept of metric differentiability for a
univariate multifunction $F$ with the aim of achieving a local approximation
of $F$ with an error of order higher than 1. We apply this concept in
constructing and analyzing a local metric linear approximant. Our approach
yields error estimates beyond first order (including second order) for these
approximants, extending classical error estimates for real-valued functions to
set-valued functions.

The notion of local metric linear approximant, introduced here may be regarded
as a special unique “realization” of the Euler approximants of M.S. Nikolskii
[17], as well as of the “directives” of Artstein [2]. They both have proposed
a definition of a set-valued derivative by a local “linear approximant” of the
given multifunction in a neghborhood of a given point with error of order
${o}(h)$ measured in the Hausdorff metric. They both have noticed that such a
definition does not provide uniqueness of the derivative.

We define here the notion of metric divided differences of first order
anchored at a point $y\in F(x_{0})$ and introduce one-sided (left and right)
metric derivatives anchored at $y$ as the one-sided limits (left and right) of
the corresponding metric divided differences. We construct the local metric
linear approximant of a metrically differentiable multifunction at a point
$x_{0}$ and obtain rate of approximation $o(h)$. To get a higher order,
$O(h^{m})$, $m>1$, we introduce a higher order smoothness property for SVFs
related to the rate of convergence of the metric divided differences to the
one-sided metric derivatives.

The paper is organized as follows: In the next section we present preliminary
material. In Section 3 we present metric divided differences of first order,
define the one-sided metric derivatives and investigate their properties. In
Section 4 we introduce and study the local metric linear approximant.

##  2 Preliminaries

We use the folowing notation. The Euclidean norm in ${\mathbb{R}}^{n}$ is
denoted by $|\cdot|$. In this paper we consider sets in ${\cal K}_{n}$, where
${\cal K}_{n}$ denotes the collection of nonempty compact subsets of
${\mathbb{R}}^{n}$. The distance from a point ${x\in{\mathbb{R}}^{n}}$ to a
set $A\subset{\mathbb{R}}^{n}$ is denoted by ${\mbox{dist}(x,A)=\min_{a\in
A}|x-a|}$.

The set of all projections of $a\in{\mathbb{R}}^{n}$ on a set $B\in{\cal
K}_{n}$ is

$\Pi_{B}{(a)}=\\{b\in B:|a-b|=\operatorname{dist}(a,B)\\}.$

We denote the collection of all metric pairs of two sets $A,B\in{\cal K}_{n}$
by

$\Pi({A},{B})=\\{(a,b)\in A\times
B:\;a\in\Pi_{A}{(b)}\;\,\mbox{or}\;\,b\in\Pi_{B}{(a)}\\}.$

######  Remark 2.1.

It follows from the compactness of $A,B\in{\cal K}_{n}$ that $\Pi({A},{B})$ is
a non-empty compact set. (see Lemma 2.4 in [11])

Using the metric pairs, the Hausdorff distance can written as

${\mbox{h}aus}(A,B)=\max\\{|a-b|:\;(a,b)\in\Pi({A},{B})\\}.$ (1)

The space ${\cal K}_{n}$ is a complete metric space with respect to this
metric [20]. In this paper we regard convergence of sets and continuity of
set-valued maps with respect to the Hausdorff metric.

We use in the paper the following property of the Hausdorff metric, which is
easy to check:

${\mbox{h}aus}(\bigcup\limits_{\lambda\in\Lambda}A_{\lambda},\bigcup\limits_{\lambda\in\Lambda}B_{\lambda})\leq\sup_{\lambda\in\Lambda}{\mbox{h}aus}(A_{\lambda},B_{\lambda}),\quad\Lambda\;\hbox{is
an index set}.$ (2)

For a given set $A$ we use the notation

$\lambda A=\left\\{\lambda a\ :\ a\in A\right\\},\quad\lambda\in{\mathbb{R}},$

and

$\|A\|=\max\left\\{|a|\ :\ a\in A\right\\}.$

A linear Minkowski combination of two sets $A,B$ is

$\lambda A+\mu B=\\{\lambda a+\mu b:a\in A,b\in
B\\},\quad\lambda,\mu\in{\mathbb{R}}.$

The following relations between operations on sets and the Hausdorff metric
are easy to verify.

${\mbox{h}aus}(\\{x\\}+A,\\{x\\}+B)={\mbox{h}aus}(A,B),\qquad{\mbox{h}aus}(\lambda
A,\lambda B)=|\lambda|{\mbox{h}aus}(A,B),$ (3)

where $\lambda\in{\mathbb{R}}$ and $\\{x\\}$ is the set consisting of the
single element $x$.

We recall here the notion of a metric chain and of a metric linear combination
[10].

Given a finite sequence of sets $A_{0},\ldots,A_{m}\in{\cal K}_{n}$, $m\geq
1$, the collection of all metric chains of ${A_{0},\ldots,A_{m}}$ is

${\mathrm{CH}}(A_{0},\ldots,A_{m}):=\left\\{(a_{0},\ldots,a_{m}):\;(a_{i},a_{i+1})\in\Pi({A_{i}},{A_{i+1}}),\
i=0,1,\ldots,m-1\right\\}.$ (4)

For $m=0$ we define ${\mathrm{CH}}(A_{0})=A_{0}$. Note that
${\mathrm{CH}}(A_{0},A_{1})=\Pi({A_{0}},{A_{1}})$.

The metric linear combination of the sets $A_{0},\ldots,A_{m}\in{\cal K}_{n}$,
is

$\bigoplus_{i=0}^{m}\lambda_{i}A_{i}:=\left\\{\sum_{i=0}^{m}\lambda_{i}a_{i}\,:(a_{0},\ldots,a_{m})\in{\mathrm{CH}}(A_{0},\ldots,A_{m})\right\\}.$
(5)

######  Remark 2.2.

  1. (1)

In the special case $m=1$ and $\lambda_{0}=1$, $\lambda_{1}=-1$ we obtain the
metric difference of two sets $A,B$ defined in [8],

$A\ominus B:=\\{\ a-b:(a,b)\in\Pi({A},{B})\ \\}.$

  2. (2)

In the special case $m=1$ and $\lambda_{0},\lambda_{1}\in[0,1]$,
$\lambda_{0}+\lambda_{1}=1$, the metric linear combination is the binary
operation between compact sets, introduced in [1].

  3. (3)

The set $\displaystyle\bigoplus_{i=0}^{m}\lambda_{i}A_{i}$ is compact (see
[11])).

In this work we consider univariate set-valued functions with values in ${\cal
K}_{n}$. The graph of a multifunction ${F:[a,b]\rightarrow{\cal K}_{n}}$ is
the set $Graph(F)=\\{(x,y):x\in[a,b],\,y\in F(x)\\}$.

##  3 First metric divided differences and one-sided metric derivative

We start with defining divided differences.

######  Definition 3.1.

Let $F:(a,b)\to{\cal K}_{n}$. Fix $x_{0},x\in(a,b)$, $x_{0}\neq x$ and
$y_{0}\in F(x_{0})$. Define the first metric divided difference of $F$ at
$x_{0},x$, anchored at $y_{0}$ , as the set of vectors

$[x_{0},x]^{M}F|_{y_{0}}:=\
\left\\{\frac{y-y_{0}}{x-x_{0}}:\,(y_{0},y)\in\Pi({F(x_{0})},{F(x)})\right\\}.$

######  Remark 3.2.

  1. (i)

By the compactness of $F(x),F(x_{0})$, and by Remark 2.1,
$[x_{0},x]^{M}F|_{y_{0}}$ is non-empty and compact.

  2. (ii)

If $(x_{0},y_{0})$ is in the interior of $Graph(F)$, then $y_{0}\in F(x)$ for
$x$ in a small enough neighborhood of $x_{0}$, thus
$(y_{0},y_{0})\in\Pi({F(x_{0})},{F(x)})$ is the only metric pair containing
$y_{0}$, and therefore $[x_{0},x]^{M}F|_{y_{0}}=\\{0\\}$.

  3. (iii)

By Definition 3.1, $F(x)=\bigcup_{y_{0}\in
F(x_{0})}\\{\\{y_{0}\\}+(x-x_{0})[x_{0},x]^{M}F|_{y_{0}}\\}$.

######  Example 3.3.

Let $F:(0,2)\to{\cal K}_{1}$ be defined by $F(x)=\\{x^{\alpha},x^{\beta}\\}$,
where $\alpha,\beta\in{\mathbb{N}}$, $\alpha\neq\beta$, and let
$x_{0}\in(0,2)$, $x=x_{0}+h$ with $|h|\neq 0$ sufficiently small. Then for
$x_{0},x_{0}+h\leq 1$ or $x_{0},x_{0}+h\geq 1$

$\Pi({F(x_{0})},{F(x_{0}+h)})=\left\\{\big{(}x_{0}^{\alpha},(x_{0}+h)^{\alpha}\big{)},\big{(}x_{0}^{\beta},(x_{0}+h)^{\beta}\big{)}\right\\}.$

For $x_{0}=1$ we get
$[1,1+h]^{M}F|_{1}=\left\\{\frac{(1+h)^{\alpha}-1}{h},\frac{(1+h)^{\beta}-1}{h}\right\\}$,
while for $x_{0}\neq 1$ we have

$[x_{0},x_{0}+h]^{M}F|_{x_{0}^{\alpha}}=\left\\{\frac{(x_{0}+h)^{\alpha}-x_{0}^{\alpha}}{h}\right\\}\quad,\quad[x_{0},x_{0}+h]^{M}F|_{x_{0}^{\beta}}=\left\\{\frac{(x_{0}+h)^{\beta}-x_{0}^{\beta}}{h}\right\\}.$

######  Definition 3.4.

Let $F:(a,b)\to{\cal K}_{n}$. Fix $x_{0},x\in(a,b)$, $x_{0}\neq x$. Define the
first metric divided difference of $F$ at $x_{0},x$ as

$[x_{0},x]^{M}F:=\bigcup\limits_{y\in F(x_{0})}[x_{0},x]^{M}F|_{y},$

######  Remark 3.5.

  1. (i)

The set $[x_{0},x]^{M}F$ is the metric divided difference of order $1$ in [11,
Remark 5.2]. It is compact as is shown in Lemma 5.5 there.

  2. (ii)

Using the metric difference of two sets (see Remark 2.2), the first metric
divided difference can be expressed, similarly to the divided difference of
real-valued functions, as

$[x_{0},x_{1}]^{M}F=\frac{1}{x_{1}-x_{0}}\left(F(x_{1})\ominus
F(x_{0})\right).$

It is clear that one can interchange the points $x_{0},x_{1}$ in this divided
difference.

  3. (iii)

In view of (1), ${\mbox{h}aus}(A,B)=\|A\ominus B\|$ and by the second property
in (3), a multifunction $F$ is Lipschitz with a constant $L$ in a domain if
and only if

$\|[x_{0},x_{1}]^{M}F\|=\frac{{\mbox{h}aus}(F(x_{1}),F(x_{0}))}{|x_{1}-x_{0}|}\leq
L$

for any $x_{0},x_{1}$ in this domain.

  4. (iv)

For a single-valued function, Definition 3.4 coincides with the classical
definition of divided difference of first order.

Now we introduce the right and the left metric derivatives. We use the
notations $\lim_{x\to x_{0}^{+}}$ and $\lim_{x\to x_{0}^{-}}$ for the right
and left limits, respectively.

######  Definition 3.6.

Let $F:(a,b)\to{\cal K}_{n},\ x_{0}\in(a,b)$.

  1. (i)

We say that $F$ is metrically differentiable from the right at $x_{0}$ if for
any $y\in F(x_{0})$ there is a non-empty set $D^{M}_{+}F(x_{0})|_{y}$ of
vectors satisfying

$D^{M}_{+}F(x_{0})|_{y}=\lim_{x\to x_{0}^{+}}[x_{0},x]^{M}F|_{y}$

uniformly in $y\in F(x_{0})$, namely, for any $\epsilon>0$ there exists
$\delta>0$ such that for any $x$ with $0<x-x_{0}<\delta$

$\sup_{y\in
F(x_{0})}{\mbox{h}aus}\left(D^{M}_{+}F(x_{0})|_{y},[x_{0},x]^{M}F|_{y}\right)<\epsilon.$

We call the set $D^{M}_{+}F(x_{0})|_{y}$ the right metric derivative of $F$ at
$x_{0}$ anchored at $y\in F(x_{0})$.
The set

${\mathfrak{D}}^{M}_{+}F=\left\\{D^{M}_{+}F(x)|_{y}:\;\,(x,y)\in
Graph(F)\right\\}$

determines a set-valued vector field of the right metric derivative of $F$ on
the graph of $F$.
We call the set

$D^{M}_{+}F(x_{0})=\bigcup_{y\in F(x_{0})}D^{M}_{+}F(x_{0})|_{y}=\lim_{x\to
x_{0}^{+}}\frac{1}{x-x_{0}}\left(F(x)\ominus F(x_{0})\right)$

the right metric derivative of $F$ at $x_{0}$.

  2. (ii)

Similarly, we say that $F$ is metrically differentiable from the left at
$x_{0}$ and define the sets $D^{M}_{-}F(x_{0})|_{y}$ ,
${\mathfrak{D}}^{M}_{-}F$ , $D^{M}_{-}F(x_{0})$.

  3. (iii)

If $F$ is both metrically differentiable from the right and from the left at
$x_{0}$, then we say that $F$ is metrically differentiable at $x_{0}$.

Below are some observations concerning Definition 3.6.

######  Remark 3.7.

  * (i)

Note that in $D^{M}_{+}F(x_{0})$, $D^{M}_{-}F(x_{0})$ the information on the
attachments to points in $Graph(F)$ is lost in contrary to
${\mathfrak{D}}^{M}_{+}F$ and ${\mathfrak{D}}^{M}_{-}F$.

  * (ii)

Although for single-valued functions the existence of derivative is defined
when the right derivative equals the left one, for set-valued functions the
equality $D^{M}_{-}F(x_{0})=D^{M}_{+}F(x_{0})$ is not satisfied in many cases,
and we do not require it for metric differentiability.

  * (iii)

Note that the sets
$D^{M}_{+}F(x_{0})|_{y},D^{M}_{-}F(x_{0})|_{y},D^{M}_{+}F(x_{0}),D^{M}_{-}F(x_{0})$
are compact as Hausdorff limits of compact sets (see Remark 3.2 (i), Remark
2.2 (ii)) and are unique as limits in a complete metric space.

  * (iv)

An equivalent definition of the right metric differentiability of $F$ at
$x_{0}$ is the existence of a continuous function
$\omega_{x_{0}}:[0,\infty)\to[0,\infty)$ with $\omega_{x_{0}}(0)=0$ and such
that for $x$ in a right neighborhood of $x_{0}$

$\sup_{y\in
F(x_{0})}{\mbox{h}aus}([x_{0},x]^{M}F|_{y},D^{M}_{+}F(x_{0})|_{y})\leq\omega_{x_{0}}(|x-x_{0}|).$

Similarly, for the left metric differentiability of $F$ at $x_{0},$ the
following condition should hold

$\sup_{y\in
F(x_{0})}{\mbox{h}aus}([x_{0},x]^{M}F|_{y},D^{M}_{-}F(x_{0})|_{y})\leq\omega_{x_{0}}(|x-x_{0}|).$

where $x$ is within a left neighborhood of $x_{0}$.

The next example demonstrates the disparity between the left and right
derivative at a boundary point of the graph.

######  Example 3.8.

Consider ${F:[-1,1]\rightarrow{\cal K}_{1}}$, $F(x)=[0,1+x]$, $x_{0}=0$.

First, we evaluate the right metric derivative $D^{M}_{+}F(0)|_{y}$, for all
$y\in F(0)$. The set $[0,h]^{M}F$, ${0<h<1}$ is derived from the following
collection of metric pairs

$\Pi({F(0)},{F(h)})=\big{\\{}(y,y),\ (1,\tilde{y}):\
y\in[0,1),\,\tilde{y}\in[1,1+h]\big{\\}}.$

For $y\in[0,1)$ we have $[0,h]^{M}F|_{y}=\\{0\\}$ and
$D^{M}_{+}F(0)|_{y}=\lim\limits_{h\to 0}[0,h]^{M}F|_{y}=\\{0\\}$, while for
$y=1$ we obtain

$[0,h]^{M}F|_{1}=\left\\{\frac{\tilde{y}-1}{h}:\,\tilde{y}\in[1,1+h]\right\\}=[0,1]\quad,\;D^{M}_{+}F(0)|_{1}=\lim_{h\to
0}[0,h]^{M}F|_{1}=[0,1].$

Now we calculate $\displaystyle D^{M}_{-}F(0)|_{y}$, for all $y\in F(0)$.
First note that for any $y\in[0,1)$ and any $0<h<1-y$ we have
$(y,y)\in\Pi({F(-h)},{F(0)})$, then we obtain ${[0,-h]^{M}F|_{y}=\\{0\\}}$ and
consequently $D^{M}_{-}F(0)|_{y}=\\{0\\}$. If $y=1$, then for any $0<h<1$ the
metric pair $(1-h,1)\in\Pi({F(-h)},{F(0)})$ leads to

$[0,-h]^{M}F|_{1}=\left\\{\frac{1-h-1}{-h}=1\right\\}\;{\mbox{a}nd}\;\,D^{M}_{-}F(0)|_{1}=\\{1\\}.$

Thus $D^{M}_{+}F(0)|_{y}=D^{M}_{-}F(0)|_{y}=\\{0\\}$ for $y\in[0,1)$, while
for the boundary point $(0,1)\in Graph(F)$ we have ${D^{M}_{+}F(0)|_{1}\neq
D^{M}_{-}F(0)|_{1}}$.
Note that this multifuction is metrically differentiable at $x=0$ according to
Definition 3.6.

In the following, we continue with Example 3.3 and show that $F$ from this
example is metrically differentiable, with the left and right derivatives
being equal.

######  Example 3.9.

Let $F:(0,2)\to{\cal K}_{1}$ be defined as in Example 3.3, namely
$F(x)=\\{x^{\alpha},x^{\beta}\\}$. Then, for $x_{0}=1$ we get
${D^{M}_{+}F(1)|_{1}=\lim\limits_{h\to
0^{+}}[1,1+h]^{M}F|_{1}=\\{\alpha,\beta\\},}$ while calculating for
${x_{0}\neq 1}$ we obtain

$D^{M}_{+}F(x_{0})|_{x_{0}^{\alpha}}=\lim\limits_{h\to
0^{+}}[x_{0},x_{0}+h]^{M}F|_{x_{0}^{\alpha}}=\\{\alpha
x_{0}^{\alpha-1}\\}\quad\hbox{and}\quad
D^{M}_{+}F(x_{0})|_{x_{0}^{\beta}}=\\{\beta x_{0}^{\beta-1}\\}.$

Thus, for each point $x\in(0,2)$, $x\neq 1$, $F(x)$ consists of two points,
and there is only one vector in $D^{M}_{+}F(x)|_{y}$ for $y\in F(x)$. In case
$x=1$ there is only one point $y=1$ in $F(1)$, yet
$D^{M}_{+}F(1)|_{1}=\\{\alpha,\beta\\}$.

It is straightforward to verify that for any ${y\in F(x)}$ and ${x\in(0,2)}$,
${D^{M}_{-}F(x)|_{y}=D^{M}_{+}F(x)|_{y}}$.

Now we extend the function of Example 3.3 to a function with values in ${\cal
K}_{2}$ and derive its metric left and right derivatives.

######  Example 3.10.

Let $F:(0,2)\to{\cal K}_{2}$ be defined by

$F(x)=\left\\{\begin{pmatrix}x^{\alpha}\\\
x^{\beta}\end{pmatrix},\begin{pmatrix}x^{\alpha+1}\\\
x^{\beta+1}\end{pmatrix}\right\\}\,\;\alpha,\beta\in{\mathbb{N}},\,\alpha\neq\beta.$

First we evaluate the metric divided differences. For $x\in(0,2)$ and $|h|$
sufficiently small, the collection of all metric pairs $\Pi({F(x)},{F(x+h)})$
for $x,x+h\leq 1$ or $x,x+h\geq 1$ is

$\displaystyle\Pi({F(x)},{F(x+h)})=\left\\{\left(\begin{pmatrix}x^{\alpha}\\\
x^{\beta}\end{pmatrix},\begin{pmatrix}(x+h)^{\alpha}\\\
(x+h)^{\beta}\end{pmatrix}\right)\ ,\ \left(\begin{pmatrix}x^{\alpha+1}\\\
x^{\beta+1}\end{pmatrix},\begin{pmatrix}(x+h)^{\alpha+1}\\\
(x+h)^{\beta+1}\end{pmatrix}\right)\right\\}.$

Denote $p_{1}=\begin{pmatrix}x^{\alpha}\\\ x^{\beta}\end{pmatrix}$,
$p_{2}=\begin{pmatrix}x^{\alpha+1}\\\ x^{\beta+1}\end{pmatrix}$. For $x\neq 1$
we have

$\displaystyle[x,x+h]^{M}F|_{p_{1}}=\left\\{\begin{pmatrix}\frac{(x+h)^{\alpha}-x^{\alpha}}{h}\\\
\frac{(x+h)^{\beta}-x^{\beta}}{h}\end{pmatrix}\right\\}\quad,\quad[x,x+h]^{M}F|_{p_{2}}=\left\\{\begin{pmatrix}\frac{(x+h)^{\alpha+1}-x^{\alpha+1}}{h}\\\
\frac{(x+h)^{\beta+1}-x^{\beta+1}}{h}\end{pmatrix}\right\\}.$

If $x=1$, the set $F(1)$ consists of a single element $p=\begin{pmatrix}1\\\
1\end{pmatrix}$, thus we obtain

$[1,1+h]^{M}F|_{p}=\left\\{\begin{pmatrix}\frac{(1+h)^{\alpha}-1}{h}\\\
\frac{(1+h)^{\beta}-1}{h}\end{pmatrix}\,,\,\begin{pmatrix}\frac{(1+h)^{\alpha+1}-1}{h}\\\
\frac{(1+h)^{\beta+1}-1}{h}\end{pmatrix}.\right\\}$

Now we find $D^{M}_{+}F(x)|_{p}$ for every $p\in F(x)$ and $x\in(0,2)$. If
$x\neq 1$ we get

$D^{M}_{+}F(x)|_{p_{1}}=\lim_{h\to
0^{+}}[x,x+h]^{M}F|_{p_{1}}=\left\\{\begin{pmatrix}\alpha x^{\alpha-1}\\\
\beta x^{\beta-1}\end{pmatrix}\right\\}\quad{\mbox{and}}\quad
D^{M}_{+}F(x)|_{p_{2}}=\left\\{\begin{pmatrix}{(\alpha+1)}x^{\alpha}\\\
{(\beta+1)}x^{\beta}\end{pmatrix}\right\\}.$

Calculations for $x=1$ result in the set
$D^{M}_{+}F(x)|_{p}=\left\\{\begin{pmatrix}\alpha\\\
\beta\end{pmatrix}\,,\,\begin{pmatrix}\alpha+1\\\
\beta+1\end{pmatrix}\right\\}$

Just as in Example 3.9 here we also have
$D^{M}_{-}F(x)|_{p}=D^{M}_{+}F(x)|_{p}$ for any ${p\in F(x)}$ and
${x\in(0,2)}$.

The following properties of the metric one-sided derivatives are derived from
their definition.

######  Theorem 3.11.

Let $F:(a,b)\to{\cal K}_{n}$ be metrically differentiable from the right
(left) at $x_{0}\in(a,b)$.

  * (i)

$F$ is right (left) continuous at $x_{0}$.

  * (ii)

If $F={\rm const}$, then $D^{M}_{+}F(x)|_{y}=\\{0\\}$
$(\,D^{M}_{-}F(x)|_{y}=\\{0\\}\,)$ for all $y\in F(x)$, $x\in(a,b)$.

  * (iii)

$D^{M}_{+}F(x_{0})|_{y_{0}}\neq\\{0\\}$
$(\,D^{M}_{-}F(x_{0})|_{y_{0}}\neq\\{0\\}\,)$ only for $(x_{0},y_{0})$ on the
boundary of $Graph(F)$.

###### Proof.

We prove the claims for the right derivative. The proofs for the left
derivative are analogous.

  * (i)

By Definition 3.6 $(i)$,

$\lim_{x\to x_{0}^{+}}\sup_{y_{0}\in
F(x_{0})}{\mbox{h}aus}\left(D^{M}_{+}F(x_{0})|_{y_{0}},\left\\{\frac{y-y_{0}}{x-x_{0}}\
:\ (y_{0},y)\in\Pi({F(x_{0})},{F(x)})\right\\}\right)=0.$

Multiplying the two sets above by $x-x_{0}$ and using the second relation in
(3), we get

$\lim_{x\to x_{0}^{+}}\sup_{y_{0}\in
F(x_{0})}\frac{1}{|x-x_{0}|}{\mbox{h}aus}\left((x-x_{0})D^{M}_{+}F(x_{0})|_{y_{0}},\\{{y-y_{0}}:\
(y_{0},y)\in\Pi({F(x_{0})},{F(x)})\\}\right)=0,$

and it follows that

$\lim_{x\to x_{0}^{+}}\sup_{y_{0}\in
F(x_{0})}{\mbox{h}aus}\left((x-x_{0})D^{M}_{+}F(x_{0})|_{y_{0}},\\{{y-y_{0}}:\
(y_{0},y)\in\Pi({F(x_{0})},{F(x)})\\}\right)=0.$

Since $D^{M}_{+}F(x_{0})$ is compact (see Remark 3.7 (iii)), we get

$\displaystyle\lim_{x\to x_{0}^{+}}\sup_{y_{0}\in
F(x_{0})}\|(x-x_{0})D^{M}_{+}F(x_{0})|_{y_{0}}\|=0,$

which leads to

$\lim_{x\to x_{0}^{+}}\sup\\{|y-y_{0}|\ :\
(y_{0},y)\in\Pi({F(x_{0})},{F(x)})\\}=0.$

Thus by (1)

$\lim_{x\to x_{0}^{+}}{\mbox{h}aus}(F(x),F(x_{0}))=0,$

which proves the continuity from the right of $F$ at $x_{0}$.

  * (ii)

This property follows from the definition of the metric divided difference and
the observation that if $F$ is constant, then for any two points
$x_{0},x_{1}$,

$\Pi({F(x_{0})},{F(x_{1})})=\\{(y,y):y\in F(x_{0})\\}$

and therefore $[x_{0},x_{1}]^{M}F|_{y}=\\{0\\}$ for any $y\in F(x_{0})$.

  * (iii)

Let $(x_{0},y_{0})$ be an interior point of $Graph(F)$. By Remark 3.2 (ii),
$[x_{0},x]^{M}F|_{y_{0}}=\\{0\\}$ for $x$ in a small enough neighborhood of
$x_{0}$, implying that $D^{M}_{+}F(x_{0})|_{y_{0}}=\\{0\\}$.

∎

##  4 The local metric linear approximant and $\alpha$-differentiability

In this section we define and study the rate of approximation of the local
metric linear approximant of the multifunction $F$ at a point $x_{0}\in(a,b)$.

######  Definition 4.1.

Let $x_{0}\in(a,b)$ and $F$ be metrically differentiable at $x_{0}$.

  1. (i)

The right local metric linear approximant of $F$ at $x_{0}$ anchored at
$y_{0}\in F(x_{0})$ is defined by

$L^{M}_{+}F|_{y_{0}}(x)=\\{y_{0}\\}+(x-x_{0})D^{M}_{+}F(x_{0})|_{y_{0}},\qquad
x\geq x_{0},\ x\in{\mathbb{R}}.$

The  right local metric linear approximant of $F$ at $x_{0}$ is the set-valued
function

$L^{M}_{+}F(x)=\bigcup\limits_{y\in F(x_{0})}L^{M}_{+}F|_{y}(x).$

  2. (ii)

Similarly, we define the left local metric linear approximant anchored at
$y_{0}\in F(x_{0})$,  $L^{M}_{-}F|_{y_{0}}(x)$, and the left local metric
linear approximant  $L^{M}_{-}F(x)$.

  3. (iii)

The (two-sided) local metric linear approximant of $F$ at $x_{0}$ anchored at
${y_{0}\in F(x_{0})}$ is

$L^{M}F|_{y_{0}}(x)=\begin{cases}L^{M}_{+}F|_{y_{0}}(x),&x\geq x_{0},\\\
L^{M}_{-}F|_{y_{0}}(x),&x<x_{0}.\end{cases}$

The (two-sided) local metric linear approximant  of $F$ at $x_{0}$ is the set-
valued function

$L^{M}F(x)=\begin{cases}L^{M}_{+}F(x),&x\geq x_{0},\\\
L^{M}_{-}F(x),&x<x_{0}.\end{cases}$

Obviously $L^{M}_{+}F(x_{0})=L^{M}_{-}F(x_{0})=L^{M}F(x_{0})=F(x_{0})$.

The following observations follow from this definition.

######  Theorem 4.2.

Let $F:{\mathbb{R}}\to{\cal K}_{n}$ be metrically differentiable at $x_{0}$.
Then for $x\in(a,b)$

  * (i)

${\mbox{h}aus}\left(L^{M}F|_{y_{0}}(x),\\{y_{0}\\}+(x-x_{0})[x_{0},x]^{M}F|_{y_{0}}\right)=o(|x-x_{0}|)$
uniformly in ${y_{0}\in F(x_{0})}$.

  * (ii)

${\mbox{h}aus}(F(x),L^{M}F(x))=o(|x-x_{0}|)$.

###### Proof.

(i) follows directly from Remark 3.7 (iv) and (3).
(ii) follows from (i), (2) and Remark 3.2 (iii). ∎

Note that (ii) in the above proposition implies that $L^{M}_{+}F(x)$ is a
right directive in the sense of Artstein [2, Definition 5.5] or local
approximation of first order in the sense of Nikolskiĭ [17, Definition 1].

To obtain error estimates of order $1+\alpha$ ($\alpha>0$) for the local
linear approximant, we introduce the notion of $\alpha$-differentiability of a
set-valued function at a point.

######  Definition 4.3.

  1. (i)

Let a multifunction $F$ be metrically differentiable from the right at
$x_{0}$. We call $F$ right metrically $\alpha$-differentiable at $x_{0}$ with
$\alpha>0$, if there is a constant $L>0$ and a right neighborhood $U$ of
$x_{0}$ such that for all $x\in U$

$\sup_{y\in
F(x_{0})}{\mbox{h}aus}([x_{0},x]^{M}F|_{y},D^{M}_{+}F(x_{0})|_{y})\leq
L|x-x_{0}|^{\alpha}.$

$F$ is called strongly right metrically differentiable at $x_{0}$ if
$\alpha=1$, namely

$\sup_{y\in
F(x_{0})}{\mbox{h}aus}([x_{0},x]^{M}F|_{y},D^{M}_{+}F(x_{0})|_{y})\leq
L|x-x_{0}|.$

  2. (ii)

Similarly we define these notions from the left.

  3. (iii)

If $F$ is right and left metrically $\alpha$-differentiable at $x_{0}$,
$\alpha>0$, then we say that $F$ is metrically $\alpha$-differentiable at
$x_{0}$ with $\alpha>0$. We say that $F$ is strongly metrically differentiable
at $x_{0}$ if it is strongly metrically differentiable from the right and from
the left at $x_{0}$.

######  Remark 4.4.

It is easy to verify that any real-valued function with Hölder continuous
derivative with exponent $\alpha\in(0,1]$ in a neighborhood of $x_{0}$ is
metrically $\alpha$-differentiable at ${x_{0}}$.

The next theorem shows that the local linear approximant of a metrically
$\alpha$-differentiable multifunction has order of approximation $1+\alpha$.

######  Theorem 4.5.

If $F$ is metrically $\alpha$-differentiable at $x_{0}$, then for $x$ from a
neighborhood of $x_{0}$

${\mbox{h}aus}\left(F(x),L^{M}F(x)\right)\leq L|x-x_{0}|^{1+\alpha}.$

###### Proof.

We consider the case $x>x_{0}$. For $x<x_{0}$ the proof is similar.

By Definition 4.1, by Remark 3.2 (iii) and by (2) we get for any $y\in
F(x_{0})$

$\displaystyle{\mbox{h}aus}(F(x),L^{M}F(x))$
$\displaystyle={\mbox{h}aus}(F(x),L^{M}_{+}F(x))$ $\displaystyle\leq\sup_{y\in
F(x_{0})}{\mbox{h}aus}\left(\\{y\\}+(x-x_{0})[x_{0},x]^{M}F|_{y}\
,L^{M}_{+}F|_{y}(x)\right).$

Then, using properties (3) of the Hausdorff distance and Definitions 4.1, 4.3,
we get

$\\\ haus(F(x),L^{M}_{+}F(x))\leq|x-x_{0}|\sup_{y\in
F(x_{0})}{\mbox{h}aus}([x_{0},x]^{M}F|_{y},D^{M}_{+}F(x_{0})|_{y})\leq
L|x-x_{0}|^{1+\alpha},$

which proves the claim for $x>x_{0}$ . ∎

A direct consequence of Theorem 4.5 is

######  Corollary 4.6.

If $F$ is strongly metrically differentiable at $x_{0}$, then

${\mbox{h}aus}(F(x),L^{M}F(x))=O(|x-x_{0}|^{2}).$

In the next example we present a strongly metrically differentiable
multifunction at a point $x=0$ and construct its local metric linear
approximant.

######  Example 4.7.

Let $F:[-1,1]\to{\cal K}_{1}$ be defined as

$F(x)=\begin{cases}[0,2-x^{2}]\,,&x\in[-1,0],\\\
[0,2+x]\cup\\{-x^{2}\\}\,,&x\in[0,1].\end{cases}$

First we show that $F$ is strongly right metrically differentiable at $x=0$.
We start by calculating the metric divided differences $[0,h]^{M}F|_{y}$,
$y\in F(0)=[0,2]$, $0<h\leq 1$

For $y\in(0,2)$, by Remark 3.2 (ii) and by Definition 3.6 we have

$[0,h]^{M}F|_{y}=\\{0\\}\quad\hbox{and}\quad D^{M}_{+}F(0)|_{y}=\\{0\\}.$

If $y=2$, similarly to the calculations in Example 3.8 we obtain

$[0,h]^{M}F|_{2}=[0,1]\quad\hbox{and}\quad D^{M}_{+}F(0)|_{2}=[0,1].$

Thus for ${y\in(0,2]}$ we have
${{\mbox{h}aus}\left([0,h]^{M}F|_{y},D^{M}_{+}F(0)|_{y}\right)=0}$.

When $y=0$, ${0<h\leq 1}$, one can easily see that the only metric pairs in
$\Pi({F(0)},{F(h)})$ containing $y=0$ are $(0,0)$ and $(0,-h^{2})$, thus
obtaining

$[0,h]^{M}F|_{0}=\left\\{0,\frac{-h^{2}-0}{h}\right\\}\quad\mbox{and}\quad
D^{M}_{+}F(0)|_{0}=\lim\limits_{h\to 0^{+}}[0,h]^{M}F|_{0}=\\{0\\}.$ (6)

This leads to ${\mbox{h}aus}([0,h]^{M}F|_{0},D^{M}_{+}F(0)|_{0})=h$.

Note that the elements in $[0,h]^{M}F|_{0}$ are the values of the first
divided differences of the real-valued functions (boundary functions)
$s_{1}(x)=0$, $s_{2}(x)=-x^{2}$ based on $x_{0}=0$, $x_{1}=h$. Accordingly,
the single element of $D^{M}_{+}F(0)|_{0}$ is the value of right derivatives
of $s_{1}(x),s_{2}(x)$ at $x=0$.

Summing up the above, we conclude that

$\sup_{y\in F(0)}{\mbox{h}aus}([0,h]^{M}F|_{y},D^{M}_{+}F(0)|_{y})\leq h,$

therefore $F$ is strongly right metrically differentiable at $x=0$.

Now we check from the left of $x=0$. For $y\in[0,2)$, we again refer to the
calculations in Example 3.8 from which one can easily get that
${[-h,0]^{M}F|_{y}=\\{0\\}}$ for any $0<h<\sqrt{2-y}$. Consequently
$D^{M}_{-}F(0)|_{y}=\\{0\\}$ and
${\mbox{h}aus}([-h,0]^{M}F|_{y},D^{M}_{-}F(0)|_{y})=0$.

When $y=2$, computations similar to those used in (6), lead to
${[-h,0]^{M}F|_{2}=\\{h\\}}$ and ${D^{M}_{-}F(0)|_{2}=\\{0\\}}$. Thus we can
derive

$\sup_{y\in F(0)}{\mbox{h}aus}([-h,0]^{M}F|_{y},D^{M}_{-}F(0)|_{y})\leq h,$

and as a result, we conclude that $F$ is strongly left metrically
differentiable at $x=0$. Thus $F$ is strongly metrically differentiable at
$x=0$.

Finally we construct the (two-sided) local metric linear approximant of $F$ at
$x=0$ in accordance with Definition 4.1. Consider first $x\geq 0$. Since
$D^{M}_{+}F(0)|_{y}=\\{0\\}$ for any $y\in[0,2)$, we get

${L^{M}_{+}F|_{y}(x)=\\{y\\}+(x-0)D^{M}_{+}F(0)|_{y}=\\{y\\}}.$

As for $y=2$, we have $D^{M}_{+}F(0)|_{2}=[0,1]$ and then
${L^{M}_{+}F|_{2}(x)=\\{2\\}+x\cdot[0,1]=[2,2+x]}$. Thus

$L^{M}_{+}F(x)=\bigcup\limits_{y\in[0,2]}L^{M}_{+}F|_{y}(x)=[0,2)\bigcup\
[2,2+x]=[0,2+x].$

If $x<0$, then $L^{M}_{-}F(x)=[0,2]$, since $D^{M}_{-}F(0)|_{y}=\\{0\\}$,
$\forall y\in[0,2]$. Summing up the results we obtain

$L^{M}F(x)=\begin{cases}[0,2]\,,&x<0,\\\ [0,2+x]\,,&x\geq 0.\end{cases}$

It is easy to see that for $|h|\leq 1$

${\mbox{h}aus}\left(F(h),L^{M}F(h)\right)=h^{2}$

in accordance with Corollary 4.6.

Acknowledgment

The authors thank Ron Goldman for initial discussions on possible definitions
of derivatives of SVFs and for suggesting the term ”anchored”.

## References

  * [1] Z. Artstein, Piecewise linear approximations of set-valued maps, Journal of Approximation Theory 56 (1989), 41–47.
  * [2] Z. Artstein, A calculus for set-valued maps and set-valued evolution equations. Set-Valued Analysis 3, (1995), 213–261.
  * [3] J.-P. Aubin and H. Frankowska, Set-valued analysis, Birkhäuser, Boston, (1990).
  * [4] R. Baier, E. Farkhi, Directed derivatives of convex compact-valued mappings, in Advances in convex analysis and global optimization, 501–514, Nonconvex Optim. Appl., 54, Kluwer Acad. Publ., 2001.
  * [5] H.T. Banks, M.Q. Jacobs, A differential calculus for multifunctions, J. of Mathematical Analysis and Applications 29, (1970), 246–272.
  * [6] E. E. Berdysheva, N. Dyn, E. Farkhi, A. Mokhov, Metric approximation of set-valued functions of bounded variation, J. Comput. Appl. Math. 349 (2019), 251–264.
  * [7] F. S. De Blasi, On the differentiability of multifunctions, Pacific J. Math. 66, (1976), 67–81.
  * [8] N. Dyn, E. Farkhi and A. Mokhov. Approximations of set-valued functions by metric linear operators. Constr. Approx. 25, (2007), 193–209.
  * [9] N. Dyn, E. Farkhi, A. Mokhov, Approximations of Set-Valued Functions: Adaptation of Classical Approximation Operators, Imperial College Press, London (2014).
  * [10] N. Dyn, E. Farkhi, A. Mokhov, The metric integral of set-valued functions, Set-Valued Var. Anal. 26 (2018), 867–885.
  * [11] N. Dyn, E. Farkhi and A. Mokhov. High order approximation of set-valued functions. Constr. Approx. 57, (2023), 521–546.
  * [12] S. Gautier, Affine and eclipsing multifunctions, Numer. Funct. Anal. Optim. 11 (1990) 679–699.
  * [13] M. Hukuhara, Integration des applications mesurables dont la valeur est un compact convexe, Funkcialaj Ekvacioj 10 (1967) 205–223.
  * [14] E. Khmaladze, W. Weil, Differentiation of sets - The general case. Journal of Mathematical Analysis and Applications 413 (2014), 291–310.
  * [15] C. Lemarechal, J. Zowe, The eclipsing concept to approximate a multi-valued mapping, Optimization 22 (1991), 3–37.
  * [16] B. Mordukhovich, Variational Analysis and Generalized Differentiation. I. Basic Theory, Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], vol. 330, Springer-Verlag (2006)
  * [17] Nikolskiĭ, M. S. Local approximation of first order to set-valued mappings. Set-valued analysis and differential inclusions, Progr. Systems Control Theory, 16, Birkhäuser, (1993), 149–-156.
  * [18] C.H.J. Pang, Generalized differentiation with positively homogeneous maps: applications in set-valued analysis and metric regularity, Mathematics of Operations Research, 36 (2011), 377–397.
  * [19] R. T. Rockafellar, R. Wets, Variational analysis , Springer, (1998).
  * [20] R. Schneider. Convex Bodies: the Brunn–Minkowski Theory. Cambridge University Press, Cambridge, 1993.
  * [21] D. Silin, On set-valued differentiation and integration, Set-Valued Anal. 5 (1997) 107–146.