\section{Introduction}
Multi-view echocardiographic sequences delineation provides important insight for clinical diagnosis.
The knowledge pattern of cardiac structures and textures associated with deforming tissues can be observed in echocardiographic sequence while in single frames the information is always missing and incomplete \cite{huang2014contour}.
Echocardiographic sequence also permits the assessment of wall motion and identification of end-diastolic (ED) and end-systolic (ES) phases.
Cardiologists usually check multi-view echocardiographic sequences in clinical decision-making \cite{madani2018fast}.
The apical-2-chamber view (A2C), A3C, and A4C are the most commonly used views for the left ventricle (LV) functional assessment.
Most clinical indexes of the LV (e.g., area, volume, and ejection fraction) are basically measured in these standard apical views.
Segmentation of the LV is generally a prerequisite for such quantitative analysis \cite{lang2015recommendations}.
In clinical routine, quantitative analysis of the LV still involves careful review and massive manual interpretation by experts, which is a tedious and time-consuming task.
Thus, automatic methods are desired to facilitate this process.
However, multi-view echocardiographic sequences segmentation remains a challenging task as illustrated in Fig. \ref{challenge}.
First, the fuzzy border, huge noise, and abounding artifacts of echocardiographic images result in local missing and incomplete of the anatomical structures;
Second, multi-view heterogeneous data varies in the anatomical structure, and image properties differ widely across vendors and centers;
Third, in the sequence, artifacts and noise are much severer, and the motion of mitral valve, trabeculation, and papillary muscles also poses additional interference;
Finally, limited labeled data restricts the performance of supervised learning based methods.

\begin{figure}[t]
    \setlength{\abovecaptionskip}{5pt}
    \centering
    \includegraphics[width=0.67\textwidth]{figure/challenge.pdf}
    \caption{Top left: multi-view samples (A2C, A3C, and A4C). Top right: A4C samples across vendors and centers. Bottom row: echocardiographic sequence.}
    \label{challenge}
\end{figure}

The application scenario of existing methods is always limited and only suitable under a specific situation.
They mostly focus on specific view \cite{carneiro2012segmentation} or single frames (i.e., without considering the sequence) \cite{chen2016iterative} or one single vendor and center \cite{leclerc2019deep}.
As for sequence segmentation, existing methods try to leverage temporal information by using a deformable model combined with the optical flow \cite{pedrosa2017fast,zhang2019deep} or fine-tuning pretrained CNN dynamically with first frame's label till the last frame \cite{yu2017segmentation}.
The major downsides of these temporal methods are that they are computational cumbersome and not an end-to-end manner.
The limited labeled data and specific application scenario confine the performance of existing methods and lead to the suboptimal solution.
\par
To achieve a unified model for multi-view echocardiographic sequences segmentation, we propose a recurrent aggregation learning method (RAL). The workflow is depicted in Fig. \ref{workflow}.
Pyramid ConvBlocks joint hierarchical ConvLSTMs are utilized to capture multi-level and multi-scale spatial-temporal information, enabling RAL the ability to harness the knowledge across heterogeneous data (multi-view, multi-center, and multi-vendor).
We further introduce a double-branch aggregation mechanism for segmentation and classification to lessen gaps across multi-view data.
Different from existing methods, RAL fully exploits the long term spatial-temporal information in an end-to-end manner and does not depend on any deformable model or optical flow or pretrained segmentation models.
RAL can accommodate heterogeneous data, not only generate accuracy segmentation results but also achieve the classification of different views at the same time and gain prominent temporal stability.

\begin{figure}[t]
    \setlength{\abovecaptionskip}{5pt}
    \centering
    \includegraphics[width=0.63\textwidth]{figure/workflow.pdf}
    \caption{Workflow overview of our method.}
    \label{workflow}
\end{figure}

\section{Method}
RAL is built as an end-to-end framework and comprised of three key components:
the feature extraction module, the segmentation branch, and the classification branch (as depicted in Fig. \ref{workflow}).
The feature extraction module consists of pyramid dilated dense convolution blocks (ConvBlocks).
The segmentation branch contains hierarchical recurrent architecture of multiple ConvLSTMs \cite{xingjian2015convolutional}.
While the classification branch involves a series of aggregation downsample and fully connected layers.
\newline
\textbf{Multi-Level and Multi-Scale Features Extraction.}
We design pyramid ConvBlocks architecture in the feature extraction module, which includes 5 ConvBlocks to extract multi-level and multi-scale features.
Multi-level information provides the global geometric characteristic of the LV, while multi-scale information can help to strengthen thin and small regions, further refine the boundaries of the LV.
They contribute to lessening the gap across views, vendors, and centers, increasing robustness to images conditions and the anatomical structure variations.
One ConvBlock contains $L$ densely connected dilated convolution layers as shown in Fig. \ref{detail}, which can expand the receptive field and meanwhile preserve the resolution of feature maps.
While the transition layer changes channels and resolution of feature maps by convolution and pooling.
The feedforward information propagation from preceding $l$ layers to $(l+1)^{th}$ layer can be formulated as
\begin{flalign}
    y_{l} = D(C(y_{1},y_{2},...,y_{l-1}))
