# Exact results for the order picking time distribution under return routing

Tim Engels  Department of Mathematics and Computer Science, Eindhoven
University of Technology  Ivo Adan  Department of Industrial Engineering and
Innovation Sciences, Eindhoven University of Technology  Onno Boxma
Department of Mathematics and Computer Science, Eindhoven University of
Technology  Jacques Resing  Department of Mathematics and Computer Science,
Eindhoven University of Technology

###### Abstract

This paper derives exact expressions for the Laplace-Stieltjes transform of
the order picking time in single- and 2-block warehouses. We consider manual
warehouses that deploy return routing and assume that order sizes follow a
Poisson distribution. The results in this paper apply to a wide range of
storage policies, including but not restricted to class-based and random
storage. Furthermore, we compare the performance of the storage policies and
warehouse lay-outs by using numerical inversion of the Laplace-Stieltjes
transforms.

##  1 Introduction

Warehouses have grown to become an essential part of today’s supply chains, as
is illustrated in Boysen et al. (2019). With this growth, the analysis and
optimization of warehouse design has also become ever more important. In
particular, the order picking process, the retrieval of items from a
warehouse, stands central in many designs due to its high cost. Tompkins et
al. (2011), for example, estimated that order picking contributes at least
$50\%$ of all operating costs.

The analysis and optimization of order picking in itself already is a
challenging subject, which is highlighted by the vast amount of literature on
the topic, see e.g. the literature review of Boysen et al. (2019). This
complexity stems from the many elements that affect the picking process, such
as: the lay-out of warehouses, routing of pickers and storage assignment. On
top of that, optimal solutions often are impractical, e.g., optimal routes
cause confusion amongst pickers (Gademann et al., 2001). To overcome such
practical issues, warehouses often apply straightforward, yet efficient,
policies.

The performance evaluation of such policies is also crucial to warehousing
design. Literature on this topic varies from simulation studies (Petersen and
Aase, 2004) to average route-length approximations (Hall, 1993). Dijkstra and
Roodbergen (2017) derived exact formulas for the average route-length, which
were then used to construct optimal class-based storage lay-outs. A different
point of view came from Chew and Tang (1999), where the authors proposed to
compare order batching policies based on the average order-lead time, i.e. the
delay experienced by a customer, rather than the average order picking time.
The authors modeled the warehouse as a multi-server queue and used a (two
moment) approximation of the expected order-lead time to find the optimal
batch size. More recently, Engels et al. (2022) provided exact results for the
first two moments of the order picking time for several routing policies.
Using similar queueing approximations the authors showed that the policy with
the shortest average route-length, but higher second moment, could possibly
result in higher average order-lead times.

In this paper we step away from the moments of the order picking time and
instead derive exact results for the Laplace-Stieltjes transform (abv. LST) of
the order picking time distribution. For this, we analyze both single- and
2-block warehouses that deploy the return routing policy. We assume that the
order size follows a Poisson distribution and exploit that the numbers of
picks in the different aisles are independent under this assumption. The
results in this paper allow for the analysis of other properties of the
picking time distribution, such as tail probabilities, and an accurate
approximation of the order-lead time distribution (see e.g. Tijms, 1986, pages
292-300).

Multi-block warehouses have been discussed extensively in literature, and have
proven to be interesting for both practical and theoretical purposes. The
substantial literature mostly focuses on the expected order picking time and
varies from exact optimization (Theys et al., 2010) to simulation studies
(Roodbergen and De Koster, 2001). Extensions of the models, as well as
queueing theoretical applications of the results, are also well-discussed.
Notwithstanding the amount of research, the actual distribution of the order
picking time was not yet determined. The design of warehouses, so far, has
also been restricted to maximization of the throughput and minimization of the
experienced delay. In this work, we extend the research by finding exact
expressions for the order picking time distribution in single- and 2-block
warehouses. These results can be used to design warehouses more carefully by
also considering other performance statistics, e.g. tail probabilities.
Additionally, distributional knowledge increases the accuracy of queueing
theoretical analyses (Gupta et al., 2010).

This paper is built up as follows. In Section 2 we provide the model
descriptions for the single- and 2-block warehouses, as well as a description
of the return routing policy. In Section 3 we derive the LST of the order
picking time distribution for single- and 2-block warehouses, under return
routing. We also show how to apply these results to two storage policies:
random and class-based. In Section 4 we use numerical inversion methods for
Laplace Transforms to compare these storage policies as well as the single-
and 2-block warehouses.

##  2 Warehousing Model

We consider single- and 2-block warehouses, consisting of $k$ equally spaced
parallel aisles of length $l$ and width $w_{a}$. We assume that the warehouses
consist of a single cross-aisle from which the pickers can enter the aisles.
The difference between the single- and 2-block warehouses comes from the
cross-aisle placement. In the single-block warehouse the cross-aisle is placed
in the front, whereas in the 2-block warehouse, the cross-aisle is placed in
the middle and splits the aisles in two sub-aisles (each of length $l/2$), see
Figure 1. In the sequel, we write sub-aisle $(i,j)$ for the sub-aisle in the
$i$-th aisle and $j$-th block, that is: sub-aisle $(i,1)$ is the upper half of
aisle $i$ in the 2-block warehouse. The depot at which the pickers receive and
deliver their orders is located on the left-hand side of the cross-aisle,
indicated as I/O in Figure 1.

We assume that the pickers walk with a fixed speed $v$ and spend a random time
$P$ to pick an arbitrary item, where $P$ has LST $\Phi_{P}(s)$. Furthermore,
we assume that the time spent picking an item is independent of the pick
locations and other picking times. The order size, $M$, is assumed to follow a
Poisson distribution with mean $\lambda$. For example, this is valid when item
demands arrive according to a Poisson process and are batched as a single
order after a fixed time.

We assume that the items in an order are independently located throughout the
warehouse according to some storage policy. Let $p_{i,j}$ denote the
probability that an arbitrary item in an order falls within sub-aisle $(i,j)$
and let $F_{i,j}(x)$ denote the probability that an arbitrary pick in this
sub-aisle is within distance $x$ (as a fraction of the sub-aisle length) of
the cross-aisle, for $i=1,...,k$ and $j=1,2$. Throughout this paper, we write
$p_{i}$ and $F_{i}(\cdot)$ in the case of a single-block warehouse. This
description allows for both continuous and discrete locations within aisles
and represents a wide range of storage policies, e.g., the random and class-
based storage policies. Under the former, the item locations in an order are
uniformly distributed amongst and within the aisles. In the single-block
warehouse this implies $p_{i}=1/k$ and $F_{i}(x)=x$. In Section 3.3 we discuss
the random and class-based policy more extensively. In particular, we show
that the latter policy results in piece-wise linear $F_{i}(\cdot)$.

(a) Single-block warehouse.

(b) 2-block warehouse.

Figure 1: Examples of the return routing policy in both single- and 2-block
warehouses.

###  2.1 Routing Policy

In this paper we focus on the return routing policy, where pickers traverse
each aisle until the furthest item that has to be picked and then traverse
back. In the single-block warehouse this entails that the pickers walk across
the front cross-aisle and enter each aisle in which they have to pick items.
Within such aisles, the pickers pick all items up until the furthest item and
then walk back to the cross-aisle. After all aisles are visited, the picker
returns to the depot. An example of this policy is given in Figure 1(a).

We introduce the following notation. Let $N_{i}$ denote the number of items to
pick in aisle $i$ and let $A_{i}$ be the location of the furthest item in
aisle $i$ (as a fraction of the aisle length $l$). Furthermore, we write
$k^{+}$ for the furthest aisle in which items have to be picked. Using this
notation, we find that the picking time contribution of aisle $i$, denoted as
$T_{i}$, is given by:

$\displaystyle
T_{i}=\sum_{j=1}^{N_{i}}P_{j}+\frac{2l}{v}A_{i}+\frac{2w_{a}}{v}\mathbbm{1}\\{i<k^{+}\\}.$
(1)

This relation consists of 3 different components: the time spent picking
items, the traveling time within the aisle and the time spent traveling
towards the next aisle. Since $k^{+}$ is the last aisle in which items have to
be picked, the last term disappears when $i=k^{+}$. For $i>k^{+}$ the aisle
will not have a picking time contribution and hence all terms will equal zero.

Consequently, the total order picking time, $T$, is given by:

$\displaystyle T=\begin{dcases}\sum_{i=1}^{k^{+}-1}T_{i}+T_{k^{+}}&\text{if
}M>0,\\\ 0&\text{else.}\end{dcases}$ (2)

Here, we explicitly distinguish between $T_{i}$ for $i<k^{+}$ and $T_{k^{+}}$
because of the third term in (1) and the fact that aisle $k^{+}$ must consist
of at least one pick, while the other aisles can still have zero picks.

Return routing in a 2-block warehouse is similar. The pickers walk across the
cross-aisle and apply the return routing policy for a single-block warehouse
on the lower block up until all items in said block are picked, hence they do
not yet return to the I/O point. The pickers than traverse to the furthest
(full) aisle with picks and apply the return routing policy on the upper block
from right to left. When all items are picked, the pickers deliver their order
at the depot-point. An illustration of this policy is given in Figure 1(b).

In line with the notation for the single-block warehouse, we write $k^{+}$ for
the furthest aisle in which items have to be picked, $N_{i,j}$ for the number
of items to pick in sub-aisle $(i,j)$, and $A_{i,j}$ for the location of the
furthest item in sub-aisle $(i,j)$ (as a fraction of the sub-aisle length
$l/2$ and with respect to the cross-aisle). Then $T_{i}$ is given by:

