# Pseudomode description of general open quantum system dynamics: non-
perturbative master equation for the spin-boson model

Graeme Pleasance<EMAIL_ADDRESS>Quantum Research Group, School of
Chemistry and Physics, University of KwaZulu-Natal, Durban, 4001, South Africa
Francesco Petruccione  Quantum Research Group, School of Chemistry and
Physics, University of KwaZulu-Natal, Durban, 4001, South Africa  National
Institute for Theoretical and Computational Sciences (NITheCS), South Africa

###### Abstract

We outline a non-perturbative approach for simulating the behavior of open
quantum systems interacting with a bosonic environment defined by a
generalized spectral density function. The method is based on replacing the
environment by a set of damped harmonic oscillators—the pseudomodes—thereby
forming an enlarged open system whose dynamics is governed by a Markovian
master equation. Each pseudomode is connected to one of the poles of the
spectral density when analytically continued to the lower-half complex
frequency plane. Here, we extend previous results to a completely generic
class of open system models, and discuss how our framework can be used as a
powerful and versatile tool for analyzing non-Markovian open system dynamics.
The effectiveness of the method is demonstrated on the spin-boson model by
accurately benchmarking its predictions against numerically exact results.

In any realistic setting a quantum system will interact with its surrounding
environment, and so the theoretical modelling of its behavior must account for
the effects of this interaction [1]. For environments which are structured or
strongly coupled to the system, these interactions can induce time-dependent
changes in the system state that are strongly characterized by memory, and for
these cases the resulting behavior is non-Markovian [2, 3, 4]. Over recent
years the study of non-Markovian systems has seen intense focus not only out
of fundamental interest [5, 6, 7, 8, 9, 10], but also due to the increasing
ability to probe quantum systems over shorter time and length scales at which
memory effects play a heightened role. For example, energy transport in light-
harvesting complexes and nanoscale devices [11, 12, 13], as well as photonic
crystals [14, 15, 16], have been shown to display distinctive features
associated with non-Markovian dynamics. Applications of non-Markovianity to
quantum metrology [17, 18] and quantum information processing [19, 20, 21]
have too been the subject of much interest.

Various techniques capable of simulating memory effects within the dynamics of
open systems have been introduced to date, including the hierarchical
equations of motion (HEOM) [22, 23, 24], path integrals [25, 26], and tensor
networks [27, 28, 29], among others [30, 31, 32, 33, 34]. While a number of
these methods rely on keeping track of the environment state explicitly [32,
33, 34], an alternative approach consists of tracing out the environmental
degrees of freedom in favor of obtaining a master equation solely for the
reduced density matrix of the system [1, 30, 31]. Such master equations can in
principle be derived non-perturbatively through the use of projection operator
techniques [35, 36, 37]. However, their complicated structure often renders an
exact description infeasible without the use of certain approximations. The
Born-Markov approximations can lead to a much more convenient description in
terms of a Markovian master equation either in Bloch-Redfield [38] or Gorini-
Kossakowski-Sudarshan-Lindblad (GKSL) form [39, 40], but this typically comes
at the expense of restricting the validity of results to weak coupling
regimes.

At the same time, another class of approaches for treating complex non-
Markovian dynamics have been developed using the idea of mapping the initial
problem onto a simpler Markovian one [41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64]. In particular, the
pseudomode method [42, 45, 46, 47, 48, 49] implements such a mapping by
expanding the system with a set of discrete bosonic modes (the pseudomodes),
whose properties are derived from certain analytical features of the
environment spectral density. The enlarged system comprising the original
system and pseudomodes obeys an exact Markovian master equation which can then
be efficiently simulated using standard techniques [65, 66, 67]. Previously,
this method has been applied to open system models valid within the rotating
wave approximation (RWA) and with zero temperature environments [42, 45], and
has recently seen extensions to non-RWA cases involving several specific types
of spectral density function [41, 43].

In this Letter, we go beyond all previous limitations to extend the validity
of the pseudomode mapping to a completely generic class of open system models.
We achieve this by adapting a general proof from [43] to show that the reduced
dynamics of an open system coupled to bosonic environment can be reproduced
exactly when coupled to a small number of pseudomodes via a possibly non-
Hermitian form of interaction. To demonstrate the method, we apply it to a
paradigmatic open system model—the spin-boson model—finding it able to capture
the exact non-Markovian dynamics at strong system-bath coupling, whilst
simultaneously being implemented at low computational costs.

Physical and auxiliary models.— We consider a microscopic model of a generic
open system $S$ coupled linearly to a bosonic environment $E$ with total
Hamiltonian $H=H_{S}+H_{E}+H_{I}$. Here, $H_{S}$ and $H_{E}$ are the free
Hamiltonians of the system and environment; the latter is written as
$H_{E}=\sum_{k}\omega_{k}a^{\dagger}_{k}a_{k}$, where $\omega_{k}$ is the
frequency of the $k$-mode of the environment and $a_{k}$ ($a^{\dagger}_{k}$)
is the corresponding annihilation (creation) operator satisfying
$[a_{k},a^{\dagger}_{k^{\prime}}]=\delta_{kk^{\prime}}$ ($\hbar=1$). We assume
an interaction of the form $H_{I}=A\otimes B$, with $A=A^{\dagger}$ and
$B=\sum_{k}g_{k}(a_{k}+a^{\dagger}_{k})$ observables of the system and
environment, respectively. Note that the more general case with $H_{I}$
replaced by $H_{I}=\sum_{\alpha}A_{\alpha}\otimes B_{\alpha}$ could also be
considered, but we choose to avoid it for simplicity.
Our main focus here will be on describing the non-Markovian dynamics of open
quantum systems interacting with environments that are both Gaussian and
stationary—namely, environments which satisfy
$\mathinner{[{H_{E}},{\rho_{E}(0)}]}=0$, where $\rho_{E}(0)$ is a Gaussian
state of $E$. Consequently, if the initial state of the total system $S+E$ is
factorized as $\rho_{SE}(0)=\rho_{S}(0)\otimes\rho_{E}(0)$, then the reduced
system dynamics obtained via

$\rho_{S}(t)=\text{Tr}_{E}\Big{\\{}e^{-iHt}\rho_{S}(0)\otimes\rho_{E}(0)e^{iHt}\Big{\\}},$
(1)

will be fully determined by the two-time correlation function (assuming
vanishing first order moments $\langle B(t)\rangle=0$)

$\displaystyle C(\tau)$ $\displaystyle=\langle B(t+\tau)B(t)\rangle$
$\displaystyle=\frac{1}{2\pi}\int^{\infty}_{-\infty}d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad\tau\geq
0,$ (2)

with $\langle B(t+\tau)B(t)\rangle=\text{Tr}_{E}[B(\tau)B\rho_{E}(0)]$ and
$B(\tau)=e^{iH_{E}\tau}Be^{-iH_{E}\tau}$. Here we have also introduced the
spectral density (SD) $\gamma(\omega)$ as the Fourier transform of the
correlation function $C(\tau)$. By definition, the SD is a positive, real-
valued function describing the frequency dependent properties of the coupling
between the system and environment [1].

In the following it will prove very useful to model the interaction using a
generalized form of spectral density. More precisely, the only assumptions we
shall impose on $\gamma(\omega)$ are for it to be a meromorphic function when
analytically continued to the complex $\omega$-plane, and for $\gamma(\omega)$
to fall of faster than $\sim O(1/|\omega|)$ in the limit
$|\omega|\rightarrow\infty$. The usefulness of these assumptions is that they
allow us to evaluate the two-time correlation function (Pseudomode description
of general open quantum system dynamics: non-perturbative master equation for
the spin-boson model) analytically in terms of the poles and residues of
$\gamma(\omega)$ in the lower-half complex $\omega$-plane. To do so, let us
first enumerate these poles according to the index $l=1,...,N$. By then
writing Eq. (Pseudomode description of general open quantum system dynamics:
non-perturbative master equation for the spin-boson model) as a contour
integral $\int d\omega\rightarrow\oint_{C}d\omega$, with $C$ a closed
semicircular contour in lower-half plane, we can straightforwardly apply the
residue theorem to obtain

$C(\tau)=-i\sum^{N}_{l=1}r_{l}e^{-iz_{l}\tau},\qquad\tau\geq 0,$ (3)

where $z_{l}=\xi_{l}-i\lambda_{l}$ are the locations of the poles of
$\gamma(\omega)$, and $r_{l}=\text{Res}_{\omega=z_{l}}[\gamma(\omega)]$ are
their corresponding residues. By setting $\tau=0$ in Eqs. (Pseudomode
description of general open quantum system dynamics: non-perturbative master
equation for the spin-boson model) and (3) we may also establish the relation
$\sum_{l}r_{l}=-(2\pi i)^{-1}\int d\omega\,\gamma(\omega)$, implying the
residues $r_{l}$ have no net real part.

Figure 1:  Schematic representation of the physical and auxiliary models. (a)
An open system $S$ couples to a structured environment $E$ with a frequency
dependent SD. (b) The same system couples to a finite number of pseudomodes
$M$ with non-unitary evolution, whose dynamics when attached to $S$ is
governed by the Markovian master equation (7).

We now proceed to outline the auxiliary model that will allow us to represent
the exact time evolution of $\rho_{S}(t)$ within an enlarged Markovian open
system (see Fig. 1). Here, the physical environment $E$ is replaced with a set
of discrete bosonic modes—the pseudomodes $M$—which are defined in terms of
the positions and residues of the poles of the SD (Pseudomode description of
general open quantum system dynamics: non-perturbative master equation for the
spin-boson model). The Hamiltonian of the enlarged system is given as
$H_{0}=H_{S}+H_{M}+H^{\prime}_{I}$, whereby

$\displaystyle H_{M}$ $\displaystyle=\sum_{l}\xi_{l}b^{\dagger}_{l}b_{l},$ (4)
$\displaystyle H^{\prime}_{I}$ $\displaystyle=A\otimes B^{\prime}.$ (5)

Above, $b_{l}$ ($b^{\dagger}_{l}$) is the annihilation (creation) operator of
the $l$-pseudomode satisfying the bosonic commutation relation
$[b_{l},b^{\dagger}_{l^{\prime}}]=\delta_{ll^{\prime}}$, and

$B^{\prime}=\sum_{l}g^{\prime}_{l}\Big{(}b_{l}+b^{\dagger}_{l}\Big{)},\qquad
g^{\prime}_{l}=\sqrt{-ir_{l}}.$ (6)

Notice here that the coupling constant $g^{\prime}_{l}$ only appears outside
the bracket rather than $g^{\prime}_{l}$ and its complex conjugate; as such,
the Hamiltonian $H_{0}$ will generally be non-Hermitian due to the residues
$r_{l}$ on which the couplings $g^{\prime}_{l}$ depend being complex, except
for cases in which the SD adopts a Lorentzian structure. The implications of
this will be discussed in more detail below once we have introduced the main
result.
For the current configuration, the combined density matrix $\rho(t)$ of the
open system and pseudomodes obeys a Markovian master equation

$\frac{d}{dt}\rho(t)=\mathcal{L}\rho(t)=-i[H_{0},\rho(t)]+\mathcal{D}\rho(t),$
(7)

with the superoperator

$\mathcal{D}\rho=2\sum_{l}\lambda_{l}\bigg{(}b_{l}\rho
b^{\dagger}_{l}-\frac{1}{2}\big{\\{}b^{\dagger}_{l}b_{l},\rho\big{\\}}\bigg{)},$
(8)

describing the local dissipation of each pseudomode occurring at rate
$2\lambda_{l}$. The reduced system state evolving under (7) can then be
obtained as

$\rho^{\prime}_{S}(t)=\text{Tr}_{M}\Big{\\{}e^{\mathcal{L}t}\rho_{S}(0)\otimes\rho_{M}(0)\Big{\\}},$
(9)

where we have assumed the pseudomodes to be initially uncorrelated with the
system.

With these details in place, we now look to establish an exact equivalence of
the reduced system dynamics expressed through Eqs. (1) and (9). We proceed by
utilizing an approach recently developed in [43], in which the authors
introduced a rigorous proof for determining the equivalence of an open system
dynamics resulting from two different types of environment; one comprising a
Gaussian bosonic reservoir with free unitary evolution, and the other a set of
discrete bosonic modes, that together with the open system satisfy a Lindblad-
type evolution. In the Supplemental Material we adapt this proof in accordance
with [42, 41] to account for the non-Hermitian form of coupling between the
system and pseudomodes 111See Supplemental Material at [URL will be inserted
by publisher] for further details on the proof.. The proof is performed in two
stages: first, by purifying $\rho(t)$ on an extended state space $S+M+R$, we
show that the reduced dynamics (9) can be exactly reproduced in an extended
auxiliary model where the pseudomodes are individually coupled to Markovian
reservoirs $R$, and the total state of $S+M+R$ follows a pseudo-unitary
evolution. Second, we show that the reduced system dynamics in the physical
model is equivalent to the reduced dynamics of $S$ in the extended auxiliary
model assuming the auxiliary environment $M+R$ is initially in a Gaussian
state. Since the free evolution of both environments $E$ and $M+R$ is
(pseudo-) unitary, this is achieved by matching the two-time correlation
functions $C(\tau)$ and $C^{\prime}(\tau)$ of the respective coupling
operators $B$ and $B^{\prime}$. Overall, we then find that

