\section{Introduction}
Systems of polynomial equations appear in many areas of mathematics, as well as in many applications in the sciences and engineering.
In physics and chemistry the geometry of molecules is often modelled with algebraic constraints on the distance or the angle between atoms.
In kinematics the relation between robot joints is defined by polynomial equations.
In systems biology the steady-state equations for many bio-chemical reaction networks are algebraic equations.
A central task in all those applications is computing the isolated zeros of a system of polynomials.

The study of zeros of polynomial systems is at the heart of algebraic geometry.
The field of \emph{computational algebraic geometry} is often associated with symbolic computations based on Gr\"obner bases.
But over the last thirty years \emph{numerical algebraic geometry} (NAG) \cite{Sommese:Wampler:2005} emerged as an alternative; enabling us to solve problems infeasible with symbolic methods. The algorithmic framework in NAG is \emph{numerical homotopy continuation}.
Several implementations of this are available: \texttt{Bertini} \cite{Bertini}, \texttt{Hom4PS-3} \cite{Hom4PSArticle}, \texttt{HomotopyContinuation.jl} \cite{HC.jl}, \texttt{NAG4M2} \cite{NumericalAlgebraicGeometryArticle} and \texttt{PHCpack} \cite{PHCpack}.
The first and the third author are the developers of \texttt{HomotopyContinuation.jl}.

Hauenstein and Sottile remark in \cite{Hauenstein:Sottile:2012} that while all of these softwares ``routinely and reliably solve systems of polynomial equations with dozens of variables having thousands of solutions'', they have the shortcoming that ``the output is not certified'' and that ``this restricts their use in some applications, including those in pure mathematics''.
To remedy this, Hauenstein and Sottile developed the software \texttt{alphaCertified} \cite{Hauenstein:Sottile:2012}.
It can rigorously certify that Newton's method, starting at a given numerical approximation, converges quadratically to a true zero by using Smale's $\alpha$-theory \cite{Smale:1986}.
Hauenstein and Sottile's contribution to numerical algebraic geometry was a milestone.
Yet, \texttt{alphaCertified} produces rigorous certificates using expensive rational arithmetic.
This turns the big advantage of numerical computations, namely that they are fast, upside-down and makes certification of large problems prohibitively expensive.

Up to this point, the majority of researchers in applied algebraic geometry were kept from using numerical methods, because certification was too expensive and because without certification numerical methods can't be used for proofs. With our article we want to initiate a paradigm shift\footnote{We do \emph{not} claim to introduce certification methods based on interval arithmetic in general. This is a well established technique. What we do claim is the introduction of interval arithmetic methods to numerical algebraic geometry with a concrete implementation.} in numerical algebraic geometry: with a fast implementation certification becomes the default and is not just an option. This enables the extensive use of numerical methods for rigorous proofs.

\subsection{Contribution}
Our contribution to the field of computational and applied algebraic geometry is an extremely fast and easy-to-use implementation of a certification method. This implementation outperforms \texttt{alphaCertified} by several orders of magnitude.
It makes the certification of solutions often a matter of seconds and not hours or days.
This leap in performance is the basis for the proposed paradigm shift in numerical algebraic geometry where certification is the default and not an option.