$\displaystyle
T_{i}=\sum_{j=1}^{N_{i,1}+N_{i,2}}P_{j}+\frac{l}{v}\Big{(}A_{i,1}+A_{i,2}\Big{)}+\frac{2w_{a}}{v}\mathbbm{1}\\{i<k^{+}\\}.$
(3)

Remark that relation (2) still holds for the 2-block warehouse.

##  3 Results

In this section we derive the LST of the order picking time distribution in
both single- (Section 3.1) and 2-block warehouses (Section 3.2). Afterwards,
in Section 3.3 we illustrate how these results can be applied to two well-
known storage policies. The analysis in this section is based upon the key
assumption that the total order size follows a Poisson distribution. Since the
number of picks in the different (sub-)aisles follows a multinomial
distribution with the order size as number of trials, we have that the numbers
of picks in the different aisles are independent under this assumption and
follow a Poisson distribution with parameter $p_{i,j}\lambda$.

###  3.1 Single-block warehouse

Recall the relation for the order picking time in Equation (2), which splits
the order picking time in the picking time contribution of the last aisle and
the sum of the contributions of the previous aisles. The latter term is a
random sum of ($k^{+}-1$) i.i.d. random variables ($T_{i}$), hence deriving
the distribution of the order picking time requires to derive the distribution
of both $k^{+}$ and $T_{i}$. We start by remarking that $k^{+}$ is the first
aisle with picks coming from aisle $k$ to aisle $1$. For the aisle
contributions $T_{i}$, we notice that $A_{i}$ is the furthest pick location in
the aisle, i.e., $A_{i}$ is the maximum of $N_{i}$ i.i.d. random variables.
Therefore we have that: $\mathbb{P}(A_{i}\leq x|N_{i}=n)=F_{i}(x)^{n}$. Using
these insights we find the following:

######  Lemma 3.1.

The following equations hold:

$\displaystyle\mathbb{P}(k^{+}=j)=\exp\Big{(}-\lambda\sum_{i=j+1}^{k}p_{i}\Big{)}\big{[}1-\exp(-\lambda
p_{j})\big{]}\quad\text{for }j=1,...,k,$ (4)
$\displaystyle\mathbb{E}[\exp(-sT_{i})|k^{+}>i]=\exp\Big{(}-\lambda
p_{i}-\frac{2w_{a}s}{v}\Big{)}\bigg{[}\begin{aligned} &\exp\Big{(}\lambda
p_{i}\Phi_{P}(s)-\frac{2ls}{v}\Big{)}\\\
&+\frac{2ls}{v}\int_{x=0}^{1}\exp\Big{(}-\frac{2ls}{v}x+\lambda
p_{i}\Phi_{P}(s)F_{i}(x)\Big{)}\mathrm{d}x\bigg{]},\end{aligned}$ (5) and
$\displaystyle\mathbb{E}[\exp(-sT_{i})|k^{+}=i]=\frac{\exp(2w_{a}s/v)\mathbb{E}[\exp(-sT_{i})|k^{+}>i]-\exp(-\lambda
p_{i})}{1-\exp(-\lambda p_{i})}.$ (6)

###### Proof.

The first statement follows immediately from the given that $k^{+}$ is the
right-most non-empty aisle in the warehouse, combined with the fact that
$\mathbb{P}(N_{i}=0)=\exp(-\lambda p_{i})$.

For (5) we distinguish between the cases $N_{i}=0$ and $N_{i}>0$:

$\displaystyle\mathbb{E}[\exp(-sT_{i})|k^{+}>i]$
$\displaystyle=\exp\Big{(}-\lambda
p_{i}-\frac{2w_{a}s}{v}\Big{)}+\mathbb{E}[\exp(-sT_{i})\mathbbm{1}\\{N_{i}>0\\}|k^{+}>i].$
(7)

For the second term we realize that, conditional on $N_{i}$, we have by
partial integration that:

$\mathbb{E}[\exp(-sA_{i})|N_{i}=n]=\exp(-s)+s\int_{x=0}^{1}\mathbb{P}(A_{i}\leq
x|N_{i}=n)\exp(-sx)\mathrm{d}x.$

Now using $\mathbb{P}(A_{i}\leq x|N_{i}=n)=F_{i}(x)^{n}$ and by summing over
all possible values of $N_{i}$ we find:

$\displaystyle\mathbb{E}[\exp(-sT_{i})\mathbbm{1}\\{N_{i}>0\\}|k^{+}>i]$
$\displaystyle=\exp\Big{(}-\lambda
p_{i}-\frac{2w_{a}s}{v}\Big{)}\sum_{n=1}^{\infty}\frac{\lambda^{n}p_{i}^{n}}{n!}\Phi_{P}(s)^{n}\bigg{[}\exp\Big{(}-\frac{2ls}{v}\Big{)}+\frac{2ls}{v}\int_{x=0}^{1}\exp(-\frac{2ls}{v}x)F_{i}(x)^{n}\mathrm{d}x\bigg{]}$
$\displaystyle=\exp\Big{(}-\lambda
p_{i}-\frac{2w_{a}s}{v}\Big{)}\bigg{[}\begin{aligned} &\exp\Big{(}\lambda
p_{i}\Phi_{P}(s)-\frac{2ls}{v}\Big{)}-\exp\Big{(}-\frac{2ls}{v}\Big{)}\\\
&+\frac{2ls}{v}\int_{x=0}^{1}\bigg{\\{}\exp\Big{(}-\frac{2ls}{v}x+\lambda
p_{i}\Phi_{P}(s)F_{i}(x)\Big{)}-\exp\Big{(}-\frac{2ls}{v}x\Big{)}\bigg{\\}}\mathrm{d}x\bigg{]}\end{aligned}$
$\displaystyle=\exp\Big{(}-\lambda
p_{i}-\frac{2w_{a}s}{v}\Big{)}\bigg{[}\begin{aligned} &\exp\Big{(}\lambda
p_{i}\Phi_{P}(s)-\frac{2ls}{v}\Big{)}-\exp\Big{(}-\frac{2ls}{v}\Big{)}\\\
&+\frac{2ls}{v}\int_{x=0}^{1}\exp\Big{(}-\frac{2ls}{v}x+\lambda
p_{i}\Phi_{P}(s)F_{i}(x)\Big{)}\mathrm{d}x+\exp\Big{(}-\frac{2ls}{v}\Big{)}-1\bigg{]}.\end{aligned}$

Remark that the terms $\exp(-2ls/v)$ cancel and that, after substitution in
(7), the $-1$ term cancels against the term for $N_{i}=0$, proving (5).

Lastly, for (6) we note that:

$\displaystyle\Big{\\{}T_{i}\Big{|}k^{+}=i\Big{\\}}\,$
$\displaystyle\overset{d}{=}\bigg{\\{}\sum_{m=1}^{N_{i}}P_{m}\bigg{|}N_{i}>0\bigg{\\}}+\frac{2l}{v}\Big{\\{}A_{i}\Big{|}N_{i}>0\Big{\\}}\overset{d}{=}\frac{-2w_{a}}{v}+\Big{\\{}T_{i}\Big{|}k^{+}>i,N_{i}>0\Big{\\}}.$

For the LST of $\Big{\\{}T_{i}\Big{|}k^{+}=i\Big{\\}}$ we can thus use the
distribution on the right-hand side. By then applying the law of total
probability on the event $N_{i}>0$ we find:

$\displaystyle\mathbb{E}\Big{[}\exp\Big{(}-sT_{i}+\frac{2w_{a}s}{v}\Big{)}\Big{|}k^{+}>i,N_{i}>0\Big{]}$
$\displaystyle=\frac{\exp(2w_{a}s/v)\mathbb{E}[\exp(-sT_{i}\Big{)}\Big{|}k^{+}>i]-\exp(-\lambda
p_{i})}{1-\exp(-\lambda p_{i})}.\qed$

The results in Lemma 3.1 provide the main ingredients for the LST of the total
order picking time. Starting from (2), we sum over the possible values of
$k^{+}$. Since the aisle contributions, conditional on $k^{+}$, are
independent (up to aisle $k^{+}$) we have that the LST is simply the product
of the transforms of each aisle contribution. In conclusion we get:

######  Theorem 3.2.

The total order picking time in a single-block warehouse, under return
routing, has the following Laplace-Stieltjes transform:

$\displaystyle\mathbb{E}[\exp(-sT)]=\exp(-\lambda)+\sum_{i=1}^{k}\exp\Big{(}-\lambda\sum_{m=i+1}^{k}p_{m}\Big{)}\Big{[}\exp\Big{(}\frac{2w_{a}s}{v}\Big{)}\psi_{T_{i}}(s)-\exp\big{(}-\lambda
p_{i}\big{)}\Big{]}\prod_{l=1}^{i-1}\psi_{T_{l}}(s),$ (8)

where $\psi_{T_{l}}(s)=\mathbb{E}[\exp(-sT_{l})|k^{+}>l]$ for $l=1,2,...k$,
given in Lemma 3.1.

###  3.2 2-block warehouse

The order picking time, $T$, in the 2-block warehouse still satisfies the
relation in (2), where $T_{i}$ is given in (3). The main idea behind Theorem
3.2 therefore still applies, resulting in the same relation between the LST of
the total order picking time and the aisle contributions. The main difference
between the transforms comes from the contribution of each aisle, which we
analyze by splitting this contribution into two independent parts: one for
each block. This results in Theorem 3.3.

