text:
\section{Introduction}

Quantum temporal correlations are analyzed with diverse methods by utilizing histories or trajectories of evolving quantum systems with more recent emphasis on mathematical formulation of the entangled superposition of quantum histories in Reference~\cite{cotler2016}, i.e., denoted with the entangled histories framework. These varying methods include Feynman's path integral (FPI) formalism~\cite{feynman} as the most fundamental of all inherently including histories, consistent histories approach defined by Griffiths~\cite{griffiths1984consistent,griffiths1993consistent,griffiths2003consistent}, and the recently formulated entangled histories framework~\cite{cotler2016} and two-state vector formalism~\cite{aharonov2008two, nowakowski2018} while all emphasizing correlations in time as standard quantum mechanical (QM) formalisms without violating Copenhagen interpretations. Multiplane diffraction (MPD) design as a simple linear optical system was recently proposed for quantum computing (QC)~\cite{bg1, gulbahar2019quantumfourier} and for modulator design in classical optical communications~\cite{gulbahar2019quantum} by exploiting the tensor product structure of quantum temporal correlations as quantum resources while utilizing only the classical light sources and conventional photon-counting intensity  detection. The MPD architecture generates interference of an exponentially increasing number of propagation trajectories along the diffraction events through multiple slits on the consecutive planes. The simplicity of source and detection in MPD setup combined with the highly important  promise of the utilization of the tensor product structure of the temporal correlations as quantum resources  motivates the definition and study of quantum trajectories or histories in MPD as novel quantum resources. These new resources denoted as \textit{quantum path entanglement} (QPE) and \textit{quantum path interference} (QPI) are defined and theoretically modeled in this article in terms of the temporal correlations and interference among the trajectories, respectively, to be exploited for future quantum computing and communications systems.

In this article, MPD design is, for the first time, proposed for defining QPE and QPI as novel quantum resources. Operator theory modeling for MPD-based resources is presented by combining the consistent histories approach of Griffiths~\cite{cotler2016, griffiths1984consistent,griffiths1993consistent,griffiths2003consistent}  and the entangled histories framework in Reference~\cite{cotler2016} with the FPI approach as the inherent structure of MPD creating Feynman paths. MPD creates quantum propagation paths through individual slits in a superposition in which the linear combinations result in evolving quantum history states.   It has low experimental complexity  with classical light sources and conventional photon-counting detection for near-future experimental verification. The theory of QPE and QPI based on MPD proposed in this article provides a set of tools to explore new structures composed of the correlations and interference among the paths for future applications in quantum computing and communications and provides QM foundational studies based on quantum histories.

The concept of the entangled histories  is defined in   References~\cite{cotler2016,cotler2017} as the quantum history which cannot be described as a definite sequence of states in time. There is a superposition of multiple timelines of sequences of events. In this article, we follow similar terminology and denote the temporal correlation among the quantum propagation paths unique to the MPD design with QPE, i.e., emphasizing the entanglement among the path histories similar to References~\cite{cotler2016,cotler2017}.  Tensor product structure among the temporal correlations of multiple time instants is utilized as a novel resource for  computing in References~\cite{bg1, gulbahar2019quantumfourier} and for communications in Reference~\cite{gulbahar2019quantum} in an analogical manner to the multiparticle spatial correlations of the conventional quantum entanglement resources. MPD provides a simple system design inherently including such states having correlations among the paths denoted with QPE. A concrete example of a history state in MPD composed of diffraction events through $N$ planes  is defined as follows:
\begin{equation}
\label{eq_0_9b3d_4dff}
 \sum_{n} \pi_n \left[ \mathbf{P}_{N, s_{n,N}} \right]  \odot \left[ \mathbf{P}_{N-1, s_{n,N-1}} \right]  \odot ...  \odot \left[\mathbf{P}_{1, s_{n,1}} \right] \odot \left[ \rho_0 \right]
\end{equation}
where $\mathbf{P}_{j, s_{n,j}}$ is the projection operator for diffraction through the slit indexed with  $s_{n,j}$ on $j$th plane and for $n$th trajectory,  $\pi_n$ as $0$ or $1$ allows to choose a compound set of trajectories, $\odot$ denotes tensor product operation, and $ \left[ \rho_0\right]$ denotes the initial  state. The quantum state of the light after diffraction through consecutive $N$ planes includes a superposition of different trajectories through the slits.  Experiments for entangled histories has just been, for the first time,  performed in Reference~\cite{cotler2017} by using the polarization states of a single photon and by creating Greenberger--Horne--Zeilinger (GHZ)-type states. MPD-based design compared with complex single photon setup allows the classicality of light sources and simple intensity detection (or photon counting) as a significantly low complexity tool to study quantum histories and QM foundations with near-future experiments. MPD  utilizes simple and widely available coherent sources such as Gaussian wave packets of standard laser output conventionally denoted as classical light.

In this article, an important property of MPD-based QPE is, for the first time, presented: Leggett--Garg Inequality (LGI) violations  as the temporal analog of Bell's inequality. One of the fundamental tools to analyze  quantum temporal correlations of a system is to check the violations of LGIs~\cite{lg0}. LGIs, as proposed by Leggett and Garg in 1985, check a system in terms of the fundamental principles of  macroscopic realism (MR) and noninvasive measurability (NIM) such that the systems obeying these rules satisfy the intuition about the classical macroscopic  world~\cite{lg1}.  QM systems violate LGIs such that MR principles implying the existence of a preexisting value of a macroscopic system and the NIM principle implying the measurement of the value without disturbing the system are both invalidated~\cite{emary2012leggett, wilde2012addressing}.  LGI   violations~\cite{lg0,lg1,emary2012leggett, lg2,katiyar2017experimental} are utilized for various purposes such as  testing temporal correlations of a single system as an indicator of  the  quantumness and analyzing  QC systems, e.g., Grover's algorithm violating temporal Bell inequality~\cite{Morikoshi}. The simple LGI inequality with three-time formulation violated with various QM setups is defined as follows:
\begin{equation}
\centering
\label{eq_d_9ec1_4bc5}
 C_{01} + C_{12} - C_{02} \leq 1
\end{equation}
where  $C_{ij} \equiv \left\langle Q_{i} Q_{j}\right \rangle $ is the expected value of the multiplication of the dichotomic observables $Q_{i}$ as the measurement outcomes at time $t_i$. The left-hand side is maximally violated by QM systems with the value of $3 \,  / \,2$. The violation analysis of LGIs  is, for the first time, performed for MPD by utilizing the recently proposed ambiguous form by Emary in Reference~\cite{lg2} with the precautions regarding the signaling-in-time (SIT) problem in order to convince a macrorealist about the noninvasive nature of measurements, i.e., to prevent signaling forward in time with measurements. This is achieved by inferring event probabilities from ambiguous measurements rather than direct measurements and by modifying the fundamental inequality in Equation (\ref{eq_d_9ec1_4bc5}) by including a signaling term and by providing a NIM-free bound as described in detail in the Results section.  The violation of LGI with no-signaling assumption reaching $>$0.2, i.e., left-hand side of $>$1.2, is numerically obtained for three-time formulation of LGI in MPD setup. The optimization study to maximize it to the calculated bounds~\cite{lg2} is left as an open issue. Besides that, a novel system design, i.e., MPD, violating LGIs with classical light sources is proposed in this article, complementing the recent experimental result in Reference~\cite{zhang2018experimental} utilizing linear polarization degree of freedom of the classical light to violate LGIs. However, MPD utilizes photon-counting intensity detection with a significantly low experimental complexity. It is also simpler compared with the LGI violating architectures utilizing single-photon sources and Mach--Zehnder interferometers~\cite{cotler2017, wang2017enhanced, xu2011}. Besides that, light sources not fully coherent in terms of spatial and temporal dimensions are theoretically modeled while the violation of LGI and QPI are numerically analyzed for specific MPD setup geometry satisfying coherence of light under Gaussian source beam assumptions.

{On the other hand,  LGIs are interpreted in a quantum contextual framework in Reference~\cite{asano2014}, where the contextuality implies the impossibility to consider a quantum measurement as revealing a preexisting property independent of the set of measurements. It is also analyzed in relation with consistent histories approach in Reference~\cite{losada2014}. Furthermore, nonlocality and contextuality are presented as important quantum resources~\cite{howard2014}. Therefore, the relation of the proposed QPE and QPI resources with quantum contextuality is an open issue to be explored.}

The other important property of quantum histories is the interference among them denoted by QPI. To the best of the author's knowledge, theoretical modeling of the interference among quantum history states leading to a counterintuitive observation to be easily verified experimentally has not been previously formulated.  Implementation of the theoretically modeled QPI setup will significantly improve our understanding about QM fundamentals regarding time. QPI is the temporal analogue of the spatial interference obtained in Young's double-slit setup. Destructive and constructive interferences among the paths are observed in the time domain for the QPI case. A special case is modeled such that decreasing the number of photons to diffract through a plane by removing a Feynman path results in an increase in the number of photons  diffracting through the next plane due to the interference between two quantum trajectories. This is proposed, for the first time, as a counterintuitive nature of the interference among the quantum histories.

The novel contributions of the article are summarized as follows:
\begin{itemize}
\item introduction and operator theory modeling of two novel quantum resources, i.e., QPE and QPI, denoting temporal correlations and the interference among quantum trajectories, respectively, in MPD while utilizing the tensor product structure for future quantum computing and communication architectures and foundational QM studies;
\item operator theory modeling of MPD-based resources QPE and QPI by combining history-based previous formulations of quantum histories~\cite{cotler2016, griffiths1984consistent, griffiths1993consistent, griffiths2003consistent} with FPI formalism;
\item theoretical modeling and numerical analysis of MPD setup for the violation of LGI, with the ambiguous and no-signaling forms recently proposed by Emary in Reference~\cite{lg2}, reaching $>1.2$ of correlation amplitude numerically obtained for three-time formulation while leaving the maximization of the violation to the boundary levels as an open issue;
\item a novel setup, i.e., MPD, violating the ambiguous form of LGI with classical light sources complementing the recent experiment utilizing linear polarization degree of freedom of the classical light~\cite{zhang2018experimental} while MPD setup with remarkably low complexity design utilizing classical light sources and photon-counting intensity detection;
\item theoretical modeling and numerical analysis of counterintuitive properties and examples of the interference among MPD-based Feynman paths denoted as QPI promising to be easily verified experimentally in future studies;
\item the modeling and numerical analysis of the coherence properties of the light sources in terms of spatial and temporal dimensions while discussing design issues for MPD setup with coherent light sources; and
\item discussion for future applications of QPE and QPI as quantum resources and  experimental implementations.
\end{itemize}

The paper is organized as follows. We firstly define MPD setup with diffractive projection and measurement operators in Sections \ref{sec2-1} and \ref{sec2-2}. It is followed by the history state modeling of QPE in Section \ref{sec2-3}. Then, we present theoretical modeling of the violation of LGI in Section \ref{sec2-4}, followed by QPI scenario in Section \ref{sec2-5}. Then, numerical analysis is presented in Section \ref{sec2-6}.  We provide the conclusions and discuss future applications of QPE and QPI based on MPD setup in Section \ref{sec3}. Finally, the methods utilized for theoretical modeling are presented in Section \ref{sec4}.

\section{Results}
\vspace{-6pt}
\subsection{MPD Setup for Quantum Temporal Correlations}
\label{sec2-1}

MPD setup is formed from $N-1$  diffraction planes of multiple slits in front of a  classical light source and the measurement of interference pattern  with $N$ sensor planes, i.e., both diffraction and sensing on the same plane, as shown in Figure \ref{Fig1}a. It is also possible to locally count the diffracted photons with the measurement planes inserted between the diffraction planes as discussed in Section \ref{sec2-6}. The utilized light source is assumed to be coherent as the closest analog of a classical light field emphasizing the absence of nonclassical states of light such as single photon generation, squeezed light, or multiple particles of entangled photons~\cite{zavatta2004quantum}. The standard laser output is almost perfectly a coherent state corresponding to the fundamental transverse modes of light field distribution producing Gaussian beams. This coherent Gaussian wave function keeps the position and momentum uncertainties stationary as emphasized by Glauber~\cite{glauber2007quantum}. It is an eigenstate of the annihilation operator $\hat{a}$ for the harmonic oscillator, i.e., $ \hat{a} \, \ket{\alpha} = \alpha \, \ket{\alpha}$, represented as follows in the complete orthonormal basis of the number  states $\ket{n}$ of the single mode oscillator~\cite{glauber2007quantum}:
\begin{equation}
\ket{\alpha} = e^{-\vert \alpha \vert^2 \,/ \, 2} \sum_{n} \frac{\alpha^n}{(n!)^{1/2}} \ket{n}
\end{equation}
where its representation in the position basis gives the Gaussian form.   Therefore, the source is assumed to have normalized Gaussian wave function  $  \Psi_0 (x_0) \equiv \mbox{exp}\big(-x_0^2/(2 \sigma_0^2)\big)/ \sqrt{\sigma_0 \, / \, \sqrt{\pi}}$ with the standard deviation term $\sigma_0$.

\begin{figure}[H]
\centering
\includegraphics[width=6.2in]{Fig1.eps}
\caption{ (\textbf{a}) System model of the free propagating light with velocity $c$ in the $z$-direction and MPD through $N$ planes, where $j$th plane includes $S_j$ slits at positions $X_{j,i}$ for $i \in [1, S_j]$ and interplane distance of $L_{j, j+1}$. (\textbf{b}) Example of three plane diffractions ($N = 3$) with two slits for the first and second planes showing all the possible seven types of histories composed of diffractions or projections $P_{1,1}$, $P_{1,2}$, $P_{2,1}$, and $P_{2,2}$ through slits and measurements $M_{1}$, $M_{2}$, and $M_{3}$ on the planes. { There are $N_p \equiv \prod_{j=1}^{N-1} S_{j} = 2 \times 2 = 4 $ paths detected on the third plane.}}
\label{Fig1}
\end{figure}

Each plane is assumed to be capable of performing measurement with photodetectors for counting the number of photons hitting the detector area. Therefore, a plane either allows  projective diffraction of  light  through slits  denoted by the operator symbol $\mathbf{P}$ or performs measurement  denoted by $\mathbf{M}$ on its sensor array positions where there are no slits.  Gaussian slits are utilized with FPI modeling for simplicity~\cite{feynman,bg1} as mathematically described in Equation (\ref{eq_a_ef97_41ec}) in the next subsection.   Light  is assumed to perform free space propagation between consecutive planes. The plane with the index $j$  has $S_j$ slits, where the central positions and widths of slits are  denoted  by $X_{j, i}$ and $W_{j,i}$, respectively, and  $j \, \in [1, N-1]$ and $i \, \in [1, S_j]$. The widths of the slits are assumed to be the same on each plane but not constrained among different planes. Distance between the $i$th and $j$th planes is denoted by $L_{i,j}$, where the distance from  the light  transmitter source to the first plane is given by $L_{0,1}$.  Light is assumed to have propagation  in the $z$-axis with the velocity given by  $c$,  while  quantum superposition interference is observed in the $x$-axis as a one-dimensional model which can be easily extended to two dimensions (2D)~\cite{bg1}.   Interplane distances and  durations  are denoted by  the vectors $\vec{L}^T =[L_{0,1} \, ... \,L_{N-1, N}]$ and $\vec{t}^T = [t_{0,1}\, ... \, t_{N-1,N}] \equiv \vec{L}^T \, / \, c$, respectively, where transpose is denoted by $(.)^T$. The value $\vec{t}^T$ is accurate with the assumption $L_{j-1,j} \gg W_{j,i}, X_{j,i}$ for  $j \, \in [1, N-1]$  and $i \in [1, S_j]$  such that  QM effects are emphasized in the $x$-axis. Nonrelativistic modeling  is assumed.  We do not consider the effects of environment dephasing or decohering of the interference pattern for double-slit setups~\cite{chen2018,divincenzo1998decoherence}. Furthermore, minor effects of exotic paths~\cite{exotic}  on the numerical results  are ignored as discussed in Reference~\cite{bg1} without affecting the main modeling.