$\rho^{\prime}_{S}(t)=\rho_{S}(t).$ (10)

Establishing an exact equivalence between the reduced dynamics of $S$ in the
two models means we have not only reaffirmed the results of [41, 42, 45], but
also now extended them to a generic class of open system models within the
restrictions above.

Application to arbitrary SDs.— Following on from this proof, we will now
discuss some details relevant to the application of the result. One key
consideration is that because the correspondence in Eq. (10) is entirely
determined by the equivalence between the two-time correlation functions
$C(\tau)$ and $C^{\prime}(\tau)$, the mapping will only be exact when
$C(\tau)$ can be written as a finite sum of complex exponentials. This is
equivalent to saying $\gamma(\omega)$ is constrained to have the same
frequency dependence as the effective spectral density
$\gamma^{\prime}(\omega)=\int^{\infty}_{-\infty}d\tau\,C^{\prime}(\tau)e^{i\omega\tau}$:

$\gamma^{\prime}(\omega)=2\sum^{N}_{l=1}\frac{r^{R}_{l}(\omega-\xi_{l})+r^{I}_{l}\lambda_{l}}{(\omega-\xi_{l})^{2}+\lambda^{2}_{l}},$
(11)

where $r_{l}\equiv r^{R}_{l}+ir^{I}_{l}$. For other forms of SD, including
those related to thermal baths (see below), the correspondence between
$\gamma(\omega)$ and $\gamma^{\prime}(\omega)$ cannot be made exact [43, 44],
and so for these cases we will instead rely on optimizing the parameters of
Eq. (11) as those most accurately representing the original SD function to
obtain

$C(\tau)\approx C^{\prime}(\tau),\quad\tau\geq 0.$ (12)

In practice this can be achieved by fitting the SD of the physical problem
with a set of $N$ basis functions
$\gamma^{\prime}_{l}(\omega)=2[r^{R}_{l}(\omega-\xi_{l})+r^{I}_{l}\lambda_{l}]/[(\omega-\xi_{l})^{2}+\lambda^{2}_{l}]$—or
with $N$ complex exponentials $C^{\prime}_{l}(\tau)=-ir_{l}\exp(-iz_{l}\tau)$
in the time domain—where each pole of $\gamma^{\prime}_{l}(\omega)$ in the
lower-half plane is associated to one of the $N$ pseudomodes in the auxiliary
model. The residues of these poles are again used to parameterize the
couplings $g^{\prime}_{l}$ within the enlarged system Hamiltonian. In fact,
since the non-Hermitian terms of $H_{0}$ are generated through the imaginary
part of these couplings $g^{\prime}_{l}$, their role in determining
$\gamma^{\prime}(\omega)$ may be conveniently illustrated by imposing a
Hermitian form of interaction, i.e., by setting $r^{R}_{l}=0$ and
$r^{I}_{l}>0$, so that $g^{\prime}_{l}$ are constrained to be real. From Eq.
(11) we see this corresponds to $\gamma^{\prime}(\omega)$ being written as a
sum of positively weighted Lorentzians, implying the effect of the non-
Hermitian terms is to generalize the form of $\gamma^{\prime}(\omega)$ from a
simple sum of Lorentzians to more complicated functions which are either (i)
non-Lorentzian, or (ii) Lorentzian but with negative weights. As one might
expect, these functions are capable of representing a much greater variety of
environments compared to when $g^{\prime}_{l}$ are real [54, 55], thus
enabling a more efficient simulation of the problem in general.

Figure 2:  Time evolution of $\langle\sigma_{z}(t)\rangle$ obtained from the
pseudomode master equation (solid blue line and solid black triangles) and the
HEOM (solid red points), (a) for $\omega_{0}=0.25\Delta$, (b)
$\omega_{0}=0.5\Delta$, and (c) $\omega_{0}=\Delta$. The black triangles
indicate the result generated by mapping the $n=1$ term of Eq. (16) onto a
single pseudomode and with a local dephasing term acting on the TLS. The other
parameters are $\beta\Delta=1$, $\epsilon=0.5\Delta$, $\Gamma=0.05\Delta$ and
$\alpha=0.25\Delta$. For all plots, the Matsubara part of the correlation
function (16) is evaluated up to a total of $1.5\times 10^{3}$ terms.

We finish this section by briefly commenting on the connection of our approach
with that of Mascherpa et. al. in Ref. [44]. There, the authors have
introduced an analogous treatment based on using network of damped harmonic
oscillators to replicate the effect of a bosonic environment on an open
system, assuming the bath correlation function can be approximated as a sum of
complex exponentials (3). Our approach differs from [44] primarily in the
ansatz used to construct the effective environment as well as the master
equation used to model the dynamics of the enlarged system; more precisely,
the auxiliary modes used to replace the environment in [44] are interacting,
and their dynamics when coupled to the system is described by a Lindblad
master equation. The main consequence of these differences is that the
couplings and decay rates of the auxiliary modes cannot be directly related to
those of the correlation function (3), and so these parameters must be
extracted by performing an additional step of inverting a system of non-linear
equations depending on the weights and exponents of $C^{\prime}(\tau)$. In
this regard, our method may offer certain numerical advantages for situations
where it is difficult to perform such an inversion, since apart from the
fitting of the two-time correlation function $C(\tau)$, no further steps are
necessary in order to simulate the physical problem.

Spin-boson model.— We now proceed to outline an example of the mapping applied
to the thermal bath case, where the dynamics of the system is described by the
spin-boson Hamiltonian [69]

$H=\frac{\epsilon}{2}\sigma_{z}+\frac{\Delta}{2}\sigma_{x}+\sum_{k}\omega_{k}a^{\dagger}_{k}a_{k}+\sigma_{z}\otimes\sum_{k}g_{k}\big{(}a_{k}+a^{\dagger}_{k}\big{)}.$
(13)

Here, the system and interaction Hamiltonians
$H_{S}=(\epsilon/2)\sigma_{z}+(\Delta/2)\sigma_{x}$ and
$H_{I}=\sigma_{z}\otimes B$ are written in terms of the Pauli
spin-$\frac{1}{2}$ operators $\sigma_{i}$ ($i=x,y,z$) in the subspace
$\\{|g\rangle,|e\rangle\\}$ of the two-level system $S$ (TLS), with
$\sigma_{z}=|e\rangle\langle e|-|g\rangle\langle g|$. The environment is
initially in a thermal equilibrium state $\rho_{E}(0)=\exp(-\beta
H_{E})/\text{Tr}\big{[}\exp(-\beta H_{E})\big{]}$ at temperature $T=1/\beta$
($k_{B}=1$), such that the corresponding two-point correlation function reads

$C(\tau)=\frac{1}{\pi}\int^{\infty}_{0}d\omega\,J(\omega)\Big{[}\coth\left(\frac{\beta\omega}{2}\right)\cos(\omega\tau)-i\sin(\omega\tau)\Big{]},$
(14)

having introduced the one-sided spectral density function
$J(\omega)=\pi\sum_{k}g^{2}_{k}\delta(\omega-\omega_{k})$ [1]. For
convenience, we will consider $J(\omega)$ to be in an underdamped Brownian
oscillator form
$J(\omega)=\alpha\omega^{2}_{0}\Gamma\omega/[(\omega^{2}_{0}-\omega^{2})^{2}+\Gamma^{2}\omega^{2}]$,
where $\Gamma$ and $\omega_{0}$ define a characteristic width and resonance
frequency, respectively, and $\alpha$ is the reorganization energy of the bath
[70].

The correlation function corresponding to $\gamma(\omega)$ can be determined
analytically as $C(\tau)=C_{0}(\tau)+M(\tau)$, where

$\displaystyle
C_{0}(\tau)=\frac{\alpha\omega^{2}_{0}}{4\Omega}e^{-\Gamma\tau/2}\left[\coth\left(\frac{\beta(\Omega-i\Gamma/2)}{2}\right)e^{i\Omega\tau}+\text{h.c.}\right]$
$\displaystyle\quad+\frac{\alpha\omega^{2}_{0}}{4\Omega}e^{-\Gamma\tau/2}\left[e^{-i\Omega\tau}-e^{i\Omega\tau}\right],$
(15)

and

$\displaystyle M(\tau)=\sum^{\infty}_{n=1}c_{n}e^{-\nu_{n}\tau},$ (16)
$\displaystyle
c_{n}=\frac{-2\alpha\Gamma\omega^{2}_{0}\nu_{n}}{\beta[\nu^{2}_{n}+(\Omega+i\Gamma/2)^{2}][\nu^{2}_{n}+(\Omega-i\Gamma/2)^{2}]},$

with $\Omega=\sqrt{\omega^{2}_{0}-(\Gamma/2)^{2}}$, and Matsubara frequencies
$\nu_{n}=2\pi n/\beta$. Since the analytic part of the correlation function
(Pseudomode description of general open quantum system dynamics: non-
perturbative master equation for the spin-boson model) is already in form of
Eq. (3), we choose to capture the effect of the remaining Matsubara
contributions by first fitting $M(\tau)$ with the sum of two real exponentials
[41], such that
$M^{\prime}(\tau)=W_{1}e^{-\gamma_{1}\tau}+W_{2}e^{-\gamma_{2}\tau}$, where
$C^{\prime}(\tau)=C_{0}(\tau)+M^{\prime}(\tau)$. We emphasize that it is the
use of complex couplings $g^{\prime}_{l}$ which allows us to represent
$C_{0}(\tau)$ exactly using only two pseudomodes. In Fig. 2, we compare the
reduced system dynamics computed from the pseudomode master equation and the
HEOM for different values of $\omega_{0}$ in the strong coupling and narrow
bath regime, assuming the TLS to be initialized in the state
$\rho_{S}(0)=|e\rangle\langle e|$. The two approaches have been implemented
numerically using the Python library QuTiP [71] and the integrated QuTiP
package BoFiN-HEOM provided in [72, 73]. As expected, we find close
correspondence between the predictions of the master equation and HEOM over
both short and long time scales. In addition to the fitting approach, we also
choose to simulate the TLS dynamics by retaining only a single Matsubara term
in the expansion (16), while treating the remaining terms as delta functions
[23]. Since $M(\tau)$ is sharply peaked around $\tau=0$, for $\tau\ll
1/\nu_{n}$ one can write it approximately as $M(\tau)\approx
c_{1}e^{-\nu_{1}\tau}+\sum^{\infty}_{n=2}\frac{c_{n}}{\nu_{n}}\delta(\tau)$.
The first term of $M(\tau)$ is then mapped to a single pseudomode following
the same procedure outlined above. On the other hand, the delta function can
be accounted for by adding a local dephasing term to the Liouvillian in Eq.
(7) of the form
$\mathcal{L}_{D}\rho=2\gamma_{D}(\sigma_{z}\rho\sigma_{z}-\rho)$, where
$\gamma_{D}=\sum_{n>1}(c_{n}/\nu_{n})$. This is done analogously to how the
‘Taninmura terminator’ is implemented within the HEOM. For the cases shown
(see Fig. 2), we were able with the approximated expansion to simulate the
master equation at an even smaller cost than using the full Matsubara
decomposition (16), which becomes particularly advantageous at strong coupling
when the local dimension of each pseudomode is generally needed to be higher
for convergence.

Summary & discussion.— In conclusion, we have presented a novel framework for
simulating the non-equilibrium behavior of open quantum systems within
structured environments. Our method relies on replicating the effect of a
bosonic environment on a general open system with a collection of pseudomodes,
i.e., discrete bosonic modes, where each pseudomode is connected to one of the
poles of the SD in the lower-half complex frequency plane. Since the mapping
is non-perturbative in the original system-bath coupling, our formalism
enables the full inclusion of memory effects within the reduced dynamics of
the simulated problem. This has been demonstrated for the spin-boson model by
benchmarking the results obtained from the pseudomode master equation against
those from the numerically exact HEOM. Besides having numerous theoretical
applications, our findings open up ways to possibly analyze the complex
dynamics occurring within experimental setups relevant to e.g., quantum
thermodynamics [53, 74], trapped ions [75] and quantum many-body physics [76].
Future developments of the method could look at extensions to fermionic
environments [62, 63, 64], as well as making improvements to its performance
based on the accuracy of the fitting algorithm.

Acknowledgements.— This work is based upon research supported by the South
African Research Chair Initiative, Grant No. 64812 of the Department of
Science and Innovation and the National Research Foundation of the Republic of
South Africa. Support from the NICIS (National Integrated Cyber Infrastructure
System) e-research grant QICSI7 is kindly acknowledged.