######  Theorem 3.3.

The total order picking time in a 2-block warehouse, under return routing, has
the following Laplace-Stieltjes transform:

$\displaystyle\mathbb{E}[\exp(-sT)]=\exp(-\lambda)+\sum_{i=1}^{k}\Big{[}$
$\displaystyle\exp\Big{(}\frac{2w_{a}s}{v}\Big{)}\psi_{T_{i,1}}(s)\psi_{T_{i,2}}(s)-\exp\big{(}-\lambda(p_{i,1}+p_{i,2})\big{)}\Big{]}$
(9)
$\displaystyle\times\exp\Big{(}-\lambda\sum_{m=i+1}^{k}(p_{m,1}+p_{m,2})\Big{)}\prod_{l=1}^{i-1}\psi_{T_{l,1}}(s)\psi_{T_{l,2}}(s),$

where:

$\displaystyle\psi_{T_{l,j}}(s)=\exp\Big{(}-\lambda
p_{l,j}-\frac{w_{a}s}{v}\Big{)}\Big{[}\exp\Big{(}\lambda
p_{l,j}\Phi_{P}(s)-\frac{ls}{v}\Big{)}+\frac{ls}{v}\int_{x=0}^{1}\exp\Big{(}-\frac{ls}{v}x+\lambda
p_{l,j}\Phi_{P}(s)F_{l,j}(x)\Big{)}\mathrm{d}x\Big{]},$ (10)

for $l=1,2,...k$ and $j=1,2$.

###### Proof.

By similar reasoning as in Lemma 3.1 we have:

$\displaystyle\mathbb{P}(k^{+}=j)=\exp\Big{(}-\lambda\sum_{i=j+1}^{k}(p_{i,1}+p_{i,2})\Big{)}\big{[}1-\exp\big{(}-\lambda(p_{j,1}+p_{j,2})\big{)}\big{]}\quad\text{for:
}j=1,...,k.$ (11)

For the aisle contributions we note that $T_{i}$ can be written as the sum of
the contributions in the upper and lower block and also split the horizontal
travel distance contribution in two, i.e.:

$\displaystyle
T_{i}=T_{i,1}+T_{i,2}=\sum_{j=1}^{N_{i,1}}P_{j}+\frac{l}{v}A_{i,1}+\frac{w_{a}}{v}\mathbbm{1}\\{i<k^{+}\\}+\sum_{j=1}^{N_{i,2}}P_{j}+\frac{l}{v}A_{i,2}+\frac{w_{a}}{v}\mathbbm{1}\\{i<k^{+}\\}.$
(12)

Observe that $T_{i,1},T_{i,2}$ are similar to $T_{i}$ in the single-block
warehouse. In fact, they have the same distribution for $i<k^{+}$ when we
take: $w_{a}/2,\lambda p_{i,j}$ and $l/2$ instead of $w_{a},\lambda p_{i}$ and
$l$. We can thus use the result in Lemma 3.1 for the contribution of a single
sub-aisle, resulting in (10). Again, these contributions are independent for
$i<k^{+}$ and hence the transform of $T_{i,1}+T_{i,2}$ is the product of the
transforms: $\psi_{T_{i}}(s)=\psi_{T_{i,1}}(s)\psi_{T_{i,2}}(s)$. For
$i=k^{+}$ we, however, have to be careful, since either $N_{i,1}$ or $N_{i,2}$
can still be zero. For this aisle, we remark:

$\displaystyle\mathbb{E}\big{[}\exp(-sT_{i})\big{|}i=k^{+}\big{]}$
$\displaystyle=\mathbb{E}\bigg{[}\exp\Big{(}-s\Big{[}-2\frac{w_{a}}{v}+T_{i,1}+T_{i,2}\Big{]}\Big{)}\bigg{|}i<k^{+},N_{i}>0\bigg{]}.$

Similar to (6) we apply the law of total probability. By independence of the
$T_{i,1},T_{i,2}$ we now find:

$\displaystyle\mathbb{E}\big{[}\exp(-sT_{i})|i=k^{+}\big{]}$
$\displaystyle=\frac{\exp(2w_{a}s/v)\psi_{T_{i,1}}(s)\psi_{T_{i,2}}(s)-\exp(-\lambda(p_{i,1}+p_{i,2}))}{1-\exp(-\lambda(p_{i,1}+p_{i,2}))}.$

The denominator cancels against the term in $\mathbb{P}(k^{+}=i)$ and we
therefore find the result as in the theorem. ∎

###  3.3 Storage policy examples

We apply these results for two different storage policies. Firstly, we discuss
the random storage policy, under which each item in an order is uniformly
located amongst and within the aisles. Under this policy, we therefore have
$p_{i}=1/k$ and $F_{i}(x)=x$ for all $i$. Substituting this in Theorem 3.2
results in the following expression:

######  Lemma 3.4.

Under the random storage policy the LST of the total order picking time, for
the single-block warehouse, is given by:

$\displaystyle\mathbb{E}[\exp(-sT)]=\exp(-\lambda)\bigg{\\{}1+$
$\displaystyle\exp\Big{(}\frac{\lambda}{k}\Big{)}\Big{[}\exp\Big{(}\frac{2w_{a}s}{v}\Big{)}\psi_{T_{i}}(s)-\exp\Big{(}-\frac{\lambda}{k}\Big{)}\Big{]}\frac{1-\exp(\lambda)\psi_{T_{i}}(s)^{k}}{1-\exp(\lambda/k)\psi_{T_{i}}(s)}\bigg{\\}},$
(13)

where:

$\displaystyle\psi_{T_{i}}(s)=\begin{aligned}
\exp&\Big{(}-\frac{\lambda}{k}-\frac{2w_{a}s}{v}\Big{)}\cdot\bigg{\\{}1+\frac{\Phi_{P}(s)\lambda
v}{\Phi_{P}(s)\lambda
v-2lks}\bigg{[}\exp\Big{(}\frac{\lambda}{k}\Phi_{P}(s)-\frac{2ls}{v}\Big{)}-1\bigg{]}\bigg{\\}}.\end{aligned}$
(14)

###### Proof.

Using $p_{i}=1/k,F_{i}(x)=x$ in Theorem 3.2 gives:

$\displaystyle\psi_{T_{i}}(s)$
$\displaystyle=\exp\Big{(}-\frac{\lambda}{k}-\frac{2w_{a}s}{v}\Big{)}\bigg{[}\exp\Big{(}\frac{\lambda}{k}\Phi_{P}(s)-\frac{2ls}{v}\Big{)}+\frac{2ls}{v}\int_{x=0}^{1}\exp\Big{(}-\frac{2ls}{v}x+\frac{\lambda}{k}\Phi_{P}(s)x\Big{)}\mathrm{d}x\bigg{]}$
$\displaystyle=\exp\Big{(}-\frac{\lambda}{k}-\frac{2w_{a}s}{v}\Big{)}\bigg{[}\exp\Big{(}\frac{\lambda}{k}\Phi_{P}(s)-\frac{2ls}{v}\Big{)}+\frac{2lks}{\Phi_{P}(s)\lambda
v-2lks}\Big{[}\exp\Big{(}\frac{\lambda}{k}\Phi_{P}(s)-\frac{2ls}{v}\Big{)}-1\Big{]}\bigg{]},$

which results in (14) by adding the terms. Next turning to (8) we use that
$\psi_{T_{i}}(\cdot)$ is the same for all $i$ and hence:

$\displaystyle\mathbb{E}[\exp(-sT)]$
$\displaystyle=\exp(-\lambda)+\sum_{j=1}^{k}\Big{[}\exp\Big{(}\frac{2w_{a}s}{v}\Big{)}\psi_{T_{i}}(s)-\exp\Big{(}-\frac{\lambda}{k}\Big{)}\Big{]}\exp\Big{(}-\frac{\lambda(k-j)}{k}\Big{)}\psi_{T_{i}}(s)^{j-1}$
$\displaystyle=\exp(-\lambda)+\Big{[}\exp\Big{(}\frac{2w_{a}s}{v}\Big{)}\psi_{T_{i}}(s)-\exp\Big{(}-\frac{\lambda}{k}\Big{)}\Big{]}\exp\Big{(}-\frac{\lambda(k-1)}{k}\Big{)}\sum_{j=0}^{k-1}\Big{[}\exp\Big{(}\frac{\lambda}{k}\Big{)}\psi_{T_{i}}(s)\Big{]}^{j}.$
(15)

Rewriting the geometric sum in (15) yields (13). ∎

Note: From this equation one can find the first moment of $T$ by evaluating
the derivative of the LST at $0$. The resulting expectation coincides with the
results in Engels et al. (2022):

$\displaystyle\mathbb{E}[T]$
$\displaystyle=\frac{2lk^{2}}{v\lambda}\Big{[}1-\exp\Big{(}-\frac{\lambda}{k}\Big{)}\Big{]}+\frac{2w_{a}}{v}\Big{[}k-\frac{1-\exp(-\lambda)}{1-\exp(-\lambda/k)}\Big{]}.$