Free-particle evolution kernel for the optical propagation paths between  time--position values $(t_j, x_j)$ and $(t_{j+1}, x_{j+1})$ is defined as follows~\cite{gulbahar2019quantum, gulbahar2019quantumfourier} with the same form for electron propagation~\cite{feynman, bg1}:
\begin{equation}
\label{eq_c_b338_4fde}
K(x_{j+1}, t_{j+1}; x_j, t_j) = \sqrt{ \frac{m}{2 \, \pi \, \imath \, \hbar \, \Delta t}} \,  \mbox{exp} \bigg( \frac{\imath \, m \,  \Delta x^2 }{2 \, \hbar \, \Delta t} \bigg)
\end{equation}
where $\Delta t = t_{j+1} - t_j$, $\Delta x = x_{j+1} - x_j$,  $m \, \equiv \hbar \, k \, / \, c$ is the virtual mass term for the photon with the wave number $k \, = 2 \, \pi \,/ \, \lambda$, and $\lambda$ is the wavelength of the light.

The validity of Fresnel diffraction formulation for quantum optical propagation  is verified based on recent  experimental~\cite{santos2018huygens} and theoretical~\cite{sawant2014nonclassical} studies, while Fourier optics~\cite{ozaktas2001fractional} extension of MPD is recently proposed in   Reference~\cite{gulbahar2019quantumfourier}. Therefore, the Fresnel diffraction integral for free space proposed in Equation (\ref{eq_c_b338_4fde}) and its consecutive application with FPI formalism are theoretically valid and highly reliable for the simple design of MPD. The proposed theoretical model significantly promises to be verified with near-future experiments due to the simplicity of the setup. Then, the propagated wave function  $\ket{\Psi_j} = \int_{ -\infty}^{\infty}  \diff x_j \ket{x_j}    \Psi_j(x_j)  $ on the $j$th plane  becomes as follows by utilizing  Equation (\ref{eq_c_b338_4fde}) consecutively in FPIs~\cite{bg1, gulbahar2019quantum, gulbahar2019quantumfourier}:
 \begin{equation}
 \label{eq_3_6538_4339}
 \Psi_j(x_j)  \, \equiv \,  \sum_{n = 0}^{N_j-1}  \psi_{j,  n  }(x_j) \, \equiv \,
  \sum_{n = 0}^{N_j-1}  \Bigg(  \Upsilon_j\,  e^{(A_{j-1} \, + \, \imath \, B_{j-1})\, x_j^2  }
 \,   e^{ \overrightarrow{x}_{n}^T \, \mathbf{H}_{j-1}  \, \overrightarrow{x}_{n } }
  \, e^{ (\overrightarrow{c}_{j-1}^T \, + \,  \imath \, \overrightarrow{d}_{j-1}^T)  \overrightarrow{x}_{n } \, x_j}          \Bigg)
\end{equation}
where $\psi_{j, n  }(x_j) $ is the contribution for each $n$th propagation path through the slits on the overall superposition, and the definitions of the notations $n$ and $N_j$ are explained next while $\Upsilon_j  = \chi_0 \, \big( \prod_{l=1}^{j-1} \sqrt{\xi_l} \big)$; the constants  $A_{j-1}$,  $B_{j-1}$, $\chi_0 $,  and  $\xi_l$ for $l  \in [1, j-1]$; $\mathbf{H}_{j-1} = \mathbf{H}_{R,j-1} \, + \,  \imath\, \mathbf{H}_{I,j-1}$;  and the vectors $\overrightarrow{c}_{j-1}$ and $\overrightarrow{d}_{j-1}$  depending on  the group of $\lbrace \hbar$,  $m$,   $\sigma_0$, $t_{l, l +1}$, and $\beta_l  \rbrace$ for $l \leq j - 1 $  are explicitly defined in Reference~\cite{bg1}. Explicit forms of the parameters required for double- and triple-plane setups are provided in Section  \ref{sec4} while formulating LGIs and QPI, respectively, in the following discussions.   The total number of paths just before diffraction on the $j$th plane is calculated by $N_j = \prod_{l=1}^{j-1} S_{l}$, while the set  of slit positions for the path indexed with $n  \in [0, N_j -1]$ is denoted by    $\overrightarrow{x}_{n } \equiv [ X_{1,s_{n,1}}  \, \, X_{2,s_{n,2}} \, \, ... \, \, X_{j-1,s_{n,j-1}} ]^T$ while each  $n$th path is indexed by the set of diffracted slits as the~following:
\begin{equation}
\label{eq_patheq}
 Path_n \equiv \lbrace s_{n,1}, \, s_{n,2}, \, ... \, s_{n,j-1}; \, s_{n,l} \in [1, S_l] {, l \in [1, j-1]} \rbrace
\end{equation}
where the specific slit on $l$th plane for $n$th trajectory  is indexed with $s_{n,l}$. The same symbol of the position vector  $\overrightarrow{x}$ is used for both the dimensions $N$ and $j$. The size of the vector is inferred from the index of the current plane analyzed throughout the text.  The position on the $j$th plane is denoted by $x_j$. In Equation (\ref{eq_3_6538_4339}) for $j =N$, each path reaching the $N$th plane is indexed by $n$ for $n \in [0, N_p-1]$ as shown in Figure \ref{Fig1}b  for a simple example of $N_p = 4$,  where total number of paths is given by multiplying the number of slits on each plane as  $N_p \equiv \prod_{j=1}^{N-1} S_{j}$.  The vector $\overrightarrow{x}_{n} \equiv$ $ [ X_{1,s_{n,1}}\, X_{2,s_{n,2}} $ $\, ...  \,$  $X_{N-1,s_{n,N-1}} ]^T$ denotes the set of slit positions ordered with respect to the plane indices for $n$th path for the case of $N$ planes.  Next, diffraction and measurement operators are theoretically defined by emphasizing the operator algebra of multiplane evolution.

\subsection{Diffractive Projection and Measurement Operators}
\label{sec2-2}

Projection operator denotes the  light  to be in the Gaussian slit  in a  coarse-grained  sense~\cite{bg1, dowker1992quantum} as~follows:
\begin{equation}
\label{eq_a_ef97_41ec}
\mathbf{P}_{j,i} \equiv   \int_{ -\infty}^{\infty}  \diff x_j \, \mbox{exp} \bigg ( - \frac{(x_j - X_{j, i})^2}{2 \,\beta_{j,i}^2 }\bigg)\ket{x_j} \bra{x_j}
\end{equation}
where $g_{j,i}(x_j) \equiv \mbox{exp}  \big( -  (x_j - X_{j, i})^2 \, / \, (2 \,\beta_{j,i}^2 ) \big)$ is the slit projection function and the effective slit width is  $W_{j,i} \equiv 2 \, \sqrt{2} \,\beta_{j,i}$, i.e., leading to a $1 \, / \,  e^2$ drop in the intensity,    where  $j \, \in [1, N-1]$ and $i \, \in [1, S_j]$. Projectors are mutually exclusive with high accuracy  such that slit distances are chosen large enough to satisfy $\mbox{exp} \big ( - (X_{j, i_1} - X_{j,  i_2 })^2 \, / \, (2 \,\beta_{j, i_1 }^2) \big) \ll 1$ for $ i_1  \neq i_2$. Total diffraction through all slits of the $j$th plane has the operator  $\mathbf{P}_{j} \equiv \sum_{i = 1}^{S_j} \mathbf{P}_{j,i}$. Measurement operators are redefined due to the proposed Gaussian slit design such that trace preserving equality is satisfied, i.e., $\mathbf{M}_{j}^{\dag} \,\mathbf{M}_{j}  +   \mathbf{P}_{j}^{\dag} \, \mathbf{P}_{j}= \mathbf{I}$, where $\mathbf{I}$ is the identity operator and $(.)^\dag$ or $(.)^H$ denotes Hermitian or conjugate transpose operation.  It is assumed that wave function at time $t = t_0$ evolves to $\ket{\Psi_j}$ and  $\ket{\Psi_j^+}$ for  just before and just after diffraction on the $j$th plane at $t_j^-$ and $t_j^+$, respectively. The state of the  light  at $t_j^+$ has experienced either $\mathbf{M}_{j}$ or $\mathbf{P}_j$. The measurement operator on the $j$th plane is defined as the following:
\begin{equation}
\label{eq_3_51b3_45ac}
\mathbf{M}_{j}^{\dag} \,\mathbf{M}_{j} \equiv \mathbf{I} -\big(\sum_{i = 1}^{S_j} \mathbf{P}_{j,i}^{\dag} \big)  \big(\sum_{i = 1}^{S_j} \mathbf{P}_{j,i} \big)  \equiv   \int_{ -\infty}^{\infty}  \diff x_j \, \Bigg(1 - \bigg(\sum_{i = 1}^{S_j}  e^{ - \frac{(x_j - X_{j, i})^2}{2 \,\beta_{j,i}^2 }} \bigg)^2    \Bigg) \ket{x_j} \bra{x_j}
\end{equation}
Therefore, if we define the measurement operator in FPI formalism as  multiplication  of the wave function with $m_j(x_j)$ reducing the probability to measure the  light  while approaching the slit center, then the following is obtained by using Equation (\ref{eq_3_51b3_45ac}):
\begin{equation}
\label{eq_7_c5a0_4d9a}
\vert m_j(x_j) \vert^2 =  1 - \bigg(\sum_{i = 1}^{S_j}  e^{ - \frac{(x_j - X_{j, i})^2}{2 \,\beta_{j,i}^2 }} \bigg)^2
\end{equation}

There are two different types of detection mechanisms in MPD design denoted by $Rec_1$ and $Rec_2$. In $Rec_1$, all of the planes for $j \in [1, N]$ have detectors measuring the incident  light  and $Rec_1$ is the model proposed in this article forming a complete set of diffractive projection $\mathbf{P}_{j}$ and measurement  $\mathbf{M}_{\widetilde{j}} $  operators until the final detector plane $N$ for $j \in [1, N-1]$ and $\widetilde{j} \in [1, N]$.  In this article, $Rec_1$ modeling is utilized to model history-based time evolution of the  light.  An example is shown in Figure \ref{Fig1}b, where there is a total of seven different sets of consecutive events forming a complete set of histories. The proposed setup is modeled compatible with  the  consistent histories approach  defined  in  Reference~\cite{griffiths1984consistent}  or the entangled histories framework in Reference~\cite{cotler2016}.  On the other hand,  the receiver type with the sensors only on the final plane is denoted by $Rec_2$. In $Rec_2$, i.e.,   the modeling utilized in Reference~\cite{bg1} for QC, only the final intensity distribution or interference pattern on the detector plane is measured.  There is either no detection at the time $t_N^+$ or the  light  is detected on the final detector plane with the index $N$. An operator denoting no detection is defined as $\mathbf{M}_{o}$ to form a complete set for $Rec_2$; then, $  \mathbf{M}_{N}^{\dag} \mathbf{M}_{N}+ \mathbf{M}_{o}^{\dag} \mathbf{M}_{o}= \mathbf{I}$.  Next, consistent histories approach is applied for MPD setup.

\subsection{History State Modeling of QPE}
\label{sec2-3}

Following the definition of consistent histories~\cite{griffiths1984consistent,griffiths1993consistent,griffiths2003consistent}  and  entangled histories~\cite{cotler2016}, a  history state is defined for MPD based on the set of projections $\mathbf{M}_{j}$ and $\mathbf{P}_{j,i}$ on each $j$th plane for $j \in [1, N]$ and $i \, \in [1, S_j]$. History Hilbert space is defined as follows:
\begin{equation}
\label{eq_3_80a4_4ee5}
\mathcal{H} \equiv \mathcal{H}_N \odot  \mathcal{H}_{N-1} \odot ... \odot \mathcal{H}_{1}\odot \mathcal{H}_{0}
\end{equation}
where $\mathcal{H}_j$ denotes the set of projections on planes and $\odot$ denotes tensor product operation.  Hilbert space until $t_j^+$ includes both projections $\mathbf{P}_j$ and $\mathbf{M}_l$ on the planes with the indices $l  \leq j$ since the  light  is detected at some plane until $t_j$ or still diffracting through the $j$th plane. A general history  state  with QPE composed of superposition of trajectories  is denoted as follows based on the notation (similar to bra-ket   but with different notations of  $( . \vert$ and $\vert . )$ for the histories corresponding to $\bra{.}$ and $\ket{.}$, respectively) in  Reference~\cite{cotler2016}:
\begin{equation}
\label{eq_0_9b3d_3dtf}
\vert \Psi_{N} ) = \sum_{n} \pi_n  \left[ \mathbf{O}_n (t_N) \right]  \odot \left[ \mathbf{O}_n (t_{N-1}) \right]  \odot ...  \odot \left[ \mathbf{O}_n (t_0) \right]
\end{equation}
where $\vert \Psi_{j} )$ is some history state between times $t_0$ and $t_j $ for $t_j  > t_0$, the projector $\left[ \mathbf{O}_n(t_j) \right]$ denotes either of $\mathbf{M}_{l}$ or $\mathbf{P}_{l,i}$ for $l \leq j$ and $i \, \in [1, S_l]$, and $\pi_n$ as $0$ or $1$ is some permutation choosing a compound set of histories indexed by $n$.  Observe that $t_j$ includes measurements $\mathbf{M}_l$ for $l \leq j$  as   possible events such that the state does not change after measurement. It also includes events with zero probability such as $j$th plane projection at times not equal to $t_j$. Some examples for $N = 4$ are as follows:
 \begin{align}
 \label{eq_9_c0ec_43a3}
\begin{split}
 \vert \Psi_{4}^a )  \equiv & \left[ \mathbf{M}_1  \right]  \odot  \left[ \mathbf{M}_1  \right]  \odot \left[ \mathbf{M}_1 \right]  \odot \left[ \mathbf{M}_1  \right] \odot \left[ \rho_0 \right] \\
 \vert \Psi_{4}^b ) \equiv &  \left[ \mathbf{M}_4  \right] \odot \left[ \mathbf{P}_{3,2}  \right]  \odot \left[ \mathbf{P}_{2,4} \right]  \odot \left[  \mathbf{P}_{1,1} \right]  \odot \left[ \rho_0 \right] \\
  \vert \Psi_{4}^c ) \equiv &  \left[ \mathbf{M}_4  \right]  \odot \left[ \mathbf{M}_2  \right]  \odot \left[ \mathbf{M}_2 \right]  \odot \left[ \mathbf{M}_1  \right]  \odot \left[ \rho_0 \right]
 \end{split}
\end{align}
The state $\vert \Psi_{4}^a)$ shows that the light is detected on the first plane at $t_1$ while not changing at consecutive time states, i.e., without diffracting even from the first plane. In $\vert \Psi_{4}^b)$, the light is diffracted from the first slit of the first plane at $t_1$, then is diffracted from the fourth slit of the second plane at $t_2$ and from the second slit of the third plane at $t_3$, and is finally measured on the fourth plane. The third example $\vert \Psi_{4}^c ) $ is a state with zero probability due to the orthogonality of the operators on different planes.  A~simple example for three planes with two slits  is shown in Figure \ref{Fig1}b with seven different history states while $N_p = 4$ of them reach the final detector  plane as consecutively diffracted trajectories. History Hilbert space summing to the identity denoted by $\overline{I}_H $ as the  family based upon an initial state and neglecting the histories with zero probability is described as follows~\cite{griffiths1984consistent}:
 \begin{align}
 \label{eq_0_18f8_421f}
\begin{split}
\overline{I}_H =   & \sum_{j =1}^{N}  \sum_{i_{j-1} =1}^{S_{j-1}} \sum_{i_{j-2} =1}^{S_{j-2}}  ... \sum_{i_1 =1}^{S_{1}}      \big( \left[ \mathbf{M}_j  \right]^{\odot \, \alpha}    \odot   \left[ \mathbf{P}_{j-1,i_{j-1}}  \right]   \odot ... \odot  \left[ \mathbf{P}_{1,i_1} \right] \odot \left[ \rho_0 \right]\big)
 \end{split}
