\section{Introduction}
\label{sec:intro}
The problem of privacy preserving machine learning has become increasingly important. Recently, there have been many works that have made rapid strides towards realizing {\em secure inference}~\cite{cryptonets,chet,secureml,minionn,gazelle,delphi,nitin,xonn,ezpc,hycc,nhe1,deepsecure,ball}.
Consider a server that holds the weights $w$ of a publicly known deep neural network (DNN), $F$, that has been trained on private data. A client holds a private input $x$; in a standard machine learning (ML) inference task, the goal is for the client to learn the prediction $F(x,w)$ of the server's model on the input $x$. In secure inference, the inference is performed with the guarantee that the server learns {\em nothing} about $x$ and the client learns nothing about the server's model $w$ beyond what can be deduced from $F(x,w)$ and $x$.

A solution for secure inference that scales to practical ML tasks would open a plethora of applications based on MLaaS (ML as a Service). Users can obtain value from ML services  without worrying about the loss of their private data, while model owners can effectively monetize their services with no fear of breaches of client data (they never observe private client data in the clear). Perhaps the most important emerging applications for secure inference  are in healthcare where
prior work~\cite{cryptflow,xonn,nitin} has explored secure inference services for privacy preserving medical diagnosis of chest diseases, diabetic retinopathy, malaria, and so on.

Secure inference is an instance of secure 2-party computation (2PC) and cryptographically secure general protocols for 2PC have been known for decades~\cite{Yao,gmw}.
However, secure inference for practical ML tasks, e.g., ImageNet scale prediction~\cite{imagenet}, is challenging for two reasons:
a) realistic DNNs use $\mathsf{ReLU}$ activations\footnote{$\mathsf{ReLU}(x)$ is defined as $\max(x,0)$.} that are expensive to compute securely; and b) preserving inference accuracy requires a faithful implementation of secure fixed-point arithmetic.
 All prior works~\cite{cryptonets,secureml,minionn,gazelle,delphi,ball} fail to provide efficient implementation of $\mathsf{ReLU}$s. Although $\mathsf{ReLU}$s can be replaced with approximations that are more tractable for 2PC~\cite{cryptonets,delphi,chet}, this approach results in significant accuracy losses that can degrade user experience. The only known approaches to evaluate $\mathsf{ReLU}$s efficiently require sacrificing security by making the untenable assumption that a non-colluding third party takes part in the protocol~\cite{chameleon,securenn,cryptflow,aby3,quantizednn} or by leaking activations~\cite{nhe2}.
Moreover, some prior works~\cite{secureml,delphi,securenn,aby3,cryptflow} even sacrifice correctness of their fixed-point implementations and the result of their secure execution can sometimes diverge from the expected result, i.e. cleartext execution, in random and unpredictable ways.
Thus, correct and efficient 2PC protocols for secure inference over realistic DNNs remain elusive.
\subsection{Our Contributions}