In the 2-block warehouse the results are very similar: relation (13) still
holds, yet now we have $p_{i,j}=1/(2k)$. Consequently, $\psi_{T_{i}}(\cdot)$
is given by the square of the expression in (14) with $k$ replaced by $2k$ and
$w_{a}$ replaced by $w_{a}/2$. The square comes from the fact that
$\psi_{T_{i,1}}(\cdot)$ and $\psi_{T_{i,2}}(\cdot)$ are equal.

In practice, warehouses often exploit knowledge on the demand of items. For
example, storing high-demand items in close locations results in lower order
picking times. One well-known example of such a storage policy is the class-
based policy, see e.g. (Dijkstra and Roodbergen, 2017).

Under the class based storage policy, each item is grouped in a certain demand
class. Each of these classes has a separate area within the warehouse in which
the corresponding items are located uniformly at random. Under this policy,
one can thus prevent long routes by placing high demand items (high class)
close to the cross-aisles.

Below we provide a mathematical definition for class-based storage. We follow
the notation and definition in Dijkstra and Roodbergen (2017) and assume that
there are a total of $Q$ classes and write $\tilde{p}_{q}$ for the probability
that an arbitrary item in an order is of class $q$. We consider the case where
each aisle consists of $Q$ different ranges, one for each class. Let
$u_{q}^{(i)}$ be the upper-bound of the locations of type $q$ items in aisle
$i$, such that all class $q$ items in aisle $i$ fall between $u_{q-1}^{(i)}$
and $u_{q}^{(i)}$. Here we use the notion that $u_{0}^{(i)}=0$ and
$u_{Q}^{(i)}=1$. We thus have that type $q$ items make up a total of $f_{q}$
space, where $f_{q}$ is given by:
$f_{q}:=\sum_{i=1}^{k}\big{(}u_{q}^{(i)}-u_{q-1}^{(i)}\big{)}$.

Although the definition above is a special case of the class-based storage, it
is known to provide optimal solutions, w.r.t. the average route-length, for
return routing when $\tilde{p}_{1}>\tilde{p}_{2}>...>\tilde{p}_{Q}$ (Dijkstra
and Roodbergen, 2017).

In case of 2 classes we can translate this to our model by noting:

$\displaystyle p_{i}$ $\displaystyle=\sum_{q=1}^{2}\mathbb{P}(\text{item of
class }q)\mathbb{P}(\text{item in aisle }i|\text{item of class
}q)=\tilde{p}_{1}\cdot\frac{u_{1}^{(i)}}{f_{1}}+\tilde{p}_{2}\cdot\frac{1-u_{1}^{(i)}}{f_{2}},$
$\displaystyle F_{i}(x)$ $\displaystyle=\sum_{q=1}^{2}\mathbb{P}(\text{item of
class }q|\text{item in aisle }i)\mathbb{P}(\text{Location in aisle $i$ before
}x|\text{item of class }q)$
$\displaystyle=\begin{dcases}\frac{\tilde{p}_{1}}{p_{i}}\frac{u_{1}^{(i)}}{f_{1}}\cdot\frac{x}{u_{1}^{(i)}}=\frac{\tilde{p}_{1}}{p_{i}}\frac{x}{f_{1}}&\text{if
}x<u_{1}^{(i)},\\\
\frac{\tilde{p}_{1}}{p_{i}}\frac{u_{1}^{(i)}}{f_{1}}+\frac{\tilde{p}_{2}}{p_{i}}\frac{x-u_{1}^{(i)}}{f_{2}}&\text{if
}x\geq u_{1}^{(i)}.\end{dcases}$

One can substitute these relations into Theorem 3.2 yielding the LST of the
total order picking time. In fact, one can analytically find the integral in
(5) by splitting the integral in two parts, one for each class. The case for a
general number of classes follows similarly.

Remark: This class-based storage model is a continuous version of the model in
Dijkstra and Roodbergen (2017). By letting the number of locations in Dijkstra
and Roodbergen (2017) go to infinity, and appropriately scaling the
probabilities, the model becomes equivalent to the one described here.

##  4 Numerical Illustration

In this section we focus on the comparison of routing and storage policies. We
do this by numeric inversion of LST’s using the method of Abate and Whitt
(Abate and Whitt, 1992).

We consider a warehouse model with $l=20m,w_{a}=2.5m,v=0.83m/s$ and let
$P\sim\mathrm{Exp}(1/5)$. For the class-based storage policy, we consider the
example in Dijkstra and Roodbergen (2017) for a $50/30/20$ distribution of
classes and a $20/30/50\%$ of space required per class, where $k=15$ and
$\lambda=10$. We take the class-based areas as graphically given in Dijkstra
and Roodbergen (2017). For the 2-block warehouse we assign the same area
division to the lower block and mirror the areas in the upper block, see the
illustration in Figure 2; the various colours indicate the regions in which
the items of the corresponding class are stored. Besides the class-based
storage, we also consider random storage. The picking time densities for these
storage policies are shown in Figure 3.

(a) Single-block warehouse.

(b) 2-block warehouse.

Figure 2: Illustrations of the class-based storage areas for the numerical
example. Figure 3: Illustration of the density of the total order picking time
under the Class-Based Storage and Random Storage policies.

We can see that the single-block warehouse is outperformed by the multi-block
warehouse and the class-based storage, in turn, outperforms the random storage
policy. Similar comparisons have previously been found in literature, yet the
quantification that is allowed by knowing the complete distribution is new.
For example the plots reveal a similarity to the normal density. This
similarity can be explained by relation (2), showing that the order total
picking time is a random sum of independent random variables. By the Central
Limit Theorem for random sums of independent random variables we thus know
that the total order picking time distribution will converge to a normal
distribution when both $k$ and $\lambda$ go to infinity at the same rate
(Shanthikumar and Sumita, 1984).

##  5 Conclusion

In this paper we have derived exact formulas for the Laplace-Stieltjes
Transform of the total order picking time under return routing for a number of
the most important storage policies, of which we discussed two in detail. The
formulas in this paper allow one to exactly determine the order picking time
distribution and can therefore be used to gain a deeper insight in the order
picking problem, as is highlighted in the numerical results.

The results in this paper exploit the independence between the aisles, which
is a result of the assumption that the order sizes follows a Poisson
distribution. While restrictive, these results can also be used to approximate
the order picking time distribution when the order size follows a mixed
Poisson distribution, by conditioning over the mixing parameter and
approximating the integral. This includes the case of a negative binomial
distribution (mixed Poisson in which the Poisson rate is Gamma distributed) or
the case where the rate of the Poisson process (for the demand) changes over
time.

This paper restricts the derivation to return routing in single- and 2-block
warehouses, but forms a foundation for the derivation of the order picking
time LST for multi-block warehouses or other routing policies. In particular,
the results can be generalized for midpoint, largest-gap and S-shaped routing
in single-block warehouses deploying the random storage policy. For more
advanced storage policies, the formulas quickly tend to become computationally
involved. Furthermore, warehouses with multiple blocks can be considered. Yet,
efficient routing policy definitions in multi-block warehouses are scarce and
often involve many steps and cases to take into account (see e.g. Roodbergen
et al., 2008).

Finally, we emphasize that the results of this paper can be used for a large
range of performance statistics. E.g., they can be used to obtain good bounds
on the tail probabilities of the order picking time, as well as approximations
of the order-lead time distribution, viewing the order picking time as service
time in a multi-server queue with $c$ pickers and applying a result of van
Hoorn and Tijms (1982). As a consequence, this paper allows for a more direct
analysis of statistics in warehouses, which enables a better comparison of
warehouse storage policies.

## Acknowledgements

The research of Tim Engels and Onno Boxma is partly funded by the NWO
Gravitation project NETWORKS, grant number 024.002.003.