\end{align}
where $\left[ \mathbf{M}_j  \right]^{\odot \,  \alpha}$ denotes $ \alpha \equiv N\,+\,1\,-\,j$ consecutive  measurements of $\left[ \mathbf{M}_j  \right]$ on the same plane. This includes all the possible history states and evolution for the light until $t = t_N$ starting from $t_0$. A  chain operator  is  presented in  Reference~\cite{cotler2016} to define the inner product between history states which maps a history state to an operator. The chain operator provides history states with positive semi-definite inner products. This operator is inherently defined in the MPD system as the free-particle evolution kernel $K(x_1, t_1; x_0, t_0)$. Assume that the free-particle evolution operator with the notation $U_{j+1, j}$ acts as the bridging operator connecting projections at times $t_j$ and $t_{j+1}$. Then, chain operator denoted by  $\chi_{t_{j+1}, t_j}$ for the time duration ($t_j$, $t_{j+1}$) is defined as follows:
\begin{eqnarray}
\label{eq_5_b64b_4c52}
\chi_{t_{j+1}, t_j} \lbrace \left[ \mathbf{P}_{j+1} \right] \odot \left[ \mathbf{P}_j\right] \rbrace & \stackrel{\mbox{1}}{=} & \mathbf{P}_{j+1} \, U_{j+1, j}   \, \mathbf{P}_{j}  \hspace{0.2in}\\
\label{eq_d_8d03_4376}
\chi_{t_{j+1}, t_j} \lbrace \left[ \mathbf{M}_{j+1} \right] \odot \left[ \mathbf{P}_j\right] \rbrace & \stackrel{\mbox{2}}{=} & \mathbf{M}_{j+1} \, U_{j+1, j}   \, \mathbf{P}_{j}  \hspace{0.2in} \\
\label{eq_c_3117_498d}
\chi_{t_{j+1}, t_j} \lbrace \left[ \mathbf{M}_{j} \right] \odot \left[ \mathbf{M}_j\right] \rbrace & \stackrel{\mbox{3}}{=} &  \mathbf{M}_{j} \, \mathit{I}  \, \mathbf{M}_{j}   \hspace{0.2in}\\
\label{eq_d_9629_4e6b}
\chi_{t_{j+1}, t_j} \lbrace \left[ \mathbf{O}_n(t_{j+1}) \right] \odot \left[ \mathbf{O}_n(t_{j}) \right] \rbrace & \stackrel{\mbox{4}}{=}  & \mathbf{O}_n(t_{j+1}) \, \mathit{I}   \, \mathbf{O}_n(t_{j})  \hspace{0.2in}
\end{eqnarray}
where $\left[\mathbf{O}_n(t_{j+1})\right]$ and $\left[ \mathbf{O}_n(t_{j}) \right]$ in $\stackrel{\mbox{4}}{=}$ denote the cases which are not presented in the first three definitions. $\mathit{I} $ is the identity operator equalizing the consecutive measurements on the same plane, i.e., $\mathbf{M}_{j}^l = \mathbf{M}_j$, for any integer $l$. Furthermore, it bridges dynamically not possible history states which have zero probability to occur as discussed in  Reference~\cite{griffiths1984consistent}. These  include  consecutive measurements on different planes such as $\left[ \mathbf{M}_{j+1}  \right]  \odot   \left[\mathbf{M}_{j}  \right]$, future projection or measurements at a previous time such as $\left[ \mathbf{M}_{j}  \right]  \odot   \left[\mathbf{P}_{j+1}  \right]$, or consecutive sets of the same projector $\mathbf{P}_{j}$ at future times such as  $\left[ \mathbf{P}_{j}  \right]  \odot   \left[\mathbf{P}_{j}  \right]$, where  free-space propagation in the $z$-axis prevents this. Then, the compound history state mapped or affected by the chain operator is defined as follows:
\begin{align}
 \begin{split}
 \label{eq_6_3602_4f02}
 \chi_{t_{N}, t_0} \vert \Psi_{N} ) \equiv
  \sum_{n} \pi_n \, \mathbf{O}_n(t_N)   \, V_{N, N-1}  \, \mathbf{O}_n(t_{N-1})    \,  ...  \, V_{1, 0} \, \mathbf{O}_n(t_0)   &
 \end{split}
\end{align}
where $V_{j+1, j}$  denotes either $U_{j+1, j}$ or $\mathit{I}$.

Besides that, MPD allows to model and explore varying kinds of superposition of history states and QPEs similar to the specific entangled states discussed in  Reference~\cite{cotler2016} resembling the temporal counterpart of Bell states. For example, entangled history states of the  GHZ type is experimentally tested in  Reference~\cite{cotler2017}. It is an open issue to utilize MPD to generate and test such states with important implications and applications based on QPE. Next, probability amplitudes of  histories are~modeled.

\subsubsection{Event Probabilities}
\label{sec2-3-1}

The probabilities characterize the statistical properties of the measurement of classical light. It is assumed that the probability is proportional to the square of the wave function with Born's postulate. It is calculated by integrating the number of photons on the detector area at a specific position for a time interval $T$ enough to obtain the statistical properties~\cite{sawant2014nonclassical, santos2018huygens}. The normalized probability is easy to calculate by measuring the number of photons for each event by forming a histogram and then by dividing the number of photons for the specific event to the total number of source photon counts. The number of photons at a particular position $x$ is frequently denoted with the integral element  $\vert \Psi(x) \vert^2 dx$, while $\vert \Psi(x) \vert^2$ is denoted as the intensity of the light at the particular position.  The probability for the particular history state is found with the positive semi-definite inner product defined as follows:
\begin{align}
 \begin{split}
\label{eq_f_79d2_4a98}
(\Phi_{N} \vert \Psi_{N} )  \equiv  \, & \tr \lbrace \big({\chi}_{t_N, t_0} \vert \Phi_{N}) \big)^H \big(\chi_{t_N, t_0} \vert \Psi_{N}) \big) \rbrace
 \end{split}
\end{align}
where $\tr \lbrace . \rbrace$ is the trace operation. Assume that two specific elementary history states corresponding  to specific diffraction paths indexed with $\widetilde{n}$ and $\widehat{n} \in [0, N_p -1]$ composing the superposition wave function  in Equation (\ref{eq_3_6538_4339})  are  denoted by $\vert \psi_{N,\widetilde{n}})$ and $ \vert \phi_{N,\widehat{n}} )$, respectively.
These paths include only the diffraction projections at the planes with the indices $j \in [1, N-1]$  denoted by $\mathbf{P}_{j, s_{\widetilde{n},j}}$ and $\mathbf{P}_{j, s_{\widehat{n},j}}$, respectively. If the initial  state  $\rho_0 = \ket{\Psi_0}\bra{\Psi_0}$  and  $\mathbf{M}_N = \mathbf{I}$ are included, then the weight of an elementary diffraction history denoted by  the inner product $ W_{\widetilde{n}} \equiv (\psi_{N,{\widetilde{n}}} \vert \psi_{N,{\widetilde{n}}} ) $ in Reference~\cite{griffiths1984consistent} becomes the following:
 \begin{eqnarray}
 \label{eq_f_b2bc_49fd}
 W_{\widetilde{n}}  & = & \tr \bigg \lbrace U_{N, N-1} \, \mathbf{P}_{N-1, s_{\widetilde{n}, N-1}} \, ... \,  \mathbf{P}_{1, s_{\widetilde{n},1}} \, U_{1,0} \, \rho_0    \, U_{1,0}^{\dag} \, \mathbf{P}_{1, s_{\widetilde{n},1}}^{\dag}  \,  ... \,  \mathbf{P}_{N-1, s_{\widetilde{n},N-1}}^{\dag} \, U_{N, N-1}^{\dag} \, \mathbf{M}_N \bigg \rbrace   \\
  & = & \tr  \bigg  \lbrace \rho_0 \, U_{1,0}^{\dag} \, \mathbf{P}_{1, s_{\widetilde{n},1}}^{\dag}   \,  ... \,  \mathbf{P}_{N-1, s_{\widetilde{n},N-1}}^{\dag} \, U_{N, N-1}^{\dag}    \, U_{N, N-1} \, \mathbf{P}_{N-1, s_{\widetilde{n}, N-1}} \, ... \,  \mathbf{P}_{1, s_{\widetilde{n},1}} \, U_{1,0} \, \rho_0 \bigg  \rbrace       \\
  & = &  \int_{x_N = -\infty }^{\infty}  \diff x_N \, \vert \psi_{N,\widetilde{n}}(x_N) \vert^2
\end{eqnarray}
where  the trace  is realized with respect to the position, $\tr\lbrace \rho_0^2\rbrace = 1$ is utilized, and $\psi_{N, \widetilde{n}}(x_N) $ in position basis of the $N$th plane is calculated by putting $j = N$ and $n = \widetilde{n}$ in the defined wave function $\psi_{j,n}(x_j)$ in Equation (\ref{eq_3_6538_4339}). Similarly,  inner product between history states is defined as follows:
\begin{align}
 \begin{split}
 \label{eq_c_2ed9_418d}
 (\psi_{N,\widehat{n}} \vert \psi_{N, \widetilde{n}} )
     = \int_{x_N = -\infty }^{\infty}  \diff x_N \, \psi_{N, \widehat{n}}^{*}(x_N) \, \psi_{N, \widetilde{n}}(x_N)  & \\
 \end{split}
\end{align}
The probability for the light to be diffracted through the $i$th slit on the $j$th plane with the projection $\mathbf{P}_{j,i}$ is  denoted by $Prob_{j,i}^P$.  Similarly, probability to be measured on  the $j$th plane with measurement projection $\mathbf{M}_j$  is denoted by $Prob_{j}^M$. $Prob_{j,i}^P$ is calculated by using the weight of the compound history $\Omega_{N,\lbrace j,i \rbrace}$ including the targeted event $\mathbf{P}_{j,i}$ as follows:
 \begin{eqnarray}
 \label{eq_7_56f7_40e5}
Prob_{j,i}^P \equiv  & \big(\Omega_{N,\lbrace j,i \rbrace} \big \vert \Omega_{N,\lbrace j,i \rbrace} \big )
 \end{eqnarray}
where $\Omega_{N,\lbrace j,i \rbrace}$ is defined as follows:
\begingroup\makeatletter\def\f@size{9}\check@mathfonts
\def\maketag@@@#1{\hbox{\m@th\fontsize{10}{10}\selectfont\normalfont#1}}%
\begin{align}
 \begin{split}
 \label{eq_f_9d11_4b99}
 \Omega_{N,\lbrace j,i \rbrace}  =  \sum_{n}  \sum_{i_{j-1} =1}^{S_{j-1}} \sum_{i_{j-2} =1}^{S_{j-2}}  ... \sum_{i_1 =1}^{S_{1}} \,
 \bigg( \left[ \mathbf{O}_n(t_N) \right]  \odot   ...  \odot \left[ \mathbf{O}_n(t_{j+1})  \right] \odot   \left[ \mathbf{P}_{j,i}  \right]
   \odot  \left[ \mathbf{P}_{j-1,i_{j-1}}  \right] \odot ... \odot  \left[ \mathbf{P}_{1,i_1}  \right] \odot \left[  \rho_0 \right]  \bigg)
  \end{split}
\end{align}
\endgroup
where elementary diffraction history states include  diffraction events $\mathbf{P}_{l,i_l}$ for $l < j$ and $i_l \in [1, S_{l}]$ until the $j$th plane and diffraction event $\left[ \mathbf{P}_{j,i}  \right]$ on the $j$th plane at $t_j$, and where the events $\left[ \mathbf{O}_{n}(t_{j+1}) \right] $ to $\left[ \mathbf{O}_{n}(t_N) \right] $ denote  any dynamically possible projector at the times between $t_{j+1}$ and $t_N$.  Probability for the events after diffraction will not have any effect on diffraction probability through $\mathbf{P}_{j,i}$, and those projections are discarded. Then, it  is easily calculated by using  Equations (\ref{eq_3_6538_4339})  and (\ref{eq_a_ef97_41ec}) and  with $ \bra{\Psi_j} \,\mathbf{P_{j,i}^{\dag}}\,\mathbf{P_{j,i}} \, \ket{\Psi_j}$ as follows:
  \begin{eqnarray}
   \label{eq_4_2570_4f7f}
Prob_{j,i}^P =
   \mathlarger{\int_{-\infty }^{\infty} } \diff x_j \, e^{ - \frac{(x_j - X_{j, i})^2}{2\, \beta_{j,i}^2 } } \, \big \vert    \Psi_j(x_j) \big \vert^2
 \end{eqnarray}
$Prob_{j}^M$ is calculated with  $Prob_{j}^M =  \mathlarger{\int_{-\infty }^{\infty} } \diff x_j \,   \big \vert  m_j(x_j)   \Psi_j(x_j) \big \vert^2$. Similarly, diffraction through one of several slits in a superposition of $s$ slits on the $j$th plane is given by the following expression:
  \begin{eqnarray}
  \label{eq_1_81ef_4942}
Prob_{j, \tilde{i}_s}^P =
   \mathlarger{\int_{-\infty }^{\infty} } \diff x_j \, \bigg( \sum_{i \in \tilde{i}_s} e^{ - \frac{(x_j - X_{j, i})^2}{2 \beta_{j,i}^2 } } \bigg)^2 \, \big \vert    \Psi_j(x_j) \big \vert^2
 \end{eqnarray}
where  $\tilde{i}_s \equiv \lbrace i_1, i_2, ..., i_s \rbrace$ and $i_l \in [1, S_j]$ for $l \in [1, s]$, $i_a \neq i_b$.

It is important to emphasize the practical meaning of the probabilities of the light diffractions or measurements on the plane. In practice, the probabilities are proportional to the number of photons for each event, e.g., the number of photons passing through a particular slit integrated over a long measurement time for calculating projection probabilities or the number of photons detected on the specific area of the planes for the measurement projection. Histogram-based modeling for counting the photons for all planes and the slits provides the normalized overall probability  for each event summing to a total of unity. Photon or particle counting with classical light is already achieved in various studies characterizing the exotic properties of the paths~\cite{sawant2014nonclassical} or Fresnel diffraction properties~\cite{santos2018huygens}.

Quantumness and temporal correlations for the MPD system  design are analyzed by explicitly providing theoretical formulation of LGIs next.

\subsection{Modeling of the Violation of LGI in MPD}
\label{sec2-4}

LGIs test the temporal correlations by measuring at different times  in analogy with spatial Bell's inequalities for the entanglement between spatially separated systems~\cite{lg0,lg1}. Three-time correlation-based inequality is defined in Equation (\ref{eq_d_9ec1_4bc5}) as $K = C_{01} + C_{12} - C_{02} \leq 1$, where the bound is violated quantum mechanically with dichotomic systems, i.e., $ Q_{j} = \pm 1$ for $j \in [0, 2]$, reaching the bound $3/2$ for a two-level system   with the maximum LGI violation of $1  \, / \,  2$. $C_{j_1, j_2} \equiv \left\langle Q_{j_1} Q_{j_2}\right \rangle $ is the expected value of the multiplication of the dichotomic observables, which is equal to $C_{j_1,j_2} = \sum_{j_1,j_2} p(j_1,j_2) Q_{j_1} Q_{j_2}$, where $p(j_1,j_2)$ is the probability for the measurement of $Q_{j_1}$  and $Q_{j_2}$ at times  $t_{j_1}$ and $t_{j_2}$, respectively, and $t_2 > t_1 > t_0$.  Noninvasiveness or non-disturbing structure of the measurements should be clearly satisfied in order to reduce the ``clumsiness loophole''~\cite{lg2, wilde2012addressing}, i.e.,  experimental limitations and disturbance of the clumsy measurements making it difficult to convince a macrorealist. In Reference~\cite{lg2}, ambiguous measurements are utilized to revise   Equation (\ref{eq_d_9ec1_4bc5}) by including the effect of signaling. In this article, the same formulation is extended for the MPD setup exploiting simple architecture of slits.

The correlation and entanglement in time are tested with the two-plane setup, where each plane includes triple slits as shown in Figure \ref{Fig2}a.  It is assumed that the  light   diffracting  through the first plane  is  taken into account while calculating probability amplitudes, i.e., utilizing negative measurement techniques. For example, if the measured state is set to $\mathbf{P}_{1,1}$, then the second and third slits are closed, forcing the light to  diffract  through only the first slit setting the measurement result. Furthermore, denote $p_1(i_{1,1})  \equiv  Prob_{1, i_{1,1} }^P$ and $p_1(\lbrace  i_{1,1},  i_{1,2} \rbrace) \equiv  Prob_{1,i_s}^P$, where   $i_s = \lbrace  i_{1,1},  i_{1,2}  \rbrace$ for $i_{1,1}, i_{1,2} \in [1, 3]$ and $i_{1,1} \neq  i_{1,2}$.  The probability $p_1(\lbrace  i_{1,1},  i_{1,2} \rbrace)$  corresponds to the measurement result for $\mathbf{P}_{1, i_{1,1}} \cup \mathbf{P}_{1,  i_{1,2}}$  being projected in one of   the  slits with the indices $ i_{1,1} $ and $ i_{1,2} $ on the first plane. Similarly, $p_1(\lbrace 1,2,3 \rbrace)$ denotes the overall projection on superposition in all three slits.   On the other hand, assume that $p_{1,2}(\lbrace  i_{1,1},  i_{1,2}  \rbrace,   \hat{i}_{2} )$ denotes the probability for the history:
\begin{align}
 \begin{split}
 \label{eq_9_9aaf_4bba}
 \left[ \mathbf{O}_{\hat{i}_2}(t_2) \right]  \odot \big(\left[ \mathbf{P}_{1,  i_{1,1} } \right] + \left[ \mathbf{P}_{1,  i_{1,2} } \right] \big) \odot      \left[  \rho_0 \right]     &
  \end{split}