\end{flalign}
where $y_{l}$ are the output of the $l^{th}$ layer, $C(\cdot)$ refers to the concatenation of previous layers' outputs. $D(\cdot)$ is a composite function of three connected operations: batch normalization (BN), rectified linear unit (ReLU), and dilated convolution.
Five ConvBlocks generate multi-level and multi-scale features $f_{t}=\{f_{t,1},f_{t,2},f_{t,3},f_{t,4},f_{t,5}\}$ for frame $t$ in the sequence.
\par
Pyramid ConvBlocks endow RAL with the superior feature extraction ability and the LV region detection capacity in multi-level and multi-scale space, further contribute to capturing the global geometric characteristic of the LV and then establishing uniform semantic features.
Thus RAL can detect and extract the LV accurately and robustly from not only ED and ES frames but also other frames in the sequence where the boundary is not clear (disturbed by noise and other tissues, see sequence samples in Fig. \ref{challenge}).
\begin{figure}[t]
    \setlength{\abovecaptionskip}{5pt}
    \centering
    \scalebox{0.93}{
    \begin{minipage}[t]{0.12\textwidth}
    \centering
    \includegraphics[width=13mm]{figure/workflow_denseblock.pdf}
    \end{minipage}
    \begin{minipage}[t]{0.87\textwidth}
    \centering
    \includegraphics[width=100mm]{figure/workflow_lstm.pdf}
    \end{minipage}}
    \caption{Left: Dilated dense convolution block. Right: Hierarchical ConvLSTMs for spatial-temporal modeling. }
    \label{detail}
\end{figure}
\newline
\textbf{Recurrent Features Fusion for Spatial-Temporal Modeling.}
\label{sequence}
For sequence segmentation, capturing the LV characteristic over time is essential for temporal stability.
Recent studies based on LSTM have shown great ability to learn sequential information.
Inspired by \cite{chen2018multiview,yang2018multiview}, we conduct hierarchical ConvLSTMs to exploit long term spatial-temporal modeling as depicted in Fig. \ref{detail}.
We add recurrence in the temporal domain to generate prediction $S_{t}$ for frame $t$ in the sequence, which carries forward the LV information from previous frames to following frames and allows the matching between consecutive frames naturally.
Additionally, we also add recurrence in the spatial domain for multi-level and multi-scale features fusion, which helps to integrate multi-level and multi-scale features efficiently.
\par
The output $y_{t,k}$ of the $k^{th}$ ConvLSTM at frame $t$ depends on the following variables:
(1) $k^{th}$ level and scale feature $f_{t,k}$ from the feature extraction module;
(2) the output $y_{t,k-1}$ of preceding $(k-1)^{th}$ ConvLSTM at the same frame $t$;
(3) the output $y_{t-1,k}$ from the $k^{th}$ ConvLSTM of previous frame $t-1$;
(4) the hidden state representation $h_{t,k-1}$ from preceding $(k-1)^{th}$ ConvLSTM at the same frame $t$, which is the spatial hidden state;
(5) the hidden state representation $h_{t-1,k}$ from the $k^{th}$ ConvLSTM of previous frame $t-1$, which is the temporal hidden state. The information flow can be formulated as
\begin{gather}
x_{input} = [f_{t,k}\ |\ B(y_{t,k-1})\ |\ y_{t-1,k}] \\
h_{state} = [h_{t,k-1}\ |\ h_{t-1,k}] \\
y_{t,k} = ConvLSTM_{k}(x_{input}, h_{state})
\end{gather}
where $B(\cdot)$ is the bilinear upsampling operator.
At each time step, every ConvLSTM accepts hidden states and encoded spatial-temporal features from previous ConvLSTMs and frame, the corresponding extracted feature from the feature extraction module, it then outputs encoded spatial-temporal features to next ConvLSTM and frame.
Finally, predictions $S_{t}$ are generated by the last ConvLSTM at every frame.
\newline
\textbf{Double-Branch Aggregation Learning.}
To further lessen the gaps across multi-view and refine multi-view segmentation results, we introduce a double-branch aggregation mechanism for simultaneous segmentation and classification of multi-view echocardiographic sequences as depicted in Fig. \ref{workflow}.
Feature from the last ConvBlock is sent to the classification branch.
Next, it goes through successive convolution and pooling operators to deeply aggregate with multi-level and multi-scale spatial-temporal features from the segmentation branch.
Finally, the classification result is produced by fully connected layers.
\par
The segmentation branch generates multi-view segmentations while the classification branch discriminates the specific view.
They are mutually promoted by deep aggregation of multi-level and multi-scale spatial-temporal features.
The segmentation branch provides multi-level and multi-scale spatial-temporal information to guide the classification while the classification branch affords multi-view discriminative regularization to refine the segmentation results and further lessen the gaps across views.
This double-branch aggregation mechanism endows RAL outstanding ability to adapt complex variations of anatomical structure.
\par
Additionally, we propose an aggregation loss to dynamically facilitate the communication between the segmentation branch and the classification branch as illustrated in Fig. \ref{workflow}.
The aggregation loss comprises the segmentation loss and classification loss.
The segmentation loss is a combination of binary cross-entropy loss and dice loss.
While the classification loss is categorical cross-entropy loss.
Thus the aggregation loss function can be formulated as
\begin{gather}
    L_{segmentation} = -[G\cdot log(P) + (1-G)\cdot log(1-P)] + \frac{2\cdot G \cdot P}{G + P} \\
    L_{classification} = - \sum_{i=1}^{3}g_{i}\cdot log(p_{i}) \\
    L_{aggregation} = \lambda_{s} \cdot L_{segmentation} + \lambda_{c} \cdot L_{classification}