## References

  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, New York, 2002).
  * de Vega and Alonso [2017] I. de Vega and D. Alonso, Rev. Mod. Phys. 89, (2017).
  * Breuer  _et al._ [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Rev. Mod. Phys. 88, (2016).
  * Li  _et al._ [2018] L. Li, M. J. W. Hall, and H. M. Wiseman, Phys. Rep. 759, 1 (2018).
  * Wolf  _et al._ [2008] M. M. Wolf, J. Eisert, T. S. Cubitt, and J. I. Cirac, Physical Review Letters 101, 150402 (2008).
  * Rivas  _et al._ [2010] Á. Rivas, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050403 (2010).
  * Rivas  _et al._ [2014] Á. Rivas, S. F. Huelga, and M. B. Plenio, Rep. Prog. Phys. 77, 094001 (2014).
  * Breuer  _et al._ [2009] H.-P. Breuer, E.-M. Laine, and J. Piilo, Phys. Rev. Lett. 103, 210401 (2009).
  * Lorenzo  _et al._ [2013] S. Lorenzo, F. Plastina, and M. Paternostro, Physical Review A 88, 020102(R) (2013).
  * Chruściński and Maniscalco [2014] D. Chruściński and S. Maniscalco, Phys. Rev. Lett. 112, 120404 (2014).
  * Ishizaki and Fleming [2012] A. Ishizaki and G. R. Fleming, Annu. Rev. Condes. Matter Phys. 3, 333 (2012).
  * Nalbach  _et al._ [2011] P. Nalbach, D. Braun, and M. Thorwart, Phys. Rev. E 84, 041926 (2011).
  * Chin  _et al._ [2013] A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio, Nat. Phys. 9, 113 (2013).
  * John and Wang [1990] S. John and J. Wang, Phys. Rev. Lett. 64, 2418 (1990).
  * John and Quang [1994] S. John and T. Quang, Phys. Rev. A 50, 1764 (1994).
  * Hoeppe  _et al._ [2012] U. Hoeppe, C. Wolff, J. Küchenmeister, J. Niegemann, M. Drescher, H. Benner, and K. Busch, Phys. Rev. Lett. 108, 043603 (2012).
  * Chin  _et al._ [2012] A. W. Chin, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 109, 233601 (2012).
  * Wu and Shi [2020] W. Wu and C. Shi, Phys. Rev. A 102, 032607 (2020).
  * Bylicka  _et al._ [2014] B. Bylicka, D. Chruściński, and S. Maniscalco, Sci. Rep. 4, (2014).
  * Huelga  _et al._ [2012] S. F. Huelga, Á. Rivas, and M. B. Plenio, Phys. Rev. Lett. 108, 160402 (2012).
  * Li  _et al._ [2020] C.-F. Li, G.-C. Guo, and J. Piilo, EPL (Europhysics Letters) 128, 30001 (2020).
  * Tanimura and Kubo [1989] Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989).
  * Ishizaki and Tanimura [2005] A. Ishizaki and Y. Tanimura, Journal of the Physical Society of Japan 74, 3131 (2005).
  * Tanimura [2020] Y. Tanimura, J. Chem. Phys. 153, 020901 (2020).
  * Makri and Makarov [1995a] N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4600 (1995a).
  * Makri and Makarov [1995b] N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4611 (1995b).
  * Rosenbach  _et al._ [2016] R. Rosenbach, J. Cerrillo, S. F. Huelga, J. Cao, and M. B. Plenio, New J. Phys. 18, 023035 (2016).
  * Strathearn  _et al._ [2018] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, Nat. Commun. 9, 3322 (2018).
  * Gribben  _et al._ [2020] D. Gribben, A. Strathearn, J. Iles-Smith, D. Kilda, A. Nazir, B. W. Lovett, and P. Kirton, Phys. Rev. Res. 2, 013265 (2020).
  * Hu  _et al._ [1992] B. L. Hu, J. P. Paz, and Y. Zhang, Phys. Rev. D 45, 2843 (1992).
  * Zhang  _et al._ [2012] W.-M. Zhang, P.-Y. Lo, H.-N. Xiong, M. W.-Y. Tu, and F. Nori, Phys. Rev. Lett. 109, 170402 (2012).
  * Chin  _et al._ [2010] A. W. Chin, Á. Rivas, S. F. Huelga, and M. B. Plenio, J. Math. Phys. 51, 092109 (2010).
  * Tamascelli  _et al._ [2019] D. Tamascelli, A. Smirne, J. Lim, S. Huelga, and M. Plenio, Phys. Rev. Lett. 123, 090402 (2019).
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Nakajima [1958] S. Nakajima, Prog. Theor. Phys. 20, 948 (1958).
  * Zwanzig [1960] R. Zwanzig, J. Chem. Phys. 33, 1338 (1960).
  * Shibata  _et al._ [1977] F. Shibata, Y. Takahashi, and N. Hashitsume, J. Stat. Phys. 17, 171 (1977).
  * Redfield [1965] A. Redfield, in  _Advances in Magnetic Resonance_ (Elsevier, 1965).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
  * Lindblad [1976] G. Lindblad, Comm. Math. Phys. 48, 119 (1976).
  * Lambert  _et al._ [2019] N. Lambert, S. Ahmed, M. Cirio, and F. Nori, Nat. Commun. 10, 3721 (2019).
  * Pleasance  _et al._ [2020] G. Pleasance, B. M. Garraway, and F. Petruccione, Phys. Rev. Res. 2, 043058 (2020).
  * Tamascelli  _et al._ [2018] D. Tamascelli, A. Smirne, S. Huelga, and M. Plenio, Phys. Rev. Lett. 120, 030402 (2018).
  * Mascherpa  _et al._ [2020] F. Mascherpa, A. Smirne, A. D. Somoza, P. Fernández-Acebal, S. Donadi, D. Tamascelli, S. F. Huelga, and M. B. Plenio, Phys. Rev. A 101, 052108 (2020).
  * Garraway [1997] B. M. Garraway, Phys. Rev. A 55, 2290 (1997).
  * Dalton  _et al._ [2001] B. J. Dalton, S. M. Barnett, and B. M. Garraway, Phys. Rev. A 64, 053813 (2001).
  * Mazzola  _et al._ [2009] L. Mazzola, S. Maniscalco, J. Piilo, K.-A. Suominen, and B. M. Garraway, Phys. Rev. A 80, 012104 (2009).
  * Dalton and Garraway [2003] B. J. Dalton and B. M. Garraway, Phys. Rev. A 68, 033809 (2003).
  * Garraway and Dalton [2006] B. M. Garraway and B. J. Dalton, J. Phys. B: At., Mol. Opt. Phys. 39, S767 (2006).
  * Iles-Smith  _et al._ [2014] J. Iles-Smith, N. Lambert, and A. Nazir, Phys. Rev. A 90, 032114 (2014).
  * Iles-Smith  _et al._ [2016] J. Iles-Smith, A. G. Dijkstra, N. Lambert, and A. Nazir, J. Chem. Phys. 144, 044110 (2016).
  * Diósi [2012] L. Diósi, Phys. Rev. A 85, 034101 (2012).
  * Strasberg  _et al._ [2016] P. Strasberg, G. Schaller, N. Lambert, and T. Brandes, New J. Phys. 18, 073007 (2016).
  * Imamoglu [1994] A. Imamoglu, Phys. Rev. A 50, 3650 (1994).
  * Stenius and Imamoglu [1996] P. Stenius and A. Imamoglu, Quantum Semiclass. Opt. 8, 283 (1996).
  * Breuer [2004] H.-P. Breuer, Phys. Rev. A 70, 012106 (2004).
  * Woods  _et al._ [2014] M. P. Woods, R. Groux, A. W. Chin, S. F. Huelga, and M. B. Plenio, J. Math. Phys. 55, 032101 (2014).
  * Martinazzo  _et al._ [2011] R. Martinazzo, B. Vacchini, K. H. Hughes, and I. Burghardt, J. Chem. Phys. 134, 011101 (2011).
  * Burghardt  _et al._ [2012] I. Burghardt, R. Martinazzo, and K. H. Hughes, J. Chem. Phys. 137, 144107 (2012).
  * Hughes  _et al._ [2009a] K. H. Hughes, C. D. Christ, and I. Burghardt, J. Chem. Phys. 131, 024109 (2009a).
  * Hughes  _et al._ [2009b] K. H. Hughes, C. D. Christ, and I. Burghardt, J. Chem. Phys. 131, 124108 (2009b).
  * Arrigoni  _et al._ [2013] E. Arrigoni, M. Knap, and W. von der Linden, Phys. Rev. Lett. 110, 086403 (2013).
  * Dorda  _et al._ [2014] A. Dorda, M. Nuss, W. von der Linden, and E. Arrigoni, Phys. Rev. B 89, 165105 (2014).
  * Chen  _et al._ [2019] F. Chen, E. Arrigoni, and M. Galperin, New J. Phys. 21, 123035 (2019).
  * Dalibard  _et al._ [1992] J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
  * Gisin and Percival [1992] N. Gisin and I. C. Percival, J. Phys. A: Math. Gen. 25, 5677 (1992).
  * Plenio and Knight [1998] M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Note [1] See Supplemental Material at [URL will be inserted by publisher] for further details on the proof.
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * Garg  _et al._ [1985] A. Garg, J. N. Onuchic, and V. Ambegaokar, J. Chem. Phys. 83, 4491 (1985).
  * Johansson  _et al._ [2013] J. Johansson, P. Nation, and F. Nori, Comput. Phys. Commun. 184, 1234 (2013).
  * Lambert  _et al._ [2020a] N. Lambert, T. Raheja, S. Ahmed, A. Pitchford, and F. Nori, arXiv:2010.10806v1 (2020a).
  * Lambert  _et al._ [2020b] N. Lambert, T. Raheja, S. Ahmed, A. Pitchford, and F. Nori, https://github.com/tehruhn/bofin (2020b).
  * Newman  _et al._ [2017] D. Newman, F. Mintert, and A. Nazir, Phys. Rev. E 95, 032139 (2017).
  * Lemmer  _et al._ [2018] A. Lemmer, C. Cormick, D. Tamascelli, T. Schaetz, S. F. Huelga, and M. B. Plenio, New J. Phys. 20, 073002 (2018).
  * Daley [2014] A. J. Daley, Adv. Phys. 63, 77 (2014).
  * Gardiner and Zoller [2005] C. W. Gardiner and P. Zoller,  _Quantum Noise_ (Springer, Berlin, 2005).
  * Carmichael [1999] H. J. Carmichael,  _Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations_ (Physics and Astronomy Online Library Springer, 1999).

## Supplemental material

In this Supplemental Material we prove the general equivalence between the
open system dynamics derived in both the physical and auxiliary models
discussed in the main text. To do this, we specially modify the proof in [43]
to account for the non-Hermitian form of coupling between the open system and
pseudomodes in line with Refs. [42, 41]. In particular, we note that while our
strategy is closely aligned to one presented by Lambert et. al. in [41], their
result was only formally derived for a specific form of SD and zero-
temperature environment. Here we not only extend this result to arbitrary
Gaussian states (including e.g., thermal states), but also to a generalized
class of SDs satisfying the analyticity constraints above.

###  .1 Auxiliary and extended auxiliary models

First, we are going to prove that the reduced dynamics stemming from the
master equation (7) in the main text is equivalent to the reduced system
dynamics of an extended auxiliary model $S+M+R$, where each of the pseudomodes
$M$ is coupled to a local Markovian reservoir. The total Hamiltonian of this
extended model reads

$H^{\prime}=H_{S}+H_{E^{\prime}}+H^{\prime}_{I}$ (17)

where $H_{E^{\prime}}=H_{M}+H_{R}+V_{MR}$,

$\displaystyle H_{R}=\sum_{l}\int^{\infty}_{-\infty}d\omega\,\omega
a^{\dagger}_{Rl}(\omega)a_{Rl}(\omega),$ (18) $\displaystyle
V_{MR}=\sum_{l}\int^{\infty}_{-\infty}\sqrt{\frac{\lambda_{l}}{\pi}}d\omega\bigg{(}b^{\dagger}_{l}a_{Rl}(\omega)+\text{h.c.}\bigg{)},$
(19)

with $a_{Rl}(\omega)$ ($a^{\dagger}_{Rl}(\omega)$) the annihilation (creation)
operator for excitation of frequency $\omega$ in the $l$-reservoir of $R$, and
$[a_{Rl}(\omega),a^{\dagger}_{Rl^{\prime}}(\omega^{\prime})]=\delta_{ll^{\prime}}\delta(\omega-\omega^{\prime})$.
Here, $H_{R}$ is the free Hamiltonian of the reservoirs, while $V_{MR}$
describes the total $M+R$ interaction. Note that the frequency independence of
the coupling constants in Eq. (19) ensures the reservoirs each have a
vanishing correlation time in accordance with [77].

For a non-unitary evolution $S+M$ in the auxiliary model,

$\rho(t)=e^{\mathcal{L}t}\Big{[}\rho_{S}(0)\otimes\rho_{M}(0)\Big{]},$ (20)

we now proceed to show the density matrix

$\rho_{SMR}(t)=e^{-iH^{\prime}t}\big{[}\rho(0)\otimes|0\rangle\langle
0|_{R}\big{]}e^{iH^{\prime}t}$ (21)

yields a valid purification of $\rho(t)$ on the extended state space $S+M+R$,
such that $\rho(t)=\text{Tr}_{R}[\rho_{SMR}(t)]$ (here,
$|0\rangle_{R}=\bigotimes^{N}_{l=1}|0_{l}\rangle$ denotes the collective
vacuum state of the reservoirs). In particular, it is worth noting that since
the Hamiltonian $H^{\prime}$ is in principle non-Hermitian, Eq. (21) defines a
pseudo-unitary evolution which neither preserves Hermiticity or positivity of
the initial state $\rho_{SMR}(0)$. The ‘pseudo’ prefix refers to the fact that
the righthand side of (21) is acted on by the inverse of the time evolution
operator $\exp(-iH^{\prime}t)$, rather than its Hermitian conjugate, as might
otherwise be employed for a time evolution involving a non-Hermitian
generator.