## References

  * Abate and Whitt (1992) J. Abate and W. Whitt.  The Fourier-series method for inverting transforms of probability distributions.  _Queueing Systems_ , 10:5–88, 1992.
  * Boysen et al. (2019) N. Boysen, R. de Koster, and F. Weidinger.  Warehousing in the e-commerce era: A survey.  _European Journal of Operational Research_ , 277(2):396–411, 2019.  ISSN 03772217.  doi: 10.1016/j.ejor.2018.08.023.  URL https://doi.org/10.1016/j.ejor.2018.08.023.
  * Chew and Tang (1999) E. P. Chew and L. C. Tang.  Travel time analysis for general item location assignment in a rectangular warehouse.  _European Journal of Operational Research_ , 112(3):582–597, 1999.  ISSN 03772217.  doi: 10.1016/S0377-2217(97)00416-5.
  * Dijkstra and Roodbergen (2017) A. S. Dijkstra and K. J. Roodbergen.  Exact route-length formulas and a storage location assignment heuristic for picker-to-parts warehouses.  _Transportation Research Part E: Logistics and Transportation Review_ , 102:38–59, jun 2017.  ISSN 13665545.  doi: 10.1016/j.tre.2017.04.003.
  * Engels et al. (2022) T. Engels, I. Adan, O. Boxma, and J. Resing.  An exact analysis and comparison of manual picker routing heuristics, 2022\.
  * Gademann et al. (2001) A. J. R. M. Gademann, J. P. van den Berg, and H. H. van der Hoff.  An order batching algorithm for wave picking in a parallel-aisle warehouse.  _IIE Transactions (Institute of Industrial Engineers)_ , 33(5):385–398, 2001.  ISSN 15458830.  doi: 10.1080/07408170108936837.
  * Gupta et al. (2010) V. Gupta, M. Harchol-Balter, J. G. Dai, and B. Zwart.  On the inapproximability of M/G/K: Why two moments of job size distribution are not enough.  _Queueing Systems_ , 64(1):5–48, 2010.
  * Hall (1993) R. W. Hall.  Distance approximations for routing manual pickers in a warehouse.  _IIE Transactions (Institute of Industrial Engineers)_ , 25(4):76–87, 1993.  ISSN 15458830.  doi: 10.1080/07408179308964306.
  * Petersen and Aase (2004) C. G. Petersen and G. Aase.  A comparison of picking, storage, and routing policies in manual order picking.  _International Journal of Production Economics_ , 92(1):11–19, 2004.  ISSN 09255273.  doi: 10.1016/j.ijpe.2003.09.006.
  * Roodbergen and De Koster (2001) K. J. Roodbergen and R. De Koster.  Routing methods for warehouses with multiple cross aisles.  _International Journal of Production Research_ , 39(9):1865–1883, 2001.  ISSN 00207543.  doi: 10.1080/00207540110028128.
  * Roodbergen et al. (2008) K. J. Roodbergen, G. P. Sharp, and I. F. Vis.  Designing the layout structure of manual order picking areas in warehouses.  _IIE Transactions (Institute of Industrial Engineers)_ , 40(11):1032–1045, 2008.  ISSN 0740817X.  doi: 10.1080/07408170802167639.
  * Shanthikumar and Sumita (1984) J. Shanthikumar and U. Sumita.  A central limit theorem for random sums of random variables.  _Operations Research Letters_ , 3(3):153–155, 1984.  ISSN 0167-6377.  doi: https://doi.org/10.1016/0167-6377(84)90008-7.  URL https://www.sciencedirect.com/science/article/pii/0167637784900087.
  * Theys et al. (2010) C. Theys, O. Bräysy, W. Dullaert, and B. Raa.  Using a TSP heuristic for routing order pickers in warehouses.  _European Journal of Operational Research_ , 200(3):755–763, 2010.  ISSN 03772217.  doi: 10.1016/j.ejor.2009.01.036.  URL http://dx.doi.org/10.1016/j.ejor.2009.01.036.
  * Tijms (1986) H. C. Tijms.  _Stochastic Models: An Algorithmic Approach_.  John Wiley & Sons, Inc., Chichester, 1986.
  * Tompkins et al. (2011) J. Tompkins, J. White, Y. Bozer, and J. Tanchoco.  Facilities Planning.  _International Journal of Production Research_ , 49(24):7519–7520, 2011.  ISSN 0020-7543.  doi: 10.1080/00207543.2011.563164.
  * van Hoorn and Tijms (1982) M. H. van Hoorn and H. C. Tijms.  Approximations for the waiting time distribution of the M/G/c queue.  _Performance Evaluation_ , 2(1):22–28, 1982.  ISSN 01665316.  doi: 10.1016/0166-5316(82)90018-9.


# Optimizing performance of quantum operations with non-Markovian decoherence:
the tortoise or the hare?

Eoin P. Butler  School of Physics, Trinity College Dublin, College Green,
Dublin 2, Ireland  Gerald Fux  SUPA, School of Physics and Astronomy,
University of St Andrews, St Andrews KY16 9SS, United Kingdom  Abdus Salam
International Center for Theoretical Physics (ICTP), Strada Costiera 11, 34151
Trieste, Italy  Brendon W. Lovett  SUPA, School of Physics and Astronomy,
University of St Andrews, St Andrews KY16 9SS, United Kingdom  Jonathan
Keeling  SUPA, School of Physics and Astronomy, University of St Andrews, St
Andrews KY16 9SS, United Kingdom  Paul R. Eastham  School of Physics, Trinity
College Dublin, College Green, Dublin 2, Ireland

(28 March 2023)

###### Abstract

The interaction between a quantum system and its environment limits our
ability to control it and perform quantum operations on it. We present an
efficient method to find optimal controls for quantum systems coupled to non-
Markovian environments, by using the process tensor to compute the gradient of
an objective function. We consider state transfer for a driven two-level
system coupled to a bosonic environment, and characterize performance in terms
of speed and fidelity. We thus determine the best achievable fidelity as a
function of process duration. We show there is a trade-off between speed and
fidelity, and that slower processes can have higher fidelity by exploiting
non-Markovian effects.

The control of quantum systems using time-dependent Hamiltonians is crucial
for quantum technologies [1], enabling the implementation of state transfer
and gate operations. An important task is to establish how optimal performance
can be achieved for such processes. In an ideal closed quantum system perfect
operations are possible given sufficient time [2]. A speed limit arises
because physical Hamiltonians are bounded, so that energy-time uncertainty
gives a maximum rate of time-evolution and hence a minimum operation time.
Beyond this ideal case, however, other considerations arise. One is a desire
for reliable operation when precise control cannot be guaranteed; this can be
achieved at the expense of speed by using adiabatic processes [3]. Another is
the impact of decoherence and dissipation. In the standard Markovian
description of such processes, this results in effects which accumulate in
time. In most cases this means the highest achievable fidelity is then
obtained by operating at the speed limit. In the present Letter, we show that
this is not the case in non-Markovian systems, and that the highest fidelity
can be improved by slower operation.

To provide a concrete demonstration of the tradeoff between speed and fidelity
we use numerical optimal control to explore the achievable performance for a
non-Markovian system, consisting of a driven two-level system interacting with
a bosonic environment. Optimal control involves determining a set of time-
dependent control fields that maximize an objective function such as the
fidelity. Here, we show that this can be done effectively in a non-Markovian
system, using an extension of our previously introduced process tensor
approach [4] to efficiently compute the gradient of the objective function.
This allows us to repeatedly optimize over many hundreds of control parameters
for different process durations, and so find the achievable fidelity as a
function of the process duration. Our results reveal a marked improvement in
fidelity for process durations $T$ longer than a value $T_{0}$, set by the
speed-limit in the closed system [2]. We further compute a widely-used measure
of non-Markovianity, based on trace-distance [5], and show that the improved
fidelity coincides with an increase in non-Markovianity. Since many types of
device exhibit regimes of non-Markovian decoherence our results could enable
performance improvements across a range of quantum technologies [6], including
superconducting circuits [7, 8], spins [9], quantum dots [4, 10], color
centers in diamond [11], and cold atoms [12, 13, 14].

Our approach applies to a generic open quantum system, comprising a
Hamiltonian $H_{S}$, interacting with an environment (bath) with Hamiltonian
$H_{B}$, through a coupling $H_{SB}$. Models of this form are typically
treated under the assumptions of weak system-bath coupling (Born) and short
bath correlation time (Markov). This implies that information about the system
only flows away into the environment, and does not return, leading to simple
time-local evolution equations for the system reduced density matrix
$\rho_{S}$ [15]. Such equations have been used to establish performance limits
[16, 17, 18] using optimal control [19] in various problems. Studying optimal
control beyond the Born-Markov approximation has so far been difficult. One
approach is to expand the system to incorporate a few modes of the
environment, which are then treated exactly, with the remaining environment
modes providing Markovian damping. This approach has been used to study
controllability [20], optimal control [21, 22] and speed limits [18]. However,
it becomes intractable when one considers more than a few modes of the
environment. Another method involves computing the time-local dissipator
describing non-Markovian dynamics using lowest-order perturbation theory [23]
or the hierarchical equations-of-motion [24]. These techniques have been used
within optimal control to demonstrate fidelity increases, and are effective
when the environment spectral density can be approximated by a small number of
Lorentzians. Some works [25, 26, 27] have obtained optimal protocols under the
assumption that the dissipation is described by a fixed time-local dissipator,
such as that for a pure dephasing channel [28, 29]. An issue, however, is that
for optimal control one must consider the effect of the time-dependent control
fields on the dissipative processes [30].

Figure 1:  (a) Tensor network representation of the propagation of an initial
state over four time steps and calculation of the fidelity. Time increases
going up the diagram. The evolved density matrix is given by this diagram
without the target state node. The fidelity is obtained by joining the legs
indicated by arrows. (b) Derivative of the fidelity with respect to the
penultimate propagator and computation of the gradient by backpropagation (see
text).

To overcome these challenges, in this Letter we extend the process tensor
method described in our previous work [4] for solving the dynamics of non-
Markovian open systems. The process tensor (PT) [31, 32] is a mathematical
object that encodes the full influence of an environment onto a system, not
assuming any Markovian approximation. It is a multi linear map from the set of
all possible system control operation sequences to the resulting output
states. It is constructed by discretizing the time evolution into a series of
time steps. Generically the PT is a high rank tensor that grows exponentially
with the number of time steps. In many cases, however, it is possible to
systematically discard negligible correlations and express the PT in matrix
product operator (MPO) form [33, 34], allowing a numerically efficient
representation. For the case of a Gaussian bosonic environment one can
efficiently compute a PT in MPO form (PT-MPO) employing the method introduced
in [35, 36, 37, 38]. However, other methods for bosonic, fermionic, and spin
environments are available too [39, 38, 40, 41, 4, 42, 43, 44, 45, 46, 47, 48,
49]. Once computed, the PT can be contracted to time-evolve $\rho_{S}$ for any
$H_{S}$, as shown in Fig. 1(a). In this tensor network diagram [34] each node
represents a tensor, each leg represents an index, and connections between
legs correspond to contractions. The diagram is in Liouville space so that
operators such as $\rho_{S}$ are rank-1 objects, and maps between operators,
such as the propagators across each time step $U_{t}=e^{\Delta
t\mathcal{L}_{S}(t)}$ with $\mathcal{L}_{S}(t)\bullet=-i[H_{S}(t),\bullet]$,
are rank-2. The PT is the region in the dashed box, which can be contracted
with an initial density matrix and the Trotterized propagators [50] to obtain
$\rho_{S}$ at later times.