\end{align}
where $\left[ \mathbf{O}_{ \hat{i}_{2}}(t_2) \right]$ is one of $\left[ \mathbf{P}_{2,1} \right]$, $\left[ \mathbf{P}_{2,2} \right]$, $\left[ \mathbf{P}_{2,3} \right]$, or $\left[ \mathbf{M}_{2} \right]$ denoted by $  \hat{i}_{2}  = 1, 2, 3$, and $4$, respectively. Similar to Equations (\ref{eq_4_2570_4f7f}) and (\ref{eq_1_81ef_4942}), $p_{1,2}(\lbrace  i_{1,1},  i_{1,2}  \rbrace,  \hat{i}_2 )$ for $ \hat{i}_2  \in [1,3]$ is found as follows:
  \begin{eqnarray}
  \label{eq_a_8c5c_4984}
p_{1,2}(\lbrace i_{1,1},  i_{1,2}  \rbrace,  \hat{i}_2 ) =
   \mathlarger{\int_{-\infty }^{\infty} } \diff x_2 \, \bigg(  \mbox{exp} \big( - (x_2 - X_{2, \hat{i}_2})^2 \, / \, (2 \, \beta_{2,  \hat{i}_2 }^2 )  \big)   \bigg)^2 \,
 \bigg \vert     \psi_{2,  i_{1,1} }(x_2) +   \psi_{2, i_{1,2} }(x_2)   \bigg \vert^2 &&
 \end{eqnarray}
where the  elementary wave function  is  found with Equation (\ref{eq_3_6538_4339}) by using $i$ as the path index as follows:
\begin{align}
 \begin{split}
 \label{eq_e_f3b2_4ec9}
\psi_{2,i}(x_2)=   \chi_{0} \sqrt{\xi_1} e^{ \Gamma_1\, x_2^2} \,
  e^{   H_1  \, X_{1,i}^2 }   \, e^{r_1  \, X_{1,i} \, x_2}
  \end{split}
\end{align}
where $\Gamma_1 = A_{1}\, + \,\imath \,B_1$, $H_1 = H_{R,1}\,+\,\imath \, H_{I,1}$, $r_1 = c_1 \, + \, \imath \,d_1$, $i \in [1,3]$, $\chi_0$, $\xi_1$, $A_1$, $B_1$, $H_{R,1}$, $H_{I,1}$, $c_1$, and $d_1$ are defined in Section \ref{sec4}. Similarly, $p_{1,2}(\lbrace i_{1,1}, i_{1,2} \rbrace, 4)$ is defined as follows:
  \begin{eqnarray}
  \label{eq_0_fe97_4605}
p_{1,2}(\lbrace i_{1,1}, i_{1,2} \rbrace, 4) =
   \mathlarger{\int_{-\infty }^{\infty} } \diff x_2 \, \big \vert  m_2(x_2) \big \vert^2 \,
  \big \vert     \psi_{2, i_{1,1}}(x_2) +   \psi_{2, i_{1,2}}(x_2)   \big \vert^2 &&
 \end{eqnarray}
where $i_{1,1}, i_{1,2}  \in [1, 3]$ and $i_{1,1} \neq i_{1,2}$.  $p_{1,2}( i_{1,1},   \hat{i}_2)$, and $p_{1,2}( \lbrace 1, 2, 3 \rbrace,  \hat{i}_2)$ denote  the probabilities for   $\left[ \mathbf{O}_{\hat{i}_2}(t_2) \right]  \odot  \left[ \mathbf{P}_{1, i_{1,1}} \right]     \odot      \left[  \rho_0 \right]$ and $\left[ \mathbf{O}_{\hat{i}_2}(t_2)\right] \odot \big( \left[ \mathbf{P}_{1,1} \right] + \left[ \mathbf{P}_{1,2} \right] + \left[ \mathbf{P}_{1,3} \right] \big) \odot \left[  \rho_0 \right]$, where $i_{1,1}  \in [1, 3]$ and $\hat{i}_2 \in [1, 4]$. The same formulation is valid also for the second plane for $p_2(i_{2,1})$, $p_2(\lbrace i_{2,1}, i_{2,2} \rbrace)$, and $p_2(\lbrace 1,2,3 \rbrace)$ for $i_{2,1}, i_{2,2} \in [1, 3]$ and $ i_{2,1} \neq i_{2,2}$.  Observe that, at time $t_2$, it is assumed that $\mathbf{M}_2$ is also included in calculations providing a complete set  $\mathbf{I} = \mathbf{M}_2 + \sum_{i=1}^{3} \mathbf{P}_{2,i} $.   Negative measurement methodology  for the first plane is utilized such that the light only diffracting  through the first plane is utilized in calculating probabilities. Therefore, all the probability calculations based on Equations~(\ref{eq_4_2570_4f7f}), (\ref{eq_1_81ef_4942}), (\ref{eq_a_8c5c_4984}), and (\ref{eq_0_fe97_4605}) are normalized by $\Gamma_c  \equiv \big( \sum_{i = 1}^3 Prob_{1,i}^P )^{-1}$. The probabilities denoted by $p_j(i_{j, 1})$, $p_j(\lbrace i_{j, 1}, i_{j, 2} \rbrace)$, $p_j(\lbrace 1,2,3 \rbrace)$ for $j \in [1,2]$, $p_{1,2}( i_{1,1},   \hat{i}_{2})$, $p_{1,2}(\lbrace i_{1,1}, i_{1,2} \rbrace, \hat{i}_{2})$, and $p_{1,2}( \lbrace 1, 2, 3 \rbrace,  \hat{i}_{2})$ are assumed to be normalized through the rest of the article. The normalized operator is defined as $\mathbf{P}_{1, i}^N  \equiv \Gamma_c \,\mathbf{P}_{1, i}$ for $i \in [1,3]$.

\begin{figure}[H]
\centering
\includegraphics[width=6.2in]{Fig2.eps} \\
\caption{ (\textbf{a}) The violation of Leggett--Garg Inequality (LGI) with the setup of two planes with triple slits  where the event set at time $t_1$ is $\left[P_{1,1}\right]$, $\left[P_{1,2}\right]$, and $\left[P_{1,3}\right]$ and, at time $t_2$, are $\left[P_{2,1}\right]$, $\left[P_{2,1}\right]$, $\left[P_{2,3}\right]$, and $\left[M_{2}\right]$ and ambiguous measurement setups by closing (\textbf{b}) the third, (\textbf{c}) the second, and (\textbf{d}) the first slits on the first plane.}
\label{Fig2}
\end{figure}

Assume that an ambiguous measurement set of three projections composed by  $\left[ \mathbf{O}^A_1(t_1) \right] \equiv \left[ \mathbf{P}_{1,1}^N \right] + \left[ \mathbf{P}_{1,2}^N \right]$, $\left[ \mathbf{O}^A_2(t_1) \right] \equiv \left[ \mathbf{P}_{1,1}^N \right] + \left[ \mathbf{P}_{1,3}^N \right]$, and $\left[ \mathbf{O}^A_3(t_1) \right] \equiv \left[ \mathbf{P}_{1,2}^N \right] + \left[ \mathbf{P}_{1,3}^N \right]$ is defined. The setups for ambiguous measurements  are shown in Figure \ref{Fig2}b--d, respectively. In addition, an assignment of  dichotomic indices for the measurement results is designed  denoted by $Q_{1, i^A_{1}} \equiv \pm 1$ and $Q_{2, \hat{i}_2} \equiv \pm 1$ for   $\left[ \mathbf{O}^A_{i^A_1}(t_1) \right]$ and  $\left[ \mathbf{O}_{\hat{i}_2}(t_2) \right]$, respectively, where $i^A_1 \in [1,3]$ and $\hat{i}_2 \in [1, 4]$. $Q_0 \equiv 1$ denotes initial condition $\left[ \rho_0 \right]$ with unity probability. These dichotomic indices can be assigned arbitrarily while they are chosen in Section \ref{sec2-6} based on the maximization of LGI violation by comparing all the possible assignment combinations. Then, utilizing a similar architecture to the ambiguous LGI, i.e., Equation (14) in Reference~\cite{lg2}, a conversion matrix $\mathbf{D}$ is defined inferring the probability $p_{1}(i_{1,1})$ from the ambiguous measurements with   $\widehat{p}_{1}(i_{1,1}) \equiv \sum_{ i^A_1} D_{i_{1,1}, i^A_1} \,p_{1}^A(i^A_1)$,  where $p_{1}^A(i^A_1)$ denotes the probability for the history $ [ \mathbf{O}^A_{i^A_1}(t_1) ] \odot \left[\rho_0 \right]$ for $ i^A_1 \in [1, 3]$,   $D_{i_{1,1}, i^A_1 }$ is the element at the $i_{1,1}$th row and the $i^A_1$th column of the conversion matrix $\mathbf{D}$, and $\widehat{p}_{1}(i_{1,1})$ denotes the inferred probability such that a macrorealist will not observe any problem. Similarly, $\widehat{p}_{1,2}(i_{1,1}, \hat{i}_2)$ becomes the following:
\begin{eqnarray}
\label{eq_2_cb6d_4864}
  \widehat{p}_{1,2}(i_{1,1}, \hat{i}_2)\equiv \sum_{ i^A_1}  D_{ i_{1,1},  i^A_1} \, p_{1,2}^A( i^A_1, \hat{i}_2)  = \sum_{ i^A_1}  D_{i_{1,1},  i^A_1} \, \Gamma_c \, p_{1,2}(\lbrace  i^A_{1,1},  i^A_{1,2} \rbrace, \hat{i}_2) &&
  \end{eqnarray}
where $p^A_{1,2}(i^A_1, \hat{i}_2)$ denotes the probability  for the history $\left[ \mathbf{O}_{\hat{i}_2}(t_2) \right] \odot [ \mathbf{O}^A_{i^A_1}(t_1) ] \odot \left[\rho_0 \right]$,  where  $i_{1,1},  i^A_1  \in [1, 3]$ and $ \hat{i}_2  \in [1, 4]$,  while it  is found in Equations (\ref{eq_a_8c5c_4984}) and (\ref{eq_0_fe97_4605}) by normalizing as follows:
\begin{equation}
\label{eq_f_563d_44a7}
p^A_{1,2}(i^A_1, \hat{i}_2) \equiv \Gamma_c \, p_{1,2}(\lbrace i^A_{1,1},  i^A_{1,2}  \rbrace,  \hat{i}_2 )
  \end{equation}
where  $  i^A_{1,1} $ and $ i^A_{1,2}$ correspond to the the event $[ \mathbf{O}^A_{i^A_1}(t_1) ] $, i.e., $ \lbrace i^A_{1,1}, i^A_{1,2} \rbrace$  equals $ \lbrace 1, 2\rbrace $, $  \lbrace 1, 3\rbrace $, and $ \lbrace 2, 3\rbrace $ for $i^A_1  \equiv 1$, $2$, and $3$, respectively.  For example, for the proposed setup, $\widehat{p}_{1}(1) =  (p_{1}^A(1) + p_{1}^A(2) - p_{1}^A(3)) \, / \, 2$ since the following probability relation holds:
\begin{equation}
\label{eq_b_bca1_4564}
Prob_{1,1}^P = \frac{Prob_{1,\lbrace 1, 2\rbrace}^P + Prob_{1,\lbrace 1, 3\rbrace}^P - Prob_{1,\lbrace 2, 3 \rbrace}^P}{2}
\end{equation}
Therefore, $D_{11} = 0.5$,  $D_{12} = 0.5$, and $D_{13} = -0.5$. Similarly, $D_{21} = D_{23} = D_{32} = D_{33} = 0.5$ and   $D_{22} = D_{31} = -0.5$.  Then, a macrorealist is convinced that the inferred probabilities are utilized for the calculations of $C_{01}$, $C_{12}$, and $C_{02}$ with  Equation (\ref{eq_d_9ec1_4bc5}) by replacing  $p_{1,2}(i_{1,1}, \hat{i}_2)$ with  $\widehat{p}_{1,2}(i_{1,1},  \hat{i}_2 )$ and by properly defining the degree of signaling level between  the  first and second planes for the ambiguous measurements increasing the required LGI bound.  Then, following the similar methodology in  Reference~\cite{lg2} (Equations (5) and (14) in Reference~\cite{lg2}), $K = C_{01} + C_{12} - C_{02}$ with $Q_0 =1$ is easily transformed into the combination of the standard LGI term free of the invasive measurement and a signaling term as follows by firstly replacing the measured probabilities with the inferred ones and then by inserting into $K$:
\begin{align}
 \begin{split}
 \label{eq_4_5389_4bf2}
K_A  \equiv & \sum_{ i^A_1  = 1}^3  \sum_{ i_{1,1}   =1}^3  \sum_{ \hat{i}_2  = 1}^4  \big(Q_{1,  i_{1,1} }  \, + \, Q_{1,  i_{1,1} } Q_{2,  \hat{i}_2 } \, - Q_{2,  \hat{i}_2 } \, \big) D_{  i_{1,1},  i^A_1 } \, p^A_{1,2}(i^A_1, \hat{i}_2)  \,
 - \, \sum_{\hat{i}_2 =1}^4 Q_{2, \hat{i}_2} \Delta_S(\hat{i}_2)
  \end{split}
\end{align}
where the first term is the standard LGI definition with inferred probabilities and the second term includes inferred signaling terms $\Delta_S(\hat{i}_2)$ between the first and second planes for the measurement $\left[ \mathbf{O}_{\hat{i}_2}(t_2) \right]$ for each $\hat{i}_2$.   It is modeled as a signaling quantifier showing the influence of the measurement at time $t_1$ to the measurement at time $t_2$ and is defined  by utilizing ambiguous measurements as~follows:

\begin{eqnarray}
\label{eq_4_b4c0_4ada}
 \Delta_S(\hat{i}_2) \, \equiv   \, p_{2}(\hat{i}_2) \, - \, \sum_{i_{1,1}=1}^3 \widehat{p}_{1,2}(i_{1,1}, \hat{i}_2) \,
 & = &  \,  p_{2}(\hat{i}_2) \, - \, \sum_{i_{1,1}=1}^3  \sum_{i^A_1 = 1}^{3}  D_{i_{1,1}, i^A_1} \, p_{1,2}^A(i^A_1, \hat{i}_2)   \, \\
 \label{eq_4_b4c0_4ada_2}
& = &   \,  p_{2}(\hat{i}_2) \, - \, \frac{1}{2} \sum_{i^A_1 = 1}^{3}   p^A_{1,2}(i^A_1, \hat{i}_2)
\end{eqnarray}
where $\big(\sum_{ i_{1,1}  = 1}^3 D_{i_{1,1}, i^A_1} \big) = 1 \, / \, 2$.  Therefore, the no-signaling-in-time (NSIT) condition for the definition with ambiguous  measurement of $\Delta_S( \hat{i}_2 ) = 0$ is expected to convince a macrorealist about the reliability of the measurement setup.   Then,  the violation   of LGI free of the invasive measurement  becomes the~following:
\begin{equation}
 \label{eq_0_e90f_4103}