Using Eq. (21), the Heisenberg equation of motion for an arbitrary operator
$O=O^{\dagger}$ of the enlarged system is given as
($O(t)=e^{iH^{\prime}t}Oe^{-iH^{\prime}t}$)

$\displaystyle\frac{d}{dt}O(t)=-i[O(t),H_{0}(t)]$
$\displaystyle\quad-i\sum_{l}\int
d\omega\sqrt{\frac{\lambda_{l}}{\pi}}\Big{\\{}[O(t),b^{\dagger}_{l}(t)]a_{Rl}(\omega,t)+a^{\dagger}_{Rl}(\omega,t)[O(t),b_{l}(t)]\Big{\\}}.$
(22)

Eliminating the dependence on the reservoir operators $a_{Rl}(\omega,t)$ and
$a^{\dagger}_{Rl}(\omega,t)$ subsequently yields the quantum Langevin equation

$\displaystyle\frac{d}{dt}O(t)$
$\displaystyle=-i[O(t),H_{0}(t)]+\sum_{l}\lambda_{l}\left\\{b^{\dagger}_{l}(t)[O(t),b_{l}(t)]-[O(t),b^{\dagger}_{l}(t)]b_{l}(t)\right\\}$
$\displaystyle-i\sum_{l}\sqrt{2\lambda_{l}}\left\\{[O(t),b^{\dagger}_{l}(t)]a^{(in)}_{Rl}(t)+a^{(in)\dagger}_{Rl}(t)[O(t),b_{l}(t)]\right\\},$
(23)

where

$a^{(in)}_{Rl}(t)=\frac{1}{\sqrt{2\pi}}\int^{\infty}_{-\infty}d\omega\,e^{-i\omega
t}a_{Rl}(\omega).$ (24)

obeys the quantum white noise relation [77]

$[a^{(in)}_{Rl}(t),a^{(in)\dagger}_{Rl^{\prime}}(s)]=\delta_{ll^{\prime}}\delta(t-s).$
(25)

Next we can write down the expectation value $\langle
O(t)\rangle=\text{Tr}\big{[}O(t)\rho_{SMR}(0)\big{]}$ of Eq. (.1) for an
initially factorized density matrix
$\rho_{SMR}(0)=\rho(0)\otimes|0\rangle\langle 0|_{R}$. In this case, the
second line of Eq. (.1) vanishes from

$a^{(in)}_{Rl}(t)\rho_{R}(0)=0=\rho_{R}(0)a^{(in)\dagger}_{Rl}(t),$ (26)

and so

$\displaystyle\frac{d}{dt}\langle
O(t)\rangle=-i\text{Tr}\Big{[}\mathinner{[{O(t)},{H_{0}(t)}]}\rho_{SMR}(0)\Big{]}$
$\displaystyle\quad+\sum_{l}\lambda_{l}\,\text{Tr}\left[\left\\{b^{\dagger}_{l}(t)[O(t),b_{l}(t)]-[O(t),b^{\dagger}_{l}(t)]b_{l}(t)\right\\}\rho_{SMR}(0)\right].$
(27)

Because Eq. (.1) only contains operators acting on the enlarged system Hilbert
space, we are then able to use the cyclic property of the trace

$\langle
O(t)\rangle\equiv\text{Tr}[O(t)\rho_{SMR}(0)]=\text{Tr}[O\rho^{\prime}(t)]$
(28)

to obtain

$\frac{d}{dt}\rho^{\prime}(t)=-i[H_{0},\rho^{\prime}(t)]+2\sum_{l}\lambda_{l}\left(b_{l}\rho^{\prime}(t)b^{\dagger}_{l}-\frac{1}{2}\big{\\{}b^{\dagger}_{l}b_{l},\rho^{\prime}(t)\big{\\}}\right),$
(29)

which, from being identical to the master equation (7) in the main text, and
with $\rho^{\prime}(t)$ subject to the same initial condition as $\rho(t)$
[Eq. (20)], implies $\rho(t)=\rho^{\prime}(t)$. Finally, tracing out the
pseudomode degrees of freedom from both sides of this expression leads to

$\rho^{\prime}_{S}(t)=\text{Tr}_{MR}[\rho_{SMR}(t)].$ (30)

###  .2 Physical and extended auxiliary models

We now move on to prove that the reduced system evolution in the physical
model $S+E$ is equivalent to reduced dynamics of $S$ in the extended auxiliary
model. As we may recall from Ref. [43] (and as a direct consequence of Wick’s
theorem), a sufficient condition to guarantee the equivalence of the open
system dynamics evolving under two physically different Gaussian environments
is for the moments of their free bath coupling operators to match to second-
order. More recently this condition was shown to hold even if coupling between
the system and one of the environments is non-Hermitian [41, 42], so long as
the total density matrix of the system-plus-auxiliary environment obeys a
pseudo-unitary evolution as per Eq. (21). As such, we are able to employ the
same arguments here to prove the reduced system density matrices $\rho_{S}(t)$
and $\text{Tr}_{MR}[\rho_{SMR}(t)]$ share an equal time dependence. To this
end, let us assume the auxiliary environment $M+R$ is initially in the vacuum
state

$\rho_{MR}(0)=|0\rangle\langle 0|_{M}\otimes|0\rangle\langle 0|_{R},$ (31)

such that

$[H_{E^{\prime}},\rho_{MR}(0)]=0.$ (32)

The reduced system dynamics of the extended auxiliary model will then be fully
determined by the second-order correlation function (assuming vanishing first
order moments $\langle B^{\prime}(t)\rangle=0$)

$C^{\prime}(\tau)=\text{Tr}\Big{[}B^{\prime}(\tau)B^{\prime}(0)\rho_{MR}(0)\Big{]},\qquad\tau\geq
0,$ (33)

where
$B^{\prime}(\tau)=e^{iH_{E^{\prime}}\tau}B^{\prime}e^{-iH_{E^{\prime}}\tau}$.

By now inserting the explicit expression for $B^{\prime}(\tau)$ into the
above, we obtain

$C^{\prime}(\tau)=-i\sum^{N}_{l,m=1}\sqrt{r_{l}r_{m}}\langle
b_{l}(\tau)b^{\dagger}_{m}(0)\rangle,$ (34)

given that all second-order moments involving normal ordered products of
operators $b_{l}$ and $b^{\dagger}_{l}$ (e.g., $\langle
b^{\dagger}_{l}(\tau)b_{m}(0)\rangle$) are zero for the initial state (31).
The time dependence of the correlation functions $\langle
b_{l}(\tau)b^{\dagger}_{m}(0)\rangle$ is contained in free evolution of the
pseudomodes, whose Heisenberg equation of motion reads

$\frac{d}{d\tau}b_{l}(\tau)=-i[b_{l}(\tau),H_{E^{\prime}}]=-i\xi_{l}b_{l}(\tau)-i\sqrt{\frac{\lambda_{l}}{\pi}}a_{Rl}(\omega,t).$
(35)

After eliminating the reservoir variables $a_{Rl}(\omega,\tau)$ this equation
becomes

$\frac{d}{d\tau}b_{l}(\tau)=-iz_{l}b_{l}(\tau)-i\sqrt{2\lambda_{l}}a^{(in)}_{Rl}(t),$
(36)

which can be formally solved to yield

$b_{l}(\tau)=e^{-iz_{l}\tau}b_{l}(0)-i\sqrt{2\lambda_{l}}\int^{\tau}_{0}ds\,e^{-iz_{l}(\tau-s)}a^{(in)}_{Rl}(s).$
(37)

In turn, we may substitute $b_{l}(\tau)$ into the above and use the initial
condition $\langle b_{l}(0)b^{\dagger}_{m}(0)\rangle=\delta_{lm}$, together
with $\Big{\langle}b_{m}(0)a^{(in)}_{Rl}(s)\Big{\rangle}=0$, to obtain

$C^{\prime}(\tau)=-i\sum^{N}_{l=1}r_{l}e^{-iz_{l}\tau},\quad\tau\geq 0.$ (38)

Since the correlation function $C^{\prime}(\tau)$ is therefore identical to
its counterpart $C(\tau)$ in the physical model, we have from the Gaussian
property of $\rho_{E}(0)$ and $\rho_{MR}(0)$ that

$\rho_{S}(t)=\text{Tr}_{MR}[\rho_{SMR}(t)].$ (39)

Finally, by direct comparison of this result with Eq. (30), it is possible to
establish the exact equivalence of the reduced dynamics occurring in both the
physical and auxiliary models,

$\rho_{S}(t)=\rho^{\prime}_{S}(t),$ (40)

as stated in Eq. (10) of the main text.

As an additional remark, it is worth noting that the same result may equally
be derived by utilizing an alternative definition of the two-time correlation
function (33) that is consistent with the quantum regression theorem [77, 78].
Since the master equation (29) is Markovian, by setting $H_{S}=0$ and
$H^{\prime}_{I}=0$, it is possible to show that $C^{\prime}(\tau)$ can be re-
written in the form (see the Supplemental Material of [43]),

$C^{\prime}(\tau)=\text{Tr}\Big{\\{}B^{\prime}e^{\mathcal{L}_{M}\tau}\Big{[}B^{\prime}\rho_{M}(0)\Big{]}\Big{\\}},\qquad\tau\geq
0,$ (41)

where $\mathcal{L}_{M}\rho=-i[H_{M},\rho]+\mathcal{D}\rho$. From here the
proof follows the same steps outlined above: using that
$\rho_{M}(0)=|0\rangle\langle 0|_{M}$, the correlation function
$C^{\prime}(\tau)$ can again be reduced to the same form as in Eq. (34) where
the time dependence of $b_{l}(\tau)$ is determined by
$b_{l}(\tau)=\exp\big{(}\mathcal{L}^{\dagger}_{M}\tau\big{)}[b_{l}(0)]$ (here,
$\mathcal{L}^{\dagger}_{M}$ is defined via the duality relation
$\text{Tr}\big{[}A\mathcal{L}_{M}(B)\big{]}=\text{Tr}\big{[}\mathcal{L}^{\dagger}_{M}(A)B\big{]}$).
Since $b_{l}$ obeys the eigenoperator relation
$\mathcal{L}_{M}b_{l}=-iz_{l}b_{l}$, it is then straightforward to demonstrate
that the above definition recovers Eq. (38), implying
$\rho^{\prime}_{S}(t)=\rho_{S}(t)$.


# Discovering Behavioral Predispositions in Data to Improve Human Activity
Recognition

Maximilian Popko
University of Rostock, Germany
<EMAIL_ADDRESS>Sebastian Bader
University of Rostock, Germany
<EMAIL_ADDRESS>Stefan Lüdtke
University of Mannheim, Germany
<EMAIL_ADDRESS>
Thomas Kirste
University of Rostock, Germany
<EMAIL_ADDRESS>
###### Abstract

The automatic, sensor-based assessment of challenging behavior of persons with
dementia is an important task to support the selection of interventions.
However, predicting behaviors like apathy and agitation is challenging due to
the large inter- and intra-patient variability. Goal of this paper is to
improve the recognition performance by making use of the observation that
patients tend to show specific behaviors at certain times of the day or week.
We propose to identify such segments of similar behavior via clustering the
distributions of annotations of the time segments. All time segments within a
cluster then consist of similar behaviors and thus indicate a behavioral
predisposition (BPD). We utilize BPDs by training a classifier for each BPD.
Empirically, we demonstrate that when the BPD per time segment is known,
activity recognition performance can be substantially improved.

_Keywords_ Human activity recognition $\cdot$ Wearable sensors $\cdot$ Machine
learning $\cdot$ Clustering

##  1 Introduction

The recognition of human behavior is important for many domains, like sports,
the analysis of manual work processes, or as a component for situation-aware
assistants. As another example, sensor-based Human Activity Recognition (HAR)
can also be used for the automatic assessment of behavioral symptoms of people
with dementia, like apathy or agitation. In contrast to questionnaires like
the Cohen-Mansfield Agitation Inventory [1], HAR allows the real-time,
objective assessment of symptoms, supporting caregivers in selecting
appropriate interventions. However, accurate HAR is challenging due to the
high inter- and intra-subject variability of movement. Specifically, the
automatic recognition of challenging behavior from wearable sensor data has
been of limited success so far [2].

Multiple approaches have been proposed to improve HAR accuracy, by making use
of additional external knowledge [3, 4, 5, 6, 7, 8]. In this paper, we
investigate how HAR performance can be improved by making use of the fact that
there is a tendency of specific subjects to show specific behaviors at certain
times of the day or week. For example, sundowning is the increase in
restlessness of some people with dementia in the evening [9]. The knowledge of
such behavioral predispositions (BPDs) could serve as a useful prior to
increase HAR accuracy.

Figure 1: Overview of the concept. We calculate the distribution of
annotations for a given segment of the day (histograms). We then cluster these
distributions. Next, we use the cluster label (BPD) to select the classifier
(color indicates the corresponding classifier). The classifier predicts the
behavior from the features of the motion sensor data.