\end{gather}
where $G$ and $P$ denote ground truth and prediction of segmentation respectively, $g$ and $p$ refer to ground truth and prediction of classification separately, $i$ indicates the type of view.
Besides, $\lambda_{s}$ and $\lambda_{c}$ are the corresponding balance coefficients, both are chosen empirically during the training process.
\begin{table}[t]
    \setlength{\belowcaptionskip}{5pt}
    \centering
    \caption{Specifications of our dataset (left) and the CAMUS dataset (right). }
    \begin{minipage}[t]{0.63\textwidth}
    \centering
    \scalebox{0.73}{
    \begin{tabular}[c]{cccccccc}
        \toprule
        Vendor Machines & Patients & Sequences & Images & ~ & A2C & A3C & A4C \\
        \midrule
        Philips EPIQ 7C & 60 & 180 & 5400 & \textbf{Sequences} & 100 & 100 & 100\\
        GE VIVID E9 & 20 & 60 & 1800 & Training & \multicolumn{3}{c}{240} \\
        Philips IE33 & 20 & 60 & 1800 & Testing & \multicolumn{3}{c}{60} \\
        \midrule
        Total & 100 & 300 & \textbf{9000} & Total & \multicolumn{3}{c}{300} \\
        \bottomrule
    \end{tabular}}
    \end{minipage}
    \begin{minipage}[t]{0.36\textwidth}
    \centering
    \scalebox{0.73}{
    \begin{tabular}[c]{cccccc}
        \toprule
        CAMUS & ~ & A2C & A4C & Vendor Machine \\
        \midrule
        \textbf{Images} & ~ & 900 & 900 & \multirow{3}*{GE VIVID E95} \\
        Training & ~ & \multicolumn{2}{c}{1600} & \\
        Testing & ~ & \multicolumn{2}{c}{200} & \\
        \midrule
        Total & ~ & \multicolumn{2}{c}{\textbf{1800}} & \\
        \bottomrule
    \end{tabular}}
    \end{minipage}
    \label{data}
\end{table}
\section{Experiments}
\textbf{Datasets.}
To validate the efficiency of RAL, we built a large multi-view echocardiographic sequences dataset, which was acquired from three centers' various vendor machines (The Second People's Hospital of Shenzhen, The Third People's Hospital of Shenzhen, and Peking University First Hospital).
We further evaluate RAL on the public CAMUS dataset \cite{leclerc2019deep}.
Our dataset contains 300 sequences from 3 views and every sequence includes 30 frames.
All 9000 frames were labeled by two experts.
Fig. \ref{segs} presents A2C, A3C, and A4C sequences samples segmented by RAL and experts.
While the CAMUS dataset only contains manual labels at ED and ES frames, which was acquired from a single vendor and center.
Table. \ref{data} shows the specifications of two datasets.
\newline
\textbf{Evaluation Metrics.}
Accuracy, Dice, Mean Absolute Distance (MAD) and Hausdorff Distance (HD) are used to measure segmentation results.
We further evaluate the segmentation performance with the ED, ES volume, and ejection fraction on the CAMUS dataset.
We utilize the output of RAL to compute clinical indices according to standard guidelines \cite{lang2015recommendations}.
\newline
\textbf{Implementation Details.}
All images are resized to $256\times256$ for computational efficiency.
We employ Adam with the learning rate of 0.001 as the optimizer.
The dilated rates of 5 ConvBlocks are $1, 1, 2, 4, 8$ respectively, and every ConvBlock contains 6 layers.
Besides, a dynamical decay mechanism is utilized to reduce the learning rate by monitoring the change of Dice.
Ten-fold cross-validation was utilized to provide an unbiased estimation.
\begin{figure}[t]
    \setlength{\abovecaptionskip}{5pt}
    \centering
    \includegraphics[width=0.8\textwidth]{figure/segs.pdf}
    \caption{The LV contours of multi-view sequences segmented by our method (red) and experts (green). Ten frames are selected from every sequence to fit the layout view. (top row: A2C; middle row: A3C; bottom row: A4C)}
    \label{segs}
\end{figure}
\begin{table}[t]
    \setlength{\belowcaptionskip}{5pt}
    \centering
    \caption{Ablation results of our methods under different configurations.}
    \scalebox{0.83}{
    \begin{tabular}[c]{ccccccc}
        \toprule
        Configurations & Accuracy & Dice & HD(mm) & MAD(mm) & Classification\\
        \midrule
        full
        & 0.987$\pm$0.005 & 0.919$\pm$0.040 & 5.87$\pm$3.46 & 2.90$\pm$1.49 & 0.933 \\
        w/o classification
        & 0.971$\pm$0.008 & 0.910$\pm$0.049 & 5.99$\pm$3.69 & 3.10$\pm$1.66 & - \\
        w/o ConvBlock
        & 0.963$\pm$0.015 & 0.907$\pm$0.057 & 6.21$\pm$4.95 & 3.27$\pm$1.85 & 0.867 \\
        w/o temporal
        & 0.955$\pm$0.019 & 0.896$\pm$0.062 & 6.64$\pm$5.04 & 3.51$\pm$1.93 & 0.917 \\
        w/o spatial
        & 0.968$\pm$0.011 & 0.911$\pm$0.054 & 6.03$\pm$4.16 & 3.08$\pm$1.71 & 0.883 \\
        \bottomrule
    \end{tabular}}
    \label{ablation}