K_A \leq K_V \equiv 1 + \sum_{ \hat{i}_2  = 1}^4 \vert \Delta_S( \hat{i}_2 ) \vert
\end{equation}
where the summation at the right side of the inequality, i.e., $K_V-1$, shows the invasiveness of the measurements and the signaling level.  Therefore, measured values $p^A_{1,2}(i^A_1, \hat{i}_2 )$ are utilized to check the violation compatible with respect to the objections of a macrorealist. If  Equations (\ref{eq_a_8c5c_4984}), (\ref{eq_0_fe97_4605}),     and (\ref{eq_f_563d_44a7}) are inserted into Equations (\ref{eq_4_5389_4bf2}), (\ref{eq_4_b4c0_4ada}), and  (\ref{eq_0_e90f_4103}), then the following is obtained:
 \begin{eqnarray}
 \label{eq_e_82f5_4bfc}
 {K_A}  & = & G_1 \, \sum_{ \hat{i}_2 = 1}^{3}  Q_{2,  \hat{i}_2} f_{2}(X_{2, \hat{i}_2}) \, - \, G_1 \, Q_{2,4} \sum_{ \hat{i}_2 = 1}^{3} f_{2}(X_{2, \hat{i}_2})
  \,  + \, G_1 \, \sum_{i_{1,1} = 1}^{3} Q_{1, i_{1,1}} \sum_{ \hat{i}_2 = 1}^{3} Q_{2,  \hat{i}_2} f_{1,2}(X_{1, i_{1,1}}, X_{2,   \hat{i}_2 }, \vec{l}_{i_{1,1}})    \nonumber  \\
  & & - \,G_1\,  Q_{2,4}  \sum_{ i_{1,1}  =1}^{3} Q_{1,  i_{1,1} }   \sum_{  \hat{i}_2  = 1}^{3} f_{1,2}(X_{1,  i_{1,1} }, X_{2,  \hat{i}_2 }, \vec{l}_{ i_{1,1}})
    \, + \, G_2 \, (1 + Q_{2,4} )\, \sum_{ i_{1,1}  =1}^{3} Q_{1, i_{1,1} }  f_{1}(X_{1, i_{1,1} }, \vec{l}_{ i_{1,1}})      \nonumber  \\
  & & - \, G_2 \,  Q_{2,4} \sum_{ i_{1,1} =1}^{3}  f_{1}(X_{1,  i_{1,1} }, \vec{l}_{ i_{1,1}}) \, - \,  \, G_2 \,Q_{2,4}  \,  f_{T} \,\,\,\,\,\,\, \,\,\,\,  \\
\label{eq_a_1595_4ac6}
 K_V & = &  1 + \sum_{  \hat{i}_2 =1}^{3} \vert f_{V}(X_{2,  \hat{i}_2 }) \vert + \big\vert f_{T} \, G_2 \, - \sum_{  \hat{i}_2 = 1}^{3} f_{V}(X_{2,   \hat{i}_2 }) \big \vert
\end{eqnarray}
where  $\vec{l}_1 = [ 1 \, \, 1 \, \, -1]^T$;   $\vec{l}_2 = [ 1 \, \, -1 \, \, 1]^T$;  $\vec{l}_3 = [ -1 \, \, 1 \, \, 1]^T$; the functions $f_{1}(.)$, $f_{2}(.)$, $f_{1,2}(.)$, and $f_{V}(.)$; and  the variables $f_{T}$,   $G_1$, and $G_2$  are defined in Section \ref{sec4}. Next, quantum interference  among the paths, i.e., QPI,  is defined for the MPD setup.

\subsection{Modeling of QPI}
\label{sec2-5}

Double-slit interference gives a clear indication of  quantumness showing wave-particle duality and  spatial interference  as emphasized by Feynman. MPD setup presents the complementary phenomenon  of the temporal  interference  among the paths   which cannot be explained in any classical way showing that paths interfere in time, destructively and constructively decreasing and increasing the probability of the consecutive events, respectively. A \textit{gedanken} experiment shown in Figure \ref{Fig3} is designed with three planes. The target is to analyze interference effects of opening both  slits on the first plane in terms of the probability of the  light  to  diffract  through first (PL-1), second (PL-2), and third (PL-3) planes. History states at times $t_1$, $t_2$, and $t_3$ with three types of projections indexed by the superscripts $a$, $b$, and $c$ are defined  with the setups shown in Figure \ref{Fig3}a--c, respectively,  as follows:
\begin{eqnarray}
  \vert \Psi_{1}^a )  \,  & \equiv & \,       \mathbf{P}_{s} \odot      \left[  \rho_0 \right]     \\
 \vert \Psi_{2}^a )    \, & \equiv & \,        \left[  \mathbf{P}_{2,1} \right]  \odot \mathbf{P}_{s} \odot      \left[  \rho_0 \right]  ;  \\
 \vert \Psi_{3}^a )  \,  & \equiv & \,     \left[  \mathbf{P}_{3,1} \right]  \odot \left[  \mathbf{P}_{2,1} \right]  \odot \mathbf{P}_{s} \odot      \left[  \rho_0 \right]    \\
\vert \Psi_{1}^b )    \,  &\equiv & \,        \left[ \mathbf{P}_{1,1} \right]    \odot      \left[  \rho_0 \right] \\
\vert \Psi_{2}^b )   \, & \equiv & \,         \left[  \mathbf{P}_{2,1} \right]  \odot  \left[ \mathbf{P}_{1,1} \right]    \odot      \left[  \rho_0 \right] \\
\vert \Psi_{3}^b )   \, & \equiv &  \,    \left[  \mathbf{P}_{3,1} \right]  \odot \left[  \mathbf{P}_{2,1} \right]  \odot  \left[ \mathbf{P}_{1,1} \right]    \odot      \left[  \rho_0 \right] \\
\vert \Psi_{1}^c )   \,&  \equiv &  \,         \left[ \mathbf{P}_{1,2} \right]   \odot      \left[  \rho_0 \right] \\
\vert \Psi_{2}^c )    \, & \equiv & \,        \left[  \mathbf{P}_{2,1} \right]  \odot  \left[ \mathbf{P}_{1,2} \right]   \odot      \left[  \rho_0 \right]  \\
 \vert \Psi_{3}^c )   \, & \equiv & \,      \left[  \mathbf{P}_{3,1} \right]  \odot \left[  \mathbf{P}_{2,1} \right]  \odot  \left[ \mathbf{P}_{1,2} \right]   \odot      \left[  \rho_0 \right]
\end{eqnarray}
where superposition event at time $t_1$ is defined as $\mathbf{P}_{s} \equiv  \left[ \mathbf{P}_{1,1} \right] + \left[ \mathbf{P}_{1,2} \right]  $ and the  event probabilities are defined as follows:
\begin{eqnarray}
&& \hspace{0.40in}  p_{1}(\lbrace 1, 2\rbrace)    \,  \equiv  \,     (\Psi_{1}^a  \vert \Psi_{1}^a ); \hspace{0.63in}
p_{1}( 1 )    \,  \equiv  \,      (\Psi_{1}^b  \vert \Psi_{1}^b ) ; \hspace{0.58in}
p_{1}(  2 )   \,  \equiv  \,    (\Psi_{1}^c  \vert \Psi_{1}^c ) \\
&& \hspace{0.2in} p_{1,2}(\lbrace 1, 2\rbrace, 1)     \,  \equiv  \,      (\Psi_{2}^a  \vert \Psi_{2}^a ) ; \hspace{0.4in}
p_{1,2}( 1, 1 )     \,  \equiv  \,    (\Psi_{2}^b  \vert \Psi_{2}^b )      ; \hspace{0.4in}
p_{1,2}(  2,1  )     \,  \equiv  \,       (\Psi_{2}^c  \vert \Psi_{2}^c )  \\
&&  p_{1,2,3}(\lbrace 1, 2\rbrace, 1, 1)    \,  \equiv  \,   (\Psi_{3}^a  \vert \Psi_{3}^a ); \hspace{0.2in}
p_{1,2,3}( 1, 1, 1 )     \,  \equiv  \,       (\Psi_{3}^b  \vert \Psi_{3}^b ) ; \hspace{0.2in}
p_{1,2,3}(  2,1, 1  )  \,  \equiv  \,      (\Psi_{3}^c  \vert \Psi_{3}^c )
\end{eqnarray}
It is observed that $p_{1}(\lbrace 1, 2\rbrace) = p_{1}(1) + p_{1}(2)$ while interference  exists among consecutive planes.  The targeted scenario   for the relation between $\vert\Psi_{3}^a)$, i.e., the superposition of $\vert\Psi_{3}^b)$ and  $\vert\Psi_{3}^c)$, and $\vert\Psi_{3}^c)$     is as~follows:
 \begin{eqnarray}
\label{eq_4_af9b_4bd3}
 p_{1}(\lbrace 1, 2\rbrace)   & \, \stackrel{\mbox{1}}{>}   \, &   p_{1}(2 )   \\
 \label{eq_4_af9b_4bd3_b}
 p_{1,2}(\lbrace 1, 2\rbrace, 1)  & \, \stackrel{\mbox{2}}{>}\, &     p_{1,2}( 2, 1 ) \\
 \label{eq_4_af9b_4bd3_c}
 p_{1,2,3}(\lbrace 1, 2\rbrace, 1, 1)  & \,  \stackrel{\mbox{3}}{<}  \, &    p_{1,2,3}( 2, 1, 1 )
\end{eqnarray}
The superposition of $\vert \Psi_{3}^b )$ and  $\vert \Psi_{3}^c )$ on PL-1 increases the probability for the  light  to diffract  at time $t_1$ in $\stackrel{\mbox{1}}{>}$. In   $\stackrel{\mbox{2}}{>}$, the superposition constructively interferes to increase  also  the probability to  diffract  through the slit on PL-2 at time $t_2$.   However,   they destructively interfere in $ \stackrel{\mbox{3}}{<} $, decreasing the probability to  diffract  through the slit on PL-3 at time $t_3$.  Assuming starting with the setup in  Figure \ref{Fig3}c, with the second slit (the one with $X_{1,2} > 0$) open on PL-1,  if the first slit is additionally opened as shown in Figure \ref{Fig3}a, then the probability for the  light  to  diffract  through PL-1 and PL-2 increases while decreasing the probability to  diffract  through PL-3.   The counterintuitive observation based on a classical logic with balls passing through the slits is described as follows.  We open  the second slit  on PL-1 and observe that    the total number of balls passed through them for a statistical experiment  increases. It becomes more probable to  pass  through PL-1 with two slits in a classically logical manner. Furthermore, the probability to pass  through the single  slit  on PL-2   or the the total number of balls passing through PL-2  somehow  increases. However, we observe that the probability to  pass  through the single  slit  on the consecutive  PL-3   counterintuitively  decreases  in spite of the fact that more balls are coming from the second plane. This is complementary to the conventional spatial interference extensively studied in double-slit  interference  experiments.

\begin{figure}[H]
\begin{center}
\includegraphics[width=6.2in]{Fig3.eps}
\caption{ Setup for constructive and destructive interferences in time for the probabilities to  diffract  through each plane showing the history states (\textbf{a})  $\vert \Psi_{3}^a )  \equiv \left[  \mathbf{P}_{3,1} \right]  \odot \left[  \mathbf{P}_{2,1} \right]  \odot \big(\left[ \mathbf{P}_{1,1} \right] + \left[ \mathbf{P}_{1,2} \right] \big) \odot      \left[  \rho_0 \right]$  as the superposition of  $\vert \Psi_{3}^b )$ and  $\vert \Psi_{3}^c )$, (\textbf{b})  $\vert \Psi_{3}^b ) \equiv \left[  \mathbf{P}_{3,1} \right]  \odot \left[  \mathbf{P}_{2,1} \right]  \odot  \left[ \mathbf{P}_{1,1} \right]    \odot      \left[  \rho_0 \right] $, and (\textbf{c}) $\vert \Psi_{3}^c )  \equiv \left[  \mathbf{P}_{3,1} \right]  \odot \left[  \mathbf{P}_{2,1} \right]  \odot  \left[ \mathbf{P}_{1,2} \right]   \odot      \left[  \rho_0 \right]$.  The targeted scenario with \textit{classically counterintuitive} nature where a specific example of interference pattern (represented as the number of lambs denoting the number of photons for a practical counting experiment) for the cases of (\textbf{d}) two slits on PL-1 both open and (\textbf{e}) only the second slit open. The operation of closing the first slit decreases the number of photons diffracted through PL-2 while counterintuitively increases the number of photons through PL-3 since we classically expect a decrease. This scenario shows the interference of histories at two different time instants for PL-2 and PL-3 with firstly constructive and then destructive effects, respectively.}
\label{Fig3}
\end{center}
\end{figure}

QPI modeling requires the calculation of $\psi_{3, j}(x_3)$  in addition to $\psi_{2, j}(x_3)$ for $j \in [0, 1]$ as modeled in Equation (\ref{eq_e_f3b2_4ec9}). The explicit parameters required for the calculation of  both $\psi_{2, j}(x_3)$  and  $\psi_{3, j}(x_3)$ are presented in Section \ref{sec4}.

\subsection{Numerical Results}
\label{sec2-6}
Two different  calculation {scenarios} are performed denoted by $Sim_1$ and $Sim_2$ as shown in Table \ref{table1}.  $Sim_1$ calculates violation of LGI while $Sim_2$ shows an example of interference in time for the numerical analysis of QPI. Physical parameters are  monochromatic laser source wavelength $\lambda = 650$ (nm) as a widely available resource,  light velocity of $3 \, \times  \,10^{8}$ (m/s) in the z-direction,  and $\hbar = 1.05  \, \times \, 10^{-34}$ (J~$\times$~s) as Planck's constant.  The wavelength allows another degree of freedom to be adapted based on experimental design requirements or the targeted system design. {The layouts used in the simulations $Sim_1$ and $Sim_2$ are shown in Figure \ref{Fig4}a,b, respectively. Furthermore,  illustrative measurement setups for practically counting the number of photons compared with the emitted photons in unit time are presented in Figure  \ref{Fig4}c,d in order to calculate the probabilities $p_{1}(\lbrace 1, 2 \rbrace )$  and $p_{1, 2}(\lbrace 1, 3 \rbrace, 2 )$, respectively.}

\setlength\tabcolsep{2 pt}
\newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
\renewcommand{\arraystretch}{1.0}
\setlength\tabcolsep{1 pt}
\renewcommand{\arraystretch}{1.5}
\begin{table}[H]
\centering
\caption{Simulation and system parameters.}
\small
\begin{tabular}{m{1cm}m{2cm}m{3.8cm}m{1cm}m{2.5cm}m{2.7cm}}
\toprule
  {\bf ID} & {\bf Property} & {\bf Value} & {\bf ID} & {\bf Property} & {\bf Value} \\
\midrule
\multirow{6}{*}{$Sim_1$} &   $\overrightarrow{X}_{1}^T$     & \makecell[{{l}}]{$ D_s \,+ \, \left[ -\Delta x \, \, \, \, 0\, \, \, \,  \Delta x \right] \times \beta_1$ }  & \multirow{6}{*}{$Sim_2$} &   $\overrightarrow{X}_{1}^T$      & \makecell[{{l}}]{$ \left[ -4 \, \, \, \,  4 \right] \times \beta_1$ }    \\ \cline{2-3} \cline{5-6}
&  $\overrightarrow{X}_{2}^T$   & $\left[ -\Delta x \, \, \, \, 0\, \, \, \,  \Delta x \right] \times \beta_2$ &  & $X_{2,1}$  ($\upmu$m)  &  $\left[0,\, 500 \right]$  \\ \cline{2-3}  \cline{5-6}
& $\Delta x$;   $D_s$   & $ \lbrace 7, 11  \rbrace$;  $ \left[ 0,  3000  \right]$ ($\upmu$m)  & &  $X_{3,1}$  ($\upmu$m)   &  $\left[-600,\, 800\right]$ \\
\cline{2-3}  \cline{5-6}
& $t_{01} = t_{12}$ (ns) &  $ \lbrace 0.1, \,0.2  \rbrace $  &  &  $t_{01}, t_{12}, t_{23}$    (ns)  &  $0.5$,  $0.2$, $0.1$ \\
\cline{2-3}  \cline{5-6}
&  $\beta_1$, $\beta_2$  ($\upmu$m)  &   $\left[  1, \,  50  \right]$, $ \left[ 1,  \, 100  \right]$  &  &  $\beta_1$, $\beta_2$, $\beta_3$   ($\upmu$m)    &  $25$, $35$, $45$  \\
\cline{2-3}  \cline{5-6}
&  $\sigma_0$   ($\upmu$m)  &   $ \left[ 10, \, 800 \right]$   & &  $\sigma_0$   ($\upmu$m)  & $200$    \\
\bottomrule
\end{tabular}
\label{table1}
\end{table}
\renewcommand{\arraystretch}{1}
\setlength\tabcolsep{6 pt}
\vspace{-6pt}
\begin{figure}[H]
\centering
\includegraphics[width=4.5in]{Fig4.eps}
\caption{ { The layouts used in (\textbf{a}) $Sim_1$ and (\textbf{b}) $Sim_2$,  where for $Sim_2$, the fixed values of the parameters are $\sigma_0 = 200$ ($\upmu$m), $t_{01} = 0.5$ (ns), $t_{12} = 0.2$ (ns), $t_{23} = 0.1$ (ns), $\beta_1 = 25 $ ($\upmu$m), $\beta_2 = 35$  ($\upmu$m), and $\beta_3  = 45 $ ($\upmu$m) in addition to the fixed values of the slit positions on the first plane. The practical measurement setups to be utilized in future experiments are illustrated for the probabilities (\textbf{c}) $p_{1} (\lbrace 1, 2\rbrace)$ and (\textbf{d})~$p_{1,2} ( \lbrace 1, 3\rbrace, 2)$. The measurement planes count the detected number of photons compared with the number of photons emitted by the source in unit time.}}
\label{Fig4}
\end{figure}