Previous work relied on external knowledge to define BPDs. For example, Lüdtke
et al. [7] partitioned warehouse order picking into distinct process steps,
the knowledge of which can inform HAR. However, in the context of dementia
care, there are no established BPDs that are externally given. The
contribution of this paper is to identify such BPDs in a data-driven way: We
calculate the distribution of annotations in time segments of a day, and then
cluster these histograms. This results in clusters of time frames in which the
subjects behaved similarly. Each cluster forms a BPD. Given the BPD, we then
train a classifier specifically for data recorded within the BPD, i.e., there
is one classifier for each cluster. This classifier is thus biased towards the
more frequent classes in the BPD. An overview of this process can be found in
figure 1. In summary, this paper provides the following contributions:

  * •

A clustering-based method for the data-driven construction of behavioral
predispositions (BPDs)

  * •

An approach for utilizing these BPDs in HAR

  * •

An evaluation of this approach for the sensor-based recognition of behavioral
symptom of persons with dementia

Empirically, our approach leads to a substantial increase of HAR performance,
when assuming that the BPD for each time segment is given.

##  2 Data and Observations

This section presents the insideDem framework, as well as the data set. We
make observations about the distribution of annotations given time frames,
which will then lead to the construction of the BPDs and the utilization of
them in the methodology section 3.

###  2.1 insideDem Project

Part of the insideDem framework was the recording of the behaviours of
dementia diagnozed patients [10]. In two nursing homes, real-time observations
and sensor data of a total of 17 residents (11 women and 6 men) with moderate
to severe dementia were recorded. The behaviour of persons in one nursing home
was also recorded on video. The age range was 73 to 94 years. All persons were
under psycho-pharmacological treatment. Dementia medication was taken by 8
residents of the nursing homes. The behaviour of the persons was recorded in
real time by an expert in the sitting room of the nursing home. This room was
accessible to all residents and their visitors as well as the staff of the
nursing home. The behaviour was annotated every five minutes following the
annotation scheme which was developed for the insideDem project. This scheme
is based on the neuropsychiatric inventory [11] and the Cohens-Mansfield
agitation inventory [1]. It consists of seven different behaviours: apathy,
general restlessness, mannerisms, pacing, aggression, trying to get to a
different place, normal behaviour.

During the day, between 8 a.m. and 6 p.m., the subjects wore two wristbands,
which were attached to the dominant hand and the ankle. During the night, a
band was worn on the ankle. The sensors recorded accelerations and rotations
as well as the sound level of the environment, the light conditions and the
sound pressure. In this work, we only use the acceleration sensor data.

###  2.2 Observations

Figure 2 shows the distribution of annotations for each subject. The behaviour
of the eight patients were observed over the course of 29 days and accumulate
around 17,000 annotated five-minute intervals. For every subject the the
behaviours occur in different portions. As example, the patient X111 showed 9%
apathy where as X126 showed nearly 50% apathy during the time span of 29 days.
However, for X111, pacing was mostly annotated but X126 nearly showed no
pacing.

Figure 2: Distribution of annotations per person over all days of recordings.

In the literature events for triggering agitated behaviors were observed. The
sunset phenomenon as described in [9] may be such a trigger for one BPD. At
sunset people with dementia often show more agitated behavior than at other
times of the day. Also, the fact of boredom can trigger agitated behavior. If
there is nothing planned for an afternoon, dementia diagnosed persons might
show agitated behavior because of the boredom [12]. As a result, it is
reasonable to look into noticeable time frames of the day which show one
person behaving agitated more often.

The figure 3 (a) now shows the average occurrences of the annotations for one-
hour segments of the day. The represented subject behaves in the morning more
apathetic than in the later course of the day. Also at noon, the patient shows
more normal behavior. Pacing is mostly happening in the afternoon. Figure 3
(b) shows the variances of the annotations between each day. This plot
represents the distribution of the proportion each behavior has on the time
segment every day. Between 10 and 11 a.m. we saw the patient behaves
apathetically to 46% in plot (a). However, in plot (b) we see that the apathy
part can differ from 0% to 100% with equal probability from day to day.

Figure 3: Percentual distribution of annotations for one subject and one-hour
intervals of the day. (a) shows the average of the annotation and (b)
describes the distribution (variance).

##  3 Methods

Behavioural predispositions (BPD) describe the tendency to show specific
behaviours. In this section, we propose a method to form such BPDs from the
annotations of behaviors. In the following, we present a simple model that can
use these newly found BPDs to simplify the classification.

###  3.1 Extracting behavioural predispositions

BPD describe the urge to show specific kinds of behaviors. Knowing such would
intuitively mean that the prediction of the current actual behavior becomes
easier since it provides an informative prior on the behavior of the subject.
Unfortunately an expert model describing BPD for patients with dementia does
not exist and thus the aim is to derive them from the data of the insideDem
project.

####  3.1.1 Time-based BPD

Following the observations of figure 3 (a), one can argue to use a similar
approach to [8] or define that every BPD is one segment of the day. This means
from 10 to 11 a.m. the subject always behaves similarly and by knowing the
current time, we can use this information to improve the classification.
However, this approach has the drawback of high variances of the annotations
between each day which lead to BPDs that do not carry a lot of information
about the actual behavior.

####  3.1.2 Clustering-based BPD

To overcome the problem of high variances, we propose to find similar
distributions globally over all days and time segments instead of grouping
them by the time of the day. Let $p(A_{t})$ be the probability distribution of
the annotation for the day and time segment $t$. Note that $t$ describes a
specific time frame e.g. 06–15 10:00 to 06–15 11:00. This categorical
distribution can be obtained by counting the relative frequencies of the
annotations in the segment. This counted distribution is called a histogram.
The BPD represents the tendency to behave in a certain way. Given a set of
histograms, the task is to divide this set into groups. These groups should
contain similar histograms. With the help of the k-means clustering algorithm,
clusters of histograms can be generated which are grouped according to their
spatial similarity (Euclidean distance). Each of these clusters describes a
behavioral disposition. The histograms of the cluster $D=d$ therefore
represent the distribution $p(A\,|\,D=d)$.

###  3.2 Classification of Sensor Data

In the scope of this work, we want to show what impact the complete knowledge
of the clustered BPD can have on the classification. Therefore, we assume the
BPD, i.e. the cluster, as given. Note that this would not be the case in a
real world environment. Given a BPD $d$, we then have a prior knowledge of the
annotation $p(A\,|\,d)$. To make use of this knowledge, we train a different
classifier for each BPD. Thus, each classifier is trained with biased data and
is therefore biased to more likely predict the classes that are most presently
in this BPD. The classifier can be seen as expert for its BPD. The model $m$
maps sensor data $s$ and the BPD $d$ to an annotation $a$ by choosing the
classifier $c_{d}$ that was specifically trained for $d$. The classifier then
maps $s$ to $a$:

$m(s,d)=c_{d}(s)=a\,.$

##  4 Experimental Evaluation

The overall goal of the experiments are the evaluation of the impact of the
BPD on the HAR performance. In the following, we describe the experimental
design in more detail.

###  4.1 Experimental Design

Table 1: Factors and levels of the experimental design Factors |  Levels |  Description
---|---|---
Clustering |  K-means, time-based |  Methods for generating the BPDs
#-clusters |  1, …, 20 |  Number of BPDs
Segment lengths |  0, …, 120 minutes (k-means) |  Length of the segment over which the histogram is calculated, only for k-means approach
Classifiers |  Majority, naive Bayes, SVM |  Classifiers to predict the annotation from sensor data
Subjects |  X110, X111, X113, X114, X121, X122, X124, X126 |  One model for each subject

To evaluate the impact of the knowledge of the BPDs on HAR performance, we use
a factorial design. The factors and levels are presented in table 1. For each
subject, a model with a different clustering method, number of clusters (i.e.
BPDs), segment length, and classifier is trained. We randomly select 70% of
the data for training and 30% for testing the model. Every combination of the
levels is evaluated 10 times. Different clusters and segment lengths are used
to test the expressiveness of the BPD. The segment length plays an important
role for the k-means clustering strategy as the length is used to define the
size over which the histogram is calculated and then clustered. At a segment
length of 30 minutes, we obtain around 400 histograms per subject. We then
cluster the histograms of a single subject, resulting in BPDs specifically for
this subject. High number of clusters and small segment lengths yield finer
BPDs. Note that one cluster means we have no further information about the BPD
and thus gives us a baseline of performance. The time-based clustering method
describes that a BPD only depends on the time. E.g., having 10 clusters and
time-based clustering this method would split the day into 10 parts and for
every part, one classifier is trained. We use three different methods for
classifying the sensor data. The majority method functions as a baseline by
always outputting the most frequent annotation of the BPD. The naive Bayes
(NB) method shows the performance of a simplistic classifier. The NB uses for
every BPD a different prior and observation model. Thirdly, a support vector
machine (SVM) is trained. For preprocessing the accelerator data, features
were calculated following the work of [13]. The original data is sampled in 50
Hz. Feature were calculated over a window size of one minute and an overlap of
50%. Statistical features, spatial relations between the three axis and
frequency features were used. As measure of accuracy, we use $F_{1}$ score.

###  4.2 Results and Discussion

Figure 4: Averaged $F_{1}$ score for varying number of clusters and segment
lengths. For time-based clustering, the segment length is implicitly given
through the number of clusters. For k-means clustering, the segment lengths of
the histograms are represented by the color. The segment length is in minutes.
Rectangles denote similar BPD structures because of the same number of
clusters and similar segment lengths.

Figure 4 shows the average $F_{1}$ scores for each classifier and the
strategies used for generating the BPDs. Firstly, we want to investigate only
the k-means approach. The SVM shows an increase of $F_{1}$ score up to a
number of about 10 clusters. From this point on, only smaller segment sizes
benefit. For the majority classifier, the $F_{1}$ score no longer increases
significantly above a number of 6 clusters. The naive Bayes classifier shows a
nearly linear increase. Thus, for all classifiers, the $F_{1}$ measure
increases with a larger number of clusters, as well as a decreasing segment
length. The fact that the majority classifier benefits from this approach
comes as no surprise. With an increasing number of clusters and finer segment
lengths, the most frequent annotation occupies a higher portion of data within
the cluster. Predicting only the most frequent annotation yields an $F_{1}$
score of around 0.1 whereas dividing the data by the BPD yields a score of
around 0.5 at 20 BPDs and a segment length of 30 minutes. Note that some BPDs
only describe one annotation and thus predicting can achieve 100% accuracy by
always outputting that one same class. According to this observation, the
classification task should become simpler. Given such BPD with only one class,
all classifiers used here can benefit and improve with knowledge of the BPD.
However, the $F_{1}$ score of the majority classifier proposes that the
knowledge of the most frequent annotation in the training data does not give
us further advantages when we have more than 6 clusters. Yet the $F_{1}$ score
of the naive Bayes classifier increases nearly linearly. Also, the SVM has a
still increasing $F_{1}$ score at 20 BPDs, which underlines the fact that
knowing the BPD simplifies the classification task by reducing the complexity
of possible behaviors or even ruling them out. Note that the knowledge
presented by the BPD converges with the BPD describing the exact underlying
histogram, i.e. we have the same number of BPDs (clusters) as histograms.

We now investigate the time-based approach. The time-based approach shows an
increase of $F_{1}$ score as the number of clusters rises. This can be
observed for every classifier evaluated. Although the increase is not high,
the knowledge of the clock time can help the classification. Thus dividing the
day into segments and training a classifier for every segment each, increases
the $F_{1}$ score. Comparing this approach with the k-means approach can be
done by comparing the points where the number of BPDs and the segment length
are approximately equal for both. E.g. splitting the day, where the subjects
were recorded, into 20 segments yields a segment length of around 30 minutes
and 20 BPDs. These points are denoted with rectangles in figure 4. The plot
yields a better performance of the k-means approach for all 4 points. The
inferiority of the time-based approach can be explained through the high
variances of annotations among the days (shown in figure 3). The k-means
approach, on the other hand, describes the histograms by their similarity and
thus achieving less variance within the cluster. Having less variance means
that the classifier has an easier task to solve.

Figure 5: Confusion matrix of the predicted behaviors using 20 BPDs and thus
20 SVMs. A segment length of 30 minutes were used to calculate the histrogams.

Overall, the $F_{1}$ score reaches a maximum of around 0.5. This score is
produced by using 20 BPDs, a segment length of 30 minutes, and an SVM as a
classifier. The same results can be achieved by using the majority method. The
SVM outperforms the majority classifier for larger segment lengths. Figure 5
shows the confusion matrix of the SVMs for the classification with 20 BPDs and
a segment length of 30 minutes. We can see that the classes mannerisms,
apathy, and normal behavior represent the majority of the annotations.
Furthermore, there is great uncertainty between these classes. This can be
explained by the similarity of these behaviors in this scenario. Intuitively,
normal behavior is similar to apathy, as dementia patients in the nursing home
live more secluded. Mannerisms are difficult to recognize when, for example,
they are performed with the non-dominant hand (sensor is worn at the dominant
hand). Another uncertainty is between the annotations pacing, getting to a
different place, and normal behavior. This can be explained by the fact that
these three annotations can be understood as the activity of walking so a
misclassification can occur between these classes. The behavior is also
annotated in five-minute intervals, but the classification takes place in
30-second intervals. This means that the person may have walked at one point
in time, but the interval was annotated differently. Another reason is the
exclusive use of the motion sensor on the wrist.

##  5 Future work