For optimal control we must define an objective function. In the following we
consider the example of a state transfer process, for which we use the
fidelity $\mathcal{F}$ [51] between the desired target state and the obtained
final state $\rho_{f}=\rho_{S}(t_{f})$. We restrict ourselves to pure target
states, with density matrix $\sigma=|\sigma\rangle\langle\sigma|$, so that
$\mathcal{F}=\langle\sigma|\rho_{f}|\sigma\rangle$. In Liouville space this is
a scalar formed from the vectors representing the final density matrix and the
transposed target density matrix. The fidelity then corresponds to the network
shown in Fig. 1(a) with the two legs marked with arrows joined.

An optimal protocol is found by numerically maximizing $\mathcal{F}$ over $N$
control parameters which determine the system propagators at each discrete
time step. Such numerical optimization becomes significantly faster if one can
efficiently calculate the gradient of $\mathcal{F}$ with respect to the
$N$-dimensional vector of control parameters, $\nabla\mathcal{F}$. A naive
calculation of $\nabla\mathcal{F}$ requires of order $N$ evaluations of
$\mathcal{F}$, strongly limiting the size of problems that can be treated.
Such gradients can however be efficiently obtained using the adjoint method
[52], which has been applied to unitary or Markovian dynamics [53]. Crucially,
one may see that the tensor network representation of $\mathcal{F}$ in Fig.
1(a) leads immediately to the form of adjoint method required for non-
Markovian dynamics. This diagram represents the fidelity as a multilinear
functional of the propagators; thus the derivative with respect to the
propagator at a given time step is the same diagram with the node for that
propagator removed, as shown in Fig. 1(b). Furthermore, all such derivatives
can be constructed from two sets of tensors found as follows. The first set
come from evaluating the network in the typical way, proceeding from the
initial state forwards in time, and then storing the two-legged tensors
obtained after each time step. The second set of tensors are found by starting
from the target state and propagating backwards in time to a given time step.
One can then contract one index in each pair of stored tensors to obtain the
derivative of $\mathcal{F}$ with respect to all the propagators—and hence the
gradient with respect to any and all choices of controls in $H_{S}$ [32]. The
key difference compared to use of the adjoint method in unitary or Markovian
dynamics is that we propagate a rank-2 tensor rather than the state
$\rho_{S}$. The additional index in this tensor encodes information about the
past and future dynamics and enables optimization of non-Markovian evolutions.

To illustrate the general principle, we consider an example of optimal state
transfer in a driven two-level system. Our system Hamiltonian is
$H_{S}=h_{x}(t)s_{x}+h_{z}(t)s_{z}$, where the fields $h_{x,z}(t)$ are the
controls. We treat them as piecewise constant with values $h_{x,z}(t_{n})$ at
the time steps $t_{n}=n\Delta t$. We seek to determine the control fields
which steer the state $\rho_{S}(t)$ from a prescribed initial state (the state
with $s_{x}=-1$) to a target state (the state $s_{x}=+1$). The bath is a set
of harmonic oscillators, and the coupling is
$H_{SB}=s_{z}\sum_{q}(g_{q}b_{q}+g_{q}^{*}b_{q}^{\dagger})$. This model is
known as the spin-boson model, and it describes many different physical
systems [54, 55]. For definiteness, we suppose it represents an optical
transition on a quantum dot driven by a laser pulse [56], in which case
$h_{z}$ and $h_{x}$ are the time-dependent detuning and Rabi frequency
respectively. In the quantum dot the bosonic modes are acoustic phonon modes,
for which we use the super-Ohmic spectral density
$J(\omega)=2\alpha\omega^{3}\omega_{c}^{-2}\exp\left[-\omega^{2}/\omega_{c}^{2}\right]$
with $\alpha=0.126$ and $\omega_{c}=3.04\,\text{ps}^{-1}$ [57, 58], and use a
bath temperature of 5 K.

Figure 2: (a) Infidelity, $1-\mathcal{F}$, for state transfer as a function of
process duration with bounded control fields. Results are shown for the
optimal protocol of the non-Markovian spin-boson model (control A) and for
driving with a constant field (control B). The latter is an optimal protocol
in the closed system. Blue/solid: infidelity of control A in the spin-boson
model. Orange/dashed: infidelity of control B in the spin-boson model.
Green/dot-dashed: infidelity of control B in the closed system. (b) Non-
Markovianity measure as a function of process duration in the spin-boson model
for controls A (blue/solid) and B (orange/dashed). The vertical dotted lines
are the speed limit time $T_{0}$ of the closed system.

We have determined optimal controls $h_{x}(t_{n}),h_{z}(t_{n})$ by numerically
maximizing $\mathcal{F}$ using the L-BFGS-B algorithm [59], with the gradient
and fidelity computed with the PT. We explore the impact of speed on fidelity
by performing this optimization for different process durations $T$, with
bounds on the fields, $|h_{x}|\leq h_{x}^{\mathrm{max}}=5\,\text{ps}^{-1}$ and
$|h_{z}|\leq h_{z}^{\mathrm{max}}=1\,\text{ps}^{-1}$, to consider the
realistic case in which the Hamiltonian, and therefore the speed of the
unitary evolution, is restricted. We note that in the unitary case, the speed
limit for state transfer has been identified from the rate of convergence of
the Krotov algorithm [16]. Here we consider a fully non-Markovian dissipative
system, and use numerically converged optimization to identify the speed limit
from the behavior of the fidelity against process duration. Figure 2(a) shows
the resulting infidelity, $1-\mathcal{F}$ for this optimized protcol, denoted
“Control A”.

It is interesting to compare the optimized results with those for the
protocols which are optimal for unitary dynamics in the closed system [60],
denoted “Control B”, which we used as the starting point for the optimization.
These protocols can be understood by noting that for unitary dynamics one can
achieve the state transfer with infidelity zero by setting $h_{x}$=0 and
requiring an area $\pi$ for the time integral of the field $h_{z}(t)$. Such a
protocol is optimal if it satisfies the constraint
$|h_{z}(t)|<h_{z}^{\mathrm{max}}$ for all times, which is possible when $T$ is
greater than the speed limit time $T_{0}=\pi/h_{z}^{\mathrm{max}}$ (which
saturates both the Mandelstam-Tamm and Margolus-Levitin bounds [2]). Among
protocols satisfying the above condition, we choose $h_{z}(t)=\pi/T$. For
$T<T_{0}$ the state transfer is not fully achievable and the optimal protocol
is that with $h_{z}(t)=\pm h_{z}^{\mathrm{max}}$. Thus, in both regimes, we
have an optimal protocol with constant fields. As can be seen in Fig. 2(a),
the infidelity of this protocol (Control B) for the unitary dynamics smoothly
decreases as the duration $T$ increases from zero, before saturating at zero
for $T>T_{0}$. Applying this same protocol in the open system gives a similar
overall behavior, with the key difference being that the saturated value of
the infidelity at large $T$ is now non-zero. This behavior differs from that
with Markovian decoherence, where slower processes would be expected to
produce higher infidelity, i.e. growing with $T$. In the “Control B” protocol,
we have $h_{z}(t)=h_{z}$ and $h_{x}=0$, so that the model is the exactly
solvable independent-boson model [61, 62], which has non-Markovian
decoherence. The constant infidelity for $T>T_{0}$ comes from the decoherence
function of the independent-boson model, which approaches a non-zero constant
for times $T\gg\tau_{b}\sim 1/\omega_{c}$.

Figure 2 shows that the optimization increases the performance for all process
durations, as one would expect, but the increase becomes marked when
$T>T_{0}$. This is because the fidelity in this regime is limited not by the
distance to the target state and the speed of unitary evolution under $H_{S}$,
but by the decoherence and dissipation produced by the bath. The improved
fidelity thus corresponds to suppressing decoherence. An approach to
suppressing decoherence one might consider is Dynamical Decoupling [62] (DD).
Such a protocol uses control fields to produce dynamics for the system that
averages away the effects of the bath, which requires that the timescale of
the system is much shorter than that of the bath $\tau_{b}$. This regime is
not accessible in our optimization due to the bounds on the control fields.
The standard DD sequence for the process we consider would be a train of $\pi$
pulses separated by a time $\tau\ll\tau_{b}$ and each with duration
$\tau_{p}\ll\tau$; this would produce rapid spin-flips and so average away
$H_{SB}\propto s_{z}$. Thus standard DD requires a field
$|h_{x}|\gg\pi/\tau_{b}$, which with $\tau_{b}\sim 1/\omega_{c}\sim
0.3\,\mathrm{ps}$ is greater than $h_{x}^{\mathrm{max}}$.

To explore the mechanism that is giving the improvement in fidelity we compute
a measure of non-Markovianity. There are many such interrelated measures [63,
55, 64] including, among others, those based on divisiblity and complete-
positivity, and those based on information backflow. We choose that introduced
by Breuer, Lane and Piilo [5], formed from the trace distance between pairs of
states $D_{12}=\Tr|\rho_{1}(t)-\rho_{2}(t)|/2$, where $\rho_{1}(t)$ and
$\rho_{2}(t)$ are obtained by time-evolving initial states $\rho_{1}$ and
$\rho_{2}$ respectively. The non-Markovianity is