\subsubsection{Violation of LGI}  There are two planes (PL-1 and PL-2) as shown in Figure \ref{Fig2} and \ref{Fig4}a. The slit positions on PL-1 and PL-2 are set to $D_s \, +  \, [ -\Delta x \, \, 0 \, \, \Delta x] \times \beta_1$ and  $[ -\Delta x \, \, 0 \, \, \Delta x] \times \beta_2$, respectively, where the positions on PL-1 are shifted with varying $D_s$ and the ratio of inter-slit distances to the slit widths is fixed in both planes.    $\Delta x $ is chosen larger than seven  to realize  independence of Gaussian slits, i.e., $ \mbox{exp}( (X_{1,i} - X_{1,i+1})^2 \, / \, 2 \, \beta_1^2) \ll 1$. The distance between the planes is set to $L$ such that the time duration is $t_{01} = t_{12}  \equiv L \, / \, c$. Gaussian source  parameter $\sigma_0$  is varied between  $10$ ($\upmu$m) and $800$ ($\upmu$m), compatible with standard laser resources including fiber lasers allowing smaller diameters reaching tens of micrometers.

The shift of the slits on PL-1 results in varying levels of  temporal correlation  for the diffraction through the slits on PL-2.  LGI violation ($K_A \, - \, K_V$)  is analyzed for varying $D_s$, $t_{01} = t_{02}$, $\beta_1$, $\beta_2$, and $\sigma_0$. In Figure \ref{Fig5}a, it is shown with the signaling  level  ($K_V \, - \, 1 $)  for varying $D_s$ for $t_{01} = 0.2$ (ns),  $t_{12} = 0.1$ (ns),  $\Delta x = 7$, $\beta_1 = 15$ ($\upmu$m), $\beta_2 = 30$ ($\upmu$m), and $\sigma_0 = 130$ ($\upmu$m).     The maximum violation is analyzed for different values of $Q_{1, i_{1,1}}$ and  $Q_{2, \hat{i}_2}$ for  $i_{1,1}  \in [1,3]$ and $ \hat{i}_2  \in [1,4]$, respectively, and the signs maximizing the violation are chosen for each $D_s$ shift. In Figure \ref{Fig5}b, distributions of the sign assignments  maximizing the violation are shown. Different setups realized with varying shift on PL-1 result in different optimized sign assignments for   maximum violation. Furthermore, violation decreases as the interplane slit distance increases, i.e., decreasing to zero with $K_A = K_V \approx 1$. It is observed that LGI is violated significantly, reaching close to  $0.3$  for the specific setup shown in Figure \ref{Fig5}a. The signaling $\sum_{k = 1}^4 \vert \Delta_S(k) \vert$ is close to zero, as shown in Figure \ref{Fig5}a for the marked region between the violation peaks. In the next simulations, the signaling is shown to decrease approximately to zero for varying setup values. It is shown  as a proof of concept  in Figure \ref{Fig6}a that the amount of signaling is  $K_V-1 \approx 2.3 \times 10^{-3}$ for a violation of $K_A - K_V \approx 0.2122$  while  satisfying NIM and signaling-in-time-related assumptions discussed in Reference~\cite{lg2} and utilizing a NIM-free violation.

\begin{figure}[H]
\centering
\includegraphics[width=6.3in]{Fig5.eps}
\caption{(\textbf{a})  LGI violation ($K_A \, - \, K_V$)   and signaling ($K_V \, - 1 \, $) for varying $D_s$, where  $t_{01} = 0.2$ (ns),  $t_{12} = 0.1$ (ns),  $\Delta x = 7$,  $\beta_1 = 15$ ($\upmu$m), $\beta_2 = 30$ ($\upmu$m), and $\sigma_0 = 130$ ($\upmu$m)  and (\textbf{b}) the corresponding   dichotomic  sign assignments for ambiguous measurements maximizing the violation for each $D_s$.}
\label{Fig5}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=5.8in]{Fig6.eps}
\caption{ (\textbf{a}) Maximum  LGI violation ($K_A - K_V$) and the corresponding amount of signaling ($K_V-1$) for varying $\sigma_0$, $\Delta x$, and $t_{01} = t_{12}$  and (\textbf{b}) the corresponding values of $\beta_1$, $\beta_2$, and $D_s$ maximizing the violation for each $\sigma_0$ assuming fully coherent sources.  Maximum violation for varying $(\beta_1, \beta_2)$ pairs for fully coherent sources where (\textbf{c}) $\Delta x = 7$ and $t_{01} = t_{12} = 0.1$ (ns) at the maximizing $\sigma_0 = 30$ ($\upmu$m), (\textbf{d})~$\Delta x = 7$ and $t_{01} = t_{12} = 0.2$ (ns) at $\sigma_0 = 230$ ($\upmu$m), and (\textbf{e}) $\Delta x = 11$ and $t_{01} = t_{12} = 0.1$ (ns) at $\sigma_0 = 150$ ($\upmu$m). It is observed that there is a large set of slit pairs and beam width resulting in LGI violation reaching $\approx 0.4$ for $\Delta x = 7$ and $\approx 0.23$ for $\Delta x = 11$, respectively, while there are local peaks for $(\beta_1, \beta_2)$ pairs for all cases. Increasing $t_{01}, t_{12}$ values expands the $(\beta_1, \beta_2)$ pairs for similar values of violations. (\textbf{f}) The comparison of the spatial coherence diameters  $D_c$  with the diffraction setup diameters  $D_1$ and $D_2$  for the first and second planes, respectively, where the targeted case is $\Delta x =11$ and $t_{01} = t_{12} = 0.1$ (ns), i.e., analyzed as the red curve in Figure \ref{Fig6}a, and (\textbf{g}) the corresponding LGI violation curve plotted again by emphasizing the coherence including the peak points.}
\label{Fig6}
\end{figure}

In Figure \ref{Fig6}, the effects of varying $\Delta x$, $\sigma_0$, and $t_{01} = t_{12}$ on  the violation of LGI  are shown for varying  $\beta_1$ and $\beta_2$ pairs.  $D_s$ and the signs of $Q_{1, i_{1,1}}$  and  $Q_{2, \hat{i}_2 }$ for  $ i_{1,1}  \in [1,3]$ and $ \hat{i}_2  \in [1,4]$, respectively, are chosen to maximize the violation for each pair and specific $\sigma_0$ value.   $\beta_1$ and $\beta_2$  are chosen in the sets $\lbrace 5, \, 10, \, 15, \hdots, 50 \rbrace$ ($\upmu$m) and $\lbrace 10, \, 20, \, 30, \hdots, 100 \rbrace$ ($\upmu$m), respectively, while choosing the maximum violation pairs. Similarly, $D_s$ is chosen  in the interval of $[1, 1000]$  ($\upmu$m) with the resolution of $1 \,$ ($\upmu$m).    There are important various observations for the specific simulation constraints in Table \ref{table1} which can be further improved by increasing the range of values and the resolution in simulations such as for the values of $\beta_1$, $\beta_2$, and $D_s$. However, the provided simple parameter set shows LGI violation under no-signaling conditions as a proof of concept. More specific observations provide more information about the nature of LGI violations for MPD setup as discussed next.

It is observed in Figure \ref{Fig6}a that violation becomes smaller as the relative distance between slits compared with  the slit width parameter  $\beta$ increases from $\Delta x = 7$ to $\Delta x = 11$ for $t_{01} = t_{12} = 0.1$ (ns). Furthermore,  $\sigma_0$ values maximizing the violation and the functional behavior with respect to varying $\sigma_0$ are approximately the same as $\Delta x$ changes for fixed $t_{01} = t_{12}$. The range of violation  reaches $\approx 0.4$ and $\approx 0.23$  for $\Delta x = 7$ and $\Delta x =11$, respectively, while decreases as $\Delta x$ increases. The signaling term is approximately zero for the $\Delta x = 11$ case with the maximum violation amplitude of $K_A - K_V \approx 0.2122$ as emphasized previously.  It is an open issue to design MPD setups maximizing the violation, e.g., similar to the boundary value of $0.5$ observed in conventional QM setups or $0.464$, as calculated in Reference~\cite{lg2} with inverted measurements for a three-level quantum optical system.   Optimum slit width maximizing violation for each $\sigma_0$ decreases as $\Delta x$ increases, as shown in Figure~\ref{Fig6}b. In other words, as $\Delta x$ increases, $\beta_1$ and $\beta_2$ are getting smaller to stabilize $\Delta x \times \beta_1$ and $\Delta x \times \beta_2$ for better violation. As a result, increasing the relative inter-slit distance results in more classical behavior  while decreasing violation.   On the other hand, increasing $\sigma_0$ further does not have any effect on the optimum $\beta_1$, $\beta_2$, and $D_s$ as the source behaves as a plane wave for the specific setup.  Violation amplitudes with respect to different $\beta_1$ and $\beta_2$ pairs for the maximizing $\sigma_0$ values  (extracted from Figure \ref{Fig6}a)  of  $30$ ($\upmu$m) and $150$ ($\upmu$m)  for $\Delta x = 7$ and $\Delta x  = 11$, respectively,  are shown in Figure \ref{Fig6}c,e, respectively, for $t_{01} = t_{12} = 0.1$ (ns).  There is a decrease  in both the range of $\beta_1$ and $\beta_2$ values and the maximum violation for $\Delta x = 11$. On the other hand, increasing interplane distance two times, i.e., making  $t_{01} = t_{12} = 0.2$  (ns), is observed not to change the maximum violation regime of  $0.4$  while  increasing  both the value of  $\sigma_0$ for the maximum violation to $230$ ($\upmu$m), as shown in Figure \ref{Fig6}a,  and the ($\beta_1$, $\beta_2$) values giving the maximum violation amplitudes, as shown in  Figure \ref{Fig6}d.  Increasing interplane distance improves the spread of the wave function on the consecutive plane while requires larger widths of source beam and slits in order to have similar violation amplitudes.

The results in Figure \ref{Fig6} assume fully coherent source both temporally and spatially. The realistic modeling of the temporal and spatial coherence of the light source is presented in Section {\ref{sec4-2}}. In our simulations, the total duration that light propagates is smaller than $1$ ns, i.e., corresponding to  $\approx 10^{-9} \times 3 \times 10^8 = 30$ (cm), which is much smaller than the temporal coherence time of the conventional single-mode lasers, i.e., $\Delta t > 10^{-6}$ (s) for single-mode  fiber lasers with $\Delta f$  of a few KHz~\cite{meng2006stable}.  On the other hand, the analysis of the spatial coherence is realized by defining the setup diameters $D_1$ and $D_2$ on the first and second planes, respectively, for the areas covering the slits. Then, these are compared with the spatial coherence diameter defined in Equation (\ref{spatialcoh}) depending on the duration of the propagation $t$ from the source to the diffraction plane and on the standard deviation of the Gaussian source denoted by $\sigma_0$. The detailed modeling and discussion are provided in Section {\ref{sec4-2}}, where the spatial coherence diameters $D_c$ and diffraction setup diameters $D_1$ and $D_2$  are described and it is targeted that $D_c$ is larger than both $D_1$ and $D_2$ as described in Equations (\ref{spatialcoh2}) and (\ref{spatialcoh3}). Then,  $D_{c}(t_{01}, \sigma_0)$ corresponds to the propagation from the source to the first plane and $D_{c}(t_{12}, \beta_1)$ is for the propagation from the first plane to the second plane. {It is assumed that the Gaussian slit modifies the diffracted wave with $\beta_1$ as the new source standard deviation while leaving the analysis with respect to the parameters of the wave function in Equation (\ref{eq_3_6538_4339}), e.g., $A_{1}$, as an open issue, as discussed in Section~\ref{sec4-2}.} Assume that the analysis providing the minimum signaling with $\Delta x = 11$ and $t_{01} = t_{12} = 0.1$ (ns) is targeted for coherence. In Figure \ref{Fig6}f, the comparisons of $D_{c}(t_{01}, \sigma_0)$ vs. $D_1$ and $D_{c}(t_{12}, \beta_1)$ vs. $D_2$ are shown. It is observed that the spatial coherence diameter covers the simulation parameters for both the peak and no-signaling cases. The LGI violation curve for $\Delta x = 11$ in Figure \ref{Fig6}a is plotted again in Figure \ref{Fig6}g by also emphasizing the violation amplitudes for the simulation parameters in Figure \ref{Fig6}f. It is an open issue to design MPD setups compatible with the coherence properties of the practical sources while satisfying quantum properties including the violation of LGI.

\subsubsection{QPI Analysis} \label{num2}

The three-plane setup (PL-1, PL-2, and PL-3) shown in Figure \ref{Fig3} is numerically analyzed for fixed  values of   $\beta_1 = 25$ ($\upmu$m), $\beta_2 = 35$ ($\upmu$m), $\beta_3 = 45$ ($\upmu$m), $\sigma_0 = 200$ ($\upmu$m), $t_{01} =  0.5$ (ns), $t_{12}  \, = \, 0.2$ (ns), $ t_{23}  \, =  \,0.1$ (ns), and $\overrightarrow{X}_{1}^T = \left[ -4 \, \, \, \,  4 \right] \times \beta_1$.  Sampling value of  $T_s = 1 $ ($\upmu$m)  is utilized in the analysis.  Constructive and destructive interferences  of two different paths at times $t_{2}$ and $t_{3}$, respectively, are performed by designing the slit positions with respect to  spatial  constructive and destructive interferences on PL-2 and PL-3,   respectively. In Figure \ref{Fig7}a, single slit position for PL-2, i.e., $X_{2,1}$, is chosen on the constructive interference regions, where $\vert \psi_{2,0}(x_2) \, + \, \psi_{2,1}(x_2) \vert$ is larger than $\vert \psi_{2,1}(x_2) \vert$  due to the superposition. Then, for each constructive $X_{2,1}$ slit position, the destructive interference regions on PL-3 such that $\vert \psi_{3,0}(x_3) \, + \, \psi_{3,1}(x_3) \vert$ is smaller than   $\vert\psi_{3,1}(x_3)\vert$  are searched  while the magnitudes $\max\limits_{x_3} \big \lbrace \vert \psi_{3,0}(x_3) \, + \, \psi_{3,1}(x_3) \vert^2 -  \vert\psi_{3,1}(x_3) \vert^2  \big \rbrace$ and the corresponding $X_{3,1} = x_3$   are shown in Figure \ref{Fig7}b,c, respectively.   It is observed that $X_{2,1} \approx 140 \, \upmu$m maximizes the destructive interference while the corresponding wave function amplitude on PL-3  is shown  in Figure \ref{Fig7}d, showing the the maximum destructive interference at $X_{3,1} \approx 143 \, \upmu$m.    Then, for each  $(X_{2,1}, X_{3,1})$  pair,  probability amplitudes  of the histories are shown in Figure \ref{Fig7}e with the marked areas where constructive interference occurs on PL-2 but with the destructive interference obtained on PL-3.  The conditions  with counterintuitive nature  in Equations (\ref{eq_4_af9b_4bd3})--(\ref{eq_4_af9b_4bd3_c}) for interference in time are satisfied.  The top two pairs of curves satisfy $p_{1}(\lbrace 1, 2\rbrace)     \, >    \,   p_{1}(2 )$ and  $p_{1,2}(\lbrace 1, 2\rbrace, 1)    \, > \,       p_{1,2}( 2, 1 ) $  due to the constructive interference while the bottom pair of the curves satisfies $p_{1,2,3}(\lbrace 1, 2\rbrace, 1, 1)  \, <   \,      p_{1,2,3}( 2, 1, 1 )$. In other words, the probability for the light to diffract through the second plane is decreased after blocking the first slit on PL-1 which counterintuitively results in an increase for the diffraction probability through the third plane. As a result, the proposed design of the setup  and utilization of spatial interference result in  interference of the quantum paths in time  for the projection histories with classically counterintuitive probabilistic~results.