The evaluation showed that the knowledge of the BPD can yield better
classification and thus further investigations are of interest. There is the
question of whether finding the BPD with k-means clustering is the best
choice. Since similar histograms and thus categorical probability
distributions are of interest, clustering methods with entropy-based
distances, such as the Jensen-Shannon divergence, are of interest.
Furthermore, Latent Dirichlet Allocation (LDA) [14] can be investigated as a
method for determining behavioral dispositions. The LDA is a generative
probabilistic model and is mainly used in applications to assign a topic to
documents based on their words. In this scenario, a segment of the day can be
regarded as a document and its words are the annotations in this segment so
that the resulting topics describe the BPD.

Currently, the BPD is only constructed by the annotations. A next step could
be to determine the BPD based on the annotation as well as the sensor data.
This would lead to BPDs that also hold information about the current motion
data and its connection to the annotation. In this work, we created BPDs for a
single person. As a next step, the evaluation of global BPDs across the
patients is needed.

We also saw that the overall accuracy of the classification was low and thus
looking into better options for classification is of interest. In this work,
we only presented the approach of training a classifier for each BPD. However,
this reduces the number of samples for the classifier. Another approach can be
the feeding of the BPD label as input to the classifier. Also, classifiers
that can make use of distributions over BPDs are of interest. Methods like
mixed-effects neural networks [15] can model BPDs as another feature for
activity recognition.

In reality, the BPD is not available and has to be found from sensor data. A
good estimation of the BPD is therefore needed to improve the classification.
With an increasing number of BPDs, the error of classification will also rise.
Thus finding a good number of BPDs and a model to predict them is a crucial
point for practical use. This can be done by using hierarchical models which
first infer the BPD from the sensor data and then make use of the relations
between behavior and BPD. Also, hidden Markov models can be of interest to
exploit possible temporal properties of the BPDs.

##  6 Related work

To improve the classification of sensor data using additional information,
firstly, the knowledge has to be modeled and perhaps has to be generated
beforehand. Depending on the domain, the kind of knowledge can differ and thus
arise different ways of incorporation into human activity recognition. Mostly
the improvement of accuracy is the main purpose. However, the simplification
of the recognition model can also be of interest. The authors of [3] are
aiming to create a model that still maintains the same level of accuracy under
less training data than state-of-the-art methods. They use region connection
calculus (RCC) [16] to describe the spatial relations between objects and
formalize the activities and their causal properties within the domain using
the planning domain definition language (PDDL) [17]. A probabilistic
concurrent constraint automaton (PCCA) [18] is generated from the PDDL model
and reasons about the predicates from the RCC model. The RCC model calculates
these predicates from high-level observations. The authors show that their
approach can yield similar results in comparison to convolutional neural
networks (CNN) but avoid the training data hungry nature of neural networks.
On the other hand, if enough training data and further information about the
domain is available, the work of Xie et al. [4] aims to improve the accuracy
of classification by incorporating prior knowledge in the training process of
a neural network. They achieve this goal by modeling rules about the
practicability of actions using linear temporal logics (LTL) [19]. To overcome
the problem of the multi-interpretable nature of LTL, the author convert them
into deterministic finite automatons such that they can be embedded as
hierarchical graph into the loss calculation of the neural network’s training
process and achieve learning consistent to the LTL rules and better accuracy
at classification. The works by Yordanaova et al. [5] and Rueda et al. [6] use
computational causal behavioral models (CCBM) [20] to model domain knowledge
and reason about the current state of the environment. CCBM also uses PPDL-
like rules to describe the domain-specific properties of the activities. Also,
the relations between user and object as well as durations of activities can
be modeled in CCBM. As an observation model, a decision tree (DT) in [5] and a
convolutional neural network in [6] are used to weight the predicted state
based on the sensor data. In contrast to these works, the prior knowledge does
not need to be embedded directly in the model. Lüdtke et al. [7] demonstrate
how providing external _context information_ impacts classification
performance. Specifically, they show that providing the current process step
of structured activities can improve HAR performance, even when that
information is noisy. Even if no external information is available, the
authors of [8] made the observation that activities feature different
quantities in times of a day. Thus segmenting the day and reasoning about the
activities in one segment yields better accuracy than reasoning about the
activities of the whole day.

##  7 Conclusion

In this paper, we proposed a data-driven method to construct behavioral
predispositions and use them as informative prior to improve the performance
of sensor-based human activity recognition. We presented a clustering-based
approach that can determine time segments that display similar behaviors of a
dementia-diagnosed person. We found such groups by first calculating the
distribution of the annotations in the time segment and then clustering these
distributions. The resulting clusters explain behavioral predispositions
(BPDs). With complete knowledge of the current BPD, we demonstrated
empirically that a substantial increase of the $F_{1}$ score can be achieved
by using classifiers that were specifically trained for each BPD. For
practical usage, the investigation of predicting the BPDs from sensor data is
of high interest, as the BPD is not given a priori and instead must be
inferred.

Furthermore, the presented approach is not only applicable to predicting the
behavior of persons with dementia, but also to any continuous HAR task.
Whenever the prior for occurrence of activities varies over time, the method
described in this paper offers a viable option to improve HAR performance.

## References

  * Cohen-Mansfield [1996] Jiska Cohen-Mansfield.  Assessment of agitation.  _International Psychogeriatrics_ , 8(2):233–245, 1996.  doi:10.1017/S104161029600261X.
  * Goerss et al. [2020] Doreen Goerss, Albert Hein, Sebastian Bader, Margareta Halek, Sven Kernebeck, Andreas Kutschke, Christina Heine, Frank Krueger, Thomas Kirste, and Stefan Teipel.  Automated sensor-based detection of challenging behaviors in advanced stages of dementia in nursing homes.  _Alzheimer’s & Dementia_, 16(4):672–680, 2020\.  doi:https://doi.org/10.1016/j.jalz.2019.08.193.  URL https://alz-journals.onlinelibrary.wiley.com/doi/abs/10.1016/j.jalz.2019.08.193.
  * Lee et al. [2019] Sang Uk Lee, Andreas Hofmann, and Brian Williams.  A Model-Based Human Activity Recognition for Human–Robot Collaboration.  In _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pages 736–743, Macau, China, November 2019. IEEE.  ISBN 978-1-72814-004-9.  doi:10.1109/IROS40897.2019.8967650.  URL https://ieeexplore.ieee.org/document/8967650/.
  * Xie et al. [2021] Yaqi Xie, Fan Zhou, and Harold Soh.  Embedding Symbolic Temporal Knowledge into Deep Sequential Models.  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 4267–4273, Xi’an, China, May 2021. IEEE.  ISBN 978-1-72819-077-8.  doi:10.1109/ICRA48506.2021.9561952.  URL https://ieeexplore.ieee.org/document/9561952/.
  * Yordanova et al. [2019] Kristina Yordanova, Stefan Lüdtke, Samuel Whitehouse, Frank Krüger, Adeline Paiement, Majid Mirmehdi, Ian Craddock, and Thomas Kirste.  Analysing Cooking Behaviour in Home Settings: Towards Health Monitoring.  _Sensors_ , 19(3):646, February 2019.  ISSN 1424-8220.  doi:10.3390/s19030646.  URL http://www.mdpi.com/1424-8220/19/3/646.
  * Rueda et al. [2019] Fernando Moya Rueda, Stefan Lüdtke, Max Schroder, Kristina Yordanova, Thomas Kirste, and Gernot A. Fink.  Combining Symbolic Reasoning and Deep Learning for Human Activity Recognition.  In _2019 IEEE International Conference on Pervasive Computing and Communications Workshops (PerCom Workshops)_ , pages 22–27, Kyoto, Japan, March 2019. IEEE.  ISBN 978-1-5386-9151-9.  doi:10.1109/PERCOMW.2019.8730792.  URL https://ieeexplore.ieee.org/document/8730792/.
  * Lüdtke et al. [2021] Stefan Lüdtke, Fernando Moya Rueda, Waqas Ahmed, Gernot A. Fink, and Thomas Kirste.  Human Activity Recognition using Attribute-Based Neural Networks and Context Information.  In _3rd International Workshop on Deep Learning for Human Activity Recognition_ , 2021.
  * Inoue et al. [2016] Sozo Inoue, Naonori Ueda, Yasunobu Nohara, and Naoki Nakashima.  Recognizing and Understanding Nursing Activities for a Whole Day with a Big Dataset.  _Journal of Information Processing_ , 24(6):853–866, 2016.  ISSN 1882-6652.  doi:10.2197/ipsjjip.24.853.  URL https://www.jstage.jst.go.jp/article/ipsjjip/24/6/24_853/_article.
  * Canevelli et al. [2016] Marco Canevelli, Martina Valletta, Alessandro Trebbastoni, Giuseppe Sarli, Fabrizia D’Antonio, Leonardo Tariciotti, Carlo de Lena, and Giuseppe Bruno.  Sundowning in Dementia: Clinical Relevance, Pathophysiological Determinants, and Therapeutic Approaches.  _Frontiers in Medicine_ , 3:73, December 2016.  ISSN 2296-858X.  doi:10.3389/fmed.2016.00073.  URL http://journal.frontiersin.org/article/10.3389/fmed.2016.00073/full.
  * Teipel et al. [2017] Stefan Teipel, Christina Heine, Albert Hein, Frank Krüger, Andreas Kutschke, Sven Kernebeck, Margareta Halek, Sebastian Bader, and Thomas Kirste.  Multidimensional assessment of challenging behaviors in advanced stages of dementia in nursing homes—The insideDEM framework.  _Alzheimer’s & Dementia: Diagnosis, Assessment & Disease Monitoring_, 8(1):36–44, January 2017.  ISSN 2352-8729, 2352-8729.  doi:10.1016/j.dadm.2017.03.006.  URL https://onlinelibrary.wiley.com/doi/abs/10.1016/j.dadm.2017.03.006.
  * Cummings [1997] Jeffrey L Cummings.  The neuropsychiatric inventory: assessing psychopathology in dementia patients.  _Neurology_ , 48(5 Suppl 6):10S–16S, 1997.
  * Cohen-Mansfield and Werner [1995] Jiska Cohen-Mansfield and Perla Werner.  Environmental influences on agitation: An integrative summary of an observational study.  _American Journal of Alzheimer’s Care and Related Disorders & Research_, 10(1):32–39, January 1995.  ISSN 0895-5336.  doi:10.1177/153331759501000108.  URL http://journals.sagepub.com/doi/10.1177/153331759501000108.
  * Alam et al. [2019] Ridwan Alam, Azziza Bankole, Martha Anderson, and John Lach.  Multiple-Instance Learning for Sparse Behavior Modeling from Wearables: Toward Dementia-Related Agitation Prediction.  In _2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , pages 1330–1333, Berlin, Germany, July 2019. IEEE.  ISBN 978-1-5386-1311-5.  doi:10.1109/EMBC.2019.8856502.  URL https://ieeexplore.ieee.org/document/8856502/.
  * Blei et al. [2003] David M Blei, Andrew Y Ng, and Michael I Jordan.  Latent dirichlet allocation.  _Journal of machine Learning research_ , 3(Jan):993–1022, 2003.
  * Xiong et al. [2019] Yunyang Xiong, Hyunwoo J. Kim, and Vikas Singh.  Mixed Effects Neural Networks (MeNets) With Applications to Gaze Estimation.  In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 7735–7744, Long Beach, CA, USA, June 2019. IEEE.  ISBN 978-1-72813-293-8.  doi:10.1109/CVPR.2019.00793.  URL https://ieeexplore.ieee.org/document/8954429/.
  * Cohn et al. [1997] Anthony G. Cohn, Brandon Bennett, John Gooday, and Nicholas Mark Gotts.  Qualitative Spatial Representation and Reasoning with the Region Connection Calculus.  _GeoInformatica_ , 1(3):275–316, October 1997\.  ISSN 1573-7624.  doi:10.1023/A:1009712514511.  URL https://doi.org/10.1023/A:1009712514511.
  * Fox and Long [2003] M. Fox and D. Long.  PDDL2.1: An Extension to PDDL for Expressing Temporal Planning Domains.  _Journal of Artificial Intelligence Research_ , 20:61–124, December 2003.  ISSN 1076-9757.  doi:10.1613/jair.1129.  URL https://jair.org/index.php/jair/article/view/10352.
  * Williams et al. [2003] B.C. Williams, M.D. Ingham, S.H. Chung, and P.H. Elliott.  Model-based programming of intelligent embedded systems and robotic space explorers.  _Proceedings of the IEEE_ , 91(1):212–237, 2003\.  doi:10.1109/JPROC.2002.805828.
  * Pnueli [1977] Amir Pnueli.  The temporal logic of programs.  In _Proceedings of the 18th Annual Symposium on Foundations of Computer Science_ , SFCS ’77, page 46–57, USA, 1977. IEEE Computer Society.  doi:10.1109/SFCS.1977.32.  URL https://doi.org/10.1109/SFCS.1977.32.
  * Krüger et al. [2012] Frank Krüger, Alexander Steiniger, Sebastian Bader, and Thomas Kirste.  Evaluating the robustness of activity recognition using computational causal behavior models.  In _Proceedings of the 2012 ACM Conference on Ubiquitous Computing_ , UbiComp ’12, page 1066–1074, New York, NY, USA, 2012. Association for Computing Machinery.  ISBN 9781450312240.  doi:10.1145/2370216.2370443.  URL https://doi.org/10.1145/2370216.2370443.