$\mathcal{N}=\max_{\rho_{1},\rho_{2}}\int\frac{dD_{12}}{dt}dt,$

where the integral extends over regions where the integrand is positive, and
the maximum is over all pairs of orthogonal states on the surface of the Bloch
sphere (which is the same as the maximum over all pairs of initial states, see
[65]). The two curves in Fig. 2(b) show the non-Markovianity as a function of
process duration for the optimal controls (solid) and the constant field
(dashed). For the constant field case, control B, the non-Markovianity
increases from zero as the process duration increases from zero, becoming
constant for $T\gtrapprox 2\,\mathrm{ps}$. This corresponds to the build-up of
correlations due to the system-bath coupling in the independent boson model,
and is unaffected by the system dynamics since $[H_{S},H_{SB}]=0$ when
$h_{x}=0$. The non-Markovianity with the optimal control A is higher, but
initially shows a similar increase and saturation as the duration $T$
increases from zero. However, there is then a marked feature at $T\approx
T_{0}$, beyond which the non-Markovianity increases rapidly with process
duration before reaching a new and much higher level. This suggests that the
improvement in fidelity for long process durations occurs because the
optimization increases the degree of non-Markovianity of the map, restoring
some of the information lost to the environment during the earlier parts of
the dynamics by the end of the process.

Figure 3:  (a) Driving fields $h_{x}(t)$ (green) and $h_{z}(t)$ (orange) for
the optimal state transfer process, control A, of duration 4.0 ps. (b)
Dynamics of the components of the Bloch vector (left axis, green/solid:
$\expectationvalue{\sigma_{x}}$ , blue/dot-dashed:
$\expectationvalue{\sigma_{y}}$, orange/dashed:
$\expectationvalue{\sigma_{z}}$), and length of the Bloch vector (right axis,
black/solid). The length of the Bloch vector for control B is also shown
(right axis, black/dashed).

Figure 3 shows the optimized control fields (top panel), and the dynamics of
the Bloch vector $\langle\sigma_{x,y,z}\rangle$, for the optimal protocol,
control A, at $T=4.0\,\mathrm{ps}$. For this duration, the optimal $h_{z}$
remains constant, while the optimized $h_{x}$ has an overall slope and
oscillations of varying magnitude, which saturate the bounds at the beginning
and end of the process. The resulting dynamics of the Bloch vector moves down
below the equator during the control and develops weak oscillations. The
clearest effect of the optimization is seen in the length of the Bloch vector,
which is shown for both protocols A and B. In the latter case (constant field)
the length of the Bloch vector is given by the decoherence function of the
independent boson model, which goes from one, for zero time, to a constant
non-zero value at late times. The optimal control A clearly avoids this decay
and furthermore shows an increase in the Bloch-vector length at the very end
of the pulse, consistent with a decoherence suppression mechanism. The control
in $h_{x}$ is very different from that for a standard DD sequence as described
above, although it does act to produce some oscillations in the Bloch vector
which may have a similar averaging effect. The oscillations in the Bloch
vector and control fields become more pronounced (and saturate the bound) as
we move to longer process durations. The behavior shown is otherwise
representative of that we find for other process durations [32]. The
consistency in the form of these solutions suggests that there is a unique
global optimum that is found by our optimization process for durations near
the speed limit. This is expected to break down for slower processes [66].

In conclusion, we have presented a method for computing optimal controls in
general non-Markovian systems. It extends the PT method [4] so as to
efficiently compute gradients, so permitting optimization over large numbers
of control parameters. The method takes full account of the changes in the
dissipator produced by the control fields, allowing it to discover protocols
which exploit non-Markovianity for improved performance. We have used it to
investigate the performance of a state transfer process in a time-dependent
spin-boson model and find that, contrary to the case of a Markovian system,
slow processes can perform better than fast ones. The identified controls
increase the degree of non-Markovianity during the state transfer, as measured
by the growth of trace-distance, and suppress the decoherence, as measured by
the decay of the Bloch vector length. These results show that performance
improvements are possible across the many qubit implementations subject to
non-Markovian decoherence, and could be identified using the methods described
here.

We acknowledge discussions with F. Binder. E. P. B. acknowledges support from
the Irish Research Council (GOIPG/2019/1871). G. E. F. acknowledges support
from EPSRC (EP/L015110/1). B. W. L. and J. K. acknowledge support from EPSRC
(EP/T014032/1). P. R. E. acknowledges support from Science Foundation Ireland
(21/FFP-P/10142).