Besides that, similar to the analysis of spatial and temporal coherence properties of the violation of LGI, the parameter set resulting in the maximum constructive  and destructive interferences on PL-2 and PL-3, respectively, is analyzed for compatibility with the coherence of practical sources. The inequalities defined in Equations (\ref{spacohinterference1})--(\ref{spacohinterference3}) in Section {\ref{sec4-2}} are targeted to be satisfied where the spatial coherence diameter $D_c$ and diffraction setup diameters $D_1$, $D_2$, and $D_3$ are described. The inequalities are calculated for the defined group of $\sigma_0$, $\beta_1$, $\beta_2$, $\beta_3$, $t_{01}$, $t_{12}$,  and $t_{13}$, where  $X_{3,1}$ is given in Figure \ref{Fig7}c and    $X_{2,1} \in [140, {170}] \, $ ($\upmu$m) is targeted. It is found that the practical coherence properties are satisfied as shown in Figure \ref{Fig7}f comparing   $D_c(t_{12}, \beta_1)$ with  $D_2$  and $D_c(t_{23}, \beta_2)$  with $D_3$ where $D_c(t_{01}, \sigma_0) \approx 607$~($\upmu$m) $ > D_1 \approx 270 $ ($\upmu$m).

\section{Discussion and Conclusions}
\label{sec3}

In this article, two novel resources for quantum technologies, i.e., QPE and QPI, are introduced based on tensor product structure of quantum history states for the simple linear optical setup of MPD. Operator theory modeling is presented by combining conventional history-based approaches of Griffith~\cite{griffiths1984consistent,griffiths1993consistent,griffiths2003consistent} and entangled histories framework of Reference~\cite{cotler2016} with FPI modeling. The inherent Feynman path generation mechanism of MPD setup is exploited for realizing quantum trajectories and for studying quantum temporal correlations.

\begin{figure}[H]
\centering
\includegraphics[width=6.2in]{Fig7.eps}
\caption{(\textbf{a}) $\vert \psi_{2,0}(x_2) \, + \, \psi_{2,1}(x_2) \vert^2$ compared with $ \vert\psi_{2,1}(x_2) \vert^2 $  and  $ \vert\psi_{2,0}(x_2) \vert^2 $ for diffraction through the layer PL-2, (\textbf{b}) $\max\limits_{x_3} \big \lbrace \vert \psi_{3,0}(x_3) \, + \, \psi_{3,1}(x_3) \vert^2 -  \vert\psi_{3,1}(x_3) \vert^2  \big \rbrace$ for varying $X_{2,1}$ on PL-3 such that destructive interference is maximized for each $X_{2,1}$ with respect to $x_3$ while $X_{2,1} \approx 140 \, \upmu$m maximizes the destructive interference, (\textbf{c}) $X_{3,1}$ maximizing the destructive interference for varying $X_{2,1}$, (\textbf{d})~the comparison of  $\vert \psi_{3,0}(x_3) \, + \, \psi_{3,1}(x_3) \vert^2$ and $\vert\psi_{3,1}(x_3) \vert^2$  on PL-3 for specific $X_{2,1} \approx 140 \, \upmu$m showing the destructive interference maximized with $X_{3,1} \approx 143 \, \upmu$m,   and (\textbf{e}) the marked regions satisfy the counterintuitive scenario in  (\ref{eq_4_af9b_4bd3})--(\ref{eq_4_af9b_4bd3_c}) for varying $X_{2,1}$  with the corresponding  $X_{3,1}$ pair in Figure~\ref{Fig7}c. Constructive and destructive interferences are observed for diffraction through PL-2 and PL-3, respectively,  with different kinds of correlation of the paths at different times as a proof-of-concept numerical simulation of  \textit{quantum path interference (QPI) in time} between the two paths.   (\textbf{f}) The comparison of setup diameters on the second and third planes, i.e., $D_2$ and $D_3$, respectively, with the spatial coherence diameters $D_c(t_{12},\beta_1)$ and $D_c(t_{23},\beta_2)$, respectively, in the targeted range of $X_{2,1} \in [140, {170}]$  ($\upmu$m) in Figure \ref{Fig7}e.}
\label{Fig7}
\end{figure}

 Following the similar terminology of entangled quantum trajectories in previous formulations,  the temporal correlation among the quantum propagation paths of MPD design is denoted with QPE as a novel quantum resource.  The state of the light after diffraction through consecutive planes is represented as a superposition of different trajectories through the slits as the main definition of QPE.  Its two fundamental properties are theoretically and numerically analyzed: violation of LGI  as the temporal analog of Bell's inequality with the ambiguous form and no-signaling recently proposed by Emary~\cite{lg2} and QPI as a counterintuitive phenomenon defining the quantum interference of histories. LGIs are violated reaching $>0.2$ with a NIM-free formulation complementing the recent experimental implementation of the violation of LGI in Reference~\cite{zhang2018experimental} utilizing linear polarization degree of freedom of the classical light with a photon-counting simple MPD setup of classical light. Besides that, it is an open issue to further design novel structures allowing the implementation of specific QPE states in analogy with GHZ  states experimentally implemented for entangled trajectories~\cite{cotler2017}.

Furthermore, QPI is numerically analyzed for a scenario where the decrease in the number of photons diffracting through a plane counterintuitively results in an increase in the number of photons diffracting through the next plane due to interference between two quantum history or trajectory states. MPD design introducing QPE and QPI is providing a test bed to understand the nature of temporal correlations improving our understanding of quantum mechanics and to design novel structures exploiting history-based formulation for practical purposes of quantum technologies. The simplicity of both the setup and proposed theoretical modeling allows varying kinds of \textit{gedanken} experiments for implementing paradoxes emphasizing the quantum nature with counterintuitive scenarios.

MPD has significantly low experimental complexity, relying on light sources of finite spatial and temporal coherence combined with conventional intensity detection or photon counting by promising near-future implementations. The experimental implementation is further simplified based on the mature science of Fourier optics with recent proposals in Reference~\cite{gulbahar2019quantumfourier}. Experimental and simulation studies for single-plane diffraction architectures verify the proposed theoretical modeling based on Fresnel diffraction of Fourier optics~\cite{santos2018huygens} and FPIs~\cite{sinha2015superposition, sawant2014nonclassical, magana2016exotic}. Therefore, proposed MPD theory simply extends previous modeling to multiple planes with FPI approach. One of the experimental challenges is the design of Gaussian slit, while any slit mask can be represented in terms of Gaussian basis as discussed in Reference~\cite{gulbahar2019quantumfourier}. Furthermore, based on the results in this article, any slit mask is potentially expected to have unique MPD setup parameters resulting in both the violation of LGI and QPI as a future work both to theoretically design and to experimentally verify. Moreover, these experimental implementations require only the calculation of total photon counts on the total area of the specific slits. This simplifies the slit design by allowing various geometries and the experimental setup regarding the detectors related to the size, position,  precision, and noise~\cite{gagnon2014effects}. The overall number of photon counts is satisfied more easily compared with detailed representations of the interference pattern in spatial basis.

QPE as an entanglement resource promises future applications for  computing~\cite{bg1, gulbahar2019quantumfourier}, communications~\cite{gulbahar2019quantum}, and varying quantum technologies while enriching the conventional resources of quantum entanglement based on correlations among multiple spatial quantum units.  In addition, MPD-based formulation of QPE and QPI allows future studies emphasizing the importance of time in QM fundamentals such as  regarding entanglement in time~\cite{brukner2004, aharonov2009, nowakowski2018},  quantum cosmology, and gravity~\cite{gell1991alternative, omnes1987, isham1994}.

\section{Methods}
\label{sec4}
\vspace{-6pt}
\subsection{Parameters for FPI Modeling of the Violation of LGI}
\label{sec4-1}
LGI and wave functions are modeled with FPI modeling~\cite{bg1}, and the resulting parameters to be utilized in Equations (\ref{eq_3_6538_4339}), (\ref{eq_e_f3b2_4ec9}), (\ref{eq_e_82f5_4bfc}), and (\ref{eq_a_1595_4ac6}) are provided in Table \ref{table2}. The functions utilized in Equations (\ref{eq_e_82f5_4bfc}) and (\ref{eq_a_1595_4ac6}) are defined as follows,  where the variables $k_{i}$ for $i \in [1,11]$ defined in Table  \ref{table2} depend on the system setup parameters  $t_{01}$, $t_{12}$, $\beta_1$, $\beta_2$, and $\sigma_0$, where
$f_{1}(x, \vec{l}) \, \equiv   \,  e^{k_{10} x^2}+ \vec{l}^T \, \vec{e}_4  $,
 $f_{2}(y)    \, \equiv   \,    -2 \,\mathbf{1}_3^T \,\vec{g}_y   -  \sum_{i=1}^3 e_1(X_{1,i}, \,y)$,
$ f_{1,2}(x, y, \vec{l})   \, \equiv   \,    \vec{l}^T  \,\vec{g}_y   + \, e_1(x,y)   $,
 $f_{V}(y)    \, \equiv   \,  G_1  \,  \mathbf{1}_3^T \,\vec{g}_y    $,
$f_{T}    \, \equiv   \,     \mathbf{1}_3^T \,\vec{e}_4  $,
$c(x_1,x_2,y)    \, \equiv   \,      \cos \big(k_{1}  (x_1^2 \,- \,x_2^2) \, + \, k_{2}  \,(x_1 \, -  \,  x_2)\,y \big) $,
$e_1(x, y)     \, \equiv   \,  e^{-2 \, k_{3} \, y \,x\,-\,k_{4} \, y^2 \, + \,k_{5} \, x^2} $,
$e_2(x_1, x_2, y)    \, \equiv   \,    e^{k_{6} \, (x_1^2 \, +  \, x_2^2) \, -\,k_{3} \,(x_1 \, + \,x_2)\, y - \,k_{4} \, y^2+ k_{7}\, x_1 \,x_2 }     $,
$ e_4(x_1, x_2)   \, \equiv   \,     e^{ k_{8} \, (x_1^2 \, +  \, x_2^2) +  k_{9}\, x_1 \,x_2}  $,
$ \vec{e}_4  \, \equiv   \, \left[ e_4(X_{1,1}, X_{1,2}) \,\, e_4(X_{1,1}, X_{1,3}) \,\, e_4(X_{1,2}, X_{1,3})\right]^T $,
$g_{i_{1}, \widetilde{i}_{1}}(y)   \, \equiv   \,    c(X_{1, i_{1 }},X_{1,   \widetilde{i}_{1}}, \,y) \, e_2( X_{1,  i_{1 } }, X_{1,  \widetilde{i}_{1} }, \,y) $,
$\vec{g}_y \, \equiv   \, \left[ g_{1,2}(y)  \,\,g_{1,3}(y) \,\,g_{2,3}(y)  \right]^T $,
$G_1   \, \equiv   \,   G_2 \, \beta_1 \, \beta_2 \, m^2 \, \sigma_0 \,  \sqrt{(\beta_1^2+d_{t,\sigma})\, / \, k_{11}} $,
$G_2   \, \equiv   \,   1 \,/ \, \big( \sum_{i=1}^3 e^{- X_{1,i}^2 \, / \, (\beta_1^2+d_{t,\sigma} )}  \big) $,
$d_{t,\sigma} \, \equiv   \,   a_{t,\sigma} /(m^2 \sigma_0^2)$,  $\Xi_i \equiv  \beta_i^2 \, + \, \sigma_0^2$, $ \vartheta_t  \equiv m \, \sigma_0^2 \,+ \,\imath \, \hbar\,  t_{0,1}$,    and $\mathbf{1}_3 = \left[ 1 \,\, 1 \,\, 1\right]^T$.
Besides that, the explicit parameters required for the calculation of $\psi_{3, {i}}(x_3)$ for ${i \in [0, 1]}$, { i.e., two different paths through two different slits on the first plane and the single slit both on the second and the third planes}, are shown in Table \ref{table3} based on the iterative modeling in Reference~\cite{bg1}.

\subsection{Temporal and Spatial Coherence of the Light Sources}
\label{sec4-2}

The coherence of the light field is characterized both temporally and spatially~\cite{mandel1965coherence}. Its temporal coherence length is inversely proportional to the the full width half maximum (FWHM) of the emission peak in the wavelength spectrum as  follows~\cite{akcay2002estimation, mandel1965coherence}:
\begin{equation}
\Delta l \equiv \Delta  t \, c =  \sqrt{ \frac{2 \, \mbox{ln} 2}{\pi \, n_r}}\frac{\lambda^2}{\Delta \lambda}
\end{equation}
where $n_r$ is the refractive index of the medium and $\Delta \lambda$ is the wavelength domain representation of the spectral width. The interference fringes in a double-slit experiment performed with this source are observed if the path length difference between different paths through the slits is smaller than $\Delta l$ as described in Reference~\cite{mandel1965coherence}.

Spatial coherence depends on the size of the light field source where it is experimentally characterized with double-slit interference experiment as clearly described in Reference~\cite{mandel1965coherence}, as shown in Figure \ref{Fig8}a. The source is assumed to have a side length $\Delta s$, and $\Delta \theta$ is the angle from the slits separated with $D_c$ to the center of the source. The distance between the source and slit planes is  $L$. Then, spatial coherence distance $D_c$ satisfies the following relation for the interference fringes to be~observed:
\begin{equation}
\Delta \theta \, \Delta s \leq \lambda
\end{equation}