# Minimal dispersion on the cube and the torus

A. Arman  A. E. Litvak

###### Abstract

We improve some upper bounds for minimal dispersion on the cube and torus. Our
new ingredient is an improvement of a probabilistic lemma used to obtain upper
bounds for dispersion in several previous works. Our new lemma combines a
random and non-random choice of points in the cube. This leads to better upper
bounds for the minimal dispersion.

AMS 2010 Classification: primary: 52B55, 52A23; secondary: 68Q25, 65Y20.
Keywords: complexity, dispersion, largest empty box, torus

##  1 Introduction

For a given set of points $P\subset Q_{d}:=[0,1]^{d}$ its dispersion is
defined as the supremum over volumes of axis-parallel boxes in $Q_{d}$ that do
not intersect $P$. Then the minimal dispersion on the cube is defined as the
infimum of the dispersions of all possible subsets $P\subset Q_{d}$ of
cardinality $n$. The dispersion on the torus is defined similarly. This notion
goes back to [15], where a notion from [10] was modified. Often it is more
convenient to work with its inverse function, which, given a positive
$\varepsilon$, measures the smallest positive integer $N=N(\varepsilon,d)$
such that there exists a configuration $X$ of $N$ points in $[0,1]^{d}$ with
the property that any axis-parallel box of volume exceeding $\varepsilon$
contains at least one point. We refer to [12] and references therein for the
history of the question and related references.

In this note we improve some upper bounds for the minimal dispersions on the
cube and on the torus and for its inverse function.

Several recent proofs on upper bounds on $N$ use the following scheme. In the
first step, approximate the set of axis-parallel boxes in $Q_{d}$ by a finite
set $\mathcal{N}$, with the property that every axis-parallel box of volume at
least $\varepsilon$ contains at least one element of $\mathcal{N}$. In the
second step, construct a set $P$ such that each element of $\mathcal{N}$
intersects $P$, in which case we say that $P$ is a piercing set for
$\mathcal{N}$. To ensure that $P$ pierces $\mathcal{N}$, the points in $P$ are
taken randomly, according to the uniform distribution on the cube, and
subsequently a union bound is applied to verify that such random $P$ with high
probability intersects every element of $\mathcal{N}$.

Our main new ingredient is in the improvement of the second step, where a
random choice of points is followed up by a deterministic phase. More
precisely, our new probabilistic Lemma 3.3 will have two phases. In the first
phase we choose a (smaller) set $Q$ of random points and estimate how many
sets in our approximation $\mathcal{N}$ have an empty intersection with $Q$
(in previous proofs one insisted that every set in $\mathcal{N}$ intersects
$Q$). In a second phase we take care of the “empty” sets in $\mathcal{N}$ (the
ones that do not intersect $Q$) by choosing a representative point in each
such set. The set $Q_{0}$ of representative points together with $Q$ forms a
piercing set for $\mathcal{N}$.

The new Lemma 3.3 subsequently leads to better upper bounds on minimal
dispersion, see (4). However in doing the second step we lose randomness, so
our new bound does not hold for a random choice of points. Moreover, it is
known that our bound cannot hold in the random setting, for instance see [8],
where the authors provide lower bounds on the inverse of minimal dispersion
for a random set of points (formula (2) below).

We would also like to mention that the idea of choosing a random set of points
in a first phase and then completing it by a deterministic choice (depending
on a realization of a random set) in the second phase is not new and was used
in literature, see e.g. [14], or [1, Chapter 3].

##  2 Notation, preliminaries and main result

Throughout the paper the following notation is used. For positive integer $d$
define the unit cube $Q_{d}=[0,1]^{d}$. We will use $|\cdot|$ to denote both
volume of a subset of $\mathbb{R}^{d}$ and a cardinality of a finite set.
Letters $C,C_{0},C_{1}$, …, $c$, $c_{0}$, $c_{1}$, … always denote absolute
positive constant (independent of all other parameters and sets), whose values
may change from line to line.

Define a set of all axis-parallel rectangles in $Q_{d}$ by

${\mathcal{R}}_{d}=\left\\{\prod_{i=1}^{d}I_{i}\,\,|\,\,I_{i}=[a_{i},b_{i})\subset[0,1]\right\\}.$

For a set $P\subset Q_{d}$ define its dispersion by

${\rm disp}(P)=\sup\left\\{|B|\,\,|\,\,B\in{\mathcal{R}}_{d},\;B\cap
P=\emptyset\right\\}.$

The minimal dispersion ${\rm disp}^{*}$ is the following function of positive
integers $n,d$,

${\rm disp}^{*}(n,d)=\inf\\{{\rm disp}(P)\,\,|\,\,P\subset Q_{d},\;|P|=n\\}.$

Finally, we define an inverse of minimal dispersion as

$N(\varepsilon,d)=\min\\{n\in\mathbb{N}\,\,|\,\,{\rm
disp}^{*}(n,d)\leq\varepsilon\\}.$

Most of the results will be stated in terms of the inverse of minimal
dispersion.

We also will be dealing with the dispersion on the torus, which is defined
similarly. Let

$\widetilde{\mathcal{R}}_{d}:=\left\\{\prod_{i=1}^{d}I_{i}(a,b)\,\,\,|\,\,\,a,b\in
Q_{d}\right\\},$

where

$I_{i}(a,b):=\begin{cases}(a_{i},b_{i}),&\mbox{ whenever }\,0\leq
a_{i}<b_{i}\leq 1,\\\ [0,1]\setminus[b_{i},a_{i}],&\mbox{ whenever }\,0\leq
b_{i}<a_{i}\leq 1.\end{cases}$

Then we follow definitions above with boxes taken from
$\widetilde{\mathcal{R}}_{d}$ instead of ${\mathcal{R}}_{d}$, that is

$\widetilde{\rm
disp}(P)=\sup\\{|B|\,\,\,|\,\,\,B\in\widetilde{\mathcal{R}}_{d},\,B\cap
P=\emptyset\\},\quad\quad\quad\widetilde{\rm
disp}^{*}(n,d)=\inf_{|P|=n}\widetilde{\rm disp}(P),$

and

$\widetilde{N}(\varepsilon,d)=\min\\{n\in{\mathbb{N}}\,|\,\,\widetilde{\rm
disp}^{*}(n,d)\leq\varepsilon\\}.$

The behavior of dispersion has been intensively studied during the last
decade, and it turns out that it behaves differently in different regimes.
When $\varepsilon$ is extremely small with respect to $d$, namely
$\varepsilon<d^{-d^{2}}$, the best known bound was obtained in [4], where the
authors proved that

$N(\varepsilon,d)\leq\frac{C\,d^{2}\ln d}{\varepsilon}.$

This improved the bound $C_{d}/\varepsilon$ with $C_{d}$ being exponential in
$d$, obtained in previous works [15, 5, 2]. Furthermore, the authors of [4]
have also shown that for $\varepsilon\leq(4d)^{-d}$ one has

$N(\varepsilon,d)\geq\frac{d}{e\varepsilon}.$

On the other hand, for relatively large $\varepsilon$, namely for
$\varepsilon>\frac{1}{d}\,\frac{\ln^{2}d}{\ln\ln d}$, the best upper bound was
obtained in [11], improving previous works [17] and [20],

$N(\varepsilon,d)\leq\frac{C(\ln d)\ln(1/\varepsilon)}{\varepsilon^{2}}.$

Very recently it was shown in [18] that this bound is sharp (up to a
logarithmic factor) whenever $1/4\geq\varepsilon\geq 1/4\sqrt{d}$, more
precisely in this regime one has

$N(\varepsilon,d)\geq\frac{C\ln d}{\varepsilon^{2}\ln(1/\varepsilon)}.$

We would like to mention here the bound from [2], which holds for all
$\varepsilon<1/4$, and which was the first (lower) bound showing that the
dispersion grows with the dimension,

$N(\varepsilon,d)\geq\frac{\log_{2}d}{8\varepsilon}.$

We would also like to note that for large $\varepsilon$, namely
$\varepsilon\in(1/4,1/2)$, the upper bound does not grow with dimension. From
[13] (see also [17]) we have

$N(\varepsilon,d)\leq\frac{\pi}{\sqrt{\varepsilon-1/4}}-3.$

Clearly, for $\varepsilon\geq 1/2$ one has $N(\varepsilon,d)=1$ by taking the
point $(1/2,...,1/2)$.

However, for $\varepsilon$ not so large and not so small with respect to $d$
the picture is different. The best known bound was obtained in [12], improving
previous results from [3, 16, 11], namely

$N(\varepsilon,d)\leq
C\,\frac{d\ln\ln(1/\varepsilon)+\ln(1/\varepsilon)}{\varepsilon}.$ (1)

The proof in [12] shows that a random choice of independent points uniformly
distributed in the cube works. Moreover, it was shown in [8] that using such a
random choice of points, one cannot expect anything better than

$\max\left\\{\frac{c}{\varepsilon}\ln\left(\frac{1}{\varepsilon}\right),\,\frac{d}{2\varepsilon}\right\\},$
(2)

that is, the bound (1) is sharp up to double logarithmic factor for a random
choice of points. In this note we improve (1) by eliminating the
$\ln(1/\varepsilon)$ summand. Of course, as the previous formula shows, the
improvement cannot be achieved for a random choice of points.

Before formulating our main result we briefly discuss the dispersion on the
torus. In [19] the lower bound was obtained valid for all
$\varepsilon\in(0,1)$,

$\widetilde{N}(\varepsilon,d)\geq\frac{d}{\varepsilon}.$

It is interesting to note that contrary to the non-periodic case, the lower
bound is at least linear in $d$ and always grows with $d$, even for
$\varepsilon>1/2$. The best known upper bound was obtained in [12], improving
previous bounds from [11] and [16], namely

$\widetilde{N}(\varepsilon,d)\leq
C\,\frac{d\ln(2d)+\ln(e/\varepsilon)}{\varepsilon}.$ (3)

Here, the random choice of points was also used. In this note we improve (3)
to $Cd(\ln\ln(e/\varepsilon)+\ln(2d))/\varepsilon$ (which is always better),
however the choice of points is, once again, not random.

The main result of this paper is the following.

######  Theorem 2.1.

Let $d\geq 2$ and $\varepsilon\in(0,1)$. Then

$N(\varepsilon,d)\leq\frac{16ed\ln\ln(8/\varepsilon)}{\varepsilon}\quad\quad\mbox{
and
}\quad\quad\widetilde{N}(\varepsilon,d)\leq\frac{8ed(\ln\ln(e/\varepsilon)+\ln(2d))}{\varepsilon}.$

In particular, there exists an absolute constant $C>1$ such that for $d\geq 2$
and $n\geq 4d$ one has

${\rm disp}^{*}(n,d)\leq C\,\frac{d\ln\ln(n/d)}{n}\quad\quad\mbox{ and
}\quad\quad\widetilde{\rm disp}^{*}(n,d)\leq C\,\frac{d(\ln\ln(n/d)+\ln
d)}{n}.$

Thus, combining bounds of Theorem 2.1 with previously known bounds, the
current state of the art for the inverse minimal dispersion on the cube is

$N(\varepsilon,d)\leq\begin{cases}\frac{C\,\ln
d}{\varepsilon^{2}}\,\ln\left(\frac{1}{\varepsilon}\right)&\mbox{ if
}\,\varepsilon\geq\frac{\ln^{2}d}{d\ln\ln(2d)},\\\\[5.81247pt]
\frac{C\,d}{\varepsilon}\,\ln\ln\left(\frac{1}{\varepsilon}\right)&\mbox{ if
}\,\frac{\ln^{2}d}{d\ln\ln(2d)}\geq\varepsilon\geq\exp{(-d^{d})},\\\\[5.81247pt]
\frac{C\,d^{2}\ln d}{\varepsilon}&\mbox{ if
}\,\varepsilon\leq\exp{(-d^{d})},\end{cases}$ (4)

which can be summarized in Figure 1.

Figure 1: Current best known upper bounds on $N(\varepsilon,d)$.

##  3 A new probabilistic lemma

For $\varepsilon>0$ consider sets of all boxes from ${\mathcal{R}}_{d}$
(resp., $\widetilde{\mathcal{R}}_{d}$) of volume at least $\varepsilon$, i.e.,

${\mathcal{B}}_{d}(\varepsilon):=\Big{\\{}B\in{\mathcal{R}}_{d}\,\,\,|\,\,\,|B|\geq\varepsilon\Big{\\}}\quad\mbox{
and
}\quad\widetilde{\mathcal{B}}_{d}(\varepsilon):=\Big{\\{}B\in\widetilde{\mathcal{R}}_{d}\,\,\,|\,\,\,|B|\geq\varepsilon\Big{\\}}.$

Notice that these collections are infinite, and we use a usual approach and
approximate them by finite collections. Following [12], we define a
$\delta$-approximations, which is a slight modification of the notion
introduced in [11]. An essentially same notion was recently considered in a
similar context by M. Gnewuch [7].

######  Definition 3.1 ($\delta$-approximation of
${\mathcal{B}}_{d}(\varepsilon)$ and
$\widetilde{\mathcal{B}}_{d}(\varepsilon)$).