Starting from version 2.1 \texttt{HomotopyContinuation.jl} has a function \texttt{certify}\footnote{The technical documentation is available at\\
\url{https://www.juliahomotopycontinuation.org/HomotopyContinuation.jl/stable/certification}
}.
The function \texttt{certify} takes as input a \emph{square polynomial system} $F$ and a numerical approximation of a complex zero $x\in\mathbb C^n$ (or a list of zeros). If the output says ``certified'', then this is a rigorous proof that a solution of $F=0$ is near $x$. If the output says ``not certified'', then this does not necessarily mean that there is no zero near $x$, just that the method couldn't find one.  Figure~\ref{fig1} shows an example of \texttt{certify}. See also the example \cite{bacillus-subtilis2020} on \url{https://www.juliahomotopycontinuation.org}.

We combine interval arithmetic and Krawczyk's method with numerical algebraic geometry to rigorously certify solutions to square systems of polynomial equations.
In technical terms, our implementation returns \emph{strong interval approximate zeros}.
We introduce this notion in Definition~\ref{def_interval_zero} below.
The strong interval approximate zero consists of a box in $\mathbb C^n$, which contains a unique true zero of the polynomial system. If the input is a list of zeros, the routine returns a list of distinct strong interval approximate zeros.

Therefore, our method can be used to \emph{prove} hard lower bounds on the number of zeros of a polynomial system.
Combined with theoretical upper bounds this can constitute rigorous mathematical proofs on the number of zeros of such systems.

In addition, if the given polynomial system is real, we give a certificate whether the certified zero is a real zero.
The returned boxes may also be used to verify if a real zero is positive real.
Therefore, our method can also be used to prove lower bounds on the number of real and positive real zeros of a polynomial system.

It is also possible to give a square system of rational functions as input to our implementation. Although this article is formulated in terms of polynomial systems, Krawczyk's method also applies to square systems of rational functions. Consequently, all statements about using our implementation for proofs are equally valid of square systems of rational functions. Nevertheless, we think that the focus on polynomial systems simplifies the exposition.

\subsection{Comparison to previous works}

There are other implementations of certification methods using Krawczyk's method and interval arithmetic, e.g., the commercial \texttt{MATLAB} package \texttt{INTLAB} \cite{Rump1999}, the \texttt{Macaulay2} package \texttt{NumericalCertification} \cite{LeeM2}, and the \texttt{Julia} package \texttt{IntervalRootFinding.jl} \cite{BS}.

Compared to \texttt{INTLAB}, the source code of our implementation is freely available and can be verified by anyone.
Additionally, \texttt{INTLAB} doesn't support the use of arbitrary precision interval arithmetic which limits its capability to certify badly conditioned solutions.
\texttt{NumericalCertification}, as of version 1.0, takes as input not the numerical approximation of a complex zero $x\in\mathbb C^n$, but instead a box $I$ in $\mathbb{C}^n$. Then, \texttt{NumericalCertification} attempts to certify that interval $I$ is a strong interval approximate.
The process of going from a numerical approximation $x$ to a good candidate interval~$I$ needs particular care, as illustrated in Section \ref{sec:implementation_details}.
\texttt{INTLAB} and \texttt{NumericalCertification} also both require manual work to obtain a list of all distinct distinct strong interval approximate zeros. The package \texttt{IntervalRootFinding.jl} can find all zeros of a multivariate function inside a given box in $\mathbb R^n$, whereas our implementation works in $\mathbb{C}^n$ and additionally certifies reality of zeros; see Section \ref{sec:reality}.

Our contribution is a significant advancement over these previous works, since it not only provides an implementation of Krawcyzk's method, but also combines it with the necessary tools and techniques to deliver an easy to use and robust certification routine.

\subsection{Acknowledgements}
We thank Pierre Lairez for a discussion that initiated this project. We also thank him for several helpful subsequent discussions on the topic.

\subsection{Outline}
The rest of this article is organized as follows: In the next section we demonstrate our implementation on three applications. This shows both the speed of the implementation, and how it can be used for proofs.
We discuss the details of our implementation in Section \ref{sec:implementation_details}.
For completeness, we include a short introduction to interval arithmetic in Section \ref{sec:IA}, and a proof of Krawczyk's method in Section \ref{sec:krawczyk}.

\section{Applications}\label{sec:experiments}
Certification methods are useful when one wants to prove statements on the number of zeros, the number of real zeros, or the number of positive real zeros of a polynomial system.
When determining these numbers
it is often most challenging to obtain lower bounds.
Methods from algebraic geometry provide upper bounds,
and applying our certification method can give a proof that the upper bound for the number of zeros is attained. A computation with our certification method always reveals lower bounds.

In the following, we discuss the application of our implementation in three different fields.
All reported timings were obtained on an desktop computer with a 3.4 GHz processor running \texttt{Julia}~1.5.2 \cite{BEKV} and \texttt{HomotopyContinuation.jl} version 2.2.2.

\subsection{3264 real conics}\label{subsec:3264}
We demonstrate how certification methods in numerical algebraic geometry allow to proof theorems in algebraic geometry.

In \cite{3264} we used \texttt{alphaCertified} to prove that a certain arrangement of five conics in the plane had 3264 real conics, which were simultaneously tangent to each of the five given conics.
Such an arrangement is called \emph{totally real}.
It was known before that such arrangements exist \cite{RTV97}, but an explicit instance was not known.
The fact that \texttt{alphaCertified} provides a proof for a totally real instance highlights the relevance of certification software in algebraic geometry.

The strategy for the computation is this.
The zeros of the system (12) in \cite{3264} give the coordinates of the 3264 conics which are tangent to five given conics.
We compute the zeros for the coordinates of the specific instance in \cite[Figure 2]{3264} using \texttt{HomotopyContinuation.jl}.
This is a numerical computation.
Therefore, it is inexact and cannot be used in a proof. Next, we take the inexact numerical zeros as starting points for our certification method.
If our implementation outputs that it has found a real certified zero, then this is an exact result and hence it is a proof that the zero is real.
This way we can prove that indeed all the 3264 conics for the instance in \cite[Figure 2]{3264} are real.
See also the proof of \cite[Proposition 1]{3264} for a more detailed discussion.

The certification with \texttt{alphaCertified} took us \emph{more than 36 hours}.
In contrast, our implementation certifies the reality and distinctness of the 3264 conics in \emph{less than three seconds}.

\begin{figure}
\begin{center}
\includegraphics[width = 0.7\textwidth]{3264_certification.png}
\end{center}
\caption{\label{fig1}
Screenshot from a \texttt{Julia} session, where we certify the 3264 real conics for the totally real arrangement from \cite{3264}. Here, \texttt{F} is the system of polynomials from (12) in \cite{3264}. The screenshot also demonstrates the simple syntax of our implementation.}
\end{figure}

\subsection{Numerical Synthesis of Six-Bar Linkages}
Now we demonstrate that the certification routine can cope with large problems. With our computation we improve a result from the literature.

We consider the kinematic synthesis of six-bar linkages that use eight prescribed accuracy points as described in \cite{Plecnik:McCarthy:2014}.
In this article, the authors derive the synthesis equations for six-bar linkages of the Watt II, Stephenson II, and Stephenson III type.
Additionally, in \cite[Eq. (35)]{Plecnik:McCarthy:2014} they construct a system of 22 polynomials in 22 unknowns and 224 parameters, that can be used as a start system in a parameter homotopy to solve the synthesis equations of all three considered six-bar linkage types.

The number of non-singular zeros of this generalized start system is reported as 92,736.
It was computed using \texttt{Bertini} and a multi-homogeneous start system.
To certify the reported count, we solved the generalized start system using the monodromy method \cite{monodromy} implementation in \texttt{HomotopyContinuation.jl}.
In our computation we obtained 92,752 non-singular zeros for a generic choice of the 224 parameters. These are sixteen \emph{more} than reported in \cite{Plecnik:McCarthy:2014}.
We certified this count using our certification routine and obtained 92,752 distinct strong interval approximate zeros.
Therefore, we have a certificate that the generalized system has in general (at least) 92,752 non-singular solution.
This establishes that the result in \cite{Plecnik:McCarthy:2014} undercounts the true number of solutions.
The certification needed only 38.34 seconds which underlines the scalability of the certification routine.

In Section \ref{sec_disjoint} below we discuss that part of the certification process is checking if the intervals are pairwise disjoint.
In this example there are over 4 billion such pairs, which underlines the need for having an efficient algorithm for comparing pairs.

\subsection{Stress response of Bacillus Subtilis}
In this section we demonstrate that our implementation certifies positivity of zeros. In many applications variables represent magnitudes, so that only positive real solutions are physically meaningful. If such zeros exist, our methods provides a rigorous proof for their existence, and in addition, it gives a certified interval in which the true zero is contained. This is of interest for researchers working at the intersection of algebraic geometry and (bio-)chemical reaction networks.

Our example is from biochemistry.
The environmental and energy stress response of the bacterium \emph{Bacillus subtilis} are modelled in \cite{NTI2016}.
The protein $\sigma_B$ is the focus of this paper.
It is responsible for activating a stress-response of the bacterium.
$\sigma_B$ belongs to the family of~$\sigma$ factors.
These are a type of so called transcription factors;
proteins which govern the expression of genes.

In \cite{NTI2016} regulatory networks are studied.
They consist of other proteins involved in feedback loops that influence the $\sigma$-factors.
Since there can be many possible reactants involved in many reactions, the resulting system of differential equations might be very complicated. The model for \emph{Bacillus subtilis} in \cite{NTI2016} is claimed to be backed up by experimental data. The activity of $\sigma_B$ is regulated by a network consisting of an anti-$\sigma$ factor $\mathrm{RsbW}$ and an anti-anti-$\sigma$ factor $\mathrm{RsbV}$.

In \cite{NTI2016} this biochemical reactions dynamical system is modelled by a system of differential equations in the 10 variables $w$, $w_\mathrm{2}$, $w_\mathrm{2v}$, $v$, $w_\mathrm{2v2}$, $v_P$, $\sigma_B$, $w_\mathrm{2\sigma B}$, $v_\mathrm{Pp}$ and $\mathrm{phos}$.
These represent the total amounts of $\sigma_B$, $\mathrm{RsbW}$, $\sigma$ factor $\mathrm{RsbV}$, and of various protein complexes formed by these components.
The variable $\mathrm{phos}$ measures the concentration of the phosphatase, which serves as a measure for the amount of stress the bacterium experiences.

With our implementation we can determine the steady states of the described dynamical system.
The vanishing of the differentials of each of the concentrations with respect to time is equivalent to the vanishing of the ten polynomials below.
{\small
\begin{align*}
(-  k_\mathrm{Deg}   w -2  k_\mathrm{bw}  \tfrac{w^2}{2} + 2  k_\mathrm{dw}  w_\mathrm{2})  (K + \sigma_B) + \lambda_W  v_0  (1 + F  \sigma_B)&=0\\
-  k_\mathrm{Deg}   w_\mathrm{2} + k_\mathrm{bw}  \tfrac{w^2}{2} -  k_\mathrm{dw}  w_\mathrm{2} -  k_\mathrm{B1}   w_\mathrm{2}  v + k_\mathrm{D1}   w_\mathrm{2v} + k_\mathrm{K1}   w_\mathrm{2v} -  k_\mathrm{B3}   w_\mathrm{2}  \sigma_B + k_\mathrm{D3}   w_\mathrm{2\sigma B}&=0\\
-  k_\mathrm{Deg}   w_\mathrm{2v} + k_\mathrm{B1}   w_\mathrm{2}  v -  k_\mathrm{D1}   w_\mathrm{2v} -  k_\mathrm{B2}   w_\mathrm{2v}  v + k_\mathrm{D2}   w_\mathrm{2v2} -  k_\mathrm{K1}   w_\mathrm{2v} + k_\mathrm{K2}   w_\mathrm{2v2} + k_\mathrm{B4}   w_\mathrm{2\sigma B}  v -  k_\mathrm{D4}   w_\mathrm{2v}  \sigma_B&=0\\
(-  k_\mathrm{Deg}   v -  k_\mathrm{B1}   w_\mathrm{2}  v + k_\mathrm{D1}   w_\mathrm{2v} -  k_\mathrm{B2}   w_\mathrm{2v}  v + k_\mathrm{D2}   w_\mathrm{2v2} -  k_\mathrm{B4}   w_\mathrm{2\sigma B}  v + k_\mathrm{D4}   w_\mathrm{2v}  \sigma_B + k_\mathrm{P}   v_\mathrm{Pp})  (K + \sigma_B)& \\+ \lambda_V  v_0  (1 + F  \sigma_B)&=0\\
-  k_\mathrm{Deg}   w_\mathrm{2v2} + k_\mathrm{B2}   w_\mathrm{2v}  v -  k_\mathrm{D2}   w_\mathrm{2v2} -  k_\mathrm{K2}   w_\mathrm{2v2}&=0\\
-  k_\mathrm{Deg}   v_P + k_\mathrm{K1}   w_\mathrm{2v} + k_\mathrm{K2}   w_\mathrm{2v2} -  k_\mathrm{B5}   v_P  \mathrm{phos} + k_\mathrm{D5}   v_\mathrm{Pp}&=0\\
(-  k_\mathrm{Deg}   \sigma_B -  k_\mathrm{B3}   w_\mathrm{2}  \sigma_B + k_\mathrm{D3}   w_\mathrm{2\sigma B} + k_\mathrm{B4}   w_\mathrm{2\sigma B}  v -  k_\mathrm{D4}   w_\mathrm{2v}  \sigma_B)  (K + \sigma_B) + v_0  (1 + F  \sigma_B)&=0\\
-  k_\mathrm{Deg}   w_\mathrm{2\sigma B} + k_\mathrm{B3}   w_\mathrm{2}  \sigma_B -  k_\mathrm{D3}   w_\mathrm{2\sigma B} -  k_\mathrm{B4}   w_\mathrm{2\sigma B}  v + k_\mathrm{D4}   w_\mathrm{2v}  \sigma_B&=0\\
-  k_\mathrm{Deg}   v_\mathrm{Pp} + k_\mathrm{B5}   v_P  \mathrm{phos} -  k_\mathrm{D5}   v_\mathrm{Pp} -  k_\mathrm{P}   v_\mathrm{Pp}&=0\\
\mathrm{phos} + v_\mathrm{Pp} - p_\mathrm{tot}&=0
\end{align*}
}
\vspace{-\baselineskip}

\noindent The 23 parameters $k_\mathrm{bw}$, $k_\mathrm{dw}$, $kD$, $k_\mathrm{B1} $, $k_\mathrm{B2} $, $k_\mathrm{B3} $, $k_\mathrm{B4} $, $k_\mathrm{B5} $, $k_\mathrm{D1} $, $k_\mathrm{D2} $, $k_\mathrm{D3} $, $k_\mathrm{D4} $, $k_\mathrm{D5} $, $k_\mathrm{K1} $, $k_\mathrm{K2} $, $k_\mathrm{P} $, $k_\mathrm{Deg} $, $v_0$, $F$, $K$, $\lambda_W$, $\lambda_V$, $p_\mathrm{tot}$ describe the speed of different reactions.
The following parameter values are derived from experimental data.
{\small
\begin{align*}
&k_\mathrm{Bw} = 3600; k_\mathrm{Dw} = 18; k_D = 18
k_\mathrm{B1}  = 3600; k_\mathrm{B2}  = 3600; k_\mathrm{B3}  = 3600; k_\mathrm{B4}  = 1800; k_\mathrm{B5}  = 3600;\\
&k_\mathrm{D1}  = 18; k_\mathrm{D2}  = 18; k_\mathrm{D3}  = 18; k_\mathrm{D4}  = 1800; k_\mathrm{D5}  = 18;
k_\mathrm{K1}  = 36; k_\mathrm{K2}  = 36;
k_\mathrm{P}  = 180; k_\mathrm{Deg}  = 0.7;\\
&v_0 = 0.4; F = 30; K = 0.2;
\lambda_W = 4; \lambda_V = 4.5;
p_\mathrm{tot} = 2;
\end{align*}}
\vspace{-\baselineskip}

\noindent As discussed above, only real positive zeros are physically meaningful. Using our implementation we can certify that there are 12 real zeros for this system and that among them there is a unique positive one. It has the following values:
{\small
\begin{alignat*}{4}
&\mathrm{phos}  &&= 0.00406661084 \pm 5.25\cdot 10^{-12},\qquad
&&v &&= 0.0557971948 \pm 4.87\cdot 10^{-12}\\
&v_P &&= 27.0899869 \pm 3.85\cdot 10^{-8},\qquad
&&v_\mathrm{Pp} &&= 1.99593338916 \pm 5.20\cdot 10^{-12}\\
&w &&=  0.10633375735 \pm 8.47\cdot 10^{-12},\qquad
&&w_\mathrm{2} &&= 0.303554095 \pm 5.47\cdot 10^{-10}\\
&w_\mathrm{2v} &&= 2.25701026 \pm 2.08\cdot 10^{-9},\qquad
&&w_\mathrm{2v2} &&= 8.288216246 \pm 9.27\cdot 10^{-10}\\
&w_\mathrm{2\sigma B} &&= 10.42034597 \pm 7.94\cdot 10^{-9},\qquad
&&\sigma_B &&= 0.240800757 \pm 5.17\cdot 10^{-10}.
 \end{alignat*}
 }
\vspace{-\baselineskip}

\noindent This is a proof that the dynamical system has a physically meaningful steady state, which is provably contained in the intervals above. The certification took 0.012 seconds. Hence, our implementation makes it possible to certify solutions for large numbers of parameters in short time.

The code for this example is available at \cite{bacillus-subtilis2020}.
We thank Torkel Loman from the Sainsbury Laboratory at the University of Cambridge for pointing out this example to us.

\section{Interval arithmetic}\label{sec:IA}

Since the 1950s researchers \cite{Moore:1966,Sunaga:1958} have worked on interval arithmetic, which allows certified computations while still using floating point arithmetic. We briefly introduce the concepts from interval arithmetic which are relevant for our article.

\subsection{Real interval arithmetic}

Real interval arithmetic concerns computing with compact real intervals. Following \cite{Mayer:2017} we denote the set of all compact real intervals by
$$\mathbb I\mathbb R:= \{[a,b]\mid a,b\in\mathbb R, a\leq b\}.$$
For $X, Y \in \mathbb{IR}$ and the binary operation $\circ \in \{ +,-, \cdot, / \}$ we define
\begin{equation}\label{real_interval_algebra}
X \circ Y= \{ x \circ y \,|\, x\in X,y\in Y\},
\end{equation}
where we assume $0 \notin Y$ in the case of division.
The interval arithmetic version of these binary operations, as well as other standard arithmetic operations, have explicit formulas. See, e.g., \cite[Sec. 2.6]{Mayer:2017} for more details.

\subsection{Complex interval arithmetic}
We define the set of \emph{rectangular complex intervals} as
$$\mathbb{IC} :=\{X + i Y\mid X,Y\in\mathbb{IR}\},$$
where $X+iY = \{x+iy\mid x\in X,y\in Y\}$ and $i = \sqrt{-1}$.
Following \cite[Ch. 9]{Mayer:2017} we define the algebraic operations for $I = X+iY, J=W+iZ\in\mathbb{IC}$ in terms of operations on the real intervals from (\ref{real_interval_algebra}):
\begin{alignat}{2}\label{complex_arithmetic}
I + J &:= (X + W) + i (Y+Z),\qquad I \cdot J &&:= (X\cdot W - Y \cdot Z) + i (X\cdot Z + Y\cdot W)\\\nonumber
I - J &:= (X - W) + i (Y-Z),\qquad \;\, \frac{I}{J} &&:= \frac{X \cdot W + Y\cdot Z}{W\cdot W + Z\cdot Z}  + i \frac{Y \cdot W - X\cdot Z}{W\cdot W + Z\cdot Z}.
\end{alignat}
It is necessary to use (\ref{real_interval_algebra}) instead of complex arithmetic for the definition of algebraic operations in $\mathbb{IC}$.
The following example from \cite{Mayer:2017} demonstrates this. Consider the intervals $I=[1,2] + i[0,0]$ and~$J = [1,1] + i[1,1]$.
Then, $\{x\cdot y| x\in I,y\in J\} = \{t(1+i)\mid 1\leq t\leq 2\}$ is not a rectangular complex interval, while $I\cdot J = [1,2] + i [1,2]$ is.

The algebraic structure of $\mathbb{IC}$ is given by following theorem; see, e.g., \cite[Theorem 9.1.4]{Mayer:2017}.
\begin{theorem}\label{thm_arithmetic}
The following holds.
\begin{enumerate}
\item $(\mathbb{IC}, +)$ is a commutative semigroup with neutral element.
\item $(\mathbb{IC}, +, \cdot)$ has no zero divisors.
\end{enumerate}
Furthermore, if $I,J,K,L\in\mathbb{IC}$, then
\begin{enumerate}[resume]
\item $I\cdot (J + K) \subseteq I\cdot J + I\cdot K$, but equality does not hold in general.
\item $I\subseteq J, K\subseteq L$, then $I \circ K\subseteq J\circ L$ for $\circ\in\{+,-,\cdot, /\}$.
\end{enumerate}
\end{theorem}
Working with interval arithmetic is challenging because of the third item from the previous theorem: distributivity does not hold in $\mathbb{IC}$.
As a consequence, in $\mathbb{IC}$ the evaluation of polynomials depends on the exact order of the evaluation steps.
Therefore, the evaluation of polynomial maps $F: \mathbb{IC}^n \to \mathbb{IC}$ is only well-defined if $F$ is defined by a straight-line program, and not just by a list of coefficients. Figure \ref{fig2} demonstrates this issue in an example. See, e.g., \cite[Sec. 4.1]{burgisser2013algebraic} for an introduction to straight-line programs.

\begin{figure}[hb]
 \begin{minipage}{0.49 \textwidth}
     \centering
     \begin{tikzpicture}
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (x) at (0,3) {x};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (y) at (1,3) {y};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (z) at (2,3) {z};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (+) at (0.5, 2) {+};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (*) at (1.5,1) {$\cdot$};
         \path [->](x) edge (+);
         \path [->](y) edge (+);
         \path [->](+) edge (*);
         \path [->](z) edge (*);
 \end{tikzpicture}
 \end{minipage}
 \begin{minipage}{0.49 \textwidth}
     \centering
     \begin{tikzpicture}
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (x) at (0,3) {x};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (z) at (1,3) {z};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (y) at (2,3) {y};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (*1) at (0.5, 2) {$\cdot$};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (*2) at (1.5, 2) {$\cdot$};
         \node[shape=circle,draw=black, minimum size=6mm,inner sep=0pt] (+) at (1,1) {+};
         \path [->](x) edge (*1);
         \path [->](y) edge (*2);
         \path [->](z) edge (*1);
         \path [->](z) edge (*2);
         \path [->](*1) edge (+);
         \path [->](*2) edge (+);
 \end{tikzpicture}
 \end{minipage}
\caption{\label{fig2}The picture shows two straight-line programs for evaluating the polynomial $f(x,y,z)=(x+y)z$.  Let $I=([-1, 0],\,[1, 1],\, [0, 1])^T$. Then, the program on the left evaluated at $I$ yields $f(I) = ( [-1, 0]  + [1, 1] ) [0, 1] = [0,1]$, while the program on the right yields $f(I) = [-1, 0]  [0, 1]+ [1, 1][0, 1]   = [-1,1]$.}
\end{figure}

Arithmetic in $\mathbb{IC}^n$ is defined in the expected way. If $I=(I_1,\ldots,I_n),J=(J_1,\ldots,J_n) \in\mathbb{IC}^n$,
$$I + J = (I_1+J_1,\ldots,I_n+J_n).$$
Scalar multiplication for $I\in\mathbb{IC}$ and $J\in\mathbb{IC}^n$ is defined as $I \cdot J = (I\cdot J_1,\ldots,I\cdot J_n)$.
The product of an interval matrix $A=(A_{i,j})\in\mathbb{IC}^{n\times n}$ and an interval vector $I\in\mathbb{IC}^n$ is
\begin{equation}\label{def_matrix_mult}A\cdot I := I_1\cdot \begin{bmatrix}
A_{1,1}  \\
\vdots\\
A_{n,1}
 \end{bmatrix} + \cdots +
 I_n\cdot \begin{bmatrix}
A_{1,n}  \\
\vdots\\
A_{n,n}
 \end{bmatrix}.
\end{equation}
Similar to the one-dimensional case $(\mathbb{IC}^n, +)$ is a commutative semigroup with neutral element.

\section{Certifying zeros with interval arithmetic}\label{sec:krawczyk}

In 1969 Krawczyk \cite{Krawczyk:1969} developed an interval arithmetic version of Newton's method. Later in 1977 Moore \cite{Moore:1977} recognized that Krawczyk's method can be used to certify the existence and uniqueness of a solution to a system of nonlinear equations.
Interval arithmetic and interval Newton's method are a prominent tool in many areas of applied mathematics; e.g., in chemical engineering \cite{Gopalan:1995}, thermodynamics \cite{Hatice:2005} and robotics \cite{Kumar:Sen:Shooe:2015}.

The results in this section are stated for square polynomial systems, but they hold equally for square systems of rational functions.
Krawcyzk's method is even valid for general square systems of analytic functions.
Nevertheless, in view of the simplicity of the exposition all statements here are only formulated for polynomial systems.

\subsection{Krawczyk's method}
In this section we recall Krawczyk's method for zeros of polynomial systems. First, we need three definitions.

\begin{definition}[Interval enclosure]\label{enclosure}
Let $F: \mathbb{C}^n \rightarrow \mathbb{C}^{n}$ be a system of polynomials. A map $\square F: \mathbb{IC}^n \rightarrow \mathbb{IC}^{n}$ is an interval enclosure of $F$ if for every $I \in \mathbb{IC}^n$ we have
$
\{F(x) \mid x \in I\} \subseteq \square F(I).$
\end{definition}
In the rest of this article we use the notation $\square F$ to denote the interval enclosure of $F$.
Also, we do not distinguish between a point $x\in\mathbb C^{n}$ and the complex interval $[\mathrm{Re}(x),\mathrm{Re}(x)] + i[\mathrm{Im}(x),\mathrm{Im}(x)]$ defined by $x$. We simply use the symbol ``$x$'' for both terms so that $\square F(x)$ is well-defined.
\begin{definition}[Interval matrix norm]\label{matrix_norm}
Let $A\in\mathbb{IC}^{n\times n}$. We define the operator norm of $A$ as
$\Vert A\Vert_\infty := \max\limits_{B\in A} \max\limits_{v\in\mathbb C^n} \tfrac{\Vert Bv\Vert_\infty}{\Vert v\Vert_\infty},$
where $\Vert (v_1,\ldots,v_n) \Vert_\infty = \max_{1\leq i\leq n}\vert v_i\vert$ is the infinity norm in $\mathbb C^n$.
\end{definition}

Next we introduce an interval version of the Newton operator, the \emph{Krawczyk operator} \cite{Krawczyk:1969}.
\begin{definition}\label{def_krawczyk}
Let $F: \mathbb{C}^n \rightarrow \mathbb{C}^n$ be a system of polynomials, and $\mathrm{J}F$ be its Jacobian matrix seen as a function $\mathbb C^n \to \mathbb C^{n\times n}$.
Let $\square F$ be an interval enclosure of $F$ and $\square \mathrm{J}F$ be an interval enclosure of $\mathrm{J}F$.
Furthermore, let $I \in \mathbb{IC}^n$ and $x \in \mathbb C^n$ and let $Y \in \mathbb C^{n\times n}$ be an invertible matrix.
We define the Krawczyk operator
\begin{align*}
K_{x,Y}(I) := x - Y \cdot \square F(x) + (\mathbf 1_n - Y \cdot \square \mathrm{J}F(I))(I-x).
\end{align*}
Here, $\mathbf 1_n$ is the ${n\times n}$-identity matrix.
\end{definition}

\begin{remark}
In the literature, $K_{x,Y}(I)$ is often defined using $F(x)$ and not $\square F(x)$.
Here, we use this definition, because in practice it is usually not feasible to evaluate $F(x)$ exactly.
Instead, $F(x)$ is replaced by an interval enclosure.
\end{remark}
\begin{remark}
The second part of Theorem \ref{thm:krawczyk} motivates to find a matrix $Y \in \mathbb{C}^{n \times n}$ such that $|| 1_n - Y \cdot \square \mathrm{J}F(I) ||_\infty$ is minimized.
 A good choice is an approximation of the inverse of $\mathrm JF(x)$.
\end{remark}

We are now ready to state the theorem behind Krawczyk's method.
The first proof for real interval arithmetic is due to Moore \cite{Moore:1977}. One of the few sources, which has this theorem in the complex setting, is \cite{BLL:2019}.
For completeness we recall their proof in this section.
Note that all the data in the theorem can be computed using interval arithmetic.

\begin{theorem}\label{thm:krawczyk}
Let $F: \mathbb{C}^n \rightarrow \mathbb{C}^n$ be a system of polynomials and $I\in\mathbb{IC}^n$.
Let $x\in I$ and $Y\in\mathbb{C}^{n\times n}$ be an invertible complex $n\times n$ matrix.
The following holds.
  \begin{enumerate}
    \item If $K_{x,Y}(I) \subset I$, there is a zero of $F$ in $I$.
    \item If additionally
    $\sqrt{2} \, \lVert \mathbf 1_n - Y \square \mathrm{J}F(I) \rVert_\infty < 1$, then
    $F$ has exactly one zero in $I$.
    \end{enumerate}
\end{theorem}

To simplify our language when talking about intervals $I \in \mathbb{IC}^n$ satisfying Theorem \ref{thm:krawczyk} we introduce the following definitions.
\begin{definition}\label{def_interval_zero}
  Let $F: \mathbb{C}^n \rightarrow \mathbb{C}^n$ be a square system of polynomials and $I \in \mathbb{IC}^n$. Let $K_{x,Y}(I)$ be the associated Krawczyk operator (see Definition \ref{def_krawczyk}).
  If there exists an invertible matrix $Y\in\mathbb C^{n\times n}$, such that $K_{x,Y}(I)\subset I$, we say that $I$ is an \emph{interval approximate zero}~$F$.
  We call $I$ a \emph{strong interval approximate zero} of $F$ if in addition $\sqrt{2} \lVert \mathbf 1_n - Y \square \mathrm{J}F(I) \rVert_\infty < 1$ .
\end{definition}
\begin{definition}
  If $I$ is an interval approximate zero, then, by Theorem \ref{thm:krawczyk}, $I$ contains a zero of~$F$. We call such a zero an \emph{associated zero} of $I$.
  If $I$ is a strong interval approximate zero then there is a unique associated zero and we refer to is as \emph{the} associated zero of $I$.
\end{definition}

The notion of strong interval approximate zero is stronger than the definition suggests at first sight. We not only certify that a unique zero of $F$ exists inside $I$, but even that we can approximate this zero with arbitrary precision. This is shown in the next proposition, which we prove at the end of this section.
\begin{proposition}\label{interesting_prop}
Let $I$ be a strong interval approximate zero of $F$ and let $x^*\in I$ be the unique zero of $F$ inside $I$. Let $x\in I$ be any point in $I$. We define $x_0:=x$ and for all $i\geq 1$ we define the iterates $x_i := x_{i-1} - Y\,F(x_{i-1})$, where $Y\in\mathbb C^{n\times n}$ is the matrix from Definition \ref{def_interval_zero}. Then, the sequence $(x_i)_{i\geq 0}$ converges (at least linearly) to $x^*$.
\end{proposition}

The idea for the proof of both Theorem \ref{thm:krawczyk} and Proposition \ref{interesting_prop} is to verify that for strong interval approximate zeros $I$ the map $G_Y(x)= x-Y\cdot F(x)$ defines a contraction on $I$.
If this is true, by Banach's Fixed Point Theorem there is exactly one fixed-point of this map in $I$.
Since $Y$ is invertible, this implies that there is exactly one zero to $F(x)$ in $I$.

Before we give the proof of Theorem \ref{thm:krawczyk}, we need a lemma.
It is a direct sequence of a complex version of the mean-value theorem which is shown implicitly in the proof of \cite[Lemma 2]{BLL:2019}.

\begin{lemma}
  \label{lemma: oracle_ext}
Fix a matrix $Y\in\mathbb C^{n\times n}$ and define $G_Y(x)= x-YF(x)$. Let $I\in\mathbb{IC}^n$ be an interval vector and $x,z\in I$. Then, we have
\begin{enumerate}
\item $G_Y(z)-G_Y(x) \in (\mathbf 1_n - Y \cdot\square \mathrm{J}F(I)) \,  \mathrm{Re}(z-x) +  (\mathbf 1_n - Y \cdot\square \mathrm{J}F(I))\, i\mathrm{Im}(z-x).$
\item $G_Y(I)\subset K_{x,Y}(I)$.
\end{enumerate}
 \end{lemma}
The following proof is adapted from \cite[Lemma 2]{BLL:2019}.
\begin{proof}[Proof of Lemma \ref{lemma: oracle_ext}]
In the proof we abbreviate $G:=G_Y$. We first show the second part assuming the first part of the lemma.
Then, we prove the first part. We fix an interval $I\in\mathbb{IC}^n$ and $x,z\in\mathbb C^n$.

For the second part, we have to show that for all $I\in\mathbb{IC}^n$ we have $G(I)\subset K_{x,Y}(I)$.
To show this we define the interval matrix $M:=(\mathbf 1_n - Y \square \mathrm{J}F(I))\in\mathbb{IC}^{n\times n}$.
By definition of $K_{x,Y}$ we have $G(x) + M(I-x) \subset K_{x,Y}(I)$.
Thus, we have to show that $G(z) - G(x) \in M(I-x)$, since $z\in I$ is arbitrary.
The first part of the lemma implies that we can find matrices $M_1,M_2\in M$ such that~$G(z) - G(x) = M_1 \mathrm{Re}(z-x) +i M_2 \mathrm{Im}(z-x)$.
Decomposing the matrices into real and imaginary part we find
$$G(z)-G(x) = \mathrm{Re}(M_1) \mathrm{Re}(z-x) - \mathrm{Im}(M_2) \mathrm{Im}(z-x) + i (\mathrm{Im}(M_1) \mathrm{Re}(z-x) + \mathrm{Re}(M_2) \mathrm{Im}(z-x)).$$
Since $z-x\in I$ and by definition of the complex interval multiplication from (\ref{complex_arithmetic}) and the interval matrix-vector-multiplication (\ref{def_matrix_mult}) we see that $G(z)-G(x)\in M (I-x)$.
This finishes the proof for the second part.

The first part of the lemma may be shown entry-wise. We will show this by combining a complex version of the mean value theorem with the following observation:  $\mathrm{J}G(x) = \mathbf 1_n - Y \cdot \mathrm{J}F(x)$, so we have the inclusion
  \begin{equation}\label{eq11}
    \mathrm{J}G(I) = \mathbf 1_n - Y \cdot\mathrm{J}F(I)   \subseteq\mathbf 1_n - Y \cdot\square \mathrm{J}F(I).
  \end{equation}
We relate $G(z)- G(x)$ to (\ref{eq11}) using the mean value theorem. First, we define $w:=\mathrm{Re}(z) + i\mathrm{Im}(x)$.
Let $1 \leq j \leq n$ and let $G_j$ denote the $j$-th entry of $G$. We define the function $h(t):=G_j(tz + (1-t)w)$. The real and imaginary part of $h(t)$ are real differentiable functions of the real variable $t$.
The mean value theorem can be applied, and we find $0<t_1,t_2<1$ such that $\mathrm{Re}(h(1)) - \mathrm{Re}(h(0)) = \tfrac{\mathrm{d}}{\mathrm d t} \mathrm{Re}(h(t_1))$ and $\mathrm{Im}(h(1)) - \mathrm{Im}(h(0))  = \frac{\mathrm{d}}{\mathrm d t} \mathrm{Im}(h(t_2))$.
Setting $c_1 = t_1z+(1-t_1)w$ and $c_2 = t_2z+(1-t_2)w$ this implies
$$G_j(w) - G_j(z) = (\nabla_{\mathrm{Re}}\mathrm{Re}(G_j(c_1)))^T (z-w) +   i\nabla_{\mathrm{Re}}\mathrm{Im}(G(c_2)))^T (z-w),$$
where  $\nabla_{\mathrm{Re}} G$ denotes the vector of partial derivatives with respect to the real variable.
Let us denote by $G_j'$ the complex derivative of $G_j$; that is, $G_j': \mathbb C^{n} \to \mathbb C^{n}$ as a function.
From the Cauchy Riemann equations it follows that $\nabla_{\mathrm{Re}}\mathrm{Re}(G_j(c_1)) = \mathrm{Re}(G_j'(c_1))$ and likewise $\nabla_{\mathrm{Re}}\mathrm{Im}(G(c_2)) = \mathrm{Im}(G'(c_2))$.
This yields  $G_j(z) - G_j(w) = (\mathrm{Re}(G_j'(c_1) +   i\mathrm{Im}(G_j'(c_2)))^T (z-w)$. Putting these equations ranging over $j$ together we find
$G(z) - G(w) = (\mathrm{Re}(\mathrm{J}G(c_1)) +   i\mathrm{Im}(\mathrm{J}G(c_2))) (z-w)$.
By construction, $c_1$ and~$c_2$ are contained in $I$, because $w$ and $z$ are contained in $I$, and $I$ is a product of rectangles and thus convex.
Combined with (\ref{eq11}) this yields
$$G(z) - G(w) \in (\mathbf 1_n - Y \cdot\square \mathrm{J}F(I))(z-w).$$
Using essentially the same arguments for the path from $x$ to $w$ we also find
$$G(w) - G(x) \in (\mathbf 1_n - Y \cdot\square \mathrm{J}F(I))(w-x).$$
By construction, $z-w = i\mathrm{Im}(z-x)$ and $w-x = \mathrm{Re}(z-x)$, which implies
$$G(z)-G(x) \in (\mathbf 1_n - Y \cdot\square \mathrm{J}F(I)) \,  \mathrm{Re}(z-x) +  (\mathbf 1_n - Y \cdot\square \mathrm{J}F(I))\, i\mathrm{Im}(z-x).$$
This finishes the proof.   \end{proof}

\begin{proof}[Proof of Theorem \ref{thm:krawczyk} and Proposition \ref{interesting_prop}]
We fix $Y\in\mathbb C^{n\times n}$. The second part of Lemma \ref{lemma: oracle_ext} implies that, if we have $K_{x,Y} (I) \subseteq I$, then $G_Y(I) \subseteq I$. Brouwer's fixed point Theorem shows that $G_Y$ has a fixed point in~$I$. Since $Y$ is assumed to be invertible, the fixed point is a zero of $F$.
This finishes the proof for the first part of Theorem \ref{thm:krawczyk}.
For the second part let $z_0,z_1\in I$.
The first part of Lemma \ref{lemma: oracle_ext} implies
\[
G_Y(z_0) - G_Y(z_1) \in (\mathbf 1_n- Y \cdot \square \mathrm{J}F(I)) \text{Re}(z_0 - z_1) +
(\mathbf 1_n- Y \cdot \square \mathrm{J}F(I)) \ i \text{Im}(z_0 - z_1)  .
\]
(Note that we can't apply the distributivity law because of Theorem \ref{thm_arithmetic} 3.).
Applying norms and using submultiplicativity yields
\[
  \lVert G_Y(z_0) - G_Y(z_1)  \rVert_\infty \leq
  \lVert (\mathbf 1_n- Y \cdot \square \mathrm{J}F(I)) \rVert_\infty\, (\lVert \text{Re}(z_0 - z_1) \rVert_\infty + \lVert\text{Im}(z_0 - z_1) \rVert_\infty).
\]
Since $ \lVert \text{Re}(z_0 - z_1) \rVert_\infty  +\lVert \text{Im}(z_0 - z_1) \rVert_\infty \leq  \sqrt{2} \lVert z_0 -z_1 \rVert_\infty$ it holds
\begin{equation}\label{speed_of_convergence}
\lVert G_Y(z_0) - G_Y(z_1)  \rVert_\infty    \leq \sqrt{2} \lVert \mathbf 1_n- Y \cdot \square \mathrm{J}F(I) \rVert_\infty
 \lVert z_0 -z_1 \rVert_\infty.
\end{equation}
By assumption $ \sqrt{2} \lVert \mathbf 1_n- Y \cdot \square \mathrm{J}F(I) \rVert_\infty$ is smaller than 1 so $G_Y$ is a contraction.
Banach's Fixed Point Theorem implies that $G_Y$ has a unique fixed point in $I$. Since $Y$ is invertible, this fixed point is the unique zero of $F$ in $I$. This shows the second part of Theorem~\ref{thm:krawczyk}.

Finally, let $(x_i)_{i\geq 0}$ be the sequence defined in Proposition \ref{interesting_prop}. By assumption, $x_0\in I$ and $x_{i+1} = G_Y(x_{i})$, $i\geq 0$. The second part of Lemma \ref{lemma: oracle_ext} together with an induction argument imply that $x_i\in I$ for all $i$. Let $x^*$ be the unique zero of $F$ in $I$.
Then, \cref{speed_of_convergence} holds if we replace $z_0$ by $x^*$ and $z_1$ by $x_i$. We get $\lVert x^* - x_{i+1}  \rVert_\infty  < c\cdot \lVert x^* -x_i \rVert_\infty$ with a constant $c:=\sqrt{2} \lVert \mathbf 1_n- Y \cdot \square \mathrm{J}F(I) \rVert_\infty<1$. This shows linear convergence and finishes the proof of Proposition \ref{interesting_prop}.
\end{proof}

\subsection{Certifying reality}\label{sec:reality}
 For many applications only the real zeros of a polynomial system are of interest.
 Since numerical homotopy continuation computes in $\mathbb C^n$, it is important to have a rigorous method to determine whether a zero is real.

 Recall from Definition \ref{def_interval_zero} the notion of \emph{strong interval approximate zero}.

\begin{lemma}\label{lemma:real_solution}
  Let $F: \mathbb{C}^n \rightarrow \mathbb{C}^n$ be a real square system of polynomials and $I \in \mathbb{IC}^n$ a strong interval approximate zero of $F$.Then there exists $x\in I$ and $Y\in \mathbb{C}^{n\times n}$ satisfying $K_{x,Y}(I) \subset I$ and $\sqrt{2} \, \lVert \mathbf 1_n - Y \square \mathrm{J}F(I) \rVert_\infty < 1$.
  If additionally $\{ \bar{z} \, | \, z \in K_{x,Y}(I) \} \subset I$, the associated zero of $I$ is real.
\end{lemma}
\begin{proof}
  Theorem \ref{thm:krawczyk} implies that $F$ has a unique zero $s \in K_{x,Y}(I) \subset I$.
  Since $F$ is a real polynomial system it follows that also the element wise complex conjugate $\bar{s}$ is a zero of $F$.
  If we have that $\bar{s} \in \{ \bar{z} \, | \, z \in K_{x,Y}(I) \} \subset I$, then $\bar{s} = s$, since otherwise
  $\bar{s}$ and $s$ would be two distinct zeros of $F$ in~$I$, contradicting the uniqueness result from Theorem \ref{thm:krawczyk}.
\end{proof}

For a wide range of applications positive real zeros are of particular interest.
\begin{corollary}
  Let $F: \mathbb{C}^n \rightarrow \mathbb{C}^n$ be a real square system of polynomials and $I \in \mathbb{IC}^n$ a strong interval approximate zero of $F$ satisfying the conditions of Lemma \ref{lemma:real_solution}.
  If $\textnormal{Re}(I) > 0$ then the associated zero of $I$ is real and positive.
\end{corollary}

If the reality test in Lemma \ref{lemma:real_solution} fails for a strong interval approximate zero $I \in \mathbb{C}^n$ then this does not necessarily mean that the associated zero of $I$ is not real.
A sufficient condition that $I$ is not real is that there is a coordinate such that the imaginary part of it does not contain zero.
\begin{lemma}\label{lemma:complex_solution}
  Let $F(x)$ be a square system of polynomials or rational functions and let $I \in \mathbb{IC}^n$ be a strong interval approximate zero of $F$.
  If there exists $k \in \{1,\ldots,n\}$ such that $0 \notin \mathrm{Im}(I_k)$ then the associated zero of $I$ is not real.
\end{lemma}
\begin{proof}
  The associated zero $x$ of $I$ is contained in $I$. Since $0 \notin \mathrm{Im}(I_k)$ follows $x_k \notin \mathbb{R}$ and $x \notin \mathbb{R}^n$.
\end{proof}

Now assume that the certification routine produced a list $\mathcal{I}$ of $m$ distinct strong interval approximate zeros for a given system $F$, and that $m$ also agrees with the theoretical upper bound on the number of isolated nonsingular zeros of $F$.
If we apply Lemma \ref{lemma:real_solution} to $I_k \in \mathcal{I}$, then we obtain only a
\emph{lower bound}, say $r$, on the number of real zeros of $F$.
However, combined with Lemma~\ref{lemma:complex_solution} we can also obtain an \emph{upper bound} of the number of real zeros.
If these two bounds agree we obtain a certificate that, among the associated zeros of the intervals in $\mathcal{I}$, there are \emph{exactly} $r$ real zeros.
An application of this is, e.g., the study of the distribution of the number of real solutions of the power flow equations \cite{Lindberg:Zachariah:Boston:Lesieutre:2020}.

\section{Implementation details}\label{sec:implementation_details}
In this section we describe the necessary considerations to implement Krawcyzk's method described in Section \ref{sec:krawczyk}, as well as the technical realization in \texttt{HomotopyContinuation.jl}.
The certification routine takes as input a square polynomial system $F: \mathbb{C}^n \rightarrow \mathbb{C}^n$ and a finite list $X \subset \mathbb{C}^n$ of  (suspected) approximations of isolated nonsingular zeros of $F$.
It is also possible to provide a square system of rational functions as input.
Similar to Section \ref{sec:krawczyk}, we restrict to the polynomial case to simplify our exposition .
It returns a list of strong interval approximate zeros $\mathcal{I}=\{I_1,\ldots,I_m\}$ in $\mathbb{IC}^n$, such that no two intervals $I_k$ and $I_\ell$, $k \ne \ell$, overlap.
If two strong interval approximate zeros don't overlap then this implies that their associated zeros are distinct.
Additionally, if $F$ is a real polynomial system then for each $I_k \in \mathcal{I}$ it is determined whether its associated zero is real.
The prototypical application of the certification routine is to take as input approximations of all isolated nonsingular solutions $X \subset \mathbb{C}^n$ of $F$ as computed by numerical homotopy continuation methods.

\subsection{Interval enclosures for polynomial systems}
As already discussed in Section \ref{sec:IA}, the fact that distributivity doesn't hold in $\mathbb{IC}$ require us to define the polynomial system $F: \mathbb{C}^n \to \mathbb{C}^n$, and its interval enclosure $\square F$, by a straight-line program, and not just by a list of coefficients.
The overestimation of the interval enclosure $\square F$ increases with the size of the straight line program.
Therefore, it is good to express $F$ and its enclosure $\square F$ by the smallest straight line program possible.
To achieve this, \texttt{HomotopyContinuation.jl} automatically applies
a multivariate version of Horner's rule to reduce the number of operations necessary to evaluate $F$ and $\square F$.

\begin{remark}
Our implementation of interval enclosures can also be used to prove that a polynomial map $F:\mathbb C^n\to\mathbb C^m$ with real coefficients, evaluated at a real point $p\in\mathbb R^n$, is positive. To verify this, one takes an interval $I\in \mathbb{IC}^n$ of the form $I=J + i[0,0]^{\times n}$ such that $p\in J$. If $\square F$ is an interval enclosure of $F$, and if $\square F(I)\subset \mathbb R_{>0}^m + i[0,0]^{\times m}$, then this is a proof that $F(p) \in \mathbb R_{>0}^m$.
\end{remark}

\subsection{Machine interval arithmetic}
In the next subsection we give a method to construct an candidate $I \in \mathbb{IC}^n$ for a strong interval approximate zero.
Before, we need to study \emph{machine interval arithmetic}; the realization of interval arithmetic with finite precision floating point arithmetic.
We assume the standard model of floating point arithmetic \cite[Section 2.3]{Higham:2002}, where the result of a floating point operation is accurate up to relative unit roundoff $u$: $\mathrm{fl}(x \circ y) = (x  \circ  y)(1 + \delta),$ where $|\delta| \le u$ and $\circ\in\{+,-,*,/\}$.
For instance, following the IEE-754 standard, the unit roundoff in double precision arithmetic is $u = 2^{-53} \approx 2.2 \cdot 10^{-16}$.
The key property in the context of interval arithmetic is that each result of a floating point operation can be rounded outwards, such that the resulting \emph{interval} contains the true (exact) result; see, e.g., \cite[Section 3.2]{Mayer:2017}.
Therefore, given $X, Y \in \mathbb{IC}$ the result of $X \circ Y$, $\circ\in\{+,-,*,/\}$, is $\mathrm{fl}(X \circ Y) := \{ (x  \circ  y)(1+\delta) \,|\, |\delta| \le u, x\in X,y\in Y \}$ in machine arithmetic. This interval contains $X\circ Y$. It is \emph{larger}.
Additionally, for a given $x \in \mathbb{IC}$, all intervals of the form
$\{x + (|\mathrm{Re}(x_j)| + i |\mathrm{Im}(x_j)|)\delta \, | \, |\delta| \le \mu \}$ with $0 < \mu \le u$ are indistinguishable when working with precision $u$.

Consequently it is possible that the Krawczyk operator $K_{\tilde{x},Y}$, see Definition \ref{def_krawczyk},
is a contraction for the interval $I$, but that machine arithmetic can't verify this, because $\mathrm{fl}(X  \circ  Y)$ is larger than $X\circ Y$.
In such a case, the unit roundoff $u$ needs to be sufficiently decreased.
For this reason our implementation uses machine interval arithmetic based on double precision arithmetic as well as, if necessary, the arbitrary precision interval arithmetic implemented in \texttt{Arb} \cite{Johansson2017arb}.

\subsection{Determining strong interval approximate zeros}\label{sec:inflation}
In a first step, the certification routine attempts to produce for a given $x \in X$ a strong interval approximate zero $I \in \mathbb{IC}^n$.
Recall that for $I \in \mathbb{IC}^n$ to be a strong interval approximate zero we need by Theorem \ref{thm:krawczyk} to have a point $\tilde{x} \in I$, and a matrix $Y \in \mathbb{C}^{n \times n}$ such that $K_{\tilde{x},Y}(I) \subset I$, and~$\sqrt{2} \, \lVert \mathbf 1_n - Y \square \mathrm{J}F(I) \rVert_\infty < 1$.

Given a point $x \in X$ and a unit roundoff $u$, the point $x$ is refined  using Newton's method to maximal accuracy. We denote this refined point $\tilde{x}$.
Here, we assume that $x$ is already in the region of quadratic convergence of Newton's method.
Next, the point $\tilde{x}$ needs to be inflated to an interval~$I$ with $\tilde{x} \in I$.
This process is called $\varepsilon$-inflation in the literature \cite[Sec. 4.3]{Mayer:2017}.
However, choosing the correct $I$ is a hard problem:
if $I$ is too small or too large, then the Krawcyzk operator is not a contraction.

In spite of these difficulties, we found the following heuristic to determine $I$ work very well.
If we assume $\tilde{x}$ to be in the region of quadratic convergence of Newton's method, it follows from the Newton-Kantorovich theorem that $||\mathrm JF(\tilde{x})^{-1}F(\tilde{x})||_\infty$ is a good estimate of the distance between~$\tilde{x}$ and the convergence limit $x^*$.
Therefore we set $Y \approx \mathrm JF(\tilde{x})^{-1}$ (computed in floating point arithmetic) and use $I = (\tilde{x}_j \pm |(Y \cdot \square F(x))_j|u^{- \frac14})_{j=1,\ldots,n}$ where the factor $u^{- \frac14}$ accounts for the overestimation by machine interval arithmetic.
If $I$ doesn't satisfy the conditions in Theorem~\ref{thm:krawczyk} the procedure is repeated with a smaller unit roundoff $u$.
This repeats until either a minimal unit roundoff is reached or the certification is successful.

\subsection{Producing distinct intervals}\label{sec_disjoint}
Assume now that the steps in Section \ref{sec:inflation} have been performed for all $x \in X$. We obtain a list of strong interval approximate zeros $I_1,\ldots,I_r \in \mathbb{IC}^n$.
In a final step we want to select a subset $M \subset \{1,\ldots,r\}$ such that for all $k, j \in M$, $k \ne j$, the intervals $I_k$ and $I_j$ do not overlap.
If two strong interval approximate zeros do not overlap then it is guaranteed that they have distinct associated zeros.
A simple approach to determine $M$ is to compare all intervals pairwise. However, this approach requires us to perform ${r \choose 2}$ interval vector comparisons.
For larger problems this becomes prohibitively expensive.

Instead, we employ the following improved scheme to determine all non-overlapping intervals.
First, we pick a random point $q \in \mathbb{C}^n$ and compute in interval arithmetic for each $I_k$, $k \in M$, the squared Euclidean distance $d_k \in \mathbb{IR}$ between $I_k$ and $q$.
Due to the guarantees of interval arithmetic we have that $d_k$ and $d_\ell$ overlap if $I_k$ and $I_\ell$ overlap (but the converse it not necessarily true).
Next, we check for all overlapping intervals $d_k, d_\ell\in \{ d_k \in \mathbb{IR} \, | \, k = 1,\ldots,r\}$, whether $I_k$ and $I_\ell$ overlap, and if so, we group them accordingly.
This allows us to construct the set $M$ by selecting those intervals which don't overlap with any other and by picking one representative of each cluster of overlapping intervals.
The worst case complexity of this procedure still requires $O(r^2)$ operations, but in the common case where no or only a small number of intervals overlap $O(r \log r)$ operations are sufficient.

 {\small
\bibliographystyle{alpha}


\section{Path integral formulation for continuous reservoirs} \label{sec:continuous}

We start with the implementation and evolution of a system subjected to a continuous harmonic reservoir. For a numerically exact description, our theoretical approach is based on a real-time path integral formulation~\cite{Feynman1965, Weiss2011, Caldeira1983, Vagov2011, Glaessl2013, Barth2016, Cosacchi2018, Makri1995, Makri1995a}. In recent breakthroughs, path integrals have been implemented in a tensor network approach based on MPS~\cite{Strathearn2018, Jorgensen2019}, allowing to solve non-Markovian dynamics and providing an efficient representation using high-dimensional tensors with restricted correlations~\cite{Strathearn2018}. In the following, we briefly recapitulate the algorithm introduced by \textit{Strathearn et al.}~\cite{Strathearn2017, Strathearn2018, Strathearn2020thesis}, which is employed as a part of our solution to multiple non-Markovian system-reservoir interactions.

Our goal is to employ path integrals for a numerically exact solution of the von-Neumann equation given a Hamiltonian $H(t)$ describing a time-dependent system-reservoir interaction~\cite{Breuer2002, Mukamel1999},
\begin{equation}
\dot{\rho} (t) = \mathcal{L}(t) \rho(t)= -i/\hbar[H(t),\rho(t)],
\label{eq:von_Neumann}
\end{equation}
with $\rho(t)$ denoting the density matrix and $\mathcal{L}(t)$ the Liouvillian superoperator.
As a main challenge, the evaluation of path integrals becomes increasingly expensive over time, since the history of all preceding paths at times $0, \ldots, t_{n-1}$ must be taken into account for the calculation of the current time step $t_n$. However, in case the system-reservoir correlations are finite in time, the \textit{augmented density tensor} scheme can be introduced for improved numerical accessibility~\cite{Makri1995, Makri1995a}.
Exploiting the finite reservoir memory length, only the last $n_c$ time steps are taken into account for the calculation of the current path. This treatment is known as the \textit{finite memory approximation} and results in the augmented density tensor representation as a solution to the system part of Eq.~\eqref{eq:von_Neumann} with traced out reservoir contributions, which reads at time $t_N=N\Delta t$
\begin{align}
\rho_{i_N i^\prime_N} (t_N) &= \prod_{n=1}^N \sum_{i_{n-1}} \sum_{i^\prime_{n-1}} M_{i_n i_{n-1}} M^*_{i^\prime_{n-1} i^\prime_n} \nonumber \\
& \times \prod_{m=n-n_c}^n \exp{ \left( S^{i_n i_m}_{i^\prime_n i^\prime_m} \right) } \rho_{i_0 i^\prime_0} (0),
\label{eq:adt}
\end{align}
and constitutes a discrete path integral formulation where indices $i_n^{(\prime)}$ contain the left (right) configuration of the system at time $t_n=n \Delta t$.
The field transformation matrix $M_{i_n i_{n-1}}$ e.g. accounts for the action of an external driving field $\Omega_0$. For the case $\Omega_0=0$ considered below, it takes the simple form $M_{i_n i_{n-1}}= \mathbb{1}\delta_{i_n,i_{n-1}}$. The influence functional is given by
\begin{equation}
S^{i_n i_m}_{i^\prime_n i^\prime_m} =  -( i_n - i^\prime_n ) [ \eta_{n-m} i_m - \eta^*_{n-m} i_m^\prime ],
\end{equation}
with
\begin{equation}
\eta_{n-m} = \int_{(n-1)\Delta t}^{n \Delta t} \mathrm{d} \tau \int_{(m-1)\Delta t}^{m \Delta t} \mathrm{d} \tau^\prime \, \phi(\tau-\tau^\prime),
\end{equation}
and $\phi(\tau-\tau^\prime)$ the reservoir autocorrelation function~\cite{Weiss2011, Caldeira1983, Leggett1987}. Using an \textit{improved} finite memory approximation, for $n-m \equiv n_c$ all former paths up to $t_{n_c}=n_c \Delta t$ are additionally incorporated in the integration, i.e. $\eta_{n_c} := \eta_{n-m} + \sum_{k=1}^{n-n_c-1} \eta_{n-k}$~\cite{Strathearn2017}.
Under this approximation, it is possible to restate the augmented density tensor and its time evolution efficiently as a tensor network~\cite{Strathearn2018}: First, Eq.~\eqref{eq:adt} is mapped to a vector $\rho_{j_n}$ in Liouville space,
\begin{align}
\rho_{j_N} (t_N) &= \prod_{n=1}^N \prod_{m=n-n_c}^n I(j_n,j_m) \rho_{j_0} (0),
\end{align}
with $I(j_n,j_m) := \sum_{j_{n-1}} \tilde{M}_{j_n j_{n-1}} \exp{ ( \tilde{S}^{j_n j_m} ) }$. Here, left and right system indices $i_k$, $i_k^\prime$ have been combined to a single index $j_k$ for each time step, resulting in Liouville space representations $\tilde{M}_{j_n j_{n-1}}$ and $\tilde{S}^{j_n j_m}$ of the field transformation matrix and the influence functional, respectively. Afterwards, the augmented density tensor is rewritten as an MPS, storing the present and up to $n_c-1$ past states in individual tensors with the oldest state located at the left end of the MPS. In this representation, tensor compression by consecutive applications of the singular value decomposition~\cite{Schollwoeck2011} reduces the memory requirements to polynomial rather than exponential scaling with respect to $n_c$~\cite{Strathearn2018}.

\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{tensornetwork_feedback_comb2}
\caption{MPS implementation of a time-discrete memory in Liouville space. (a)-(d) show the operations performed during the first time step, updating the system state (red square) interacting with present and past reservoir bins (orange and blue circles) by application of the Liouvillian $\mathcal{L}(t)$. (e) shows the second time step with analogous operations. (f) Occupation dynamics of a two-level emitter subjected to feedback [Eq.~\eqref{eq:H_D}], calculated at $\tau=3\,$ps and $\Gamma=31.6\,\mathrm{ps}^{-1}$. Numerical results obtained from the MPS algorithm (solid blue line) are in excellent agreement with the analytical solution for the single excitation case (dashed grey line). (g) MPS feedback dynamics for an additional, analytically not accessible phenomenological dephasing rate $\gamma=0.5\,\mathrm{ps}^{-1}$ (orange line) and without dephasing (blue line).}
\label{fig:mps_feedback}
\end{figure*}

The time evolution is carried out by a network of matrix product operators (MPOs), shown schematically in Fig.~\ref{fig:mps_harmonic}(a) (dark grey shapes). During the first time step, the initial system state $\rho_{j_0} (0)$ (red shape) is contracted with the first MPO in the network [dashed frame in Fig.~\ref{fig:mps_harmonic}(a)]. As a result, the system state is updated and the preceding path is stored to its left, increasing the length of the MPS by one. Fig.~\ref{fig:mps_harmonic}(b) shows the MPS after completion of the first time step. Once step $n=n_c$ is reached, the oldest path in the MPS is summed over by application of a delta tensor [semicircular shape in Fig.~\ref{fig:mps_harmonic}(a)], corresponding to the improved finite memory approximation~\cite{Strathearn2017, Strathearn2018}. At this stage, the MPS length is fixed for the rest of the time evolution. Moreover, for time-independent problems, apart from the index nomenclature the structure of the MPO remains unchanged for all time steps $n \geq n_c$, resulting in an additional performance gain.

To provide an example of a continuous reservoir of non-interacting harmonic oscillators, we consider the Hamiltonian~\cite{Mahan2000}
\begin{align}
H_C/\hbar&= \int \! \mathrm{d}^3 q \, \left[ \omega_q b^\dagger_{\bm{q}} b_{\bm{q}} + g_{\bm{q}} \sigma_{11} \left( b^\dagger_{\bm{q}} e^{i\omega_q t} + \mathrm{H.c.} \right) \right],
\label{eq:H_C}
\end{align}
corresponding to diagonal system coupling without inducing system transitions, with system operators $\sigma_{ij}=\ket{i}\bra{j}$, bosonic annihilation (creation) operators $b_{\bm{q}}^{(\dagger)}$ of reservoir modes with frequency $\omega_q=c_s |\bm{q}|$, $c_s$ the sound velocity, and a mode $\bm{q}$-dependent system-reservoir coupling amplitude $g_{\bm{q}}$. The corresponding correlation function reads
\begin{align}
\phi({\tau-\tau^\prime}) = \int \mathrm{d}^3 q \ g_{\bm{q}}^2 \bigg\{ &\coth \left( \dfrac{\hbar \omega_{\bm{q}}}{2 k_B T} \right) \cos [\omega_q (\tau-\tau^\prime)] \nonumber \\
&- i \sin [\omega_q (\tau-\tau^\prime) ] \bigg\},
\end{align}
with temperature $T$ and $k_B$ the Boltzmann constant.
In the following, we choose a generic coupling element e.g. describing acoustic bulk phonons interacting with a quantum emitter~\footnote{As generic coupling, we choose the acoustic bulk phonon coupling element of GaAs, given by $g_{\bm{q}}^{ii} = \sqrt{ \hbar q/(2 \rho c_s )} D_i \exp [ - \hbar q^2/(4 m_i \omega_i) ]$, resulting in a transition coupling element $g_{\bm{q}} := g_{\bm{q}}^{22} - g_{\bm{q}}^{11}$~\cite{Carmele2019, Foerstner2003, *Foerstner2003pssb, CarmeleMilde2013}. Here, $D_i$ are deformation potentials, $m_i$ denote effective masses, $\hbar \omega_i$ refer to the confinement energies and $\rho$ is the mass density of GaAs, respectively.}.
As a benchmark of the tensor network-based path integral implementation, we first calculate the analytically solvable \textit{Independent Boson model}, consisting of a single two-level emitter subjected to pure dephasing by a structured harmonic reservoir, as described by Eq.~\eqref{eq:H_C}. In Fig.~\ref{fig:mps_harmonic}(c), we prepare the initial polarization at $\rho_{01}(0)=0.5i$ and calculate the resulting dynamics at varying temperatures $T$ (solid lines). The corresponding analytical solution (dashed grey lines) is given by~\cite{Carmele2019}
\begin{align}
\rho_{01}(t) &= \exp \Bigg\{ \int \mathrm{d}^3 q \ \Bigg[ -\dfrac{ig_{\bm{q}}^2}{\omega_q} t + \dfrac{ig_{\bm{q}}^2}{\omega_q^2} \sin (\omega_q t) \nonumber \\
&- \dfrac{g_{\bm{q}}^2}{\omega_q^2} \coth \left( \dfrac{\hbar \omega_q}{2 k_B T} \right) \left[ 1- \cos (\omega_q t) \right] \Bigg] \Bigg\}  \rho_{01} (0),
\end{align}
exhibiting excellent agreement with the numerical results at all considered temperatures. In addition, the employed method features very high performance, enabling reservoir memory depths of $n_c=100$ and beyond. To exemplify the capabilities of the tensor network implementation, we calculate the time evolution dynamics of the \textit{Spin-Boson model}, corresponding to Eq.~\eqref{eq:H_C} with an additional continuous driving field term at amplitude $\Omega_0$, $H_{SBM}= H_C + \Omega_0 (\sigma_{01}+\sigma_{10})$.  Fig.~\ref{fig:mps_harmonic}(d) shows the resulting dynamics at parameters $n_c=100$, $\Omega_0= 0.5\,\mathrm{ps}^{-1}$ and $T=77\,$K for the memory, driving field and temperature, respectively.

\section{Time-discrete memory in Liouville space} \label{sec:discrete}

As a second non-Markovian reservoir, we consider a discrete time-bin based quantum memory. Recently established implementations rely on an MPS-based wave function ansatz to compute the quantum stochastic Schrödinger equation~\cite{Pichler2016, Guimond2016, Guimond2017}. However, this formulation is inherently incompatible with the previously introduced path integral formulation. As a solution to this problem, we present an MPS implementation of a time-discrete quantum memory \textit{in Liouville space}.
Here, the dynamics of the system density matrix is prescribed by a Liouvillian superoperator [see Eq.~\eqref{eq:von_Neumann}], with a Hamiltonian $H_D$ containing the time-delayed system-reservoir coupling, such that interactions occurring at time $t$ couple back into the system and affect its state at a subsequent time $t+\tau$, with $\tau$ the retardation time.
Such a time-discrete coupling e.g. arises in a two-level emitter with states $\ket{0}$, $\ket{1}$ at an energy difference $\hbar \omega_0$, placed in front of a mirror with a round trip time $\tau$. The corresponding Hamiltonian reads
\begin{align}
&H_D / \hbar = \omega_0 \sigma_{11} \nonumber \\
&+ \sqrt{\dfrac{2\Gamma}{\pi}} \int \! \mathrm{d} k \, \sin \left( \dfrac{\omega_k \tau}{2} \right) \left[ \sigma_{10} r_{k} e^{i(\omega_0-\omega_k)t} + \mathrm{H.c.} \right],
\label{eq:H_D}
\end{align}
describing off-diagonal system coupling leading to energy exchange between system and reservoir and system phase relaxation, with system operators $\sigma_{ij}=\ket{i}\bra{j}$, bosonic annihilation (creation) operators $r_k^{(\dagger)}$ of photon modes with frequency $\omega_k=ck$, $c$ the speed of light, and a constant electron-photon coupling amplitude $\Gamma$.
For an efficient evaluation, the dynamics imposed by $\mathcal{L}(t)$ is translated in a time bin-based MPS formalism~\cite{Pichler2016, Guimond2016, Guimond2017, Finsterhoelzl2020} which maintains the relevant system-environment correlations, scaling with $\tau$. In case of additional phenomenological dissipative channels, the Liouvillian can be extended by the standard Lindblad operator~\cite{Breuer2002, Mukamel1999}.

For the MPS implementation of the time-ordered Liouvillian, we start from the formal solution of the system part of Eq.~\eqref{eq:von_Neumann} for the density matrix,
\begin{equation}
\rho(t) = T \exp \left[ \int_0^t \mathrm{d} t^\prime \ \mathcal{L}(t^\prime)  \right] \rho(0),
\label{eq:rho_evo1}
\end{equation}
with $T$ the time-ordering operator. For an MPS-based approach and in analogy to the time-discrete path integral formulation, we restate Eq.~\eqref{eq:rho_evo1} in a time-discrete basis, which reads at time $t_N=N \Delta t$
\begin{equation}
\rho(t_N) = L(N,N-1)L(N-1,N-2)\ldots L(1,0),
\end{equation}
with at time discretization $\Delta t$ and with time-bin normalized operators
\begin{equation}
L(n,n-1) = \exp \left[\sqrt{\Delta t} \int_{(n-1)\Delta t}^{n \Delta t} \mathrm{d} t^\prime \ \mathcal{L}(t^\prime)  \right].
\end{equation}
For the MPS evolution of the density matrix during each time step $n$, the discrete Liouvillian time step operator $L(n,n-1)$ is approximated as a \textit{tenth order} series expansion, i.e.,
\begin{equation}
L(n,n-1) = \sum_{m=0}^{10} \dfrac{\sqrt{\Delta t}^m}{m!} \left[ \int_{(n-1)\Delta t}^{n \Delta t} \mathrm{d} t^\prime \ \mathcal{L}(t^\prime) \right]^m.
\label{eq:liou_evo}
\end{equation}

Figs.~\ref{fig:mps_feedback}(a)-(e) show the tensor network scheme for the implementation of the time-discrete memory.
The square red tensor in Fig.~\ref{fig:mps_feedback}(a) contains the system density matrix at the initial time $t=0$. To consider a time-discrete memory, here $n_{d}=4$ circular tensors to its \textit{left} store the reservoir state in Liouville space at preceding times, with the oldest state located on the left end of the MPS (blue).
The reservoir states for all future time steps are initialized to the \textit{right} of the system bin, containing full reservoir entanglement e.g. at finite temperature. In the following we assume an initial vacuum state. Therefore, during each time step a new empty reservoir bin (orange) is added to the MPS from the right, representing the \textit{present} reservoir state [see Fig.~\ref{fig:mps_feedback}(a)].
The memory loop realization explained in detail below introduces the retardation time $\tau = n_d \Delta t$ by the number of initial memory bins $n_{d}$.

The first step of the time evolution is carried out as follows: By applications of the singular value decomposition algorithm~\cite{Schollwoeck2011}, the first memory bin (blue) is pushed to the left of the system (red) while maintaining relevant entanglement information in the swapping procedure [see Fig.~\ref{fig:mps_feedback}(a)]. The Liouvillian operator $L(1,0)$ for the first time step is then applied to the system bin, current memory bin and present reservoir bin, as shown in Fig.~\ref{fig:mps_feedback}(b).
Afterwards, the processed memory bin (grey) is swapped back to its original position and stored for the rest of the time evolution. The updated present reservoir bin (green) is pushed to the left, taking the role of a memory bin [see Fig.~\ref{fig:mps_feedback}(c),(d)]. Fig.~\ref{fig:mps_feedback}(d) shows the MPS after completion of the first time step. The second time step is carried out in the same fashion, as shown in Fig.~\ref{fig:mps_feedback}(e).
After completion of $n_{d}$ time steps, all initial memory bins have been processed. At step $n_{d}+1$, the reservoir bin modified during the first time step [green bin in Fig.~\ref{fig:mps_feedback}(d)] becomes the current memory bin, containing information of a previous system state and setting off reservoir-induced memory effects in the system in complete agreement with the time-ordered problem.

As a first benchmark for the presented time-discrete quantum memory in Liouville space, we calculate the system dynamics imposed by Eq.~\eqref{eq:H_D}. Fig.~\ref{fig:mps_feedback}(f) shows the unfolding emitter population dynamics at a feedback time $\tau=3.0\,$ps and $\Gamma=31.6\,\mathrm{ps}^{-1}$ calculated using the MPS implementation (solid blue line) and compared to its analytical solution up to $t=3\tau$ (dashed grey line). The latter is given by $\left\langle \sigma_{11}(t)\right\rangle=|\left\langle \sigma_{10}(t)\right\rangle|^2$ only valid in the single-excitation regime with
\begin{equation}
\left\langle \sigma_{01}(t)\right\rangle = \sum_{n=0}^\infty \frac{e^{-\Gamma t}}{n!}
\!
\left[ \Gamma e^{(\Gamma-i\omega_0)\tau} (t \! - \! n \tau) \right]^n \Theta(t \! - \! n \tau),
\label{eq:fb_ana}
\end{equation}
as calculated in \cite{Dorner2002,Kabuss2015,Carmele2020} and shows excellent agreement with the numerical result. In this regime, i.e. $\Gamma\tau\gg1$, the delay in the amplitude governs the dynamics and leads to re-excitations at multiples of the round trip time $\tau$. The phase of the amplitude $\phi=\omega_0\tau$, however, loses importance in the first $\tau$-intervals due to a stronger decay of the mixing terms in the absolute square of Eq.~\eqref{eq:fb_ana}.
The advantage of our implementation of quantum feedback in Liouville space becomes evident if the impact of phase destroying processes is in question. Up until now, this impact has only been investigated for a special case, finding the emergence of an Ornstein-Uhlenbeck process during the first $\tau$-intervals~\cite{Carmele2020}. These results have been obtained via analytical calculations, limiting the investigation to a small number of feedback intervals. Steady-state scenarios, however, are out of reach in this case as the evaluation of the phase-noise kernels must be done analytically. In our method, these limits have been overcome. Due to the here presented Liouville architecture, additional Lindblad-based dissipation can be easily implemented without increased numerical expense.

In Fig.~\ref{fig:mps_feedback}(g), we present the dynamics of a decaying, intially excited two-level emitter under the influence of quantum coherent feedback and additional phenomenological dephasing at rate $\gamma$, realized by adding a Lindblad dissipator to the Liouvillian [Eq.~\eqref{eq:von_Neumann}] and $H=H_D$,~\cite{Breuer2002, Mukamel1999}
\begin{equation}
\mathcal{D}[ \sqrt{\gamma/2} \tilde{\sigma}_{11}] \rho(t) = \dfrac{\gamma}{2} [ 2 \tilde{\sigma}_{11} \rho(t) \tilde{\sigma}_{11} - \{ \rho(t), \tilde{\sigma}_{11} \}],
\label{eq:lindblad}
\end{equation}
with a redefined system operator in full configuration space, $\tilde{\sigma}_{11} = \mathbb{1}_{D} \sigma_{11} \mathbb{1}_{D}$, including the time-discrete reservoir basis $\mathbb{1}_{D} = \int \mathrm{d} k \sum_{n=0}^\infty \ket{\{ n_k\}} \bra{\{ n_k\}}$.
The time trace shows long time calculations of the emitter population, comparing the cases $\gamma=0$ (blue line) and $\gamma=0.5\,\mathrm{ps}^{-1}$ (orange line). Most importantly, we see that the pure dephasing process becomes important only after the feedback signal re-excitates the emitter and the stabilization of the incoming and outgoing phase comes into play. In the presence of an additional pure dephasing $\gamma \neq 0$, the initial decay process is unchanged but the re-excitation becomes less efficient until only incoherent re-excitation takes place, leading to a faster decay to zero without population trapping, regardless of the choice of $\phi$.

This important result sheds light on the robustness of quantum feedback processes in the presence of additional Markovian dissipation channels. As expected, additional Markovian decoherence leads to a washing out of the signal since a loss of quantum feedback-induced coherence is inevitable. However, this does not have to be the case in the presence of an additional non-Markovian dissipation channel, which we discuss in the following.

\section{Quasi-2D tensor network} \label{sec:2dnetwork}

As a next step, we expand upon the MPS architecture for time-discrete quantum memory in Liouville space by combining it with the previously discussed tensor network-based path integral implementation for continuous harmonic reservoirs, resulting in a \textit{quasi-2D tensor network}.
The technical connection of the networks via link indices which store arising entanglement information enables the numerically exact description of correlation buildup \textit{in between the reservoirs}. When considering a scenario involving two non-Markovian reservoirs not isolated from each other, such inter-reservoir correlations may have fundamental impact on the system dynamics, therefore prohibiting a strict truncation of the arising inter-reservoir entanglement, e.g. in the form of a low Schmidt value cutoff precision $d_{cut}$. As a result of not only two system-reservoir interactions but additionally arising reservoir-reservoir entanglement, the overall grade of entanglement in the system rises \textit{intensively} with respect to the twofold single reservoir case. On the other hand, in setups where two non-Markovian reservoirs are present but do not crucially interact, e.g. via dynamical decoupling, a much more restrictive truncation is possible without cost of accuracy. For the presented results, we have employed a high Schmidt value cutoff precision $d_{cut}=10^{-12}$, such that no relevant entanglement information is lost during the time evolution. The dynamical interplay between the two reservoirs with the system and with each other poses an immense numerical challenge and strongly limits the accessible memory depths in the here considered system: While the two presented algorithms by themselves enable simulations of a single reservoir with deep memories, their combination is accompanied by limitations due to the arising inter-reservoir entanglement. As a result, the combined number of memory bins in the quasi-2D network is limited to $n_c+n_d < 20$ for our model of choice, as is the case in traditional path integral implementations for a single reservoir~\cite{Caldeira1983, Leggett1987, Vagov2011, Glaessl2013, Barth2016, Cosacchi2018}.
However, we stress once more that this limitation is a natural consequence of the high grade of entanglement in between the two reservoirs and the system. With the presented quasi-2D network architecture, we take first steps to unravel the mostly unexplored field of multiple interacting non-Markovian reservoirs by explicitly considering memory-enabled information backflow in between them.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{tensornetworks_connect6}
\caption{Quasi-2D tensor network implementation. (a) Tensor network-based path integral implementation for continuous reservoirs. (b) System MPS containing current (red) and preceding states (grey) after the first network contraction. To realize the quasi-2D network, the current state tensor now features a link connection to the discrete memory MPS (blue diagonal line). (c) Time-discrete memory MPS at the beginning of the first time step, with the common current system state acting as a junction between the reservoirs (dashed circles).}
\label{fig:mps_connect}
\end{figure}

The construction of the quasi-2D network is sketched in Fig.~\ref{fig:mps_connect}: The two tensor networks for the continuous [Figs.~\ref{fig:mps_connect}(a),(b)] and time-discrete reservoirs [Fig.~\ref{fig:mps_connect}(c)] are connected to each other via the common tensor representing the current system state in both MPS algorithms (red shape). The system state tensor is employed to act as a junction connecting the two reservoirs [dashed circles in Figs.~\ref{fig:mps_connect}(b),(c)] and thereby enables the buildup and storage of inter-reservoir correlations in the connecting link indices.
The time evolution of the quasi-2D network is carried out as follows: During each time step, the system is first evolved under the influence of the continuous reservoir by a single contraction of the network, as shown in Fig.~\ref{fig:mps_connect}(a) (dashed frame). Afterwards, the new current system state [red shape in Fig.~\ref{fig:mps_connect}(b)] is subjected to the second tensor network algorithm accounting for the time-discrete reservoir [Fig.~\ref{fig:mps_connect}(c)]. In consequence, the quasi-2D network stores the history of both interactions, maintaining crucial entanglement information and enabling the calculation of two dynamically interacting time-delayed processes.
As an example application for the quasi-2D network, in the following we investigate the interplay of off-diagonal coherent quantum feedback and a diagonal reservoir of independent oscillators. As illustrated below, for certain memory depths and initial states, this results in dynamical protection of coherent quantum feedback properties in the open system.

\section{Memory-induced dynamical population trapping} \label{sec:results}

Coherent quantum feedback mechanisms exhibit a rich variety of non-Markovian phenomena~\cite{Grimsmo2015, Wilson2003, Hetet2011, Hoi2015, Kabuss2015, Svidzinsky2018, Barkemeyer2019}, e.g. enabling coherent population trapping~\cite{Dubin2007, Glaetzle2010, Carmele2013, Guimond2016, Nemet2019, CarmeleNemet2020}, Ornstein-Uhlenbeck-type events in the presence of white noise~\cite{Carmele2020}, and formation of large entangled photon states~\cite{Pichler2017}. However, so far these effects have not been explored in the presence of additional non-Markovian decoherence or dissipative channels.
To investigate the impact of dephasing on feedback-induced decoherence, we consider a two-level emitter placed in front of a mirror with a round trip time $\tau$, taking the role of a time-discrete reservoir (see Fig.~\ref{fig:system}) described by Eq.~\eqref{eq:H_D}.
The photon-induced feedback imprints a time-delayed coherence in the form of a feedback phase $\varphi=\omega_0 \tau / (2 \pi)$ on the system, critically influencing its dynamics. It is given by the delay time $\tau$ and the transition frequency $\omega_0$ of the electronic coherence operator $\sigma_{12}$.
To study a pronounced quantum optical effect, we consider the case of an initially excited two-level system where coherent population trapping occurs as a result of a bound state in continuum at feedback phases $\varphi \in \mathbb{Z}$~\cite{Dubin2007, Glaetzle2010, Carmele2013, Guimond2016, Nemet2019, CarmeleNemet2020, Carmele2020}. The unfolding dynamics are evaluated for parameters $\Gamma=0.9\, \mathrm{ps}^{-1}$, $\tau = 1.2\,$ps and $n_{d}=4$, resulting in a time discretization $\Delta t = 0.3\,$ps, e.g. typical for semiconductor quantum dot based devices~\cite{Carmele2013, Nemet2019, Barkemeyer2019}. Moreover, we get $\Gamma \tau \approx 1.1$, corresponding to the strong non-Markovian regime~\cite{baranger_2019, Pichler2016, Carmele2020}. The employed series expansion of the Liouvillian up to tenth order [see Eq.~\eqref{eq:liou_evo}] justifies this coarse time discretization, making $n_d=4$ time bins sufficient for our investigation while resulting in convergent results (see Appendix~\ref{app:convergence}).

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{FBDeph_gam003_phi117_lindblad2}
\caption{Excited state population dynamics of the two-level emitter under time-discrete photon feedback at varying feedback phases $\varphi=\omega_0 \tau/(2\pi)$ and constant dephasing rates $\gamma$.}
\label{fig:lindblad}
\end{figure}

To illustrate the power of our method, we compare the cases of Markovian and non-Markovian dephasing introduced by an additional diagonal system-reservoir coupling, representing the continuous reservoir in Fig.~\ref{fig:system}. As a first step, we calculate the system dynamics given by Eq.~\eqref{eq:H_D} in the presence of phenomenological dephasing at rate $\gamma$, introduced by the Lindblad dissipator stated in Eq.~\eqref{eq:lindblad}.
Fig.~\ref{fig:lindblad} shows resulting excited state population dynamics at varying $\gamma$ and feedback phases $\varphi$. At $\gamma=0$ there exists a periodic ideal feedback phase $\varphi \in \mathbb{Z}$ such that the system decouples from its environment by constructive interference, resulting in coherent population trapping (solid blue line in Fig.~\ref{fig:lindblad}).
Choosing a nonzero dephasing $\gamma=0.001\,\mathrm{ps}^{-1}$ has no impact on the population dynamics until feedback sets in, since the radiative decay is frequency-independent until $t=\tau$ (dashed blue line). Thereafter, phenomenological dephasing destroys the phase interference and with it the trapping mechanism, resulting in an asymptotic decline of occupation to zero.
At a nonideal feedback phase, here $\varphi=1.17$, and no dephasing, destructive interference leads to an asymptotic decline towards zero as well (solid orange line). The decay is further accelerated by setting $\gamma > 0$ (dashed orange line), since any phenomenological decoherence attacking the phase relation $\varphi$ results in faster decay. In conclusion, choosing a feedback phase $\varphi \notin \mathbb{Z}$ without a structured phonon reservoir inevitably results in asymptotic population decay via thermalization, and a Lindblad formulation of decoherence never preserves quantum correlations between the reservoir and system states.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{FBDeph_gam003_phi117_path2}
\caption{Excited state population dynamics of the two-level emitter subject to photon and phonon reservoir interactions at $\varphi \neq 1$ and varying temperatures $T$. Inset: Dimension of the link index connecting the discrete memory MPS to the system in the presence of phonons (red), without phonon coupling (blue) and for phenomenological dephasing (orange).}
\label{fig:path_dynamics}
\end{figure}

As a next step, we show this is not necessarily the case if the decoherence process itself is the result of a non-Markovian reservoir interaction. Using the quasi-2D tensor network, we calculate the emitter dynamics imposed by $H_D+H_C$ [Eqs.~\eqref{eq:H_C},~\eqref{eq:H_D}] at $n_c=4$ and leaving all remaining parameters unchanged. Fig.~\ref{fig:path_dynamics} shows resulting population dynamics at different temperatures. For the chosen parameters, we find population trapping for the non-ideal feedback phase $\varphi=1.17$, i.e. $\varphi \notin \mathbb{Z}$, at $T=4\,$K (solid line).
Time-delayed excitation backflow from the continuous reservoir of oscillators to the system enables a decoupling from destructive interference with the time-discrete photon environment. This information backflow results in correlation buildup and information exchange between the reservoirs, dynamically protecting feedback-induced coherence in the system for long times.
As long as diffusion processes at finite temperature take place on a comparable time scale as the coherent feedback dynamics, we always find dynamical population trapping by tuning of $\varphi$ after a typical excitation backflow time. At higher temperatures, it is to be expected that thermal properties of the phonon reservoir start to dominate the dynamics: Dashed and dotted lines in Fig.~\ref{fig:path_dynamics} show corresponding thermalization dynamics at $T=30\,$K and $T=77\,$K, respectively, exhibiting population decay. The temperature dependence clearly shows that correlation lengths within the full system-reservoir dynamics are of importance, and the observed effect allows to probe these otherwise inaccessible microscopic environmental properties.
This formation of self-stabilizing dissipative structures is closely related to a localized phase stabilization in the coherent driving case $\Omega_0\neq0$ at Ohmic spectral densities and without photons~\cite{Leggett1987, Strathearn2018}. There, above a critical coupling strength the system transitions into a localized phase with nonzero steady state population rather than decaying to zero. In our case, a similiar phenomenon is established with incoherent feedback instead of coherent external driving, addressing the localized phase stabilization process from a dissipative non-Markovian side.
The complexity of this phenomenon is illustrated in the inset in Fig.~\ref{fig:path_dynamics}, showing the dimension of the link index connecting the current open system bin to the discrete memory MPS over time (see Fig.~\ref{fig:mps_connect}). The red line shows the case including the continuous reservoir. After slowly increasing during times $t < \tau$, it exhibits a vast increase once feedback sets in and quickly reaches a maximum due to finite memory. The high grade of entanglement between the two reservoirs even at arguably low finite memory sizes $n_d=4$, $n_c=4$ underlines the crucial role of the interplay between the two non-Markovian processes for the observed protection of coherence.
Switching off the phonon coupling, $g_{\bm{q}}=0$, results in a much lower maximum link dimension (blue line), as no entanglement between the reservoirs arises. For phenomenological dephasing (orange line), the link connecting system and memory bins has an even lower dimension due to the highly decreased complexity of the then 1D network.

\section{Conclusions} \label{sec:conclusions}

We have presented an MPS algorithm for the description of a time-discrete quantum memory in Liouville space. By combining this technique with a path integral tensor network implementation for continuous non-interacting harmonic reservoirs, we have established a quasi-2D tensor network, allowing for simulations of quantum systems subject to two non-Markovian environments while maintaining crucial entanglement information in the coupled system with both diagonal and off-diagonal system-reservoir interactions.
Due to arising reservoir-reservoir correlations, system correlations scale intensively with respect to the twofold single reservoir case. In consequence, the achievable memory depth is limited to $n_c+n_d < 20$ in our study. However, appropriate tuning of the relevant system and reservoir time scales via the employed parameters still opens up a wide array of accessible systems and scenarios where numerical convergence can be achieved. The next step will be to trace out the time-discrete feedback bins as well after their interaction with the systems' degrees of freedom to further improve numerical efficiency. This will allow for longer delay times and therefore increased time discretizations without changing the qualitative results. Hence, the presented quasi-2D tensor architecture is a first step towards unraveling the mostly unexplored field of multiple interacting non-Markovian reservoirs in a numerically exact fashion.

As an example application, we have demonstrated that the interplay of a structured phonon reservoir and photon feedback can dynamically protect the system from destructive interference by time-delayed backflow of coherence, resulting in dynamical population trapping. Tuning the non-Markovian interactions with respect to each other allows for the formation of inter-reservoir correlations, dynamically preservering feedback-induced coherence in the system.
These findings have implications for the fields of quantum thermodynamics and nonequilibrium physics, as well as dynamical quantum phase transitions with ergodic, entropic or negentropic information exchange, where taking account of such dissipative structures may unravel new phenomena. Future works will aim to advance our architecture to a full 2D representation via projected entangled pair states with combined memory bins, potentially allowing for simulations of multi-level systems at improved time resolutions.