## References

  * Glaser  _et al._ [2015] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Training Schrödinger’s cat: Quantum optimal control, Eur. Phys. J. D 69, 279 (2015).
  * Deffner and Campbell [2017] S. Deffner and S. Campbell, Quantum speed limits: From Heisenberg’s uncertainty principle to optimal quantum control, J. Phys. A: Math. Theor. 50, 453001 (2017).
  * Torrontegui  _et al._ [2013] E. Torrontegui, S. Ibáñez, S. Martínez-Garaot, M. Modugno, A. del Campo, D. Guéry-Odelin, A. Ruschhaupt, X. Chen, and J. G. Muga, Shortcuts to adiabaticity, in  _Advances in Atomic, Molecular, and Optical Physics_, Vol. 62, edited by E. Arimondo, P. R. Berman, and C. C. Lin (Academic Press, 2013) pp. 117–169.
  * Fux  _et al._ [2021] G. E. Fux, E. P. Butler, P. R. Eastham, B. W. Lovett, and J. Keeling, Efficient Exploration of Hamiltonian Parameter Space for Optimal Control of Non-Markovian Open Quantum Systems, Phys. Rev. Lett. 126, 200401 (2021).
  * Laine  _et al._ [2010] E.-M. Laine, J. Piilo, and H.-P. Breuer, Measure for the Non-Markovianity of Quantum Processes, Phys. Rev. A 81, 062115 (2010).
  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, The quantum technologies roadmap: A European community view, New J. Phys. 20, 080201 (2018).
  * Fischer  _et al._ [2019] J. Fischer, D. Basilewitsch, C. P. Koch, and D. Sugny, Time-optimal control of the purification of a qubit in contact with a structured environment, Phys. Rev. A 99, 033410 (2019).
  * Goerz  _et al._ [2017] M. H. Goerz, F. Motzoi, K. B. Whaley, and C. P. Koch, Charting the circuit QED design landscape using optimal control theory, npj Quantum Inf 3, 37 (2017).
  * Dolde  _et al._ [2014] F. Dolde, V. Bergholm, Y. Wang, I. Jakobi, B. Naydenov, S. Pezzagna, J. Meijer, F. Jelezko, P. Neumann, T. Schulte-Herbrüggen, J. Biamonte, and J. Wrachtrup, High-fidelity spin entanglement using optimal control, Nat. Commun. 5, 3371 (2014).
  * Wilson-Rae and Imamoğlu [2002] I. Wilson-Rae and A. Imamoğlu, Quantum dot cavity-QED in the presence of strong electron-phonon interactions, Phys. Rev. B 65, 235311 (2002).
  * Norambuena  _et al._ [2020] A. Norambuena, J. R. Maze, P. Rabl, and R. Coto, Quantifying phonon-induced non-Markovianity in color centers in diamond, Phys. Rev. A 101, 022110 (2020).
  * Choi  _et al._ [2014] T. Choi, S. Debnath, T. A. Manning, C. Figgatt, Z.-X. Gong, L.-M. Duan, and C. Monroe, Optimal Quantum Control of Multimode Couplings between Trapped Ion Qubits for Scalable Entanglement, Phys. Rev. Lett. 112, 190502 (2014).
  * Bason  _et al._ [2012] M. G. Bason, M. Viteau, N. Malossi, P. Huillery, E. Arimondo, D. Ciampini, R. Fazio, V. Giovannetti, R. Mannella, and O. Morsch, High-fidelity quantum driving, Nat. Phys 8, 147 (2012).
  * Treutlein  _et al._ [2006] P. Treutlein, T. W. Hänsch, J. Reichel, A. Negretti, M. A. Cirone, and T. Calarco, Microwave potentials and optimal control for robust quantum gates on an atom chip, Phys. Rev. A 74, 022312 (2006).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Caneva  _et al._ [2009] T. Caneva, M. Murphy, T. Calarco, R. Fazio, S. Montangero, V. Giovannetti, and G. E. Santoro, Optimal Control at the Quantum Speed Limit, Phys. Rev. Lett. 103, 240501 (2009).
  * del Campo  _et al._ [2013] A. del Campo, I. L. Egusquiza, M. B. Plenio, and S. F. Huelga, Quantum Speed Limits in Open System Dynamics, Phys. Rev. Lett. 110, 050403 (2013).
  * Deffner and Lutz [2013] S. Deffner and E. Lutz, Quantum Speed Limit for Non-Markovian Dynamics, Phys. Rev. Lett. 111, 010402 (2013).
  * Koch [2016] C. P. Koch, Controlling open quantum systems: Tools, achievements, and limitations, J. Phys.: Condens. Matter 28, 213001 (2016).
  * Reich  _et al._ [2015] D. M. Reich, N. Katz, and C. P. Koch, Exploiting Non-Markovianity for Quantum Control, Sci. Rep. 5, 12430 (2015).
  * Rebentrost  _et al._ [2009] P. Rebentrost, I. Serban, T. Schulte-Herbrüggen, and F. K. Wilhelm, Optimal Control of a Qubit Coupled to a Non-Markovian Environment, Phys. Rev. Lett. 102, 090401 (2009).
  * Floether  _et al._ [2012] F. F. Floether, P. de Fouquieres, and S. G. Schirmer, Robust quantum gates for open systems via optimal control: Markovian versus non-Markovian dynamics, New J. Phys. 14, 073023 (2012).
  * Hwang and Goan [2012] B. Hwang and H.-S. Goan, Optimal control for non-Markovian open quantum systems, Phys. Rev. A 85, 032321 (2012).
  * Mangaud  _et al._ [2018] E. Mangaud, R. Puthumpally-Joseph, D. Sugny, C. Meier, O. Atabek, and M. Desouter-Lecomte, Non-Markovianity in the optimal control of an open quantum system described by hierarchical equations of motion, New J. Phys. 20, 043050 (2018).
  * Ohtsuki [2003] Y. Ohtsuki, Non-Markovian effects on quantum optimal control of dissipative wave packet dynamics, J. Chem. Phys. 119, 661 (2003).
  * Ying-Hua  _et al._ [2016] J. Ying-Hua, H. Ju-ju, H. Jian-Hua, and K. Qiang, Optimal control-based states transfer for non-Markovian quantum system, Physica E Low Dimens. Syst. Nanostruct. 81, 77 (2016).
  * Cui  _et al._ [2008] W. Cui, Z. R. Xi, and Y. Pan, Optimal decoherence control in non-Markovian open dissipative quantum systems, Phys. Rev. A 77, 032117 (2008).
  * Chin  _et al._ [2012] A. W. Chin, S. F. Huelga, and M. B. Plenio, Quantum Metrology in Non-Markovian Environments, Phys. Rev. Lett. 109, 233601 (2012).
  * Bylicka  _et al._ [2014] B. Bylicka, D. Chruściński, and S. Maniscalco, Non-Markovianity and reservoir memory of quantum channels: A quantum information theory perspective, Sci. Rep. 4, 5720 (2014).
  * Addis  _et al._ [2016] C. Addis, E.-M. Laine, C. Gneiting, and S. Maniscalco, Problem of coherent control in non-Markovian open quantum systems, Phys. Rev. A 94, 052117 (2016).
  * Pollock  _et al._ [2018a] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Non-Markovian quantum processes: Complete framework and efficient characterization, Phys. Rev. A 97, 012127 (2018a).
  * [32] See supplemental information at [$\ldots$] for additional details of process tensors, the method for computation of gradients, and results for different process durations.
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Matrix product density operators: Simulation of finite-temperature and dissipative systems, Phys. Rev. Lett. 93, 207204 (2004).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys. (N. Y.) 349, 117 (2014).
  * Makri and Makarov [1995a] N. Makri and D. E. Makarov, Tensor propagator for iterative quantum time evolution of reduced density matrices. II. Numerical methodology, J. Chem. Phys. 102, 4611 (1995a).
  * Makri and Makarov [1995b] N. Makri and D. E. Makarov, Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory, J. Chem. Phys. 102, 4600 (1995b).
  * Strathearn  _et al._ [2018] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, Efficient non-Markovian quantum dynamics using time-evolving matrix product operators, Nat. Commun. 9, 3322 (2018).
  * Jørgensen and Pollock [2019] M. R. Jørgensen and F. A. Pollock, Exploiting the Causal Tensor Network Structure of Quantum Processes to Efficiently Simulate Non-Markovian Path Integrals, Phys. Rev. Lett. 123, 240602 (2019).
  * Bañuls  _et al._ [2009] M. C. Bañuls, M. B. Hastings, F. Verstraete, and J. I. Cirac, Matrix product states for dynamical simulation of infinite chains, Phys. Rev. Lett. 102, 240603 (2009).
  * Strathearn [2020] A. Strathearn,  _Modelling Non-Markovian Quantum Systems Using Tensor Networks_, Springer Theses (Springer International Publishing, Cham, 2020).
  * Lerose  _et al._ [2021] A. Lerose, M. Sonner, and D. A. Abanin, Influence matrix approach to many-body floquet dynamics, Phys. Rev. X 11, 021040 (2021).
  * Cygorek  _et al._ [2022] M. Cygorek, M. Cosacchi, A. Vagov, V. M. Axt, B. W. Lovett, J. Keeling, and E. M. Gauger, Simulation of open quantum systems by automated compression of arbitrary environments, Nat. Phys. 18, 662 (2022).
  * Sonner  _et al._ [2021] M. Sonner, A. Lerose, and D. A. Abanin, Influence functional of many-body systems: Temporal entanglement and matrix-product state representation, Ann. Phys. (N. Y.) 435, 168677 (2021), special issue on Philip W. Anderson.
  * Ye and Chan [2021] E. Ye and G. K.-L. Chan, Constructing tensor network influence functionals for general quantum dynamics, J. Chem. Phys. 155, 044104 (2021).
  * White  _et al._ [2022] G. White, F. Pollock, L. Hollenberg, K. Modi, and C. Hill, Non-markovian quantum process tomography, PRX Quantum 3, 020344 (2022).
  * Thoenniss  _et al._ [2022a] J. Thoenniss, A. Lerose, and D. A. Abanin, Non-equilibrium quantum impurity problems via matrix-product states in the temporal domain, arXiv:2205.04995 [cond-mat.str-el] (2022a).
  * The TEMPO collaboration [2022] The TEMPO collaboration, OQuPy: A Python 3 package to efficiently compute non-Markovian open quantum systems. (2022).
  * Thoenniss  _et al._ [2022b] J. Thoenniss, M. Sonner, A. Lerose, and D. A. Abanin, An efficient method for quantum impurity problems out of equilibrium, arXiv:2211.10272 [cond-mat.str-el] (2022b).
  * Ng  _et al._ [2022] N. Ng, G. Park, A. J. Millis, G. K. Chan, and D. R. Reichman, Real time evolution of anderson impurity models via tensor network influence functionals, arXiv:2211.10430 [cond-mat.str-el] (2022).
  * Trotter [1959] H. F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. 10, 545 (1959).
  * Benenti  _et al._ [2018] G. Benenti, G. Casati, D. Rossini, and G. Strini,  _Principles Of Quantum Computation And Information: A Comprehensive Textbook_ (World Scientific, 2018).
  * Plessix [2006] R.-E. Plessix, A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, Geophys. J. Int. 167, 495 (2006).
  * Mishra  _et al._ [2021] S. D. Mishra, R. Trivedi, A. H. Safavi-Naeini, and J. Vučković, Control Design for Inhomogeneous-Broadening Compensation in Single-Photon Transducers, Phys. Rev. Applied 16, 044025 (2021).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).
  * de Vega and Alonso [2017] I. de Vega and D. Alonso, Dynamics of non-Markovian open quantum systems, Rev. Mod. Phys. 89, 015001 (2017).
  * Ramsay [2010] A. J. Ramsay, A review of the coherent optical control of the exciton and spin states of semiconductor quantum dots, Semicond. Sci. Technol. 25, 103001 (2010).
  * Ramsay  _et al._ [2010] A. J. Ramsay, T. M. Godden, S. J. Boyle, E. M. Gauger, A. Nazir, B. W. Lovett, A. M. Fox, and M. S. Skolnick, Phonon-Induced Rabi-Frequency Renormalization of Optically Driven Single InGaAs/GaAs Quantum Dots, Phys. Rev. Lett. 105, 177402 (2010).
  * Lüker  _et al._ [2012] S. Lüker, K. Gawarecki, D. E. Reiter, A. Grodecka-Grad, V. M. Axt, P. Machnikowski, and T. Kuhn, Influence of acoustic phonons on the optical control of quantum dots driven by adiabatic rapid passage, Phys. Rev. B 85, 121302 (2012).
  * Virtanen et al. [SciPy 1.0 Contributors] P. Virtanen et al. (SciPy 1.0 Contributors), SciPy 1.0: Fundamental algorithms for scientific computing in python, Nat. Methods 17, 261 (2020).
  * Hegerfeldt [2013] G. C. Hegerfeldt, Driving at the Quantum Speed Limit: Optimal Control of a Two-Level System, Phys. Rev. Lett. 111, 260501 (2013).
  * Mahan [2013] G. D. Mahan,  _Many-Particle Physics_ (Springer Science & Business Media, 2013).
  * Viola  _et al._ [1999] L. Viola, E. Knill, and S. Lloyd, Dynamical Decoupling of Open Quantum Systems, Phys. Rev. Lett. 82, 2417 (1999).
  * Pollock  _et al._ [2018b] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Operational Markov Condition for Quantum Processes, Phys. Rev. Lett. 120, 040405 (2018b).
  * Li  _et al._ [2018] L. Li, M. J. W. Hall, and H. M. Wiseman, Concepts of quantum non-Markovianity: A hierarchy, Phys. Rep. 759, 1 (2018).
  * Wißmann  _et al._ [2012] S. Wißmann, A. Karlsson, E.-M. Laine, J. Piilo, and H.-P. Breuer, Optimal state pairs for non-Markovian quantum dynamics, Phys. Rev. A 86, 062108 (2012).
  * Larocca  _et al._ [2018] M. Larocca, P. M. Poggi, and D. A. Wisniacki, Quantum control landscape for a two-level system near the quantum speed limit, J. Phys. A: Math. Theor. 51, 385305 (2018).