\setlength\tabcolsep{2 pt}
\newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
\renewcommand{\arraystretch}{2.0}
\begin{table}[H]
\caption{Parameters for modeling {LGI} and  path integrals  $\big(\psi_{2,i}(x_2)$ for  $i \in {[0, 2]:  \mbox{ three  paths} } \big)$.}
\centering
\scriptsize
\begin{tabular}{m{0.3cm}m{4.85cm}m{0.65cm}m{3.4cm}m{0.65cm}m{4.8cm}}
\toprule
&  {\bf Formula}  & & {\bf Formula}   & &  {\bf Formula} \\
\midrule
$k_{1}$  & $  \frac{  - \hbar\, m \,t_{12}  (a_{t,\sigma} \,+ \,\hbar^2  \, t_{01} \, t_{12} )   }{ 2 \, k_{11} }     $ &   $k_{8}$ & \makecell[{{l}}]{ $  - \frac{1}{4}   \,  \big( \frac{1}{\beta_1^2   }   \,+ \, \frac{1}{  \beta_1^2\,+\, d_{t,\sigma}}     \big)$} & $ \xi_{1}$ & $\frac{ \beta_1^2 \,m \,   \vartheta_t   }{ \big( \beta_1^2 \, m \, (m \,\sigma_0^2 \, + \, \imath \, \sqrt{b_t} ) \, + \,   \imath \,  \hbar \, t_{1,2} \, \vartheta_t  \big)}$   \\
\hline
$k_{2}$  & $  \hbar\, m^3 \,\sigma_0^2 t_{12} \left(\beta_1^2 \,+ \,  d_{t,\sigma} \right) \, / \, k_{11}  $ &  $k_{9}$ &  $  \frac{1}{2}  \,   \big( \frac{1}{\beta_1^2}  \,  - \, \frac{1}{\beta_1^2\,+\, d_{t,\sigma}}   \big)$  &  $A_{1}$  & $\frac{ - \beta_1^2 \, m^2 \left(\hbar^2  \, t_{0,1}^2 \, + \, m^2  \, \sigma_0^2 \,  \Xi_1 \right) }{( 2 \, \alpha_{t, \sigma, \beta} )  }$   \\
\hline
$k_{3}$  & $ \big(-\beta_1^2 m^2  ( a_{t,\sigma}   \,+ \,\hbar^2  \,t_{01}  \,t_{12} ) \big) \, / \,  k_{11}   $ &  $k_{10}$  &   \makecell[{{l}}]{$ - 1 \, / \, ( \beta_1^2\,+\,  d_{t,\sigma}  )   $} & $B_{1}$ & $\frac{ \big( \beta_1^4 \,m^3 \,\hbar\,  t_{0,1} \, + \,m  \, \hbar  \, t_{1,2}  (\hbar^2 \,t_{0,1}^2 \, + \, m^2 \,  \Xi_1^2 ) \big )}{ (2 \, \alpha_{t, \sigma, \beta})}$  \\
\hline
$k_{4}$  & $  \beta_1^2 \,m^4\, \sigma_0^2  (\beta_1^2 \,+ \, d_{t,\sigma}   ) \, / \,  k_{11}   $ &  $k_{11}$ &  \makecell[{{l}}]{$\beta_1^4 \,m^2\, \left(m^2 \,\sigma_0^2\,  \Xi_2 \,+\,
b_t \,\right) $ \\ $ + \,\beta_1^2 m^2 \left(\beta_2^2 \, a_{t,\sigma} \, + \,2 \,c_{t, \sigma}\right)$ \\ $+\,\hbar^2 \,t_{12}^2  \,a_{t,\sigma}  $}  & $\mathbf{H}_{R,1}$ & $\frac{- m^2 \big(\beta_1^2  (  b_{t}  \,+\,m^2 \sigma_0^4 )\,+\,  c_{t,\sigma}  \big)}{(2 \, \alpha_{t, \sigma, \beta})}$   \\
\hline
$k_{5}$ & $ \frac{- m^2 \big(\beta_1^2  (m^2 \, \sigma_0^2  \, \Xi_2  \,+ \,b_{t}) \,+ \,c_{t, \sigma}\big)}{ k_{11}}     $ &  $a_{t,\sigma}$, $b_{t}$, $c_{t, \sigma}$  &  \makecell[{{l}}]{ $\hbar^2  \,t_{01}^2 \,+ \,m^2 \sigma_0^4$, \\ $\hbar^2  \,(t_{01} \,+ \,t_{12})^2$, \\ $\hbar^2  \,\sigma_0^2  \,t_{12}^2$ }   &  $\mathbf{H}_{I,1}$ &  $\frac{  m \,\hbar \, t_{1,2}  \left(  a_{t, \sigma} \, + \, \hbar^2 \,t_{0,1} \, t_{1,2}     \right)}{(2 \, \alpha_{t, \sigma, \beta})}$   \\
\hline
$k_{6}$  & \makecell[{{l}}]{ $  -m^2 \big(2 \, \beta_1^2\, (m^2 \,\sigma_0^2   \,  \Xi_2  \, +  b_{t}) \big)  \, / \, (4 \,  k_{11})   $  \\ $  \, + \,  \, m^2  (- \beta_2^2  \,a_{t,\sigma}   \, - \, 2 \,  c_{t, \sigma} )  \, / \, (4 \,  k_{11})   $ } &  $\alpha_{t, \sigma, \beta}$ &   \makecell[{{l}}]{$\beta_1^4 m^2 \left(  b_{t}  \,+ \,m^2 \sigma_0^4\right)  \, $ \\ $ +\, 2 \, \beta_1^2\, m^2 \,  c_{t, \sigma}  $ \\ $+ \,\hbar^2  \, t_{1,2}^2  \,  a_{t,\sigma} $}  & $c_{1}$ &   \makecell[{{l}}]{$  \beta_1^2  \,m^2 \left(  a_{t, \sigma} \, + \, \hbar^2 \,t_{0,1} \, t_{1,2}     \right) \, / \,  \alpha_{t, \sigma, \beta}  $} \\
\hline
$k_{7}$ & $  \beta_2^2 \, m^2  \, a_{t,\sigma} \, / \, (2 \,  k_{11})  $ &   $\chi_0$ & \makecell[{{l}}]{$ \pi^{-1/4}  \sqrt{ \frac{m \,\sigma_0}{ m \, \sigma_0^2 \, + \,\imath \,\hbar\,  t_{0,1} }    } $}  & $d_{1}$ & $\frac{- m \,\hbar \, t_{1,2} \left(\hbar^2 t_{0,1}^2 \, + \, m^2 \, \sigma_0^2 \,  \Xi_1 \right)}{\alpha_{t, \sigma, \beta}  }$ \\
\bottomrule
\end{tabular}
\label{table2}
\end{table}

\setlength\tabcolsep{2 pt}
\newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
\renewcommand{\arraystretch}{2.0}
\begin{table}[H]
\centering
\caption{ Parameters for modeling the path integrals of QPI $ \big( \psi_{3,i}(x_3)$ for   $i \in {[0, 1]: \mbox{ two  paths} } \big)$.}
\scalebox{0.92}[0.92]{
\scriptsize
\begin{tabular}{m{0.58cm}m{5.23cm}m{0.6cm}m{2.5cm}m{1.5cm}m{5.0cm}}
\toprule
&  {\bf Formula}  & & {\bf Formula}  &  \boldmath{\makecell[{{c}}]{$j  \in [1, 2]$}} & {\bf Formula}     \\
\midrule
\multirow{2}{*}{$\mathbf{H}_{2}$}  & \multirow{2}{*}{$\left(
\begin{array}{cc}
 {\nu}_{2,2} \left(\zeta_{1,c} \, + \imath \, \zeta_{1,d}\right)^2 \, + \, {\nu}_{1,1} & 0 \\
 {\nu}_{3,2} \left(\zeta_{1,c} \, + \imath \, \zeta_{1,d}\right) & {\nu}_{1,2} \\
\end{array}
\right) $ }   & ${\nu}_{2,2}$   &  $   - \frac{\beta_2^2 \, \hbar \, t_{2, 3}}{2 \, \imath \, \varsigma_2}   $   &  ${\nu}_{1, j}$ &  $- \frac{ 2 \, \hbar \, t_{j, j+1} (A_{j-1} \, + \,\imath \, B_{j-1}) \, +  \, \imath \, m  }{ 2 \, \imath \, \varsigma_j }$     \\
\cline{3-6}
  &  &    ${\nu}_{3,2}$    &  $ -\frac{  \hbar \, t_{2, 3}}{ \imath \, \varsigma_2 }  $ &  $\zeta_j$ & \makecell[{{l}}]{ $   4 \, B_{j-1} \, \beta_j^4 \, \hbar \, m \, t_{j, j+1} \, + \, \beta_j^4 \, m^2 \, $ \\ $+ \, \hbar^2  \, t_{j, j+1}^2  \, \varrho_{j}$} \\
\hline
\multirow{2}{*}{ \makecell[{{l}}]{$\overrightarrow{c}_{2}$, \\ $\overrightarrow{d}_{2}$} }  & \multirow{2}{*}{ \makecell[{{l}}]{ $ \left( \begin{array}{cc} \, {\nu}_{4,2} \, \zeta _{1,c} \, + \, {\nu}_{5,2}  \, \zeta _{1,d} \\
\zeta _{2,c} \end{array} \right)$, \\ $ \left( \begin{array}{cc} {\nu}_{4,2} \, \zeta _{1,d} \, - \, {\nu}_{5,2} \, \zeta _{1,c}\\
\zeta _{2,d}\end{array} \right)$}} & ${\nu}_{4,2}$   &   $ \beta_2^2 \, \zeta_{2,c}$ &  $\zeta_{j,c}$  & $  (2 \, B_{j-1}  \, \hbar  \, m  \, t_{j, j+1} \,\beta_j^2      + \beta_j^2 \,  m^2) \, / \, \zeta_j$  \\
\cline{3-6}
  &   &    ${\nu}_{5,2}$  & $ - \frac{ 2 \,\hbar \, t_{2, 3}  \,  A_2 }{m}   $   &  $\zeta_{j,d}$ &   $   \hbar \,  m  \, t_{j, j+1} \,  \left(2  \, A_{j-1} \,  \beta_j^2 \, - \, 1\right) \, / \, \zeta_j$    \\
\hline
$A_0$  & $ - m^2 \sigma_0^2 \, / \, (2\, \hbar^2 t_{0,1}^2\,+\,2\,m^2 \sigma_0^4)$    &      $A_2$  &  $\frac{ \beta_2^2 \, m^2 \,  \left(2 \, A_{1} \, \beta_2^2 \, -\, 1\right)}{ 2 \, \zeta_2 }$   & \makecell[{{l}}]{$\varrho_j$, \\$\xi_{j}$} & \makecell[{{l}}]{$   4  \, \beta_j^4    \left(A_{j-1}^2   +   B_{j-1}^2 \right) \, - \, 4 \, A_{j-1}   \beta_j^2   +   1$,\\ $  \beta_j^2 \, m \, / \, \varsigma_j$} \\
\hline
 $B_0$ & $ \hbar \, m \,t_{0,1} \,/ \,(2 \,\hbar^2 t_{0,1}^2 \, + \, 2 m^2 \sigma_0^4) $ &   $B_2$ & $\frac{ 2 \, B_{1} \, \beta_2^4 \, m^2  +   \hbar \, m \, t_{2, 3} \, \varrho_2 }{2 \, \zeta_2}$    & $\varsigma_j$  & \makecell[{{l}}]{$  \hbar \, t_{j, j+1} \,  \big(2  \, \beta_j^2 \, (B_{j-1}   - \, \imath  \, A_{j-1}) \, + \, \imath \big)$ \\ $\, +\beta_j^2 \, m  $ }       \\
\bottomrule
\end{tabular}}
\label{table3}
\end{table}

In Reference~\cite{latychevskaia2017spatial}, more detailed calculation of the spatial coherence diameter of the 2D Gaussian intensity source $I(x_0, y_0) = I_0 \, \mbox{exp}\big( -(x_0^2 + y_0^2) \, /  \, (2 \, \sigma_0^2) \big)$ is achieved by using van Cittert--Zernike theorem~\cite{zernike1938concept}. It is  observed that an incoherent source of uniform intensity with circular diameter $\Delta s$ results in $D_c \propto \lambda \, L  \, / \, (\pi \, \Delta s)$ while the coherence diameter of an incoherent source with Gaussian-distributed
intensity and the standard deviation of a propagating coherent Gaussian source are both some multiple of  $  \lambda \, L \,/ \, (\pi \, \sigma_0)$. Therefore, in this article, it is assumed that $D_c$ is given by the beam width of the propagating Gaussian wave, which is approximated as $2 \, \sqrt{2} \, \sigma_D$ for the far field, as shown in Figure \ref{Fig8}b, where the intensity drops to $1/e^2$ at $D_c \, / \, 2$. The accuracy of the estimation is further improved to include near field due to the diverse parameter ranges of the LGI and time domain interference setups as follows:
\begin{equation}
\label{spatialcoh}
D_c(t, \sigma_0)  \approx 2 \, \sqrt{2} \, \sigma_D = \frac{2}{\sqrt{-A_0}}
\end{equation}
where $\sigma_D \equiv 1 \, / \, \sqrt{-2 \, A_0}$ { for $A_0 < 0$} since the free-space propagating Gaussian wave function on the detector plane $\Psi_{D}(x)$ is proportional to  $ \mbox{exp}( A_0 \,x^2 \, + \, \imath \, B_0 \, x^2 )$, $A_0 = -  m^2 \, \sigma_0^2  \, / \, (2 \, \hbar^2  \, L^2 \, / \, c^2 \, + \, 2 \,  m^2 \, \sigma_0^4)$, and $t \equiv L \, / \,  c$ based on trivial application of FPI kernel in Equation (\ref{eq_c_b338_4fde}). The slit positions on the first plane should be inside the coherence area for the proposed numerical simulations to be more compatible with the future experimental implementations.

In numerical simulations of $Sim_1$ for the violation of LGI with two planes of triple slits, the position of the furthest edge of a slit on PL-1 is found by $   D_s + \, \Delta x \, \beta_1 + \sqrt{2} \, \beta_1$ {for $D_s > 0$} by assuming that the Gaussian slit has a width of $2 \, \sqrt{2} \, \beta_1$ similar to the Gaussian source formulation. Then,  for the specific set of $(\sigma_0, D_s, \beta_1, \Delta x)$ shown in Figure  \ref{Fig6}a,b, the simulation results are reliable if the slits are in the spatial coherence area formulated as follows:
\begin{equation}
\label{spatialcoh2}
D_c(t_{01}, \sigma_0) \geq D_{1} \equiv 2 \, \big( D_s + \, (\Delta x  \, + \, \sqrt{2})\, \beta_1 \big)
\end{equation}
where $D_{1}$ corresponds to the diameter of the area on PL-1, where the slits reside as shown in Figure \ref{Fig8}c.

Similarly, the maximum difference between the positions of the slits on PL-1 and PL-2, i.e., the upper slit on PL-1 and the lower slit on PL-2, is calculated by using  $D_2  \equiv  2 \, \big( D_s + \, (\Delta x  \, + \, \sqrt{2})\, \beta_1 \big) +  2\, (\Delta x  \, + \, \sqrt{2})\, \beta_2$.
{ It is assumed that the projected light through the slit on the $(j-1)$th plane has spatial coherence larger than $D_c(t_{j-1,j}, \beta_{j-1})$ by assuming that the slit masking modifies the beam width to $\approx 2 \, \sqrt{2} \, \beta_{j-1}$.

\begin{figure}[H]
\centering
\includegraphics[width=5.0in]{Fig8.eps}
\caption{(\textbf{a})  The conventional modeling for the spatial coherence of light sources based on double-slit diffraction~\cite{mandel1965coherence}, where $\Delta \theta \, \Delta s \leq \lambda$ is required for the fringes to be observed determining the spatial coherence diameter ($D_c$); (\textbf{b}) free-space propagation of Gaussian beam, where $D_c$ is approximated as the $1/e^2$ intensity beamwidth of $2 \, \sqrt{2} \, \sigma_0$ with the standard deviation of $\sigma_D$. The descriptions of the calculation of the setup diameters on the planes to include the slits are denoted by $D_j$ for $j \in [1, 3]$ with respect to the location and the standard deviation of the source on the previous plane ($\sigma_0$ for the first plane and {$\beta_{j-1}$} for the $j$th plane) for (\textbf{c}) LGI violation numerical analysis $Sim_1$ with two planes of triple slits on each plane and (\textbf{d}) interference in time scenario $Sim_2$ with three planes. }
\label{Fig8}
\end{figure}

More accurate analysis left as an open issue requires formulation of the coherence for the wave function $\psi_{j,n}(x_j)$ for each $n$th path in Equation (\ref{eq_3_6538_4339}) by also considering the parameter $A_{j}$ for calculating the beam width on the $j$th plane  instead of utilizing $- m^2 \, \beta_{j-1}^2  \, / \, (2 \, \hbar^2  \, L^2 \, / \, c^2 \, + \, 2 \,  m^2 \, \beta_{j-1}^4)$ calculated by replacing $\sigma_0$ with $\beta_{j-1}$ in the expression of $A_0$. Then,} the validity of the simulation requires the following inequality between PL-1 and PL-2:
\begin{equation}
\label{spatialcoh3}
D_c(t_{12}, \beta_1) \geq D_{2}
\end{equation}
Similar formulations for $Sim_2$ with three planes as shown in Figure \ref{Fig8}{d} for the time domain interference result in the following requirements for the results in Figure \ref{Fig7} to be compatible with the sources practically not fully coherent:
\begin{eqnarray}
\label{spacohinterference1}
D_c(t_{01}, \sigma_0) & \geq & D_1 \equiv 2 \, \big(\Delta x  \, + \, \sqrt{2}\big)  \, \beta_1 \\
\label{spacohinterference2}
D_c(t_{12}, \beta_1) & \geq & D_2 \equiv 2 \,   \Delta x   \, \beta_1 \, + 2 \, X_{2,1} \, + \, 2 \, \sqrt{2} \, \beta_2 \\
\label{spacohinterference3}
D_c(t_{23}, \beta_2) & \geq & D_3 \equiv 2 \, \big( {\vert} X_{3,1} - X_{2,1}  {\vert} \, + \, \sqrt{2} \beta_3\big)
 \end{eqnarray}

\vspace{6pt}

\funding{This research was funded by Ozyegin University Research Grant.}

\abbreviations{
\noindent
\setstretch{0.2}
\begin{tabular}{@{}ll}

MPD  & Multiplane diffraction   \\
QC & Quantum computing\\
QM & Quantum mechanical \\
QPE & Quantum path entanglement\\
QPI & Quantum path interference  \\
FPI & Feynman's path integral \\
LGI & Leggett-Garg Inequality \\
MR &  Macroscopic realism \\
NIM & Non-invasive measurability \\
SIT & Signaling-in-time \\
GHZ & Greenberger-Horne-Zeilinger \\
FWHM & Full width half maximum\\

\end{tabular}}

\reftitle{References}

redpajama_set_name: RedPajamaArXiv