In this work, we address the above two challenges and build new semi-honest secure 2-party cryptographic protocols for secure computation of DNN inference.
Our new efficient protocols enable the first secure implementations of ImageNet scale inference that complete in under a minute!
We make three main  contributions:
\begin{tiret}
\item First, we give new protocols for millionaires' and $\mathsf{DReLU}$\footnote{$\mathsf{DReLU}$ is the derivative of $\mathsf{ReLU}$, i.e., $\mathsf{DReLU}(x)$ is $1$ if $x \geq 0$ and $0$ otherwise.} that enable us to securely and efficiently evaluate the non-linear layers of DNNs such as $\mathsf{ReLU}$, $\mathsf{Maxpool}$ and $\mathsf{Argmax}$. \item Second, we provide new protocols for division. Together with new theorems that we prove on fixed-point arithmetic over shares, we show how to evaluate linear layers, such as convolutions, average pool and fully connected layers, faithfully.
\item Finally, by providing protocols that can work on a variety of input domains, we build a system\footnote{Implementation is available at \url{https://github.com/mpc-msri/EzPC}.}
    ${\textsc{CrypTFlow2}}\xspace$ that supports two different types of Secure and Correct Inference ({\textsc{SCI}}\xspace) protocols where linear layers can be evaluated using either homomorphic encryption ($\textsc{SCI}_{\mathsf{HE}}$) or through oblivious transfer ($\textsc{SCI}_{\mathsf{OT}}$).
\end{tiret}
We now provide more details of our main contributions.
\\\\
\noindent{\textbf{New millionaires' and $\mathsf{DReLU}$ protocols.}}
Our first main technical contribution is a novel protocol for the well-known {\em millionaires'} problem~\cite{Yao}, where parties \ensuremath{P_0}\xspace and \ensuremath{P_1}\xspace hold  $\ell-$bit integers $x$ and $y$, respectively, and want to securely compute $x < y$ (or, secret shares of $x<y$).
The theoretical communication complexity of our protocol is $\approx 3\times$ better than the most communication efficient prior millionaires' protocol~\cite{Cou18,Yao,gmw,emp-toolkit,GSV07}.
In terms of round complexity, our protocol executes in $\log \ell$ rounds (e.g. $5$ rounds for $\ell=32$ bits);
see \tableref{comp-mill} for a detailed comparison and \cite{Cou18} for a detailed overview of the costs of other comparison protocols.

Using our protocol for millionaires' problem, we build new and efficient protocols for computing $\mathsf{DReLU}$ for both $\ell-$bit integers (i.e., $\mathbb{Z}_L$, $L = {2^\ell}$) and general rings $\mathbb{Z}_n$.
Our protocol for $\mathsf{DReLU}$ serves as one of the main building blocks for non-linear activations such as $\mathsf{ReLU}$ and $\mathsf{Maxpool}$, as well as division over both input domains.
Providing support for  $\ell-$bit integers $\mathbb{Z}_{L}$ as well as arbitrary rings $\mathbb{Z}_{n}$, allows us to securely evaluate the linear layers (such as matrix multiplication and convolutions) using  the approaches of Oblivious Transfer (OT)~\cite{beaver,secureml} as well as Homomorphic Encryption (HE)~\cite{gentryfhe,gazelle,delphi}, respectively.
This provides our protocols great flexibility when executing over different network configurations.
Since all prior work~\cite{secureml,minionn,gazelle,delphi} for securely computing these activations rely on Yao's garbled circuits~\cite{Yao}, our protocols are much more efficient in both settings.
Asymptotically, our $\mathsf{ReLU}$ protocol over $\mathbb{Z}_L$ and $\mathbb{Z}_n$ communicate $\approx 8\times$ and $\approx 12\times$ less bits than prior works~\cite{Yao,secureml,minionn,gazelle,delphi,emp-toolkit} (see \tableref{comp-relu} for a detailed comparison).
Experimentally, our protocols are at least an order of magnitude more performant than prior protocols when computing $\mathsf{ReLU}$ activations at the scale of ML applications.
\\\\
\noindent{\textbf{Fixed-point arithmetic.}}
The ML models used by all prior works on secure inference are expressed
using fixed-point arithmetic; such models can be obtained from~\cite{tflite,bnn,cryptflow,dfq}. A faithful implementation of fixed-point arithmetic is quintessential to ensure that the secure computation is {\em correct}, i.e., it is {\it equivalent} to the cleartext computation for all possible inputs.
Given a secure inference task $F(x,w)$, some prior works~\cite{secureml,delphi,securenn,aby3,cryptflow} give up on correctness
when implementing division operations and instead compute an approximation $F'(x,w)$. In fixed-point arithmetic, each multiplication requires a division by a power-of-2 and multiplications are used pervasively in linear-layers of DNNs.
 Moreover, layers like average-pool require division for computing means.
Loss in correctness is worrisome as the errors can accumulate and $F'(x,w)$ can be arbitrarily far from $F(x,w)$.
 Recent work~\cite{delphi} has shown that even in practice the approximations can lead  to significant losses in classification accuracy.

As our next contribution, we provide novel protocols to compute division by power-of-2 as well as division by arbitrary integers that are both correct and efficient. The inputs to these protocols can be encoded over both $\ell-$bit integers $\mathbb{Z}_{L}$ as well as $\mathbb{Z}_{n}$, for arbitrary $n$.
To the best of our knowledge, the only known approach to compute division correctly is via garbled circuits which we compare with in Table~\ref{tab:comp-avgpool}.
While garbled circuits based protocols require communication which is quadratic in $\ell$ or $\log n$, our protocols are asymptotically better and incur only linear communication. Concretely, for average pool with $7\times 7$ filters and 32-bit integers, our protocols have $\approx 54\times$ less communication.
\\\\
\noindent{\textbf{Scaling to practical DNNs.}} These efficient protocols, help us securely evaluate practical DNNs like \textsf{SqueezeNet}\xspace\   on ImageNet scale classification tasks in under a minute. In sharp contrast, all prior works on secure  2-party inference (\cite{cryptonets,chet,secureml,minionn,gazelle,delphi,nitin,xonn,ezpc,hycc,nhe1,deepsecure,ball}) has been limited to small DNNs on tiny datasets like MNIST and CIFAR.
While MNIST deals with the task of classifying black and white handwritten digits given as $28\times 28$ images into the classes 0 to 9, ImageNet tasks are much more complex: typically $224\times 224$ colored images need to be classified into thousand classes (e.g., agaric, gyromitra, ptarmigan, etc.) that even humans can find challenging .
 Additionally, our work is the first to securely evaluate practical {\em convolutional neural networks} (CNNs) like \textsf{ResNet50}\xspace\ and \textsf{DenseNet121}\xspace; these DNNs are at least an order of magnitude larger than the DNNs considered in prior work, provide over $90\%$ Top-5 accuracy on ImageNet, and have also been shown to predict lung diseases from chest X-ray images~\cite{cryptflow, chestxray2018}.
Thus, our work provides the first implementations of practical ML inference tasks running securely.  Even on the smaller CIFAR scale DNNs, our protocols require an order of magnitude less communication and $20\times$-$30\times$ less time than the state-of-the-art~\cite{delphi} (see  Section~\ref{sec:exp-comparison-prior-work}).
\\\\
\noindent{\textbf{OT vs HE.}} Through our evaluation, we also resolve the OT vs HE conundrum:
although the initial works on secure inference~\cite{secureml,minionn} used OT-based protocols for evaluating convolutions, the state-of-the-art protocols~\cite{gazelle,delphi}, which currently provide the best published inference latency, use HE-based convolutions. HE-based secure inference has much less communication than OT but HE requires more computation.
Hence, at the onset of this work, it was not clear to us whether HE-based convolutions would provide us the best latency for ImageNet-scale benchmarks.

To resolve this empirical question, we implement two classes of protocols, $\textsc{SCI}_{\mathsf{OT}}$ and $\textsc{SCI}_{\mathsf{HE}}$, in {\textsc{CrypTFlow2}}\xspace.
In $\textsc{SCI}_{\mathsf{OT}}$, inputs are in $\mathbb{Z}_{L}$ ($L = 2^\ell$, for a suitable choice of $\ell$). Linear layers such as matrix multiplication and convolution are performed using OT-based techniques~\cite{beaver,secureml}, while the activations such as $\mathsf{ReLU}$, $\mathsf{Maxpool}$ and $\mathsf{Avgpool}$  are implemented using our new protocols over $\mathbb{Z}_{L}$.
In $\textsc{SCI}_{\mathsf{HE}}$, inputs are encoded in an appropriate prime field $\mathbb{Z}_{n}$ (similar to~\cite{gazelle,delphi}).
Here, we compute linear layers using homomorphic encryption and the activations  using our protocols over $\mathbb{Z}_{n}$.
In both $\textsc{SCI}_{\mathsf{OT}}$ and $\textsc{SCI}_{\mathsf{HE}}$ faithful divisions after linear layers are performed using our new protocols over corresponding rings. Next, we evaluate ImageNet-scale inference tasks with both $\textsc{SCI}_{\mathsf{OT}}$ and $\textsc{SCI}_{\mathsf{HE}}$ . We observe that in a WAN setting, where communication is a bottleneck, HE-based inference is always faster and in a LAN setting OT and HE are incomparable.
\begin{table}
  \centering
      \begin{tabular}{|c|c|c|c|}
    \hline
				Layer & Protocol & Comm. (bits) & Rounds \\ \hline \hline
				\multirowcell{4}{Millionaires' \\ on $\{0,1\}^\ell$}
				& GC \cite{Yao, emp-toolkit} & $4\lambda\ell$ & 2 \\ \cline{2-4}
				& GMW$\footnotemark/$GSV~\cite{gmw,GSV07}& $\approx 6\lambda\ell$ & $\log \ell+3$ \\ \cline{2-4}
				& SC3\footnotemark\cite{Cou18} & $> 3\lambda\ell$ & $\approx 4\log^* \lambda$ \\ \cline{2-4}
				& This work ($m = 4$) & $< \lambda \ell + 14\ell $ & $\log {\ell}$\\ \hline \hline
				\multirowcell{5}{Millionaires' \\example \\$\ell = 32$}
				& GC \cite{Yao, emp-toolkit} & 16384 & 2 \\ \cline{2-4}
				& GMW$/$GSV \cite{gmw, GSV07} & 23140 & 8 \\ \cline{2-4}
				& SC3 \cite{Cou18} & 13016 & 15 \\ \cline{2-4}
				& This work ($m = 7$) & 2930 & 5 \\ \cline{2-4}
				& This work ($m = 4$) & 3844 & 5 \\ \hline
			\end{tabular}
 \caption{Comparison of communication with prior work for millionaires' problem. For our protocol, $m$ is a parameter. For concrete bits of communication we use $\lambda=128$. }
\label{tab:comp-mill}
\end{table}
\addtocounter{footnote}{-1}
\footnotetext{Here we state the communication numbers for GMW~\cite{gmw} for a depth-optimized circuit. The circuit that would give the best communication would still have a complexity of $> 2\lambda\ell$ and would additionally pay an inordinate cost in terms of rounds, namely $\ell$.}
\addtocounter{footnote}{1}
\footnotetext{Couteau~\cite{Cou18} presented multiple protocols; we pick the one that has the best communication complexity.}
 \begin{table}
  \centering
      \begin{tabular}{|c|c|c|c|}
    \hline
				Layer & Protocol & Comm. (bits) & Rounds \\ \hline \hline
				\multirowcell{2}{$\mathsf{ReLU}$ for \\$\mathbb{Z}_{2^\ell}$}

				& GC \cite{Yao, emp-toolkit} & $ 8\lambda\ell-4\lambda$ & 2 \\ \cline{2-4}
				& This work & \multirow{1}{*}{ $< \lambda\ell+18\ell$ } & $\log \ell + 2$\\
				 \hline

				\multirowcell{2}{$\mathsf{ReLU}$ for \\general $\mathbb{Z}_n$}

				& GC \cite{Yao, emp-toolkit} & $18\lambda\eta-6\lambda$ & 2 \\ \cline{2-4}
				& This work & $< \frac{3}{2}\lambda (\eta+1) + 31\eta$ & $\log \eta + 4$\\
				\hline \hline

				\multirowcell{2}{$\mathsf{ReLU}$ for \\ $\mathbb{Z}_{2^\ell}$, $\ell = 32$}

				& GC \cite{Yao, emp-toolkit} & 32256 & 2 \\ \cline{2-4}
				& This work  & 3298 & 7 \\ \cline{2-4}
				\hline

				\multirowcell{2}{$\mathsf{ReLU}$ for \\$\mathbb{Z}_n$, $ \eta = 32$}

				& GC \cite{Yao, emp-toolkit} & 72960 & 2 \\ \cline{2-4}
				& This work  & 5288 & 9 \\ \cline{2-4}
				\hline
			\end{tabular}
 \caption{Comparison of communication with garbled circuits for $\mathsf{ReLU}$. We define $\eta = \lceil \log n \rceil$. For concrete bits of communication we use $\lambda=128$. }
\label{tab:comp-relu}
\end{table}
 \begin{table}
  \centering
      \begin{tabular}{|c|c|c|c|}
    \hline
				Layer & Protocol & Comm. (bits) & Rounds \\ \hline \hline
				\multirowcell{2}{$\mathsf{Avgpool}_d$ \\$\mathbb{Z}_{2^\ell}$}

				& GC \cite{Yao, emp-toolkit} & $2\lambda(\ell^2 + 5\ell -3)$ & 2 \\ \cline{2-4}
				& This work &  $< (\lambda + 21) \cdot ( \ell + 3 \delta )$   & $\log (\ell\delta ) + 4$\\
              \hline

				\multirowcell{2}{$\mathsf{Avgpool}_d$  \\ $\mathbb{Z}_n$}

				& GC \cite{Yao, emp-toolkit} & $2\lambda(\eta^2 + 9\eta-3)$ & 2 \\ \cline{2-4}
				& This work & $< (\textstyle\frac{3}{2} \lambda + 34) \cdot (\eta + 2 \delta)$ & $\log (\eta\delta) + 6$\\
             \hline \hline

                \multirowcell{2}{$\mathsf{Avgpool}_{49}$ \\  $\mathbb{Z}_{2^\ell}, \ell = 32$}

				& GC \cite{Yao, emp-toolkit} & 302336 & 2 \\ \cline{2-4}
				& This work  & 5570 & 10 \\ \cline{2-4}
				\hline

                \multirowcell{2}{$\mathsf{Avgpool}_{49}$  \\ $\mathbb{Z}_n, \eta = 32$}

				& GC \cite{Yao, emp-toolkit} & 335104 & 2 \\ \cline{2-4}
				& This work & 7796 & 14 \\ \cline{2-4}
				\hline
			\end{tabular}
 \caption{Comparison of communication with garbled circuits for $\mathsf{Avgpool}_d$. We define $\eta = \lceil \log n \rceil$ and $\delta = \lceil \log (6 \cdot d) \rceil$. For concrete bits of communication we use $\lambda=128$. Choice of $d =49$ corresponds to average pool filter of size $7\times 7$. } \label{tab:comp-avgpool}
\end{table}
 \subsection{Our Techniques}

\noindent\textbf{Millionaires'.} Our protocol for securely computing the millionaires' problem (the bit $x<y$) is based on the following observation (first made in~\cite{GSV07}).
Let $x = x_1||x_0$ and $y = y_1||y_0$ (where $||$ denotes concatenation and $x_1, y_1$ are strings of the same length).
Then, $x<y$ is the same as checking if either $x_1<y_1$ or $x_1=y_1$ and $x_0<y_0$.
Now, the original problem is reduced to computing two millionaires' instances over smaller length strings ($x_1<y_1$ and $x_0<y_0$) and one equality test ($x_1 = y_1$).
By continuing recursively, one could build a tree all the way where the leaves are individual bits, at which point one could use 1-out-of-2 OT-based protocols to perform the comparison/equality.
However, the communication complexity of this protocol is still quite large.
We make several important modifications to this approach.
First, we modify the tree so that the recursion is done $\log (\ell/m)$ times to obtain leaves with strings of size $m$, for a parameter $m$.
We then use 1-out-of-$2^m$ OT to compute the comparison/equality at the leaves, employing the lookup-table based approach of \cite{DKSSZZ17}.
Second, we observe that by carefully setting up the receiver's and sender's messages in the OT protocols for leaf comparisons and equality, multiple 1-out-of-$2^m$ OT instances can be combined to reduce communication.
Next, recursing up from the leaves to the root, requires securely computing the $\mathsf{AND}$ functionality\footnote{This functionality takes as input shares of bits $x, y$ from the two parties and outputs shares of $x~\mathsf{AND}~y$ to both parties.}  that uses Beaver bit triples~\cite{beaver}.
We observe that the same secret value is used in 2 $\mathsf{AND}$ instances.
Hence, we construct correlated pairs of bit triples using 1-out-of-8 OT protocols~\cite{kkot} to reduce this cost to $\lambda+8$ bits (amortized) per triple, where $\lambda$ is the security parameter and typically $128$.
Some more work is needed for the above technique to work efficiently for the general case when $m$ does not divide $\ell$ or $\ell/m$ is not a power of $2$.
Finally, by picking $m$ appropriately, we obtain a protocol for millionaires' whose concrete communication (in bits) is nearly $5$ times better than prior work.
\\

\noindent\textbf{$\mathsf{DReLU}$.}
Let $a$ be additively secret shared as $a_0, a_1$ over the appropriate ring.
$\mathsf{DReLU}(a)$ is $1$ if $a \geq 0$ and $0$ otherwise; note that $a\geq0$ is defined differently for $\ell-$bit integers and general rings.
Over $\mathbb{Z}_L$, where values are encoded using 2's complement notation, $\mathsf{DReLU}(a) = 1\oplus \mathsf{MSB}(a)$, where $\mathsf{MSB}(a)$ is the most significant bit of $a$.
Moreover,  $\mathsf{MSB}(a) = \mathsf{MSB}(a_0) \oplus \mathsf{MSB}(a_1) \oplus \mathsf{carry}$.
Here, $\mathsf{carry} = 1$ if $a_0'+a_1' \geq 2^{\ell-1}$, where $a'_0, a'_1$ denotes the integer represented by the lower $\ell-1$ bits of $a_0, a_1$.
We compute this $\mathsf{carry}$ bit using a call to our millionaires' protocol.
Over $\mathbb{Z}_n$, $\mathsf{DReLU}(a) = 1$ if $a \in [0,\lceil n/2 \rceil)$. Given the secret shares $a_0, a_1$, this is equivalent to $(a_0+a_1) \in [0,\lceil n/2 \rceil) \cup [n, \lceil 3n/2 \rceil)$ over integers.
While this can be na\"{\i}vely computed by making 3 calls to the millionaires' protocol, we show that by carefully selecting the inputs to the millionaires' protocol, one can do this with only 2 calls.
Finally, we set things up so that the two calls to millionaires' have correlated inputs that reduces the overall cost to $\approx$ 1.5 instances of millionaires' over $\mathbb{Z}_n$.
\\

\noindent\textbf{Division and Truncation.} As a technical result, we  provide a {\em correct} decomposition of division of a secret ring element in $\mathbb{Z}_L$ or $\mathbb{Z}_n$ by a public integer into division of secret shares by the same public integer and  correction terms (\theoremref{general-division}).
These correction terms consist of multiple inequalities on secret values.
As a corollary, we also get a much simpler expression for the special case of {\em truncation}, i.e., dividing $\ell$-bit integers by a power-of-$2$ (\corolref{truncate-int}).
We believe that the general theorem as well as the corollary can be of independent interest.
Next, we give efficient protocols for both general division (used for $\mathsf{Avgpool}$, \tableref{comp-avgpool}) as well as division by a power-of-$2$ (used for multiplication in fixed-point arithmetic). The inequalities in the correction term are computed using our new protocol for millionaires' and the division of shares can be done locally by the respective parties.
Our technical theorem is the key to obtaining secure implementation of DNN inference tasks that are bitwise equivalent to cleartext fixed-point execution.

\subsection{Other Related Work}

Perhaps the first work to consider the secure computation of machine learning inference algorithms was that of~\cite{bost}.
SecureML~\cite{secureml} was the first to consider secure neural network inference and training. Apart from the works mentioned earlier, other works include those that considered malicious adversaries~\cite{helen,mlspdz,leviosa} (for simpler ML models like linear models, regression, and polynomials) as well as specialized DNNs with 1 or 2 bit weights~\cite{xonn,deepsecure,nitin}.
Recently, \cite{edabits} gave protocols for faithful truncation (but not division) over $\ell$-bit integers and prime fields in various adversarial settings. For 2-party semi-honest setting, our protocols have up to $20\times$ less communication for the truncations required in our evaluation.
\cite{RSS19} proposed an HE-based triple generation protocol over $\mathbb{Z}_{2^{\ell}}$, which requires less communication than generating triples using OT.

\subsection{Organisation}
We begin with the details on security and cryptographic primitives used in Section~\ref{sec:prelims} on preliminaries.
In Section~\ref{sec:comp-protocols} we provide our protocols for millionaires' (Section~\ref{sec:prot-mill}) and $\mathsf{DReLU}$ (Section~\ref{sec:prot-drelu-int},~\ref{sec:prot-drelu-ring}), over both $\mathbb{Z}_L$ and general ring $\mathbb{Z}_n$. In Section~\ref{sec:division}, we present our protocols for division and truncation. We describe the various components of DNN inference  in Section~\ref{sec:secureinference} and show how to construct secure protocols for all these components given our protocols from Sections~\ref{sec:comp-protocols} and ~\ref{sec:division}. We present our implementation details in Section~\ref{sec:impl} and our experiments in Section~\ref{sec:experiments}. Finally, we conclude and discuss future work in \sectionref{conclusion}.

 \section{Preliminaries}
\label{sec:prelims}
\paragraph{Notation.}
For a set $\mathcal{W}$, $w \getsr \mathcal{W}$ denotes sampling an element $w$, uniformly at random from $\mathcal{W}$. $[\ell]$ denotes the set of integers $\{0, \cdots, \ell-1\}$. Let $\ind{b}$ denote the indicator function that is $1$ when $b$ is {\em true} and $0$ when $b$ is {\em false}.

\subsection{Threat Model and Security}
We provide security in the simulation paradigm~\cite{gmw, canetti00, lindellsim} against a {\em static semi-honest} probabilistic polynomial time (PPT) adversary $\mathcal{A}$. That is, a computationally bounded adversary $\mathcal{A}$ corrupts either \ensuremath{P_0}\xspace or \ensuremath{P_1}\xspace at the beginning of the protocol and follows the protocol specification honestly.
Security is modeled by defining two interactions: a real interaction where \ensuremath{P_0}\xspace and \ensuremath{P_1}\xspace execute the protocol in the presence of $\mathcal{A}$ and the environment $\mathcal{Z}$ and an ideal interaction where the parties send their inputs to a trusted functionality that performs the computation faithfully.
Security requires that for every adversary $\mathcal{A}$ in the real interaction, there is an adversary
$\mathcal{S}$ (called the simulator) in the ideal interaction, such that no environment $\mathcal{Z}$ can distinguish between real and ideal interactions. Many of our protocols invoke multiple sub-protocols and we describe these using the \emph{hybrid model}. This is similar to a real interaction, except that sub-protocols are replaced by the invocations of  instances of  corresponding functionalities. A protocol invoking a functionality $\mathcal{F}$ is said to be in ``$\mathcal{F}$-hybrid model.''

\subsection{Cryptographic Primitives}

\subsubsection{Secret Sharing Schemes}
\label{sec:ss}

Throughout this work, we use 2-out-of-2 additive secret sharing schemes over different rings~\cite{shamir,blakeley}. The 3 specific rings that we consider are the field $\mathbb{Z}_2$, the ring $\mathbb{Z}_L$, where $L = 2^\ell$ ($\ell = 32$, typically), and the ring $\mathbb{Z}_n$,  for a positive integer $n$ (this last ring includes the special case of prime fields used in the works of \cite{gazelle,delphi}). We let $\genshare{L}{x}$ denote the algorithm that takes as input an element $x$ in $\mathbb{Z}_L$ and outputs shares over $\mathbb{Z}_L$, denoted by $\share{x}{L}{0}$ and $\share{x}{L}{1}$.
Shares are generated by sampling random ring elements $\share{x}{L}{0}$ and $\share{x}{L}{1}$, with the only constraint that $\share{x}{L}{0}+\share{x}{L}{1} = x$ (where $+$ denotes addition in $\mathbb{Z}_L$). Additive secret sharing schemes are perfectly hiding, i.e., given a share $\share{x}{L}{0}$ or $\share{x}{L}{1}$, the value $x$ is completely hidden. The reconstruction algorithm $\reconst{L}{\share{x}{L}{0}, \share{x}{L}{1}}$ takes as input the two shares and outputs $x = \share{x}{L}{0}+\share{x}{L}{1}$. Shares (along with their corresponding $\genshare{}{}$ and $\reconst{}{}$ algorithms) are defined in a similar manner for $\mathbb{Z}_2$ and $\mathbb{Z}_n$ with superscripts $B$ and $n$, respectively. We sometimes refer to shares over $\mathbb{Z}_L$ and $\mathbb{Z}_n$ as arithmetic shares and shares over $\mathbb{Z}_2$ as boolean shares.

\subsubsection{Oblivious Transfer}

Let $\kkot{k}{\ell}$ denote the 1-out-of-$k$ Oblivious Transfer (OT) functionality~\cite{bcr1outofnot} (which generalizes 1-out-of-2 OT~\cite{rabinot, eglot}). The sender's inputs to the functionality are the $k$ strings $m_0, \cdots, m_{k-1}$, each of length $\ell$ and the receiver's input is a value $i \in [k]$.
The receiver obtains $m_i$ from the functionality and the sender receives no output.
 We use the protocols from~\cite{kkot}, which are an optimized and generalized version of the OT extension framework proposed in~\cite{beaverotextension,iknp}.
This framework allows the sender and receiver,  to ``reduce''  $\lambda^c$ number of oblivious transfers to $\lambda$ ``base'' OTs. We also use the notion of correlated 1-out-of-2 OT~\cite{ALSZ13}, denoted by $\iknpcot{\ell}$. In our context, this is a functionality where the sender's input is a ring element $x$ and the receiver's input is a choice bit $b$. The sender receives a random ring element $r$ as output and the receiver obtains either $r$ or $x+r$ as output depending on $b$.
The protocols for $\kkot{k}{\ell}$ \cite{kkot} and $\iknpcot{\ell}$ \cite{ALSZ13} execute in $2$ rounds and have total communication\footnote{The protocol of $\kkot{k}{\ell}$ \cite{kkot} incurs a communication cost of $\lambda + k\ell$. However, to achieve the same level of security, their security parameter needs to be twice that of $\iknpcot{\ell}$. In concrete terms, therefore, we write the cost as $2\lambda+k\ell$.} of $2\lambda + k\ell$ and $\lambda + \ell$, respectively. Moreover, simpler $\kkot{2}{\ell}$ has a communication of $\lambda+2\ell$ bits \cite{iknp,ALSZ13}.

\subsubsection{Multiplexer and B2A conversion}
\label{sec:mux}

The functionality $\fmux{n}$ takes as input arithmetic shares of $a$ over $n$ and boolean shares of choice bit $c$ from $\party{0}, \party{1}$, and returns shares of $a$ if $c =1$, else returns shares of $0$ over the same ring.
A protocol for $\fmux{n}$ can easily be implemented by 2 simultaneous calls to $\kkot{2}{\eta}$ and communication complexity is $2(\lambda + 2\eta)$, where $\eta = \lceil \log n \rceil$.

The functionality $\fBtoA{n}$ (for boolean to arithmetic conversion) takes boolean (i.e., over $\mathbb{Z}_2$) shares as input and gives out arithmetic (i.e., over $\mathbb{Z}_{n}$) shares of the same value as output. It can be realized via one call to $\iknpcot{\eta}$ and hence, its communication is $\lambda+\eta$. For completeness, we provide the protocols realizing $\fmux{n}$ as well as $\fBtoA{n}$ formally in


\section{Introduction}\label{sec:int}

We consider a quantum dynamics, which is
a quantum counterpart of turbulent diffusion
in comb geometry.
A possible realization of the turbulent transport
can be discussed
in the framework of a Langevin equation in a
so-called Matheron - de Marsily form \cite{MaMa80},
\begin{equation}\label{MadeMa}
\dot{x_1}=\delta(x_2)vx_1\, \quad \dot{x_2}=\eta(t)\, .
\end{equation}
Here $\eta(t)$ is a random Gaussian delta correlated
process
$\langle\eta(t)\eta(t')\rangle=2D\delta(t-t')$,
where $D$ is a diffusion coefficient.
In this case, white noise in the $x_2$ axis
affects the velocity of inhomogeneous advection along
the $x_1$ axis at $x_2=0$.
Introducing a probability distribution function (PDF)
$P(x,y,t)= \langle\delta(x_1-x)\delta(x_2-y)\rangle$, we
obtain a Fokker-Planck equation (FPE) in the form of a 2D comb
model \cite{BaIo04}
\begin{equation}\label{comb-eq-1}
\partial_tP(x,y,t)=-v\delta(y)\partial_x xP(x,y,t)+
D\partial_y^2P(x,y,t)\, .
\end{equation}

The comb Eq. \eqref{comb-eq-1} has been considered in
Ref.~\cite{BaIo04} as a realistic mechanism of
superdiffusion of heavy particles, related to various
applications.
Important role in this turbulent acceleration
of the comb transport and superdiffusion
plays inhomogeneous advection in the form of a
dilatation operator. In the Hamiltonian form,
the dilatation (contraction) operator
\begin{equation}\label{int-01}
\hat{\clH}_0=-ix\hbar\frac{d}{dx}-i\hbar/2\equiv
x\hat{p}_x-i\hbar/2
\end{equation}
attracts much attention
in quantum mechanics and number theory in connection with
the Riemann zeros. This operator has been suggested by Berry
and Keating \cite{BeKe99a,BeKe99b,Connes} in connection to the Riemann
Hypothesis, and $\hclH_0$ is known as a Berry-Keating-Connes
Hamiltonian \cite{Sierra}. In a more general content, the operator has been
considered as a generator of space-time conformal
transformations \cite{AlFuFu76,Jackiw} and entropy problem
in quantum mechanics \cite{vanW}, a random walk approximation to the Riemann
Hypothesis \cite{Armitage}, and the quantum Mellin transform \cite{TwMi06}.
Also it corresponds to quantum dynamics in inverted
potentials \cite{BhKhLa95}, and near hyperbolic points \cite{NoVo97},
including singular behavior \cite{BeVi03} and exponential
spreading of the phase space \cite{Io2013}.
In comb geometry, this operator leads to the exponential
spread of the phase space as well \cite{Io2009}.
Recently, it has been considered as a turbulent diffusion search
optimisation \cite{SaIoKo2020}.

Due to the comb geometry, the backbone inhomogeneous advection
along the $x$ axis acts\footnote{Note that also due to this comb geometry,
inhomogeneous advection is replaced by a random
motion \cite{BaIo04,SaIoKo2019}, see Sec. \ref{sec:Vc}. }
only at the $y=0$.
This geometrical constraint of transport coefficients
leads (at least) to two possible scenarios of quantum mechanical
transport related to turbulent diffusion
described by Eq. \eqref{comb-eq-1}.
In particular,
the comb FPE \eqref{comb-eq-1} has straightforward
relation to quantum mechanics. It is just a one
dimensional quantum mechanical problem
at $y=0$ with the backbone Schr\"odinger equation
$i\hbar\partial_t\psi(x,t)=-iv\hbar x\partial_x \psi(x,t)$.
The classical PDF is mapped on the wave/Green's functions:
$P(x,y=0,t)\rightarrow \psi(x,t)/G_x(x,t)$.
In this case, one can consider
the quantum dynamics governed by the non-Hermitian Hamiltonian
according to the Schr\"odinger equation
\begin{equation}\label{in-sche-1}
i\hbar\partial_t\psi(x,y,t)=\hat{\clH}\psi(x,y,t)
=\left[v\delta(y)\hclH_0 -\tfrac{iD}{\hbar}\hat{p}_y^2\right]\psi(x,y,t)\, .
\end{equation}
In this case, the quantum (non-unitary) dynamics inside the backbone
along the $x$ axis at $y=0$
is accompanied by relaxation due to diffusion inside fingers
along the $y$ axis.

The second scenario is according to the Wick rotation of time
$t\rightarrow it$ with $D=\hbar/2m$, where $m=1$ is a unit mass of
a quantum particle. Then the non-Hermitian Hamiltonian is
\begin{equation}\label{int-02}
\hclH=\hat{p}_y^2/2+i\hclH_0\delta(y)\equiv \hat{A}+i\hat{B}\, .
\end{equation}
In this case, the backbone is
a quantum trap for the $y$ axis transport, where relaxation is due to
the dilatation operator $\hclH_0$, which works now as a
scattering.
This approach is a generalization of the dynamics of quantum particles
in the presence of traps due to a complex optical potential
\cite{KrLuMa14,DeKhRo2003}.

\section{Turbulent diffusion}\label{sec:Vc}

In this section, we consider the comb Eq. \eqref{comb-eq-1}
as a realistic mechanism of exponential (or turbulent)
superdiffusion and extend a picture of possible
solutions of the problem, described in Ref.~\cite{BaIo04}.
For the completeness of the presentation
we first consider the logarithmic scale of possible solutions.

\subsection{Logarithmic scale}\label{sec:lsv}

Let us consider dynamics for $x(x_1)>0$.
Following to  Ref.~\cite{TwMi06}, we introduce new variable
$w=\ln(x_1)$, and Eq. \eqref{MadeMa}  reads
\begin{equation}\label{lsv-1}
\dot{w}=v\delta(x_2), \quad \dot{x_2}=\eta(t)\, .
\end{equation}
The comb Fokker-Planck equation \eqref{comb-eq-1} in the
log scale for the PDF $P=P(w,y,t)$ now reads
\begin{equation}\label{lsv2}
\partial_tP=-v\delta(y)\partial_wP+D\partial_y^2P
\end{equation}
with the initial condition\footnote{To relate this initial
condition to the initial condition of
Eq. \eqref{comb-eq-1} in the form of Dirac
$\delta$ function, we introduce the multiplier
$1/x_0$ to keep $P_0(x,y)=P_0(w(x),y)=
\frac{1}{x}\delta(\log(x)-\log(x_0))=
\delta(x-x_0)$.}
$P_0(w,y)=\frac{1}{x}\delta(w-w_0)\delta(y)=
e^{-w}\delta(w-w_0)\delta(y)$ and $w\in(-\infty\, ,\infty)$.
After the Laplace transform with respect to $t$, we  look for the
solution in the form
\begin{equation}\label{td-022}
P(w,y,s)=\clL\left[P(w,y,t)\right](s)=f(w,s)g(y,s)=
f(w,s)\exp\left(-|y|\sqrt{s/D}\right)\, .
\end{equation}
This yields
\begin{equation}\label{lsv-3a}
v\partial_wf(w,s)+2\sqrt{Ds}f(w,s)=\frac{1}{x}\delta(w-w_0)\, .
\end{equation}
We obtain the PDF of the backbone diffusion as follows
\begin{multline}\label{lsv-3}
f(w,t)\equiv P(w,y=0,t)= \\
=\frac{(w-w_0)D}{\sqrt{\pi t^3v^2}}
\exp\left(-\frac{D(w-w_0)^2}{tv^2}\right)e^{-w}\theta(w-w_0)\, .
\end{multline}
In the $x$ space, it corresponds to the
log-normal distribution\footnote{In Eq. \eqref{lsv-4}, we
replace $\theta(w-w_0)$ by $\theta(x-x_0)$ using
a chain of transformations as follows $\frac{d}{dx}\theta(x-x_0)=\delta(x-x_0)=
\frac{1}{x}\delta(w-w_0)=\frac{1}{x}\frac{d}{dw}\theta(w-w_0)
=\frac{d}{dx}\theta(w-w_0)$.}
\begin{equation}\label{lsv-4}
P(x,y=0,t)=\frac{D\ln(x/x_0)}{x\sqrt{\pi t^3v^2}}
\times\exp\left(-\frac{D\ln^2(x/x_0)}{tv^2}\right)\theta(x-x_0)\, .
\end{equation}
Solution in Eq. \eqref{lsv-4} can be extended on $x<0$.
Introducing symmetric initial condition
$\frac{1}{2}[\delta(x-a)+\delta(x+a)]$ at the backbone,
one obtains
\begin{equation}\label{lsv-4a}
P(x,y=0,t)=\frac{\ln(|x|/a)}{2|x|\sqrt{\pi t^3v^2/D}}
\times\exp\left(-\frac{\ln^2(|x|/a)}{tv^2/D}\right)\theta(|x|-a)\, .
\end{equation}

\subsection{Symmetric solution}\label{sec:Ss}

Let us consider the comb equation \eqref{comb-eq-1}
with a symmetrical initial condition at $(x,y)=(\pm a,0)$:
$P_0(x,y)=\delta(y)[\delta(x-a)+\delta(x+a)]/2$, where $a>0$.
Taking into account the Laplace transformation with respect to $t$
and ansatz $P(x,y,s)=f(x,s)\exp\left(-|y|\sqrt{s/D}\right)$,
we obtain
\begin{equation}\label{ss-1}
vx\partial_xf+(2\sqrt{Ds}+v)f=\delta(|x|-a)/2\, ,
\end{equation}
where $x\in (-\infty\, ,\infty)$
and $\delta(|x|-a)\equiv\delta(x-a)+\delta(x+a)$.
The solution of Eq. \eqref{ss-1} is
$f(x,s)=C|x/a|^{-\mu(s)}\theta(|x|-a)$,
where $\mu(s)=2\sqrt{Ds}/v+1$. We also take into account that
$\partial_x|x|=\rm{sgn}(x)=2\theta(x)-1$, and
$\rm{sgn(x)}|x|=x$. An arbitrary constant $C$ is
determined from the normalization condition of $P(x,y,s)$.
Therefore, we obtain a ``dead zone'' in $x\in (-a,a)$ with
$f(x,t)=0$ and log normal distribution according
Eq. \eqref{lsv-4a} for $|x|\ge a$. The ``dead zone'' is due
to the initial condition at $x=\pm a\neq 0$ and
directional outward streaming
from $x=\pm a$ to $x\rightarrow \pm\infty$.
For $a=0$, the solution of Eq. \eqref{ss-1}
reads\footnote{Here we use that $\partial_xx\delta(x)=
\delta(x)+x\delta'(x)=\delta(x)+x(-\delta(x)/x)=0$.}
\cite{BaIo04}
$f(x,s)=C|x|^{1-\mu(s)}+\delta(x)/(\mu(s)-1)$.
In this case, the PDF consists of two parts: the log-normal
distribution, which contributes to the exponentially fast
spreading and the second part is due to the pining delta
function. The latter does not contribute to the mean
squared displacement (MSD).

Correspondingly, the solution for the PDF $P(x,y,t)$ is a
convolution integral
\begin{equation}\label{ss-2}
P(x,y,t)=P(x,y=0,t)\ast P_y(y,t)
=\int_0^tP(x,y=0,t-t')P_y(y,t')dt'\, ,
\end{equation}
where $P(x,y=0,t)=f(x,t)$ and $P_y(y,t)$ is
the L\'evy-Smirnov density defined by the inverse Laplace
transform,
\begin{equation}\label{ss-3}
P_y(y,t)=\clL^{-1}\left\{\exp[-|y|\sqrt{s/D}]\right\}\, .
\end{equation}

\section{Quantum dynamics: Eigenvalue's expansion}\label{sec:Ee}

We consider two scenarios of the quantum dynamics
described by non-Hermitian operators.
The first one is the quantum dynamics along the backbone
with relaxation in the form of classical diffusion along the fingers.
Contrary to that, the second scenario is the quantum dynamics in fingers
while relaxation in the backbone is due to the dilatation operator.
In both cases, the dilatation operator $H_0$ plays
important role for the quantum tasks.
Therefore we start from the eigenvalue problem of the dilatation operator.

The dilatation operator \eqref{int-01}
$$\hat{\clH}_0=\hbar [-ix\partial_x-i/2]=\hat{p}_xx+i\hbar /2$$
determines the complete set of eigenfunctions $\chi_{\omega}(x)$ with
the eigenvalues $\omega$ according to the eigenvalue problem
$\hat{\clH}\chi_{\omega}(x)=\hbar\omega\chi_{\omega}(x)$,
where $\omega$ is the continuous spectrum and the
eigenfunctions are \cite{BeKe99a,Sierra}
\begin{equation}\label{Ee1}
\chi_{\omega}(x)=\frac{1}{\sqrt{N|x|}}
\exp[i\omega\ln|x|] \,,
\end{equation}
which satisfies the boundary conditions
$\chi_{\omega}(x=\pm\infty)=0$ and $ N=4\pi$.
For the continuous spectrum, the normalization
condition is
\begin{equation}\label{Ee2}
\int_{-\infty}^{\infty}\chi_{\omega'}^*(x)\chi_{\omega}(x)dx
=\delta(\omega-\omega')\, ,
\end{equation}
while the completeness relation is
$\int \chi_{\omega}^*(x')\chi_{\omega}(x)d\omega=\delta(x-x')$
(see {\em e.g.} \cite{LL}). Note also that for $x>0$ the
normalization constant is half as larger, $N=2\pi$
\cite{Sierra}.

It is worth mentioning that a mathematically rigorous calculation of the
normalization constant for the wave function $\chi_{\omega}(x)$
can be presented by following\footnote{We present these arguments
from Ref. \cite{Io2013} with necessary corrections.} the presentation due to
the monograph by V.A. Fock \cite{fok}. Since the operator $ x\hat{p}$ has
continuous spectrum  $\omega$, the eigenfunctions
$\chi(\omega,x)\equiv\chi_{\omega}(x)$ are not square
integrable. Therefore, the normalization condition exists not for
the eigenfunction but for the ``eigendifferential'' \cite{fok}
$\Delta\chi(\omega,x)$, which reads
$ \Delta\chi(\omega,x)=\int_{\omega}^{\omega+\Delta\omega}
\chi(\omega',x)d\omega'$.
Substituting here Eq. \eqref{Ee1}, one obtains
$$ \Delta\chi(\omega,x)=\frac{2}{\sqrt{N|x|}\ln|x|}
\exp[i(\omega+\Delta\omega)\ln|x|]\sin\frac{\Delta\omega\ln|x|}{2}\, .$$
This solution is already square integrable and has the
normalization form
\begin{equation}\label{Ee3}
\lim_{\Delta\omega\to 0}\frac{1}{\Delta\omega}
\int_{-\infty}^{\infty}dx\left|\Delta\chi(\omega,x)\right|^2=1\, .
\end{equation}
To take the limit in Eq. \eqref{Ee3}, the integrand can be presented as follows
$$\left|\Delta\chi(\omega,x)\right|^2=\frac{4}{N|x|\ln^2|x|}
\sin\frac{\Delta\omega\ln|x|}{2}\sin\frac{\Delta\Omega\ln|x|}{2}\, ,$$
where $\Delta\Omega=\Delta\omega+\Delta\omega_1$.
We can perform this trick, since
the ``eigendifferentials'' for not overlapping spectral regions $\Delta\omega$ and $\Delta\omega_1$ are orthogonal. \cite{fok}.
Then taking the limit
$\Delta\omega=0$ and carrying out the variable change
$z=(\Delta\omega_1/2)\ln |x|$ and taking into account that
$\frac{4}{N}\int_{-\infty}^{\infty}\frac{\sin z}{z}dz=\frac{4\pi}{N}$,
we obtain $N=4\pi$, which coincides exactly with the dimensionless
normalization constant $N$ in Eq. \eqref{Ee1}.

\section{Quantum dynamics in backbone}\label{sec:bQD}

The Schr\"odinger equation \eqref{in-sche-1} for the Green's
function, described the dynamics of the first scenario, is
\begin{equation}\label{bQD-1}
i\hbar\partial_tG(x,y,t)
=v\delta(y)(x\hat{p}_x-i\hbar/2)G(x,y,t)
-\tfrac{iD}{\hbar}\hat{p}_y^2G(x,y,t)
\end{equation}
with the initial condition $G(x,y,t=0)=\delta(x-x_0)\delta(y)$
and $D=\hbar/2$.
In the Laplace domain $\clL[\psi(t)](s)=\tilde{\psi}(s)$,
Eq. \eqref{bQD-1} reads
\begin{multline}\label{bQD-2}
is\hbar \tilde{G}(x,y,s)=
v\delta(y)(x\hat{p}_x-i\hbar/2)\tilde{G}(x,y,s) \\
-\tfrac{iD}{\hbar}\hat{p}_y^2\tilde{G}(x,y,s)
-i\hbar\delta(x-x_0)\delta(y)\, .
\end{multline}
The Laplace image of the Green's function
has the multiplicative form of both backbone
and finger's dynamics
\begin{equation}\label{bQD-3}
\tilde{G}(x,y,s)=\tilde{G}(x,y=0,s)g(y,s)=
f(x,s)e^{-|y|\sqrt{s/D}}\, ,
\end{equation}
where $g(y,s)=e^{-\sqrt{s/D}|y|}$ describes relaxation of
quantum dynamics in the form of classical
diffusion in fingers, while $f(x,s)$ describes
the quantum dynamics in the backbone.
Substituting solution \eqref{bQD-3} in Eq. \eqref{bQD-2},
we obtain equation for the backbone quantum dynamics
\begin{equation}\label{bQD-4}
-iv(x\partial_x+\tfrac{1}{2})f(x,s)-i\sqrt{2\hbar s}f+i\delta(x-x_0)=0.
\end{equation}

Formally, it corresponds to the PDF of classical diffusion,
described in Sec. \ref{sec:Ss}, with
$\mu(s)=1+\sqrt{2\hbar s}/v$. When $x_0=0$, the solution reads
\begin{equation}\label{bQD-5}
f(x,s)=|x|^{\tfrac{1}{2}-\mu(s)} +\delta(x)/[\mu(s)-1/2]\, .
\end{equation}

Further quantum mechanical analysis is carried out
in the framework of the spectral decomposition according to
the eigenvalue problem of Sec. \ref{sec:Ee}.
Since $\chi_{\omega}(x)$ is a complete set
of eigenfunction, $f(x,s)$ can be considered as a superposition
\begin{equation}\label{bQD-6}
f(x,s)=\int d\omega\chi_{\omega}(x)b_{\omega}(s)\, .
\end{equation}
From Eq. \eqref{bQD-4}, this yields equation for
$b_{\omega}(s)\equiv b(s) $ as follows
\begin{equation}\label{bQD-7}
v\omega b(s)-i\sqrt{2\hbar s}b(s) =-i\chi^*_{\omega}(x_0)\, .
\end{equation}
Resolving this equation with respect to $b(s)$ and substituting it in Eq.
\eqref{bQD-6}, we obtain
\begin{multline}\label{bQD-8}
f(x,s)=\frac{1}{iv}\int\frac{ d\omega\chi_{\omega}(x)\chi^*_{\omega}(x_0)}
{\omega-i\sqrt{2\hbar s}/v} =\\
=\frac{1}{4\pi iv|x|}\int\frac{d\omega e^{i\omega\ln|x/x_0|}}
{\omega-i\sqrt{2\hbar s}/v}= \\
=\frac{1}{2v|x|}e^{-\sqrt{2\hbar s}\ln|x/x_0|/v} \, ,
\end{multline}
where $x>x_0$ for $\rm{Re}\sqrt{s} >0$ and $x<x_0$ if $\rm{Re}\sqrt{s} <0$.
We note that there is no anymore the ``dead zone'', and
this quantum situation changes drastically from classical
turbulent diffusion discussed in Sec. \ref{sec:Vc}.
We call this effect a
\textit{quantum swimming upstream}\footnote{Note that
here we have only an indication on this phenomenon.
A more accurate and elegant way to obtain this result
will be carried out in Sec. \ref{sec:FF}.}.

Substituting the result of Eq. \eqref{bQD-8} in Eq. \eqref{bQD-3},
we obtain the Green's function in the form of the Laplace inversion
\begin{multline}\label{bQD-9}
G(x,y,t)=\frac{1}{2v|x|}
\clL^{-1}\left[e^{-(v|y|+\ln|x/x_0|)\sqrt{2s\hbar}/v}\right]= \\
=\frac{(v|y|+\ln|x/x_0|)}{4v^2|x|\sqrt{2\pi t^3/\hbar}}
\exp\left[-\hbar\frac{(v|y|+\ln|x/x_0|)^2}{2v^2t}\right]\, ,
\end{multline}
which is the quantum Green's function in the form of the
combination of the log-normal and L\'evy-Smirnov
densities.

\subsection{Fractional dynamics in backbone}\label{sec:Dis1}

Let us show that behind this classical
L\'evy-Smirnov form of the Green's function
there is a fractional quantum dynamics along the backbone.
To this end, let us perform the Laplace inversion
before performing integration with respect to the
spectrum in Eq. \eqref{bQD-8}.
Then the expansion coefficient
$b_{\omega}(s)$ in Eq. \eqref{bQD-7} reads
\begin{equation}\label{Dis1-1}
 b_{\omega}(s) = \frac{A_{\omega}(x_0)}{\sqrt{s}+i\tilde{\omega}}\, ,
\end{equation}
where $\tilde{\omega}=v\omega/\sqrt{2\hbar}$ and
$A_{\omega}(x_0)=\chi^*_{\omega}(x_0)/\sqrt{2\hbar}$.
Performing the Laplace inversion, we obtain the
temporal behavior of the coefficients in the form of the
Mittag-Leffler functions $E_{\alpha,\beta}(z)$ \cite{BaEr55}:
\begin{equation}\label{Dis1-2}
B_{\omega}(t)=\clL^{-1}[b_{\omega}(s)](t)=A_{\omega}(x_0)
t^{-\frac{1}{2}}
E_{\frac{1}{2},\frac{1}{2}}\left(-i\tilde{\omega}t^{\frac{1}{2}}\right)\, .
\end{equation}
Therefore, the backbone component of the Green's function reads
\begin{multline}\label{Dis1-3}
f(x,t)=\int B_{\omega}(t)\chi_{\omega}(x)d\omega = \\
=\frac{1}{\sqrt{2\hbar t}}\int d\omega \chi^*_{\omega}(x_0)
E_{\frac{1}{2},\frac{1}{2}}\left(-i\tilde{\omega}t^{\frac{1}{2}}\right)
\chi_{\omega}(x) = \\
=\frac{1}{\sqrt{2\hbar t}}
E_{\frac{1}{2},\frac{1}{2}}\left(-i\hat{H}t^{\frac{1}{2}}\right)
\delta(x-x_0)\equiv G_{x}\left(i\hat{\clH}t^{\frac{1}{2}}\right)\, ,
\end{multline}
where the ``Hamiltonian''
is $\hat{\clH}= \tfrac{v\hclH_0}{\sqrt{2\hbar}}$.

The Green's function is the convolution of two non-unitary
dynamics, which are the quantum fractional dynamics along the backbone
with the Green's function $G_x(x,t)=f(x,t)$
and the finger's dynamics according to the L\'evy-Smirnov density
\begin{equation}\label{Levy-Smirnov}
G_y(y,t)=
\clL^{-1}\left[\exp\left(-|y|\sqrt{2s\hbar}\right)\right](t)
=\frac{\hbar|y|}{\sqrt{\pi t^3}}
\exp\left(-\hbar\frac{y^2}{2t}\right)\, .
\end{equation}
That is
\begin{equation}\label{Dis1-4}
G(x,y,t)=G_x\left(\hat{\clH}t^{\frac{1}{2}}\right)\ast G_y(y,t),
\end{equation}
which yields the result in Eq. \eqref{bQD-9}.
Therefore, the quantum dynamics at the backbone is not unitary.

\section{Fox H function for quantum swimming upstream}\label{sec:FF}

The Mittag-Leffler function can be presented in the form of the Fox
$H$ function  \cite{MaHa08} (see also Exampl. 2.7
in Ref. \cite{IMH2018}) defined in terms of the Mellin-Barnes integral,
\begin{equation}\label{FF-1}
E_{\alpha,\beta}(-z)=H^{1,1}_{2,1}(z)=\frac{1}{2\pi i}
\int_C\frac{\Gamma(\xi)\Gamma(1-\xi)}{\Gamma(\beta-\alpha\xi)}
(z)^{-\xi}d\xi\, ,
\end{equation}
where $\Gamma(\xi)$ is a gamma function and the contour $C$ starts
at $c-i\infty$, ends at $c+i\infty$ and separates the
poles of the gamma functions $\Gamma(\xi)$ and $\Gamma(1-\xi)$, such that
$0<\rm{Re}(\xi)<1$. In Eq. \eqref{FF-1},
$z=i\tilde{\omega}t^{\frac{1}{2}}$ is the imaginary variable
and $\alpha=\beta=1/2$.
Integration with respect to the spectrum reads
\begin{equation}\label{FF-2}
I(x/x_0)=\frac{1}{4\pi |x|}\int_{-\infty}^{\infty}
(i\omega)^{-\xi} e^{i\omega\ln|x/x_0|} d\omega\, .
\end{equation}
To perform this integration, the variable change
$i\omega=\tau$ yields the standard Laplace inversion
\begin{equation}\label{FF-3}
I(x/x_0)=\frac{1}{2|x|}\cdot\frac{1}{2\pi i}
\int_{-i\infty}^{i\infty}\tau^{-\xi}e^{w\tau}d\tau
=\frac{w^{\xi-1}}{2|x|\Gamma(\xi)}\, ,
\end{equation}
where $w\equiv w(x/x_0)=\ln|x/x_0|>0$ and $0<\rm{Re}(\xi)<1$.
When $w<0$, it means that a quantum particle is in
the dead zone $|x|<|x_0|$, then the variable change is
$\tau=-i\omega$. The integral in Eq. \eqref{FF-3} reads
\begin{equation}\label{FF-4}
I(x/x_0)=\frac{1}{2|x|}\frac{1}{2\pi i}\int_{-i\infty}^{i\infty}
(-\tau)^{-\xi} e^{\tau|w|} d\tau
=\frac{e^{-i\pi\xi}|w|^{\xi-1}}{2|x|\Gamma(\xi)}
=-\frac{w^{\xi-1}}{2|x|\Gamma(\xi)}\, .
\end{equation}
Therefore, taking into account results in Eqs.
\eqref{FF-1} - \eqref{FF-4} for the integration \eqref{Dis1-3}.
we obtain the backbone Green's function as follows
\begin{multline}\label{FF-5}
G_x(x,t)\equiv f(x,t)=
\int B_{\omega}(t)\chi_{\omega}(x)d\omega = \\
=\pm
\frac{w^{-1}}{2|x|\sqrt{2\hbar t}}
\frac{1}{2\pi i} \int_C \frac{\Gamma(1-\xi)}
{\Gamma(\frac{1}{2}-\frac{1}{2}\xi)}
\left(\frac{vt^{\frac{1}{2}}}{w\sqrt{2\hbar}}
\right)^{-\xi}d\xi \, ,
\end{multline}
where $(+)$ sign is for $w>0~(|x/x_0|>1)$ and
$(-)$ sign is for $w<0~(|x/x_0|<1)$.
Now, using the Legendre duplication formula,
$\Gamma(z)\Gamma(z+1/2)=2^{1-2z}\sqrt{\pi}\Gamma(2z)$,
which yields $\tfrac{\Gamma(1-\xi)}{\Gamma(1/2-\xi/2)}
=2^{-\xi}\Gamma(1-\xi/2)/\sqrt{\pi}$.
Defining $p=1-\xi/2$, we obtain Eq. \eqref{FF-5}
as follows (see \textit{e.g.}, Eq. (2.43) in Ref. \cite{IMH2018})
\begin{multline}\label{FF-6}
G_x(x,t)=
\pm \frac{w^{-1}}{2|x|\sqrt{2\pi\hbar t}} \\
\times \frac{\hbar w^2}{2v^2 t}
\frac{1}{2\pi i} \int_C\left(\frac{\hbar w^2}{2v^2t}\right)^{-p}
\Gamma(p)dp = \\
=\pm \frac{\ln\left|\tfrac{x}{x_0}\right|}{4|x|v^2\sqrt{2\pi t^3 /\hbar}}
\exp\left[-\hbar\frac{\ln^2\left|\tfrac{x}{x_0}\right|}{2v^2t}\right]\, .
\end{multline}
The Green's function in Eq. \eqref{FF-6} is a more accurate result on the quantum dynamics inside the ``dead zone''  for $ |x|<|x_0|$. It
contains the phase multiplier $e^{-i\pi}$, which corresponds to the
quantum swimming upstream.

\section{Quantum dynamics in fingers}\label{sec:QDF}

In this section we consider the second scenario,
which is performed according to the Wick rotation of
time $t\rightarrow it$ in Eq. \eqref{comb-eq-1}
with $D=\hbar/2m$, and $m=1$ is a unit mass of
a quantum particle. Then the Schr\"odinger equation
reads
\begin{equation}\label{QDF-1}
i\hbar\partial_t\psi(x,y,t)=\hclH\psi(x,y,t)\, ,
\end{equation}
where the non-Hermitian Hamiltonian is
$\hclH=\hat{p}_y^2/2+iv\hclH_0\delta(y)$.
An important point of the analysis should be
admitted: it is the imaginary optical potential,
which acts as a trap \cite{KrLuMa14,DeKhRo2003,AgMuBl10}.
In the finger's quantum dynamics the trap relates
to the dilatation operator and leads to the imaginary delta
function, see also \ref{sec:App-A}.
To understand this phenomenon, the dynamics of both the wave and Green's functions is analysed.

Again, we consider the wave function as a decomposition over
the complete set $\chi_{\omega}(x)$,
\begin{equation}\label{QDF-2}
\psi(x,y,t)=\int \phi_{\omega}(y,t)\chi_{\omega}(x)d\omega\, .
\end{equation}
Then taking into account properties of the dilatation/contraction
operator $\hclH_0$ and owing to Eq. \eqref{QDF-2}, we present the initial
condition as the spectral decomposition as well with the
Gaussian weight
$\rho_a(\omega)=\left[\frac{2a}{\pi}\right]^{\frac{1}{4}}\exp(-a\omega^2)$
with real $a>0$  \cite{Io2013}.
This yields the initial condition for the backbone
wave function in the form of a log-normal distribution.
Eventually, we obtain the initial condition for $\phi_{\omega}(y,t=0)$
as follows
\begin{equation}\label{QDF-2a}
\phi_{\omega}(y,t=0)=\rho_a(\omega)\Phi_0(y,\omega)\, ,
\end{equation}
where the choice of $\Phi_0(y,\omega)$ depends on the choice of
the boundary conditions for every $\omega$ mode.
In particular, free boundary conditions is considered in what follows.

\subsection{Wave function}\label{sec: FreeBC}

In the case of free conditions at the boundaries of the
fingers, we arrive at the $y$ transport in the form
of a scattering problem with the imaginary potential
$iv\hbar\omega\delta(y)$ as follows
\begin{equation}\label{QDF-3}
i\hbar\partial_t\phi_{\omega}=-\frac{\hbar^2}{2}\partial_y^2\phi_{\omega}
+iv\hbar\omega\delta(y)\phi_{\omega}\, .
\end{equation}
Looking for the eigenvalue solutions
$\phi_{\omega}(y,t)=e^{-iet/\hbar}\Phi_e(y)$,
we arrive at the eigenvalue problem
\begin{equation}\label{QDF-4}
-\frac{\hbar^2}{2}\partial_y^2\Phi_e
+i\lambda\delta(y)\Phi_e=e\Phi_e\, ,
\end{equation}
with free boundary conditions at infinities
and $\lambda=\hbar v\omega$.
Contrary to the real (attractive) delta potential \cite{Bl88},
there are no bound states.
However, still we can consider two classes of the
spectrum \cite{IMH2018,Bl88}.
The first class corresponds to the scattering problem and contains the
continuous spectrum solutions, which are standing waves. The solutions
can be considered as the result of
the scattering by the $i\delta$ potential of the left and right
incident waves \cite{Bl88},
\begin{equation}\label{QDF-5}
\Phi^{\pm}_k(y)=\frac{1}{\sqrt{2\pi}}\left[e^{\pm iky}+
\mathcal{B}(k)e^{ik\lvert y\rvert}\right],
\end{equation}
where $\mathcal{B}(k)={\lambda}/{(\hbar^2k-\lambda)}$ is the
``scattering amplitude'', determined by Eq.\ \eqref{QDF-4}
and related to imaginary potential\footnote{In contrast to
the scattering amplitude of a real potential, here the scattering amplitude
does not correspond to the optical theorem,
see \ref{sec:App-A}.},  and $e=e(k)=\hbar^2k^2/2$.

In further analysis we are interesting in the second
class of the solutions, which is but one (not bound) state
\begin{equation}\label{QDF-6}
\Phi_0(y)=(-i\lambda)^{1/2}e^{i\lambda|y|}\, , \quad e_0=\lambda^2/2\, .
\end{equation}
We choose $\Phi_0(y,\omega)=\Phi_0(y)$  to be the initial condition
for the finger's quantum dynamics.
This immediately yields the
quantum evolution of the finger's wave function
\begin{equation}\label{QDF-7}
\phi_{\omega}(y,t)=(-i\hbar\omega v)^{\frac{1}{2}}
e^{-i\frac{\hbar}{2}(\omega v)^2t}
\cdot e^{i\hbar\omega v|y|}\,
\end{equation}
and the dynamics (including the initial condition)
of the comb wave function reads from Eqs. \eqref{QDF-2} and \eqref{QDF-2a},
and \eqref{QDF-7} as follows
\begin{equation}\label{QDF-8}
\psi(x,y,t)
=\frac{(-i\hbar v)^{\frac{1}{2}}}{\sqrt{4\pi |x|}}
\left[\frac{2a}{\pi}\right]^{\frac{1}{4}}
\times\int_{-\infty}^{\infty}
\frac{e^{-A(t)\omega^2}}{\omega^{\frac{1}{2}}}
e^{i\omega (\hbar v|y|+\ln|x|)}d\omega\, ,
\end{equation}
where $A(t)=a+i\hbar v^2t$. When $t=0$ it corresponds to the
initial condition of the comb wave function.
The exact solution of the wave function is expressed in the form
of the hypergeometric and Bessel functions.
However it is instructive to present the result in the form of
elementary functions\footnote{We consider
the table integral \cite{BaEr54} in Eq. \eqref{QDF-8} as follows
$\int z^{-\frac{1}{2}}e^{-Az^2}e^{iBz}dz=
(1-i)2^{-1}A^{\frac{1}{4}}\Gamma(\tfrac{1}{4}){}_1F_1[\frac{1}{4};
\frac{1}{2};-\frac{B^2}{4A}]+(i-1)2^{\frac{3}{2}}\pi B^{\frac{1}{2}}
e^{-\frac{B^2}{8A}} I_{\frac{1}{4}}\left(\frac{B^2}{8A}\right)$.
The stationary phase approximation for the integral yields
$\left[\frac{2\pi}{iB}\right]^{\frac{1}{2}}e^{-\frac{B^2}{2A}}$.}.
To this end, the integral in Eq. \eqref{QDF-8} is evaluated by the
stationary phase approximation, which yields
\begin{equation}\label{QDF-9}
\psi(x,y,t)
=\frac{i(\hbar v)^{\frac{1}{2}}}{\sqrt{|x|(\hbar v|y|+\ln|x|)}}
\cdot\left[\frac{a}{2\pi}\right]^{\frac{1}{4}}
\times\exp\left[-\frac{(\hbar v|y|+\ln|x|)^2}{(a+i\hbar v^2 t)}\right]\, .
\end{equation}

\subsection{Green's function}\label{sec:Dis2}

\begin{figure}[htbp]
\includegraphics[width=1.0\hsize]{./Count1.eps}
\caption{Contour presentation of the real part of
Green's function \eqref{Dis2-3}
as a function of $t$ and $x$ for $t\in[10^{-5}\, ,  2\cdot 10^{-4}]$ and
$x\in [10^{-1}\cdot x_0\, ,5\cdot x_0] $, where the initial condition $x_0=10^{-3}$.
The imaginary part exhibits the same $x-t$ structure.
All parameter and constants are taken in dimensionless form
such that $\sqrt{\hbar}=2v$. The contour
is plotted on $x-t$ grid with $\Delta x=10^{-5}$ and $\Delta t=10^{-7}$.
In the case the axes are measured the number of steps on the grid.
At the initial stage, for small
$x$ and $t$ its behavior is random oscillations. Then it reflexes a
complicated quasi-periodic structure of the oscillations.
The white color corresponds to the positive values of $\mathrm{Re}[G_x]$,
while the turquoise color corresponds to the negative values.
However, the turquoise color regions are more colorful as zoom shows in
Fig.~\ref{fig:fig2}. }
\label{fig:fig1}
\end{figure}

Let us analyse the Green's function $G(x,y,t)$,
in particular the Green's function
of the backbone $G_x(x,t)$ in spirit of Sec. \ref{sec:bQD}.
Rewriting Eq. \eqref{QDF-1} for the Laplace image of the
Green's function  $\clL[G(x,y,t)](s)=\tilde{G}(x,y,s)=
\tilde{G}_x(x,s)\tilde{G}_y(y,s)$, we have
\begin{equation}\label{Dis2-1}
i\hbar s \tilde{G}=-\frac{\hbar^2}{2}\partial_y^2\tilde{G}
+iv\hclH_0\delta(y)\tilde{G} +i\hbar\delta(x-x_0)\delta(y)\, .
\end{equation}
Now, the ansatz for the finger Green's
function is $\tilde{G}_y(y,s)=e^{i\sqrt{2si/\hbar}}$,
while the backbone Green's function is a superposition
\eqref{bQD-6}, namely
$\tilde{G}_x(x,s)=\int b_{\omega}(s)\chi_{\omega}(x)$.
Taking these expressions into account in Eq. \eqref{Dis2-1},
we obtain the expression for the Green's function $G_x(x,t)$
in the form of Eq. \eqref{bQD-8} as follows
\begin{multline}\label{Dis2-2}
\tilde{G}_x(x,s)=\frac{-1}{v}\int
\frac{d\omega \chi_{\omega}(x) \chi^*_{\omega}(x_0)}
{\omega-i(1-i)\sqrt{\hbar s}/v} =\\
=\frac{-1}{4\pi v|x|}\int\frac{d\omega e^{i\omega\ln|x/x_0|}}
{\omega-i(1-i)\sqrt{\hbar s}/v}= \\
=\frac{-i}{2v|x|}e^{-(1-i)\sqrt{\hbar s}\ln|x/x_0|/v} \, .
\end{multline}
Performing the Laplace inversion we obtain
\begin{equation}\label{Dis2-3}
G_x(x,t)=\frac{\sqrt{i\hbar}\ln\left|\frac{x}{x_0}\right|}
{2v|x|\sqrt{2\pi t^3}}
\exp\left[\frac{i\hbar\ln^2\left|\frac{x}{x_0}\right|}{2v^2t}\right]\, .
\end{equation}
The behaviour of this solution is extremely complicated: it locks
completely random at the initial time and space and then changes to
complicated quasi-periodic behavior. The contour plot of $G_x(x,t)$
versus $(x,t)$ is shown in Figs. \ref{fig:fig1} and \ref{fig:fig2},
where Fig. \ref{fig:fig2} is a zoom of part of the $(x,t)$ space.

\begin{figure}[htbp]
\includegraphics[width=1.0\hsize]{./ZoomCount.eps}
\caption{Zoom of a part of the region in  Fig.~\ref{fig:fig1},
where $x/x_0\in[1.5\, , 2]$ and $t\in [4\cdot 10^{-5}\, , 7\cdot 10^{-5}]$.
More details are shown with smaller cells of the grid:
$\Delta x=2\cdot 10^{-6}$ and $\Delta t=5\cdot 10^{-8}$}
\label{fig:fig2}
\end{figure}

\section{Conclusion}

Two examples of quantum dynamics is considered in
the form of quantum counterparts
of turbulent diffusion in comb geometry. The latter is described by
either Langevin equation in the Matheron - de Marsily
form of Eq. \eqref{MadeMa} or equivalently by the Fokker-Planck equation
(FPE) in the form of the 2D comb model \eqref{comb-eq-1}.
The solution of the latter form corresponds to the
log-normal-L\'evy-Smirnov density
of the random variable $x$ with the variance growing exponentially:
$\mathrm{var}(x)=\langle (x-\langle x\rangle)^2\rangle\sim e^t$,
It describes relative diffusion of a pair of particles,
namely the distance between them. Therefore one finds
that the relative diffusivity of any two particles grows
with the inter-particle distance; and due to this growth,
this process has a property of turbulent diffusion.

Having concerned with the quantum dynamics, we have constructed
two variants of continuous time quantum walks related
to non-Hermitian operators of the form
$\hat{\clH}=\hat{A}+i\hat{B}$, where $\hat{A}$ is responsible for
the unitary transformation, while $i\hat{B}$ is responsible
for quantum/classical relaxation.
The first model in Eq. \eqref{bQD-1} is obtained by multiplication
of the comb model by $\hbar\sqrt{-1}$. As a rule, such simple
multiplication of an equation by a constant does not lead to a
new solution. However, due to the forms of the classical advection
$x\partial_x$ and the dilatation operator $\hat{\clH}_0=
-i\hbar(x\partial_x+1/2)=\hat{A}$, one immediately arrives at
non-Hamiltonian quantum mechanics inside the backbone,
which leads to quantum effects. In particular,
it is quantum tunneling inside a dead zone, which is a
quantum swimming upstream, that is a quantum initial
wave packet moves against
the classical streaming. The fingers play a role of classical environment
with classical relaxation according to normal diffusion.

The second model corresponds to a standard quantum counterpart
of a classical FPE by the Wick rotation of time
$t\rightarrow it$. The unitary dynamics inside fingers is due to
free motion, $\hat{A}=\frac{\hbar^2}{2}\partial_y^2$. The
relaxation is due to the dilatation operator $\hat{B}=
\hat{\clH}_0\delta(y)$, which is an imaginary optic
scattering potential for the incident plain wave coming from the fingers.
It is also a trap, which violets the optical theorem.

In conclusion, we comment symmetrical-asymmetrical positioning of the initial conditions, which is a generic task related to the comb geometry for both quantum mechanics and classical diffusion.
To be specific, let us consider the initial condition for Eq. \eqref{comb-eq-1} in a finger with $y=y_0\neq 0$  with the initial PDF $P_0(x,y)=\delta(x-x_0)\delta(y-y_0)$.
Then, following Eq. \eqref{td-022} with $g(y,s)=e^{-|y|\sqrt{s/D}}$, the finger transport is started first
with the L\'evy-Smirnov density \eqref{ss-3},
$P_y(y-y_0,t)=\tfrac{|y-y_0|}{\sqrt{4D\pi t^3}}
e^{-\frac{(y-y_0)^2}{4Dt}}$,
which is also
a first arrival time distribution to the backbone \cite{SaIoKo2019}.
Then, in the backbone's transport equation,
the L\'evy-Smirnov distribution $P_y(y_0, t)$ plays a role of a
strength of a $\delta$ source - sink term at $x = x_0$.
This situation differs from both scenarios considered above
and has been discussed for anomalous diffusion and random
search in two and three dimensional combs \cite{IoSa_MDPI,LSRJIK2020}.