\end{table}
\newline
\textbf{Ablation Study.}
We evaluate our method under different configurations to corroborate the necessity of every component in RAL.
The classification branch, ConvBlock, spatial modeling and temporal modeling are removed respectively.
Table. \ref{ablation} shows the ablation results, we can see that full RAL achieves higher mean values of Accuracy and Dice, lower mean values of HD and MAD, and lower standard deviations of all metrics compared against other configurations.
RAL also achieves the best classification accuracy (0.933).
Every single component brings important improvement for the LV segmentation, especially when adding recurrence in the temporal domain.
\newline
\textbf{Comparison Study I: Geometrical.}
We compare RAL with U-Net, ACNN, and U-Net++ on our multi-view sequences dataset.
As shown in Table. \ref{compare}, RAL outperforms other methods on all metrics, achieving the highest mean values of Accuracy (0.987) and Dice (0.919), the lowest mean values of HD (5.87mm) and MAD (2.90mm), and significantly lower standard deviations of all metrics.
These strongly prove that RAL is able to accomplish the best region coverage, the highest contour accuracy, and the minimum distance error when processing multi-view echocardiographic sequences across multi-vendor and multi-center.
\begin{table}[t]
    \setlength{\belowcaptionskip}{5pt}
    \centering
    \caption{Geometrical comparison results on our multi-view echocardiographic sequences dataset.}
    \scalebox{0.8}{
    \begin{tabular}[c]{cccccc}
        \toprule
        Methods & Accuracy & Dice & HD(mm) & MAD(mm) \\
        \midrule
        RAL & 0.987$\pm$0.005 & 0.919$\pm$0.040 & 5.87$\pm$3.46 & 2.90$\pm$1.49 \\
        U-Net    & 0.942$\pm$0.030 & 0.883$\pm$0.068 & 8.94$\pm$6.87 & 3.72$\pm$1.87 \\
        ACNN     & 0.959$\pm$0.013 & 0.893$\pm$0.061 & 7.70$\pm$6.58 & 3.40$\pm$1.57 \\
        U-Net++  & 0.937$\pm$0.032 & 0.880$\pm$0.072 & 9.01$\pm$7.14 & 3.86$\pm$2.01 \\
        \bottomrule
    \end{tabular}}
    \label{compare}
\end{table}
\begin{figure}[t]
    \setlength{\abovecaptionskip}{2pt}
    \begin{minipage}[t]{0.49\textwidth}
    \centering
    \includegraphics[width=42mm]{figure/frame.pdf}
    \end{minipage}
    \begin{minipage}[t]{0.49\textwidth}
    \centering
    \includegraphics[width=43mm]{figure/bland.pdf}
    \end{minipage}
    \caption{Left: Mean of Accuracy, Dice, HD, and MAD at different frames of the cardiac cycle. Right: Bland–Altman analysis (EF$_a$ and EF$_m$: ejection fraction calculated from automatic segmentations and manual labels) }
    \label{bland_frames}
\end{figure}
\newline
\textbf{Comparison Study II: Clinical.}
We compare RAL with U-Net, ACNN, and U-Net++ on the CAMUS dataset to calculate clinical indices.
As shown in Table. \ref{camus},
RAL obtained high correlation scores (0.952 for EDV, 0.960 for ESV, and 0.839 for EF), reasonably small biases and standard deviations, and relatively low mae (8.8ml for EDV, 7.1ml for ESV, and 5.0\% for EF).
Fig. \ref{bland_frames} presents a more intuitional result by Bland–Altman plot.
94\% of the measurements locate in the $\pm$1.96 standard deviation in Bland–Altman plot.
These results reveal the clinical potential of RAL.
\newline
\textbf{Temporal Stability.}
We compute the mean of Accuracy, Dice, HD, and MAD at different frames of all echocardiographic sequences and then observe the volatility of each metric to assess the temporal stability.
As shown in Fig. \ref{bland_frames}, RAL achieves stable mean values of all four metrics in the cardiac cycle, only exists moderate fluctuating in the middle of the sequence.
This means spatial-temporal modeling of RAL is efficient.
RAL achieves not only superior segmentation accuracy but also a good coherence of consecutive frames in the sequence.
\newline
\textbf{Limitation.}
In Fig. \ref{bland_frames}, from ED to ES frames, we observe that Accuracy and Dice decay slightly while HD and MAD increase mildly, and all metrics keep relatively stable in the diastole but show feeblish recoverability.
The sequential process carries errors forward resulting in accumulation of temporal errors in the cardiac cycle.
Fortunately, the fluctuating rate is moderate and the worst results are still fairly good.
This limitation could be alleviated via Bi-direction LSTM.