For $0<\delta\leq\varepsilon\leq 1$ we say that a collection
${\mathcal{N}}\subseteq{\mathcal{R}}_{d}$ is a $\delta$-approximation for
${\mathcal{B}}_{d}(\varepsilon)$ iff for all
$B\in{\mathcal{B}}_{d}(\varepsilon)$ there is $B_{0}\in{\mathcal{N}}$ such
that $B_{0}\subseteq B$ and $|B_{0}|\geq\delta$. We define a
$\delta$-approximation for $\widetilde{\mathcal{B}}_{d}(\varepsilon)$ in a
similar way.

In several works on minimal dispersion a variant of the following
probabilistic lemma was a key ingredient (Theorem 1 in [16], Lemma 2.3 with
Remark 2.4 in [11], Lemma 2.2 in [12]).

######  Lemma 3.2.

Let $d\geq 1$ and $\varepsilon,\delta\in(0,1)$. Let ${\mathcal{N}}$ be a
$\delta$-approximation for ${\mathcal{B}}_{d}(\varepsilon)$ and let
$\widetilde{\mathcal{N}}$ be a $\delta$-approximation for
$\widetilde{\mathcal{B}}_{d}(\varepsilon)$. Assume both $|{\mathcal{N}}|\geq
3$ and $|\widetilde{\mathcal{N}}|\geq 3$. Then

$N(\varepsilon,d)\leq\frac{3\ln|{\mathcal{N}}|}{\delta}\quad\quad\mbox{ and
}\quad\quad\widetilde{N}(\varepsilon,d)\leq\frac{3\ln|\widetilde{\mathcal{N}}|}{\delta}.$

Moreover, the random choice of independent points (with respect to the uniform
distribution on $Q_{d}$) gives the result with probability at least
$1-1/|{\mathcal{N}}|$.

Our next lemma improves the bounds of Lemma 3.2. We would like to emphasize
that our proof uses two phases: the first random phase similar to the initial
proof of Lemma 3.2 followed up by the second non-random phase. Notice that
unlike Lemma 3.2, the bounds in Lemma 3.3 will not hold for a random choice of
independent points in $Q^{d}$, in view of (2).

######  Lemma 3.3.

Let $d\geq 1$, $0<\delta\leq 1/3$, $\delta\leq\varepsilon\leq 1$. Let
${\mathcal{N}}$ be a $\delta$-approximation for
${\mathcal{B}}_{d}(\varepsilon)$ and let $\widetilde{\mathcal{N}}$ be a
$\delta$-approximation for $\widetilde{\mathcal{B}}_{d}(\varepsilon)$.
Provided $\delta|{\mathcal{N}}|\geq e$ and
$\delta|\widetilde{\mathcal{N}}|\geq e$ we have

$N(\varepsilon,d)\leq\frac{\ln(4\delta|{\mathcal{N}}|)}{\delta}\quad\quad\mbox{
and
}\quad\quad\widetilde{N}(\varepsilon,d)\leq\frac{\ln(4\delta|\widetilde{\mathcal{N}}|)}{\delta}.$

###### Proof.

We prove the bound for $N(\varepsilon,d)$, the proof for
$\widetilde{N}(\varepsilon,d)$ is similar. Let $|{\mathcal{N}}|=N$. It is
sufficient to show that there is $P\subset Q_{d}$ with
$|P|\leq(\ln(4\delta|{\mathcal{N}}|)/\delta$ such that ${\rm
disp}(P)\leq\varepsilon$. To show that ${\rm disp}(P)\leq\varepsilon$ it is
sufficient to show that for all $A\in{\mathcal{N}}$, we have $A\cap
P\neq\emptyset$.

Let $M$ be a positive integer which will be specified later. Consider a
collection $\mathcal{X}=\\{x_{1},\ldots,x_{M}\\}$ of points in $Q_{d}$ chosen
independently and uniformly at random from $Q_{d}$.

For each $A\in{\mathcal{N}}$ consider the “bad” event
$B_{A}=\\{A\cap\mathcal{X}=\emptyset\\}$. By assumptions on the volume of sets
in ${\mathcal{N}}$ and by independence of $x_{i}$’s, for all
$A\in{\mathcal{N}}$ we have $\mathbb{P}(B_{A})\leq(1-\delta)^{M}$.

Next let $b=b(\mathcal{X})$ be the random variable that counts the number of
bad events, i.e., $b$ counts the number of sets in $\mathcal{N}$ that do not
intersect $\mathcal{X}$. Clearly,

$b=\sum_{A\in{\mathcal{N}}}\chi_{B_{A}},$

where $\chi_{E}$ denotes the indicator of the event $E$. Then

$\mathbb{E}\,b=\sum_{A\in{\mathcal{N}}}\mathbb{E}\chi_{B_{A}}\leq
N(1-\delta)^{M}.$

Therefore there is an instance $Q=\\{x_{1},\ldots,x_{M}\\}\subset Q_{d}$ for
which $b(Q)\leq N(1-\delta)^{M}$. In other words, the set
${\mathcal{N}}_{0}:=\\{A\in{\mathcal{N}}\;|\;A\cap Q=\emptyset\\}$ satisfies
$|{\mathcal{N}}_{0}|\leq N(1-\delta)^{M}.$

Finally, let $Q_{0}$ be a minimal collection of points which has a non-empty
intersection with every set in ${\mathcal{N}}_{0}$ (we always have
$|Q_{0}|\leq|{\mathcal{N}}_{0}|$), and define $P=Q\cup Q_{0}$. Then, by
construction, for all $A\in{\mathcal{N}}$, we have $A\cap P\neq\emptyset$.

Note that

$|P|\leq M+N(1-\delta)^{M}\leq M+Ne^{-\delta M},$

and choose

$M=\left\lceil\frac{\ln(\delta N)}{\delta}\right\rceil.$

Then

$|P|\leq\frac{1}{\delta}+\frac{\ln(\delta
N)}{\delta}+1\leq\frac{\ln(e^{1+\delta}\delta N)}{\delta},$

which completes the proof. ∎

Remark.  Note that in order to have a random choice of points, i.e., in order
to show that for all $A\in{\mathcal{N}}$ one has
$A\cap\mathcal{X}\neq\emptyset$, one needs to work with the event $b=0$. It is
easy to see that for $M=\lfloor(3\ln|{\mathcal{N}}|)/\delta\rfloor$ one has

$\mathbb{P}(b=0)\geq 1-\frac{1}{|{\mathcal{N}}|}.$

This proves Lemma 3.2 which was used in [16, 11, 12].

##  4 Completing the proof of Theorem 2.1

To complete the proof we use the following result from [12] (Propositions 3.2
and 3.3).

######  Proposition 4.1.

Let $d\geq 2$ be an integer, $\varepsilon\in(0,1)$, and $\gamma>0$. Let
$\delta=\varepsilon^{1+\gamma}/4$. There exists a $\delta$-approximations
${\mathcal{N}}$ for ${\mathcal{B}}_{d}(\varepsilon)$ and
$\widetilde{\mathcal{N}}$ for $\widetilde{\mathcal{B}}_{d}(\varepsilon)$ of
cardinalities

$|{\mathcal{N}}|\leq 7d\ln
d\,\frac{(1+1/\gamma)^{d}(\ln(e/\varepsilon^{1+\gamma}))^{d}}{\varepsilon^{1+\gamma}}\quad\mbox{
and }\quad|\widetilde{\mathcal{N}}|\leq 7d\ln
d\,\frac{(1+1/\gamma)^{d}(2d)^{d}}{\varepsilon^{1+\gamma}}.$

The proof of Theorem 2.1 is a straightforward application of Lemma 3.3 and
Proposition 4.1. We provide it here for the sake of completeness.

###### Proof of Theorem 2.1..

Let $\gamma=1/\ln(1/\varepsilon)$. Then
$\varepsilon^{1+\gamma}=\varepsilon/e$. Let
$\delta=\varepsilon^{1+\gamma}/4=\varepsilon/(4e)$. Let ${\mathcal{N}}$ and
$\widetilde{\mathcal{N}}$ be $\delta$-approximations constructed in
Proposition 4.1 with cardinalities

$|{\mathcal{N}}|\leq 7d\ln
d\,\frac{(1+1/\gamma)^{d}\,(\ln(e/\varepsilon^{1+\gamma}))^{d}}{\varepsilon^{1+\gamma}}\leq
7ed\ln
d\,\frac{(\ln(e/\varepsilon))^{d}\,(\ln(e^{2}/\varepsilon))^{d}}{\varepsilon}$

and

$|\widetilde{\mathcal{N}}|\leq 7d\ln
d\,\frac{(1+1/\gamma)^{d}(2d)^{d}}{\varepsilon^{1+\gamma}}\leq 7ed\ln
d\,\frac{(\ln(e/\varepsilon))^{d}\,(2d)^{d}}{\varepsilon}.$

Then

$\ln(4\delta|{\mathcal{N}}|)\leq 2d\ln\ln(e^{2}/\varepsilon)+\ln(7d\ln d)\leq
4d\ln\ln(8/\varepsilon)$

and

$\displaystyle\ln(4\delta|\widetilde{\mathcal{N}}|)$ $\displaystyle\leq
d\ln\ln(e/\varepsilon)+d\ln(2d)+\ln(7d\ln d)\leq
2d\left(\ln\ln(e/\varepsilon)+\ln(2d)\right).$

Hence, Lemma 3.3 now implies the result. ∎

### Acknowledgments

The first named author was supported by the University of Manitoba Research
Grant Program project “Covering problems with applications in Euclidean Ramsey
theory and convex geometry” and PIMS PDF. The second named author gratefully
acknowledge the support of the Leibniz Center for Informatics, where several
discussions about this research were held during the Dagstuhl Seminar
“Algorithms and Complexity for Continuous Problems” (Seminar ID 23351). Part
of this work was conducted during the first named author’s visit to University
of Alberta.

## References

  * [1] N. Alon, J.H. Spencer, The probabilistic method, 4th edition, Wiley Series in Discrete Mathematics and Optimization, Hoboken, NJ: John Wiley & Sons, 2016\.
  * [2] C. Aistleitner, A. Hinrichs, D. Rudolf, On the size of the largest empty box amidst a point set, Discrete Appl. Math. 230 (2017), 146–150.
  * [3] A. Blumer, A. Ehrenfeucht, D. Haussler, M. Warmuth, Learnability and the Vapnik–Chervonenkis dimension, J. Assoc. Comput. Mach. 36 (1989), 929–965.
  * [4] B. Bukh, T. Chao, Empty axis-parallel boxes, IMRN 18 (2022), 13811–13828.
  * [5] A. Dumitrescu, M. Jiang, On the largest empty axis-parallel box amidst $n$ points, Algorithmica 66 (2013), 225–248.
  * [6] A. Dumitrescu, M. Jiang, Computational geometry column 60, ACM SIGACT News 45 (2014), 76–82.
  * [7] M. Gnewuch, private communications.
  * [8] A. Hinrichs, D. Krieg, R.J. Kunsch, D. Rudolf, Expected dispersion of uniformly distributed points, J. Complexity 61 (2020), Paper No. 101483, 9pp.
  * [9] A. Hinrichs, J. Prochno, M. Ullrich, J. Vybíral, The minimal k-dispersion of point sets in high-dimensions, J. Complexity 51 (2019), 68–78.
  * [10] E. Hlawka, Abschätzung von trigonometrischen Summen mittels diophantischer Approximationen, Österreich. Akad. Wiss. Math.-Naturwiss. Kl. S.-B. II 185 (1976), 43–50.
  * [11] A.E. Litvak, A remark on the minimal dispersion, Commun. Contemp. Math. 23 (2021), Paper No. 2050060, 18pp.
  * [12] A.E. Litvak, G.V. Livshyts, New bounds on the minimal dispersion, J. Complexity 72 (2022), Paper No. 101648, 13pp.
  * [13] K. S. MacKay, Minimal dispersion of large volume boxes in the cube, J. Complexity 72 (2022), Paper No. 101650, 23pp.
  * [14] C.A. Rogers, A note on coverings, Mathematika 4 (1957), 1–6.
  * [15] G. Rote, R.F. Tichy, Quasi-Monte Carlo methods and the dispersion of point sequences, Math. Comput. Modelling 23 (1996), 9–23.
  * [16] D. Rudolf, An upper bound of the minimal dispersion via delta covers, Contemporary Computational Mathematics - A Celebration of the 80th Birthday of Ian Sloan, Springer-Verlag (2018), 1099–1108.
  * [17] J. Sosnovec, A note on the minimal dispersion of point sets in the unit cube, European J. of Comb. 69 (2018), 255–259.
  * [18] M. Trödler, J. Volec, J. Vybíral, A tight lower bound on the minimal dispersion, European J. Combin. 120 (2024), Paper No. 103945, 5 pp.
  * [19] M. Ullrich, A lower bound for the dispersion on the torus, Mathematics and Computers in Simulation 143 (2018), 186–190.
  * [20] M. Ullrich, J. Vybíral, An upper bound on the minimal dispersion, Journal of Complexity 45 (2018), 120–126.

Andrii Arman
Department of Mathematics,
University of Manitoba,
Winnipeg, MB, R3T 2N2, Canada
e-mail<EMAIL_ADDRESS>

Alexander E. Litvak
Dept. of Math. and Stat. Sciences,
University of Alberta,
Edmonton, AB, Canada, T6G 2G1.
e-mail<EMAIL_ADDRESS>