\begin{table}[t]
    \setlength{\belowcaptionskip}{5pt}
    \centering
    \caption{Clinical comparison results on CAMUS dataset. (EDV: ED volume; ESV: ES volume; EF: ejection fraction; corr: Pearson correlation; mae: mean absolute error)}
    \scalebox{0.73}{
    \begin{tabular}[c]{cccccccccccc}
        \toprule
        \multirow{2}*{Methods} & \multicolumn{3}{c}{EDV} & ~ & \multicolumn{3}{c}{ESV} & ~ & \multicolumn{3}{c}{EF} \\
        ~ & corr & bias(ml) & mae(ml) & ~
        & corr & bias(ml) & mae(ml) & ~
        & corr & bias(\%) & mae(\%) \\
        \midrule
        RAL
        & 0.952 & -7.5$\pm$11.0 & 8.8 & ~
        & 0.960 & -3.8$\pm$9.2 & 7.1 & ~
        & 0.839 & -0.9$\pm$6.8 & 5.0 \\
        U-Net
        & 0.954 & -6.9$\pm$11.8 & 9.8 & ~
        & 0.964 & -3.7$\pm$9.0 & 6.8 & ~
        & 0.823 & -1.0$\pm$7.1 & 5.3 \\
        ACNN
        & 0.945 & -6.7$\pm$12.9 & 10.8 & ~
        & 0.947 & -4.0$\pm$10.8 & 8.3 & ~
        & 0.799 & -0.8$\pm$7.5 & 5.7 \\
        U-Net++
        & 0.946 & -11.4$\pm$12.9 & 13.2 & ~
        & 0.952 & -5.7$\pm$10.7 & 8.6 & ~
        & 0.789 & -1.8$\pm$7.7 & 5.6 \\
        \bottomrule
    \end{tabular}}
    \label{camus}
\end{table}
\section{Conclusion}
In this paper, we present a recurrent aggregation learning method to exploit long term spatial-temporal information for simultaneous segmentation and classification of multi-view echocardiographic sequences.
Multi-level and multi-scale features are recurrently aggregated on both spatial domain and temporal domain for effective spatial-temporal modeling.
A double-branch aggregation mechanism further brings multi-view discriminative regularization to refine the segmentation results.
Adequate experiments of geometrical and clinical evaluation demonstrate that RAL achieves not only superior segmentation and classification accuracy, prominent temporal stability, but also high correlations on clinical indices.
\newline
\\
\textbf{Acknowledgment.}
This work is funded by the Shenzhen Basic Research Program (JCYJ20170818164343304, JCYJ20180507182432303).


\section{\label{sec:introduction}Introduction}
\textit{Introduction} - The study of energy transport phenomena in heterogeneous systems, like synthethic or biological media, ranging from tissues to computation nodes, enlightens the role of quantum coherence \cite{engel2007evidence,lambert2013quantum} which is believed to enhance the rate of these processes \cite{mohseni2008environment,hoyer2010limits}. In such conditions, a superdiffusive dynamic settles in the transport or propagation processes.
This is the case of heat excitation transport in particular condensed matter systems such as one dimensional (1D) metal lattices, described by the Luttinger liquids theory \cite{bulchandani2019superdiffusive}. Studies performed in this context, have shown a violation in the dichotomy between ballistic and diffusive transport regimes. There are several cases in which superdiffusion occurs, for example in classical one-dimensional systems \cite{van2012exact}, when both disorder and non linear effects are present \cite{kim2018direct}, and  even in quantum systems experiencing Many-Body Localization phase transitions \cite{vasseur2016nonequilibrium}.\\
Quantum walks (QWs) have long been found to efficiently describe coherent energy transport \cite{mohseni2008environment, hoyer2010limits}.
By suitably introducing disorder in a QW, it is possible to modify its spreading behaviour. A paradigmatic example is given by Anderson localization \cite{anderson1958absence}, originally formulated for condensed matter systems, and extensively demonstrated in the case of all-optical systems \cite{segev2013anderson,schwartz2007transport,crespianderson2013}. The theoretical and experimental studies carried out on this topic have covered a wide range of phenomena, including "hyper-transport" of light \cite{levi2012hyper, schreiberecoherence2011}, and disordered QWs preserving the time dependence of spreading \cite{brun2003quantum}, and decoherence effects determining the transition from Quantum to Classical Random Walks (CRWs) \cite{schreiberecoherence2011, broomedecoherence2010}.
In these studies the role of disorder source was played by decoherence or various types of unitary evolutions.
The present work deals with the description of the transition between quantum and classical regimes, while the coherence of the evolution is preserved. This has been done experimentally by implementing proper configurations of disorder in a discrete Quantum Walk (QW).\\
\textit{Theoretical Model} - The theoretical description of a discrete QW is based on a walker, a coin and the operators acting on them \cite{venegasQW2012}\cite{kempeQW2003}. In a 1D QW, the walker is given by a quantum particle which is in a superposition of the position states, described by the set of states $\{\ket{i}_p\}$, each of them corresponding to a particular site of the line. The coin is represented by any internal binary degree of freedom of the particle, given by the basis kets $\{\ket{0}_c,\ket{1}_c\}$. The evolution of the quantum particle is controlled by the coin operator $\hat{C}$, acting on the coin state, and by the shift operator $\hat{S}$, which moves the walker according to the coin state. The evolution of the walker is then described by the repeated action of coin and shift operators on the general state of the particle $\ket{\psi}=\sum_{i,j}P_{i,j}\ket{i}_p\otimes\ket{j}_c$, $\ket{\psi'}=\hat{S}\cdot(\hat{I}\otimes \hat{C})\ket{\psi}$
where the coin operator is expressed as $\hat{C}=\frac{1}{\sqrt{2}}(\ket{0}_c\bra{0}_c+\ket{0}_c\bra{1}_c+\ket{1}_c\bra{0}_c-\ket{1}_c\bra{1}_c)$ while the shift operator reads
$\hat{S}=\sum_{i}\ket{i-1}_p\bra{i}_p\otimes\ket{0}_c\bra{0}_c+\ket{i+1}_p\bra{i}_p\otimes\ket{1}_c\bra{1}_c.$
These equations describe the evolution of a walker in a completely ordered QW, in which the coin operator is uniform both in space and time. Other structures of the coin operator, based on suitable maps of the QW phases, are possible \cite{ahlbrechtasymptotic2011}.\\
Here, by randomly shifting the phase difference imposed by the coin operator, we are able to simulate the superdiffusive behaviour of a quantum walker and to investigate the transition from the ballistic regime of an ordered QW to the diffusive one of a Classical Random Walk (CRW).\\
The operator describing the phase shift at a certain step $n$ of the evolution can be written as $ \hat{P}_{n}=\sum_{i}\ket{i}_p\bra{i}_p\otimes(\ket{0}_c\bra{0}_c+\ket{1}_c\bra{1}_ce^{i\Delta\phi^i(n)})$, where $\Delta\phi^i(n)$ corresponds to the phase difference imposed by the coin operator in the site $i$. We define the degree of disorder $p$ as the percentage of random phases that the walker experiences during the evolution. As a consequence, the case $p=0$ corresponds to a standard ordered QW while for $p=1$ the QW becomes completely disordered (case of a CRW). We call this kind of Quantum Walk $p$\emph{-diluted} QW, because of the dilution of a given degree of disorder $p$ during the QW evolution. With no further assumptions this model describes a completely space- and time-uncorrelated disorder.\\
We limit our analysis to the case in which the phases difference $\Delta\phi^i(n)$ can only be $0$ or $\pi$. We are interested to study how the quantum to classical transition depends on the parameter $p$. For a given value of $p$ a number of different phase maps can be realized, in principle, corresponding to different evolutions of the probability distribution. Thus, the values of the relevant parameters identifying a given value of $p$ are obtained by averaging over many different phase map configurations characterized by the same degree of disorder.\\
\begin{figure}[!h]
 	\centering
 	\includegraphics[scale=0.55]{var_vs_p_7_20_step_1000_tries_teo_normalized.eps}%
 	\caption{\emph{Simulation of the variance behaviour for 7 and 20 steps as a function of $p$. The average values of the variance have been calculated on a sample of 1000 different phase maps for each value of $p$. In the vertical axes $Var_1(p)$ for 7 and 20 steps have been normalized to their maximum values. Error bands correspond to the standard deviation of the variance distributions.}}
 	\label{fig:var_simulated_vs_p}
 \end{figure}
A useful quantity in the study of the single particle evolution is the variance of the position probability distribution, defined as:
\begin{equation}
 Var_{1}(n)=\sum_{i=-n}^ni^2\cdot P_i(n)-(\sum_{i=-n}^ni\cdot P_i(n))^2
\end{equation}
where $P_{i}(n)$ is the probability to find the walker on site $i$ at step $n$.
We simulated the behaviour of the variance for different values of $p$ up to 20 steps, in the range $p=0$ to $1$, with a progressive increment of $0.01$. For each value of $p$ we simulated the evolution of $1000$ different phase maps, then we computed the mean value of the variance. In order to show how good this model of disorder is to study the superdiffusion, we computed the standard deviation of the variance distribution for each value of $p$. Results of simulation for 7 and 20 steps are shown in Fig. \ref{fig:var_simulated_vs_p}. In both cases the mean value of variance decreases with increasing of $p$, since a constant value around $p\sim 0.40$ is approached for both 7 and 20 steps. This indicates that the output distribution rapidly merges in a classical one for $p>0.40$, while a genuine superdiffusive behaviour occurs between $0<p\le 0.20$. Focusing on this region we observe that, for a given value of $Var_1(p)$, the uncertainty on the corresponding value of $p$ decreases with increasing of $n$.\\
A further useful quantity is the similarity which is a measure of how two distributions G and G' are similar: $S(G,G')=\frac{(\sum_{i}\sqrt{G_iG'_i})^2}{\sum_{i}G_i\sum_{j}G_j}$. In the graph of Fig. \ref{fig:KL_vs_p} we show the behaviour of the similarity between the average distribution for a given value of disorder $p$ and the ordered or disordered distribution at varying of $p$ (respectively, solid and dashed lines). The similarity with the ordered (disordered) distribution decreases (increases) with $p$ since the higher the disorder the more destructive interference tend to reproduce the disordered distribution. In the inset, the crossing points between the two quantities identify the transition between quantum and classical regimes for a given number of step $n$. It is worth noting that the transition point decreases with $n$.
\begin{figure}[!h]
 	\centering
 	\includegraphics[scale=0.55]{Similarity_vs_P_modified.eps}

 	\caption{\emph{Computational simulations for similarity between ordered (solid line) and disordered (dashed line) distributions (respectively $G(p=0)\equiv G(0)$ and $G(p=1)\equiv G(1)$) and the average distribution for a given value $p$ as a function of $p$. Here $S(0,p)$ $(S(1,p))$ stands for $S(G(0),G(p))$ $(S(G(1),G(p)))$. In the inset a zoom on the intersection region is reported. The crossing point shifts towards lower values of $p$ when the step number increases. Different colors stand for different step numbers of the evolution.}}
 	\label{fig:KL_vs_p}
 \end{figure}
This continuous transition occurs without loosing completely the quantum features of the evolution and can be interpreted as an evidence of the fact that the QW modifies its behaviour because of the total amount of disorder experienced along its spreading.\\
\textit{Experimental Implementation} - A 1D Quantum Walk can be realized by a network of Beam Splitters (BSs), each of them representing a particular site of the line and acting both as $\hat{C}$ and $\hat{S}$ operators. Several implementations of optical QW have been realized, based on bulk optics schemes, bulk-fiber circuits and femtosecond laser written photonic circuits \cite{bulkQW2005, crespianderson2013, silberhornQW2016, flaminireview2018, sansoni2012two}.
The bulk optic setup used in the present work consists of two displaced multi-pass Sagnac Interferometers ($SIs$) linked by a common Beam Splitter ($BS_1$ in Fig. \ref{fig:setup}).

\begin{figure}[!h]
 	\centering
 	\includegraphics[scale=0.23]{setup.eps}%
 	\caption{\emph{Sketch of the setup used in the experiment. BS: beam splitter, BD: beam displacer, RP: rotating glass plates, MM: moving mirror, SMF: single mode fiber. Each side of the square $SI$ is $50$cm long. Blue and red beams circulate in opposite directions and impinge on the $BS_1$ in the same horizontal point but at different heights along the z direction, due to the effect of BD.}}
 	\label{fig:setup}
 \end{figure}

Beam displacement is obtained by translating one mirror of $SI_1$ along the $x$ direction. This particular configuration allows to create in the $x-y$ horizontal plane of Fig. \ref{fig:setup} a (quasi) infinite loop, which is equivalent to a chain of phase-stable and independently tunable Mach-Zehnder interferometers \cite{cuevasmarkovianity2019}. The entire BS network necessary for a QW is realized by exploiting also the $z$ direction. On this purpose, suitable Beam Displacers (BDs) are inserted along clockwise directions in $SI_1$ and counterclockwise ones in $SI_2$ allowing to increase the number of $x-y$ planes on which the photons can travel (see Supplementary Material (SM)). Phases are independently addressed in each unit mesh of the QW by using Rotating Glass Plates (RP). The output modes of each step can be extracted for measurement by a set of Moving Mirros (MMs), intersecting and extracting from the setup only the modes at the selected step $N_j$, without intersecting the steps $N<N_j$. The extracted radiation is then coupled on a single mode fiber and measured (for further details on the setup see \cite{cuevasmarkovianity2019} and \cite{geraldiQW2019}). A slight modification of the setup allows to inject two distinguishable or indistinguishable photons through the two different input ports of $BS_1$. After the extraction of the modes, a further BS ($BS_2$) allows to separate photons travelling along the same mode. Adding a further coupler intercepting the extracted modes, coincidences between all possible modes at a given step can be measured and the whole two photon probability distribution can be experimentally reconstructed (see SM for further details).\\
\textit{Experimental Results} - It is well known that in a superdiffusive process the spread of the walker position follows a power law, expressed in our case by $Var_1(n)\propto n^{\beta}$, with $n$ corresponding to the number of steps and $1<\beta<2$ \cite{havlin2002diffusion}. Here, the case $\beta=2$ ($\beta=1$) is typical of a QW (CRW).
We experimentally implemented five different values of $p$, namely $0, 0.05, 0.10, 0.20$ and $1$. For each value of $p>0$ we selected three different phase maps and experimentally reconstructed the output probability distribution for each step of the evolution. In Fig. \ref{fig:var_vs_step_1photon} the average variance of the walker is reported as a function of the step number $n$.

\begin{figure}[!h]
 	\centering
 	\includegraphics[scale=0.60]{var_1_sper_teo_teoreal_new.eps}%
 	\caption{\emph{Experimental values of the variance as a function of the number of steps for different values of $p$ ($0,0.05,0.1,0.2,1$). Dots correspond to experimental data. For a given value of $p$, dashed line shows the expected average behaviour (Exp) of variance for the three selected phase maps. They were computed by taking into account the actual optical parameters of the setup. Dotted lines show the ideal behaviour (I) of a perfectly symmetric $BS_1$ in the two limit cases, namely the ordered and disordered QW. Note that, up to the third step, the variance is not affected by the phase thus the red and blue dotted lines are coincident; the same effect occurs for all the dashed lines.}}
 	\label{fig:var_vs_step_1photon}
 \end{figure}
The experimental setup allows to start the QW only from the central position of the 1D-lattice and by using any of the two coin modes (input ports of $BS_1$). Accordingly, we studied the evolution for the input state $\ket{\psi}=\ket{0}_{p}\otimes\ket{0}_{c}$.
Experimental data are in excellent agreement with the theoretical prediction for the three chosen phase maps (dashed lines). The whole superdiffusive region, namely the one between the diffusive (red lines) and the ballistic (blue lines) behaviours, can be exploited by varying the values of $p$. A small amount of phase instability was present due to the large size of the apparatus. This effect has been included in the computation of errors bars. We also show in Fig. \ref{fig:var_vs_step_1photon} the expected behaviours simulated by considering ideal optical elements for the two limit cases of ordered and completely disordered QW (respectively blue and red dotted lines). They differ from the experimental data because of the non perfectly symmetric behaviour of $BS_1$ (see SM for details) \cite{geraldiQW2019}.\\
The actual superdiffusive behaviour of the experimental evolutions is shown in Tab. \ref{tab:fit_results}, where we report the values of $\beta$ obtained by fitting the experimental data ($\beta_{Fit}$) and comparing them with the values obtained through theoretical simulations ($\beta_{Theo}$). The observed discrepancies can be explained as follows:

\begin{itemize}
    \item The $\beta_{Theo}$ values have been obtained without considering the experimental components imperfections, e. g., BS asymmetry and losses.
    \item Numerical simulations for $p>0$ have been obtained, as said, by averaging over $1000$ disorder configurations, while three phase maps were used in the experiment. Thus, we expect the simulation results to be far more accurate, while the experimental ones are more sensitive to stochastic deviations.
\end{itemize}
It is worth noting that, for the ordered case, the values of $\beta_{Theo}$ and $\beta_{Fit}$ are lower than the expected value of $2$. Indeed the quadratic growth is an asymptotic behaviour, expected to be achieved for long evolution times.
We can conclude that a superdiffusive behaviour, uniquely identified by the condition $1<\beta<2$, can be simulated using $p$-diluted QWs, making them a useful tool for further investigations on diffusion processes.
 \begin{center}
 \begin{table}[!h]
 \color{black}{
 {\small{}
  \hfill{}
  \begin{tabular}{c|c|c}
  \em $p$ & \em $\beta_{Theo}$ & \em $\beta_{Fit}$ \\
  \hline
  $0$ & 1.69 & 1.64 $\pm$ 0.10 \\
  \hline
  $0.05$ & 1.540 & 1.433 $\pm$ 0.067 \\
  \hline
  $0.10$  & 1.414 & 1.277 $\pm$ 0.061 \\
  \hline
  $0.20$ & 1.198 & 1.160 $\pm$ 0.060 \\
  \hline
  $1$ & 0.921 & 0.961 $\pm$ 0.022
 \end{tabular}}
 }
 \hfill{}
  \caption{\emph{Expected ($\beta_{Theo}$) and experimental ($\beta_{Fit}$) values of the exponent $\beta$ for different values of $p$. The values between $1$ and $2$ for $p>0$ confirm the superdiffusive behaviour of the evolution.}}
 \label{tab:fit_results}
 \end{table}
 \end{center}
We also investigated the case of two photon QW, with photon pairs entering into the setup through the two $BS_1$'s input ports (see Fig.3). The variance of the mean position of two indistinguishable walkers, whose expression is
\begin{equation}
\label{vartwophotons}
Var_{2}(n)=\sum_{i,j=-n}^n(\frac{i+j}{2})^2\cdot P_{i,j}(n)-(\sum_{i,j=-n}^n\frac{i+j}{2}\cdot P_{i,j}(n))^2
\end{equation}
is given in Fig. \ref{fig:var_vs_step_2photons} for $p=0.10$. In this case the number of steps is limited to $n=5$ mostly for two reasons: the losses present in the apparatus, limiting the number of detectable coincidences, and the growth of the number of modes on which the photons can travel. Experimental data are in very good agreement with theoretical simulation taking into account the real setup parameters. Simulations of the behaviour of $Var_2$ in the ordered and disordered case are also reported. Error bars are larger than the ones in the single photon case due to the fact that in this case a little variation of the phase must be considered twice, because it is experienced by both photons (see SM for further details on the two photon case).
\begin{figure}[!h]
 	\centering
 	\includegraphics[scale=0.60]{var_2_sper_teo_teoreal.eps}%
 	\caption{\emph{Experimental results for the two photon QW. Dots represent experimental data, dashed lines indicate the expected behaviour obtained by taking into account the real parameters of the setup (Exp), dotted lines correspond to the limit cases of the ideal ordered and disordered two photon QW (I). Note that, up to the third step, the variance is not affected by the phase thus the red and blue dotted lines are coincident.}}
 	\label{fig:var_vs_step_2photons}
 \end{figure}

\textit{Conclusion} - We have experimentally investigated the superdiffusion process, intermediate between the quantum and classical transport regimes, by introducing tunable disorder, through suitable phase map configurations, in a discrete QW. The intrinsically stable bulk optic system used in the experiment allows to operate with any kind of phase map. By increasing the number of steps exploited by the optical setup, the subdiffusive regime, localized between the diffusive and the localized one, can also be experimentally explored.
The $p$-diluted QW approach described in this work makes possible to simulate a decoherent process without altering the coherent evolution of the walker.
We are currently investigating on the relationship existing between the two models.

\date{\today}

\begin{acknowledgments}
We acknowledge support from the European Commission grants FP7-ICT-2011-9-600838 (QWAD - Quantum Waveguides Application and Development). We thank Mauro Paternostro, Beatrice Polacchi, Camilla Sarra and Federico Pegoraro for useful help and discussions.
\end{acknowledgments}