Young’s Experiment with Entangled Bipartite Systems: The Role of Underlying Quantum Velocity Fields

Abstract

We consider the concept of velocity fields, taken from Bohmian mechanics, to investigate the dynamical effects of entanglement in bipartite realizations of Young’s two-slit experiment. In particular, by comparing the behavior exhibited by factorizable two-slit states (cat-type state analogs in the position representation) with the dynamics exhibited by a continuous-variable Bell-type maximally entangled state, we find that, while the velocity fields associated with each particle in the separable scenario are well-defined and act separately on each subspace, in the entangled case there is a strong deformation in the total space that prevents this behavior. Consequently, the trajectories for each subsystem are not constrained any longer to remain confined within the corresponding subspace; rather, they exhibit seemingly wandering behavior across the total space. In this way, within the subspace associated with each particle (that is, when we trace over the other subsystem), not only interference features are washed out, but also the so-called Bohmian non-crossing rule (i.e., particle trajectories are allowed to get across the same point at the same time).

1. Introduction

The 2022 Nobel Prize in Physics was awarded to Alain Aspect, John F. Clauser, and Anton Zeilinger “for experiments with entangled photons, establishing the violation of Bell inequalities and pioneering quantum information science” [1]. All these ideas started in the mid-1930s as a heated debate about a “spooky action at a distance” [2], the separability of quantum systems [3], the impossibility of hidden variables in quantum mechanics [4], and, in sum, the problem of quantum measurement [5,6,7]. However, entanglement has become, together with coherence, the cornerstone of modern quantum information theory [8] and the quantum technologies that have emerged after the so-called second quantum revolution [9].

The analysis of entanglement and its properties is often based on rather abstract (algebraic) descriptions, where the fact that quantum systems go from somewhere to somewhere else in a given time is typically not considered. Young’s two-slit experiment is an ideal candidate to render some light on the phenomenology of how this happens, as it constitutes a paradigm in the investigation of properties related to quantum coherence, a key element strongly connected to entanglement. This experiment not only makes evident the reality and implications of the quantum phase through the appearance of interference fringes, but also the gradual or total suppression of the latter as soon as the system, diffracted by the slits, interacts with, and becomes entangled with, other surrounding systems. This coherence loss mechanism is what we call decoherence [10] (to be distinguished from the decrease of fringe visibility due to the presence of incoherence factors [11]), which leads to behaviors analogous to those observed in the ordinary (incoherent) classical world [12,13].

Among the several pictures of quantum mechanics, the results of the Bohmian one [14,15,16] are particularly appealing to investigate the above type of problem, as it puts the emphasis on the role of the phase in quantum dynamics [17,18]. Different questions related to entanglement have already been considered from this point of view in the literature [19,20,21,22], including those aimed at proving, or disproving, flaws in the approach [23,24,25,26,27,28,29,30]. Leaving aside interpretive issues, from the position that Bohmian trajectories are simply (quantum) flux streamlines describing the spatial dispersion of the probability density along time, we find that Bohmian mechanics becomes a suitable tool to study and understand the evolution of quantum systems in an alternative way to that of considering the global information supplied by the probability density. This relies on a well-defined velocity field, determined from the relationship between the quantum flux and the probability density (the so-called guidance condition in canonical Bohmian mechanics [16]). As was shown in [18,31], for instance, this velocity field is very convenient in the analysis of the dynamical consequences intrinsically associated with the concept of quantum coherence.

In this work, we focus on the relationship between entanglement and decoherence in investigating the continuous coherence flow between two spatially separated subsystems with the aid of both the (Bohmian) velocity field and the trajectories obtained from it. To this end, we consider a realization of Young’s two-slit experiment with two identical particles, X and Y, with a mass m, given that this paradigm has been used in the past to test experimentally different aspects of Bohmian mechanics [25,32,33]. Without any loss of generality, particle X is taken as the reference system, which describes the coherent superposition of the two diffracted beams in the experiment. As for Y, it is taken as the “environmental” partner, which leads to the suppression of the interference traits present in X when both are entangled. From the most elementary aspects of the theory of open quantum systems [34], it is known that such suppression is an effect of Y on the projection of the state vector for the full system $XY$ onto the subspace for X. In this paper, the objective is to determine the specific dynamical aspects that make the entangled state different with respect to the factorizable one, and which can be extrapolated to more general frameworks, regardless of whether the two subsystems are spatially separated or they are interacting [35].

Accordingly, the work has been organized as follows. In Section 2, we introduce and discuss basic theoretical aspects involved in the description of continuous variable systems, including a brief overview on both localized and delocalized states, described in terms of Gaussian wave packets, mimicking diffraction by Gaussian slits [36,37]. The results from numerical simulations illustrating the main dynamical features for both factorizable and entangled states are shown and discussed in Section 3. Finally, the main conclusions from this work are summarized in Section 4.

2. Theoretical Aspects

2.1. Basic Formal Grounds

To be self-contained, let us first introduce some general, basic considerations. Thus, let us consider two coherently emitting sources (e.g., two slits, two coupled fibers, etc.), with the transmitted waves traveling ahead faster than along the transverse direction, i.e., the direction along which we expect to detect the effects of coherence (interference) or to analyze their suppression (decoherence). This simplified phenomenological physical system has been considered previously with excellent results; for instance, in the analysis of grating diffraction [38]. Accordingly, if we consider two identical particles, X and Y, their descriptions rely on their transverse coordinates; henceforth, denoted by x and y, respectively. Additionally, to simplify the analysis, let us consider that X and Y are both nonrelativistic, spinless particles of mass m.

In a first instance, following the standard approach, the analysis of the dynamical implications of the coupling between X and Y requires a description of the problem based on a density matrix approach. If the joint state for X and Y is represented by the state vector $|\Psi (t)\rangle$, then the density matrix in the position representation (i.e., the space described by the transverse coordinates x and y) for the joint system is

$$\widehat{\rho }(x,y;x^{\prime },y^{\prime }|t)=\rho (x,y;x^{\prime },y^{\prime }|t)|x,y\rangle \langle x^{\prime },y^{\prime }|,$$

(1)

where

$$\rho (x,y;x^{\prime },y^{\prime }|t)=\langle x,y|\Psi \left(t\right)\rangle \langle \Psi \left(t\right)|x^{\prime },y^{\prime }\rangle ={\Psi }^{*}(x^{\prime },y^{\prime }|t)\Psi (x,y|t).$$

(2)

denotes the (density matrix) element in the position representation. The diagonal elements provide us with the usual probability density, while the off-diagonal elements provide an account of the preservation of the full-system ($XY$) quantum coherence. For simplicity, but without loss of generality, we always work under the hypothesis that the full system is pure and, hence, the (full-system) density matrix takes the usual form $\widehat{\rho }=|\Psi \left(t\right)\rangle \langle \Psi \left(t\right)|$, which gives rise to Equation (1) after considering the corresponding projection operators. This is justified by the fact that, here, we are investigating the joint dynamics. However, the current analysis can readily be generalized to mixed states.

To quantify the effects of the coupling for pure bipartite systems, it is also common to trace the density matrix over one of the subsystems and, then, to analyze the properties displayed by the reduced density matrix associated with the remaining subsystem [35]. Proceeding in this way, tracing over Y, for instance (the same description can be achieved the other way around, tracing over X instead of Y, without any detrimental effect regarding the conceptual aspects that are here investigated), leads to the reduced density matrix for X:

$$\begin{array}{ccc}\hfill {\widehat{\tilde{\rho }}}_X(x,x^{\prime }|t)={\mathrm{Tr}}_Y\left[\widehat{\rho }(x,y;x^{\prime },y^{\prime }|t)\right]& =& \int \int \left(\int \langle y^{″}|\rho (x,y,x^{\prime },y^{\prime })|x,y\rangle \langle x^{\prime },y^{\prime }|y^{″}\rangle d{y}^{″}\right)dyd{y}^{\prime }\hfill \\ & =& \int \int \rho (x,y;x^{\prime },y^{\prime }|t)\left(\int \langle y^{″}|x,y\rangle \langle x^{\prime },y^{\prime }|y^{″}\rangle d{y}^{″}\right)dyd{y}^{\prime }\hfill \\ & =& {\tilde{\rho }}_X(x,x^{\prime }|t)|x\rangle \langle x^{\prime }|,\hfill \end{array}$$

(3)

where

$${\tilde{\rho }}_X(x,x^{\prime }|t)=\int \rho (x,y;x^{\prime },y|t)dy$$

(4)

denotes the (density matrix) element for the subsystem X referred to its own subspace (described by the transverse coordinate x). If the two subsystems are uncorrelated (separable), the reduced density matrix element reduces to the trivial result

$${\tilde{\rho }}_X(x,x^{\prime }|t)\equiv \langle x|{\Psi }^X\left(t\right)\rangle \langle {\Psi }^X\left(t\right)|x^{\prime }\rangle ={\Psi }^{X,*}(x^{\prime },t){\Psi }^X(x,t)={\rho }_X(x,x^{\prime }|t),$$

(5)

where $|{\Psi }^X\left(t\right)\rangle$ refers to the separate state vector for X. Otherwise, the effects of Y manifest in a more or less relevant suppression of the off-diagonal terms, $(x,x^{\prime })$, in Equation (4), depending on the nature and strength of the coupling between X and Y. This phenomenon is what we regard as decoherence, which leads to seemingly classical density distributions [12], even though the full or joint state remains totally quantum mechanical.

In order to obtain a hydrodynamical picture [18] of the correlation effects, there are two ways to proceed. First, one could consider the reduced version of the problem, that is, the description enabled by Equation (3), and, then, compute the associated reduced quantum trajectories, as is done in Ref. [39]:

$$v^X(x,t)=\frac{1}{m}\frac{\mathrm{Re}\left[{\widehat{p}}_x{\tilde{\rho }}_X(x,x^{\prime }|t)\right]}{\mathrm{Re}\left[{\tilde{\rho }}_X(x,x^{\prime }|t)\right]}{|}_{x^{\prime }=x},$$

(6)

where ${\widehat{p}}_x=-i\hslash \partial /\partial x$ is the usual momentum operator. Note that the above transport equation describes the relationship between the quantum flux within the subspace for X, namely, $\tilde{j}(x,t)=\mathrm{Re}{\left[{\widehat{p}}_x{\tilde{\rho }}_X(x,x^{\prime }|t)\right]}_{x^{\prime }=x}$, and its corresponding reduced density distribution, ${\tilde{\rho }}_X(x,t)$. The use of Equation (6) thus serves only to analyze the dynamics constrained to the X-subspace, neglecting any direct or explicit contribution from the subsystem Y, since such a guidance equation only takes into account the (reduced) density matrix elements for X, Equation (4). In this way, we can only deal with reduced dynamics, which lack precise global information on the coherence swapping between the two subsystems. Only the overall effects felt by X are described (e.g., the gradual suppression of interference features), which does not remove other aspects, such as the non-crossing displayed by Bohmian trajectories, directly related to the suppression of mutual coherence between the waves passing through each slit (alternative mechanisms have to be incorporated in order to mimic this effect [40]).

To overcome the lacks and flaws involved in the above-mentioned procedure, particularly in those cases where the focus lies on investigating how the coherence swapping takes place, it is also possible to proceed by directly tackling the full problem through Equation (1), and, then, computing the trajectories for both X and Y, bearing in mind that each set of trajectories renders a dynamical picture on the corresponding subspace, and, hence, any interpretation must rely on this. In this case, though, it is not necessary to consider the density matrix approach (unless the joint state describes a statistical mixture, which is not the case here); the usual wave vector approach suffices for the purpose. Therefore, if the joint probability amplitude in the position space is denoted as $\Psi (x,y|t)$, then the corresponding equations of motion (the so-called guidance equations in Bohmian mechanics) are

$$\begin{array}{ccc}\hfill \dot{x}(x,y|t)& =& v^X(x,y|t)=\frac{1}{m}\mathrm{Re}\left[\frac{{\widehat{p}}_x\Psi (x,y|t)}{\Psi (x,y|t)}\right]=\frac{1}{m}\frac{\partial S(x,y|t)}{\partial x},\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \dot{y}(x,y|t)& =& v^Y(x,y|t)=\frac{1}{m}\mathrm{Re}\left[\frac{{\widehat{p}}_y\Psi (x,y|t)}{\Psi (x,y|t)}\right]=\frac{1}{m}\frac{\partial S(x,y|t)}{\partial y},\hfill \end{array}$$

(8)

which arise from the transport equation

$$\dot{\mathbf{r}}(x,y|t)=\mathbf{v}(x,y|t)=\frac{1}{m}\mathrm{Re}\left[\frac{\widehat{\mathbf{p}}\Psi (x,y|t)}{\Psi (x,y|t)}\right]=\frac{1}{m}\frac{\partial S(x,y|t)}{\partial y},$$

(9)

where $\widehat{\mathbf{p}}=-i\hslash \nabla$ denotes the usual momentum vector operator, with components $({\widehat{p}}_x,{\widehat{p}}_y)$$=-i\hslash (\partial /\partial x,\partial /\partial y)$, and

$$S(x,y|t)=\frac{\hslash }{2i}ln\left[\frac{\Psi (x,y|t)}{{\Psi }^{*}(x,y|t)}\right]$$

(10)

is a field that specifies the local phase associated with the wave function $\Psi (x,y|t)$. As can be noted, proceeding this way, the coupling between both subsystems becomes more apparent through the dependence of both velocities on the two coordinates. Moreover, also note that, in this manner, the information provided by the two velocity components allows us to not only better understand the reduced dynamics displayed by X later on, but also to better understand phenomena such as quantum erasure [24], or the so-called coherence swapping [41].

2.2. Dynamical Behavior for Uncorrelated Systems

The case of uncorrelated bipartite systems is trivial, because the quantum state of each party evolves separately. Yet, a brief analysis of the separate dynamics exhibited by one of the parties is important to better understand the behavior of both the full uncorrelated system and, of course, the entangled system. Thus, consider that the joint state vector for the uncorrelated bipartite system is

$$|\Psi \left(t\right)\rangle =|{\Psi }^X\left(t\right)\rangle \otimes |{\Psi }^Y\left(t\right)\rangle .$$

(11)

In the position representation, the corresponding wave function reads as

$$\Psi (x,y|t)={\Psi }^X(x,t){\Psi }^Y(y,t),$$

(12)

which leads to independent probability densities and transverse velocity fields for each subsystem. In other words, if these quantities are plotted in the corresponding subspaces (X or Y), no influence of one of the susbsystems on the other, and vice versa, is observed. Therefore, if we consider a joint representation, we can easily identify the features that belong to one or the other subsystem by simply inspecting what happens along the coordinate (x or y) that represents the corresponding subspace.

Below, we briefly analyze the main features associated with the dynamics exhibited by the time evolution of both a single Gaussian wave packet and a coherent superposition of two Gaussian wave packets, since they are used later to construct the bipartite states (see Section 2.3). Further details on the dynamical behaviors exhibited by these particular types of quantum states can be found in [17,31,42]. Physically, we here consider this particular type of wave packet, because it facilitates direct comparison with realistic physical scenarios of atom or molecular diffraction experiments, where short-range interactions between diffracting mask and diffracted object lead to transmission amplitudes describable in terms of Gaussian functions in a good approximation (e.g., see [38]). These Gaussian amplitudes are, thus, connected to the concept of the Gaussian slit introduced by Feynman and Hibbs [36,37], differing from the usual semiclassical connotations associated with this type of wave packet. Analogously, in the case of light, Gaussian beams can readily be associated with coherent states produced by laser sources coupled to optical fibers, which also model in a good approximation of real experiments (e.g., see [43]).

2.2.1. Single Gaussian Wave-Packet Dynamics

In this section, as well as in Section 2.2.2, we analyze the dynamics of a single system. Consequently, to simplify notations, all unnecessary superscripts and subscripts are removed, considering the x-coordinate to be generic. Thus, consider that the wave function at $t=0$, $\Psi (x,0)$, has the general functional form of a Gaussian wave packet with a phase (momentum) term

$${\mathcal{G}}_0(x,0)={\left(\frac{1}{2\pi {\sigma }_0}\right)}^{1/4}{e}^{-{(x-x_0)}^2/4{\sigma }_0^2+i{p}_0(x-x_0)/\hslash },$$

(13)

where $(x_0,p_0)$ specifies the phase–space position of its centroid and ${\sigma }_0$ its initial width. In free space, the time-evolved form reads as [42]

$${\mathcal{G}}_0(x,t)={\left(\frac{1}{2\pi {\tilde{\sigma }}_t^2}\right)}^{1/4}{e}^{-{(x-x_t)}^2/4{\sigma }_0{\tilde{\sigma }}_t+i{p}_0(x-x_t)/\hslash +i{E}_{0}t/\hslash },$$

(14)

which is also a Gaussian wave packet that propagates following the classical trajectory $x_t=x_0+v_{0}t$, with $v_0=p_0/m$ and $E_0=p_0^2/2m$, and where its spatial dispersion or spreading is described by the time-dependent complex function

$${\tilde{\sigma }}_t={\sigma }_0\left(1+\frac{i\hslash t}{2m{\sigma }_0^2}\right).$$

(15)

Rewriting the spreading function (15) in polar form, as ${\tilde{\sigma }}_t={\sigma }_t{e}^{i{\phi }_t}$, we note that the dispersion of the wave packet (14) is given by the function

$${\sigma }_t={\sigma }_0\sqrt{1+{\left(\frac{\hslash t}{2m{\sigma }_0^2}\right)}^2},$$

(16)

while the phase factor that develops as the wave packet evolves in time is

$${\phi }_t={\left(tan\right)}^{-1}\left(\frac{\hslash t}{2m{\sigma }_0^2}\right),$$

(17)

which can be related to the so-called geometric or Pancharatnam–Berry phase (or the Gouy phase in the context of optical Gaussian beams [44]). This latter term is responsible for the underlying dynamical behavior exhibited by the wave packet, as can readily be noticed if the Gaussian wave packet (14) is recast as

$$\begin{array}{ccc}\hfill {\mathcal{G}}_0(x,t)& =& {\left(\frac{1}{2\pi {\sigma }_t^2}\right)}^{1/4}{e}^{-{(x-x_t)}^2{e}^{-i{\phi }_t}/4{\sigma }_0{\sigma }_t}{e}^{(i/\hslash )\left[p_0(x-x_t)+E_{0}t-{\phi }_t/2\right]}\hfill \\ & =& {\left(\frac{1}{2\pi {\sigma }_t^2}\right)}^{1/4}{e}^{-{(x-x_t)}^2/4{\sigma }_t^2}{e}^{(i/\hslash )\left[\left(\right({\hslash }^{2}t)/2m{\sigma }_0^2){(x-x_t)}^2/4{\sigma }_t^2+p_0(x-x_t)+E_{0}t-{\phi }_t/2\right]}.\hfill \end{array}$$

(18)

This specific functional form for the wave packet makes the development of time-dependent phase factors more apparent and how this leads, in general, to dispersive effects, while the translation (described by the average momentum $p_0$ in this case) is only concerned with an intrinsically classical (Newtonian-like) law of motion (here, a uniform motion).

In any of the standard quantum pictures, the presence of the above-mentioned phase factor only manifests indirectly if the wave packet is made to interfere with itself (e.g., scattering off of an impenetrable barrier) or in a coherent superposition with another wave packet (see Section 2.2.2, since it is not a quantum observable, unlike the probability density, which is directly measurable with a detector, even on an event-by-event basis [45]. However, the hydrodynamical or Bohmian picture of quantum mechanics precisely highlights the role played by the spatial variations of the phase [17,18], and, hence, renders a different perspective on this phase factor. More specifically, such variations determine the local value of the transverse velocity field that acts on the probability density, which manifests in the form of a quantum flux [46]. This allows us to introduce the transport Equation (9), which, in the case of the single Gaussian wave packet (14), acquires the specific functional form

$$\dot{x}(x,t)=v(x,t)=v_0+\frac{{\hslash }^{2}t}{4m^2{\sigma }_0^2}\frac{(x-x_t)}{{\sigma }_t^2}.$$

(19)

This equation of motion contains information on both the translational motion, through the classical velocity, $v_0$, and the dispersive effect, connected to the phase factor, ${\phi }_t$. As can be noted in Equation (19), although the transverse velocity field $v(x,t)$ has the same value ($v_0$) at any position at $t=0$, at later times it becomes proportional to the ratio $x/t$, which explains the linear asymptotic spreading exhibited by the wave packet.

The integration in time of the equation of motion (19) is analytical and renders

$$x\left(t\right)=x_t+\frac{{\sigma }_t}{{\sigma }_0}\left[x\left(0\right)-x_0\right].$$

(20)

Following this equation, we find that, in the long-time limit, the Bohmian trajectories also describe a classical-like (Newtonian) uniform motion,

$$x\left(t\right)\sim v_{\infty }t,$$

(21)

with the asymptotic constant velocity being

$$v_{\infty }=v_0+\frac{\hslash \left[x\left(0\right)-x_0\right]}{2m{\sigma }_0^2}.$$

(22)

This constant velocity, which only depends on how far from the centroid the Gaussian wave packet trajectories are launched, is a direct consequence of the asymptotic linear behavior displayed by the transverse velocity field (19). Note that Equation (22) corresponds to the long-time limit of Equation (19) when it is evaluated along a given trajectory $x\left(t\right)$:

$$v\left[x\left(t\right)\right]=v_0+\frac{{\hslash }^{2}t}{4m^2{\sigma }_0^2}\frac{[x\left(0\right)-x_0]}{{\sigma }_0{\sigma }_t}.$$

(23)

2.2.2. Young-Type Superposition Dynamics

Young’s interference experiment belongs to a class of quantum scenarios where the state vector can be recast as a linear superposition (S) of different spatially separated alternatives. In particular, we describe the full state of a particle diffracted by two slits as

$$|{\Psi }_S\left(t\right)\rangle =\frac{1}{\sqrt{2}}\left[|{\Psi }_A\left(t\right)\rangle +|{\Psi }_B\left(t\right)\rangle \right],$$

(24)

which physically represents a coherent superposition having an incident wave diffracted by two slits, A and B, spatially separated (ensuring orthogonality) and with the same amplitude (probability). In the position representation, we can represent this situation by means of the coherent superposition of two Gaussian wave packets describing the diffracted wave coming up from each slit (e.g., the outcomes in an experiment such as the one reported in [32]):

$${\Psi }_S(x,t)=\mathcal{N}\left[{\mathcal{G}}_A(x,t)+{\mathcal{G}}_B(x,t)\right],$$

(25)

where

$$\begin{array}{ccc}\hfill {\mathcal{G}}_A(x,t)& =& {\left(\frac{1}{2\pi {\tilde{\sigma }}_t^2}\right)}^{1/4}{e}^{-{(x-x_{A,t})}^2/4{\sigma }_0{\tilde{\sigma }}_t+i{p}_A(x-x_{A,t})/\hslash +i{E}_{A}t/\hslash },\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill {\mathcal{G}}_B(x,t)& =& {\left(\frac{1}{2\pi {\tilde{\sigma }}_t^2}\right)}^{1/4}{e}^{-{(x-x_{B,t})}^2/4{\sigma }_0{\tilde{\sigma }}_t+i{p}_B(x-x_{B,t})/\hslash +i{E}_{B}t/\hslash },\hfill \end{array}$$

(27)

where $X_{A/B,t}$ has the same meaning as $X_t$ in the previous case, but with A and B denoting the corresponding centroid position, namely, $x_{A,0}=d/2=-x_{B,0}$, with d being the separation between both centroids (which is the same as the center-to-center distance between the two slits). Regarding the respective momenta, for simplicity, we consider no transverse drift, as this is not relevant to the discussion here, so $p_A=p_B=0$ (and, hence, $E_A=E_B=0$), On the other hand, the two Gaussian wave packets might not have a vanishing overlapping, since

$$\mathcal{N}=\frac{1}{\sqrt{2\left(1+e^{-d^2/8{\sigma }_0^2}\right)}},$$

(28)

thus, invalidating the condition of mutual orthogonality between the states representing the diffraction through the slits. To overcome this inconvenience, from now on we assume the condition $d/2\gg {\sigma }_0$, so that $\mathcal{N}\approx 1/\sqrt{2}$.

From the density matrix coefficients associated with the wave function (25),

$$\rho (x,x^{\prime }|t)\approx \frac{1}{2}\left[{\mathcal{G}}_A(x,t){\mathcal{G}}_A^{*}(x^{\prime },t)+{\mathcal{G}}_B(x,t){\mathcal{G}}_B^{*}(x^{\prime },t)+{\mathcal{G}}_A(x,t){\mathcal{G}}_B^{*}(x^{\prime },t)+{\mathcal{G}}_B(x,t){\mathcal{G}}_A^{*}(x^{\prime },t)\right],$$

(29)

we readily obtain the probability density, which reads as

$$\rho (x,t)\approx \frac{1}{2}{\left(\frac{1}{2\pi {\sigma }_t^2}\right)}^{1/2}\left[e^{-{(x-d/2)}^2/2{\sigma }_t^2}+e^{-{(x+d/2)}^2/2{\sigma }_t^2}+2e^{-[x^2+{(d/2)}^2]/2{\sigma }_t^2}cos{k}_{t}x\right],$$

(30)

where

$$k_t=\frac{\hslash td}{4m{\sigma }_0^2{\sigma }_t^2}.$$

(31)

For relatively short times, compared to $\tau =2m{\sigma }_0^2/\hslash$ (i.e., $t\ll \tau$), we have ${\sigma }_t\approx {\sigma }_0$, and Equation (42) essentially consists of two separated Gaussian wave packets centered around $x_A$ and $x_B$:

$$\rho (x,t)\approx \frac{1}{2}{\left(\frac{1}{2\pi {\sigma }_0^2}\right)}^{1/2}\left[e^{-{(x-d/2)}^2/2{\sigma }_0^2}+e^{-{(x+d/2)}^2/2{\sigma }_0^2}\right].$$

(32)

On the other hand, for $t\gg \tau$, the width of each wave packet increases linearly with time, since ${\sigma }_t\approx \hslash t/2m{\sigma }_0$ (see Section 2.2.1), and a Young-type interference pattern arises with a Gaussian envelope:

$$\rho (x,t)\approx 2{\left(\frac{2m^2{\sigma }_0^2}{\pi {\hslash }^2{t}^2}\right)}^{1/2}{e}^{-2m^2{\sigma }_0^2{x}^2/{\hslash }^2{t}^2}{cos}^2(k_{\infty }x/2),$$

(33)

with

$$k_{\infty }=\frac{md}{\hslash t}.$$

(34)

As is inferred from (33), at a given time, the condition to observe an interference minimum is

$$x_n=\left(n+\frac{1}{2}\right)\frac{2\pi \hslash t}{md},$$

(35)

where $n=0,\pm 1,\pm 2,\dots$. Therefore, the size of the interference fringes (distance between two consecutive minima) increase linearly with time, as

$$\Delta x=\frac{2\pi \hslash t}{md},$$

(36)

while the number of interference fringes within the region covered by the Gaussian envelope increase rapidly as d decreases.

In order to understand the appearance of fringes from the dynamical point of view rendered by the phase variations of the wave function (25), let us now examine the expression for $v(x,t)$, which is readily obtained after substituting this wave function into Equation (7) in the factorizable case:

$$\begin{array}{ccc}\hfill \dot{x}(x,t)& \approx & \frac{1}{2\rho (x,t)}{\left(\frac{1}{2\pi {\sigma }_t^2}\right)}^{1/2}\frac{{\hslash }^{2}t}{4m^2{\sigma }_0^2{\sigma }_t^2}\{(x-d/2)e^{-{(x-d/2)}^2/2{\sigma }_t^2}+(x+d/2)e^{-{(x+d/2)}^2/2{\sigma }_t^2}\hfill \\ & & +2x{e}^{-[x^2+{(d/2)}^2]/2{\sigma }_t^2}cos{k}_{t}x+\frac{2m{\sigma }_0^2{x}_0}{\hslash t}{e}^{-[x^2+{(d/2)}^2]/2{\sigma }_t^2}sin{k}_{t}x\}.\hfill \end{array}$$

(37)

This highly nonlinear expression already gives us an idea of the complexity of the dynamics. Yet, some conclusions can be extracted from it even without solving it numerically. Since this general expression was already discussed in detail in [31], here we focus on the above two limits. Note that this expression cancels out at $x=0$ (which is at the origin of the so-called Bohmian non-crossing rule [31]), although it remains finite on either side of this boundary. Thus, in the short-time limit, although there is a non-negligible contribution, the density mainly concentrates around the slit A or the slit B (i.e., around $x_A=d/2$ or $x_B=-d/2$). Accordingly, the trajectories launched from a vicinity of these points describes a motion analogous to that exhibited by the trajectories associated with a free Gaussian wave packet, because Equation (37) approximates

$${\dot{x}}_c(x,t)\approx \frac{{\hslash }^{2}t}{4m^2{\sigma }_0^4}(x-x_0),$$

(38)

with $c=A,B$. Thus, the larger d, the more the trajectories resemble those for a free Gaussian (centered either on site A or site B).

In the long-time limit, though, we find rather different behavior with respect to the single-Gaussian case, as Equation (37) approximates

$$\dot{x}(x,t)\approx \frac{x}{t}\left[1-\frac{m{\sigma }_0^{2}d}{2\hslash tx}tan(k_{\infty }x/2)\right].$$

(39)

Accordingly, at a given time, we observe a nearly linear increase with x, except at the regions around an interference minimum, specified by the condition (35), where the velocity undergoes a sudden change because the argument of the tangent equals $\pi /2$. For positive x, these sudden changes in the velocity are in the form of negative dips or spikes, while for negative x they have the form of positive spikes. As a consequence, at relatively large t, trajectories tend to mainly accumulate within the region determined by two consecutive neighboring spikes, where they display a nearly uniform motion, with a (quantized) mean value for the momentum given, in a good approximation, by the expression

$$\hslash {\kappa }_n=mv(\overline{x},t)=mv[(x_{n+1/2}+x_{n-1/2})/2,t]\approx \frac{2\pi \hslash }{d}n,$$

(40)

with $n=0,\pm 1,\pm 2,\dots$, and where $x_{n+1/2}$ and $x_{n-1/2}$ denote the position of the two spikes surrounding the region with momentum ${\kappa }_n$. This momentum value, quantized in terms of the unit $2\pi \hslash /d$, determines the long-time position (on average) of interference maxima: $x_n\approx \hslash {\kappa }_{n}t/m$. In this way, swarms of trajectories with nearly the same constant velocity constitute the dynamical consequence of the organized phase variations arising after the overlapping of the two mutually (phase) coherent Gaussian wave packets. The probability density, statistically specified in terms of these swarms of trajectories, thus tends to accumulate in the regions where the transverse velocity field is nearly constant, between any two consecutive velocity “spikes”, while it is expelled by means of a rather fast flow when it is in the vicinity of these “spikes”, always from the inner regions to the outer ones. This is analogous to the accumulation of sediments on the sandbanks of a river, where the watercourse is slow, while they are abruptly expelled in the rapids because of the fast streams.

2.3. Entangled Bipartite Systems

Let us now focus on the case of a Bell-type maximally entangled (E) state:

$$|{\Psi }_E\left(t\right)\rangle =\frac{1}{\sqrt{2}}\left[|{\psi }_A^X\left(t\right)\rangle |{\psi }_B^Y\left(t\right)\rangle +|{\psi }_B^X\left(t\right)\rangle |{\Psi }_A^Y\left(t\right)\rangle \right],$$

(41)

which describes a situation with two possible slits, A and B, and two correlated identical particles, X and Y, that can get diffracted through either of the slits, but not both, as in the standard superposition case seen in the previous section. More specifically, unlike the (factorizable) bipartite vector state for both particles undergoing diffraction through both slits,

$$\begin{array}{ccc}\hfill |{\Psi }_S\left(t\right)\rangle & =& |{\Psi }_S^X\left(t\right)\rangle \otimes |{\Psi }_S^Y\left(t\right)\rangle \hfill \\ & =& \frac{1}{2}\left[|{\psi }_A^X\left(t\right)\rangle |{\psi }_A^Y\left(t\right)\rangle +|{\psi }_A^X\left(t\right)\rangle |{\psi }_B^Y\left(t\right)\rangle +|{\psi }_B^X\left(t\right)\rangle |{\psi }_A^Y\left(t\right)\rangle +|{\psi }_B^X\left(t\right)\rangle |{\psi }_B^Y\left(t\right)\rangle \right],\hfill \end{array}$$

(42)

the entangled state (41) lacks the possibility to have both particles diffracted through the same slit [first and fourth terms in (42)], thus, describing a highly delocalized two-particle system. To further investigate the associated dynamics, let us consider the wave function in the position representation for (41),

$${\Psi }_E(x,y|t)={\mathcal{N}}_E\left[{\mathcal{G}}_A(x,t){\mathcal{G}}_B(y,t)+{\mathcal{G}}_B(x,t){\mathcal{G}}_A(y,t)\right],$$

(43)

where the wave packets are as given by (26) and (27), although substituting x by y wherever it corresponds, and the normalizing prefactor is

$${\mathcal{N}}_E=\frac{1}{\sqrt{2\left(1+e^{-d^2/4{\sigma }_0^2}\right)}}.$$

(44)

As before, the choice of d makes ${\mathcal{N}}_E\approx 1/\sqrt{2}$.

From the above considerations, we obtain the expression for the density matrix elements,

$$\begin{array}{ccc}\hfill \rho (x,y;x^{\prime },y^{\prime }|t)& \approx & \frac{1}{2}[{\mathcal{G}}_A(x,t){\mathcal{G}}_A^{*}(x^{\prime },t){\mathcal{G}}_B(y,t){\mathcal{G}}_B^{*}(y^{\prime },t)+{\mathcal{G}}_B(x,t){\mathcal{G}}_B^{*}(x^{\prime },t){\mathcal{G}}_A(y,t){\mathcal{G}}_A^{*}(y^{\prime },t)\hfill \\ & & +{\mathcal{G}}_A(x,t){\mathcal{G}}_B^{*}(x^{\prime },t){\mathcal{G}}_B(y,t){\mathcal{G}}_A^{*}(y^{\prime },t)+{\mathcal{G}}_B(x,t){\mathcal{G}}_A^{*}(x^{\prime },t){\mathcal{G}}_A(y,t){\mathcal{G}}_B^{*}(y^{\prime },t)],\hfill \end{array}$$

(45)

from which the diagonal elements provide us with the bipartite probability density,

$$\begin{array}{ccc}\hfill \rho (x,y|t)& \approx & \frac{1}{2}\left(\frac{1}{2\pi {\sigma }_t^2}\right)\{e^{-{(x-d/2)}^2/2{\sigma }_t^2-{(y+d/2)}^2/2{\sigma }_t^2}+e^{-{(x+d/2)}^2/2{\sigma }_t^2-{(y-d/2)}^2/2{\sigma }_t^2}\hfill \\ & & +2e^{-[x^2+y^2+2{(d/2)}^2]/2{\sigma }_t^2}cos\left[k_t(x-y)\right]\}.\hfill \end{array}$$

(46)

Let us now proceed as before and investigate the short-time and long-time behaviors exhibited by this probability density. In the short-time limit ($t\ll \tau$), it basically describes two highly delocalized systems with maximal uncertainty about both [35]:

$$\rho (x,y|t)\approx \frac{1}{2}\left(\frac{1}{2\pi {\sigma }_0^2}\right)\left[e^{-{(x-d/2)}^2/2{\sigma }_0^2-{(y+d/2)}^2/2{\sigma }_0^2}+e^{-{(x+d/2)}^2/2{\sigma }_0^2-{(y-d/2)}^2/2{\sigma }_0^2}\right].$$

(47)

On the other hand, in the long-time limit ($t\gg \tau$), we have

$$\rho (x,y|t)\approx 2\left(\frac{2m^2{\sigma }_0^2}{\pi {\hslash }^2{t}^2}\right)e^{-(x^2+y^2)/2{\sigma }_t^2}{cos}^2\left[k_{\infty }(x-y)/2\right].$$

(48)

According to this expression, the density distribution exhibits full fringe visibility along the direction $y=-x$, since

$${\rho }_{\infty }(x,-x|t)\approx 2\left(\frac{2m^2{\sigma }_0^2}{\pi {\hslash }^2{t}^2}\right)e^{-4m^2{\sigma }_0^2{x}^2/{\hslash }^2{t}^2}{cos}^2\left(k_{\infty }x\right),$$

(49)

while it is totally suppressed along the $y=x$, where

$${\rho }_{\infty }(x,x|t)\approx 2\left(\frac{2m^2{\sigma }_0^2}{\pi {\hslash }^2{t}^2}\right)e^{-4m^2{\sigma }_0^2{x}^2/{\hslash }^2{t}^2}.$$

(50)

The presence of the oscillatory term in Equation (46) can be interpreted as a signature of the coherence exchange, due to correlations between the X and Y subsystems, which is maximal, as we have seen, along the direction that joins them directly. It disappears gradually as the perpendicular direction is approached (as should be expected, since there was no probability to find both subsystems at the same position initially, i.e., around $x=y=\pm d/2$). This strong correlation has a direct implication when we investigate the behavior of the subsystem X. Tracing over the coordinates for the Y subsystem in (45), we obtain the following expression for the reduced density matrix elements of the subsystem X:

$$\tilde{\rho }(x,x^{\prime }|t)\approx \frac{1}{2}\left[{\mathcal{G}}_A(x,t){\mathcal{G}}_A^{*}(x^{\prime },t)+{\mathcal{G}}_B(x,t){\mathcal{G}}_B^{*}(x^{\prime },t)+{\Lambda }_{AB}^{*}{\mathcal{G}}_A(x,t){\mathcal{G}}_B^{*}(x^{\prime },t)+{\Lambda }_{AB}{\mathcal{G}}_B(x,t){\mathcal{G}}_A^{*}(x^{\prime },t)\right],$$

(51)

where

$${\Lambda }_{AB}=\int {\mathcal{G}}_A(y,t){\mathcal{G}}_B^{*}(y,t)dy=e^{-d^2/8{\sigma }_0^2}.$$

(52)

Therefore, the reduced probability density reads as

$$\tilde{\rho }(x,t)\approx \frac{1}{2}{\left(\frac{1}{2\pi {\sigma }_t^2}\right)}^{1/2}\left[e^{-{(x-d/2)}^2/2{\sigma }_t^2}+e^{-{(x+d/2)}^2/2{\sigma }_t^2}+2{\Lambda }_{AB}{e}^{-[x^2+{(d/2)}^2]/2{\sigma }_t^2}cos{k}_{t}x\right],$$

(53)

where the interference term is partially, or even totally, suppressed by the prefactor ${\Lambda }_{AB}$, depending on the balance between the width of the wave packets and their distances. Since we typically consider $d/2$ to be larger than ${\sigma }_0$ to avoid the initial overlapping of the wave packets, it is clear that (53) represents in both the short-time and the long-time limits an incoherent sum of two independent wave packets. It is precisely because of this incoherence effect in the X-subspace, in the long-time limit, that the density distribution is basically described by a seemingly single Gaussian distribution:

$$\tilde{\rho }(x,t)\approx {\left(\frac{2m{\sigma }_0^2}{\pi \hslash t}\right)}^{1/2}{e}^{-x^2/2{\sigma }_t^2}.$$

(54)

This arises as a consequence of the screening produced by the correlation with the Y subsystem, which leads to a statistical mixture, without phase relations, in the X subsystem, according to Schrödinger’s picture of entanglement [47].

The above behavior is what we commonly regard as decoherence [6,48,49,50], that is, the loss of phase correlations by interaction with other quantum systems. However, even after having lost such phase correlations and observing seemingly classical behavior for the X subsystem, the point is that, if we analyze the corresponding velocity fields, the correlations survives due to the nonlinear nature of these fields, both in the full space and also in the reduced one. More specifically, if both subsystems are considered, we obtain the following equations of motion:

$$\begin{array}{ccc}\hfill \dot{x}(x,y|t)& \approx & \frac{1}{2\rho (x,y|t)}\left(\frac{1}{2\pi {\sigma }_t^2}\right)\frac{{\hslash }^{2}t}{4m^2{\sigma }_0^2{\sigma }_t^2}\hfill \\ & & \times \{(x-d/2)e^{-{(x-d/2)}^2/2{\sigma }_t^2-{(y+d/2)}^2/2{\sigma }_t^2}+(x+d/2)e^{-{(x+d/2)}^2/2{\sigma }_t^2-{(y-d/2)}^2/2{\sigma }_t^2}\hfill \\ & & +2x{e}^{-[x^2+y^2+2{(d/2)}^2]/2{\sigma }_t^2}cos\left[k_t(x-y)\right]+\frac{m{\sigma }_0^{2}d}{\hslash t}{e}^{-[x^2+y^2+2{(d/2)}^2]/2{\sigma }_t^2}sin\left[k_t(x-y)\right]\},\hfill \end{array}$$

(55)

$$\begin{array}{ccc}\hfill \dot{y}(x,y|t)& \approx & \frac{1}{2\rho (x,y|t)}\left(\frac{1}{2\pi {\sigma }_t^2}\right)\frac{{\hslash }^{2}t}{4m^2{\sigma }_0^2{\sigma }_t^2}\hfill \\ & & \times \{(y+d/2)e^{-{(x-d/2)}^2/2{\sigma }_t^2-{(y+d/2)}^2/2{\sigma }_t^2}+(y-d/2)e^{-{(x+d/2)}^2/2{\sigma }_t^2-{(y-d/2)}^2/2{\sigma }_t^2}\hfill \\ & & +2y{e}^{-[x^2+y^2+2{(d/2)}^2]/2{\sigma }_t^2}cos\left[k_t(x-y)\right]-\frac{m{\sigma }_0^{2}d}{\hslash t}{e}^{-[x^2+y^2+2{(d/2)}^2]/2{\sigma }_t^2}sin\left[k_t(x-y)\right]\},\hfill \end{array}$$

(56)

which provide us with relevant information on the strong correlation (and its effects) between the two subsystems through their nonseparability. According to these equations, the evolution of X is strongly connected to the evolution of Y, and vice versa.

On the other hand, in the reduced case, although the explicit presence of the y-coordinate has been removed, its effects are still felt through the prefactor ${\Lambda }_{AB}$, and we have

$$\begin{array}{ccc}\hfill \dot{\tilde{x}}(x,t)& \approx & \frac{1}{2\tilde{\rho }(x,t)}{\left(\frac{1}{2\pi {\sigma }_t^2}\right)}^{1/2}\frac{{\hslash }^{2}t}{4m^2{\sigma }_0^2{\sigma }_t^2}\{(x-d/2)e^{-{(x-d/2)}^2/2{\sigma }_t^2}+(x+d/2)e^{-{(x+d/2)}^2/2{\sigma }_t^2}\hfill \\ & & +2x{\Lambda }_{AB}{e}^{-[x^2+{(d/2)}^2]/2{\sigma }_t^2}cos{k}_{t}x+{\Lambda }_{AB}\frac{m{\sigma }_0^{2}d}{\hslash t}{e}^{-[x^2+{(d/2)}^2]/2{\sigma }_t^2}sin{k}_{t}x\},\hfill \end{array}$$

(57)

which, given the value chosen for d, can be simplified to

$$\dot{\tilde{x}}(x,t)\approx \left(\frac{{\hslash }^{2}t}{4m^2{\sigma }_0^2{\sigma }_t^2}\right)\left\{\frac{(x-d/2)e^{-{(x-d/2)}^2/2{\sigma }_t^2}+(x+d/2)e^{-{(x+d/2)}^2/2{\sigma }_t^2}}{e^{-{(x-d/2)}^2/2{\sigma }_t^2}+e^{-{(x+d/2)}^2/2{\sigma }_t^2}}\right\}.$$

(58)

Observe that, in the short-time limit, this expression renders two swarms of trajectories that seem to be independent of one another, each one being associated with a single wave packet. On the other hand, in the long-time limit, we have a nearly homogeneous distribution of trajectories that seem to be reproducing the behavior of a single motionless ($v_0=0$) wave packet centered at $x=0$, with their effective equation of motion being

$$\dot{\tilde{x}}(x,t)\approx \left(\frac{{\hslash }^{2}tx}{4m^2{\sigma }_0^2{\sigma }_t^2}\right).$$

(59)

Although we have seen that the full and the reduced dynamics are very different from the point of view of the flux (trajectories), the former still has to provide limiting behaviors analogous to those exhibited by the latter. In other words, although the full dynamics does not need to preserve the Bohmian non-crossing property on the corresponding subspace, because it only remains valid in the full-dimensional space, the trajectories have to show the same asymptotic behavior in both cases (i.e., for $t\gg \tau$).

3. Results and Discussion

3.1. Single-System Dynamical Behaviors

In order to better understand the behaviors exhibited by the separable and entangled bipartite systems, let us first briefly revisit the main dynamical aspects of single localized and delocalized systems, respectively, described by the wave functions (14) and (25), and represented in the upper and lower rows of Figure 1. For the single Gaussian wave packet (14) (see Section 2.2.1), the time evolution of the corresponding probability density is represented in the form of a contour plot in Figure 1a. As is seen, it exhibits a slow dispersion initially until it undergoes a fast increase around $t\approx \tau$ (with the parameters used in the simulations, $\tau =0.5$), which eventually leads to a linearly increasing dispersion regime. The associated quantum flux is illustrated with the depiction of 21 Bohmian trajectories, which follow Equation (20) and have equidistant initial positions. We clearly observe a transition from a swarm of nearly parallel trajectories to a linear separation among them, with increasingly faster speeds as the initial condition gets further away from the center of the distribution.

The dynamics described by the probability density are a manifestation of the underlying flux, ruled by the local variations of the phase term that appears in the wave function along its time evolution, as can clearly be seen in Equation (18). These phase variations, as mentioned before, generate an underlying transverse velocity field, Equation (19), which gradually changes in time, as can be observed in Figure 1b in the form of a contour plot, and in Figure 1c in the form of cross sections at specific times. As seen in Figure 1c, at a given time, this velocity field increases linearly with x, pivoting around $x=0$; thus, provoking trajectories with $x\left(0\right)>0$ to acquire positive outwards momentum and trajectories with $x\left(0\right)<0$ acquiring negative outwards momentum. This eventually leads to the spread of the whole swarm as time increases. Indeed, as is also seen in Figure 1c, there is an almost linear increase of the slope of the velocity field from zero, until it reaches a maximum value at $t=\tau$ [as can be inferred from Equation (19)]; then, the slope starts decreasing again, with nearly $t^{-1}$, until it approaches zero again.

Such behaviors become more apparent by inspecting the diffusive time-dependent prefactors in Equations (19) and (23), for $v(x,t)$ and $v\left[x\right(t\left)\right]$, respectively, which are displayed, in Figure 2, as functions of time. In the case of the overall velocity field (black line), and taking into account the above-mentioned transition at $t=\tau$, we effectively appreciate two regimes, namely, a nearly uniform acceleration for $t<\tau$ and, then, an asymptotic fall with $t^{-1}$ for $t>\tau$. On the other hand, if we consider the velocity along a given trajectory (red line), the first stage is characterized by nearly uniform acceleration, which provokes fast diffusion of the trajectories, and, then, beyond $t=\tau$, the trajectories acquire a constant speed, proper for uniform motion, as in classical mechanics.

The behaviors that explain the dynamics of a single Gaussian wave packet also hold for a coherent superposition of two of them, like the one described by Equation (25), provided they are far enough from each other, as is seen at early stages in Figure 1d,e, compared to Figure 1a,b, respectively. There is a major difference, though, between the picture offered by the probability density here and the plot of the velocity field. Although the velocity field is zero at $t=0$, because of the existence of phase coherence, since the very early stages a sharp spatial division arises, from a dynamical point of view, with two independent velocity fields that do not overlap, but that undergo a sharp shear along the symmetry line $x=0$. This behavior is observed until $t\approx 2$, when the dispersion of both wave packets is enough to provoke important overlapping between them, analogous to that exhibited by the single wave packet, although with the slope of the velocity field pivoting around the respective centroids. This behavior is better appreciated with the aid of the corresponding trajectories, which show the specific flux evolution at a more local level. These trajectories, superimposed to the density plots, are obtained after numerical integration of Equation (37) with initial conditions equally distributed within the regions covered by each wave packet. As can be seen, the two swarms of trajectories show a seemingly independent evolution (diffraction) at short times.

At later times, an increasingly prominent overlapping of the two wave packets leads to the appearance of interference traits in the probability density, which are associated with the appearance of regions with relatively constant velocity, interrupted by regions wherein the velocity undergoes sudden changes. The first regions, more stable from a dynamical point of view, enable the accumulation of trajectories along them with nearly the same velocity (on average), which explains the interference maxima. On the other hand, the latter imprint an important shake on trajectories passing by too close to them, which leads them to quickly change from one stable region to the neighboring one, in an inwards staggered motion. This promotion of trajectories from the outer regions of the wave function to the innermost ones eventually leads to a higher population of the corresponding (innermost) interference maxima with respect to the more marginal ones.

Another interesting feature in the long-time regime is that, from two spatially divided dynamical regions, the velocity tends to exhibit a single average slope (interrupted by a series of sudden kicks). However, although this average slope seems to slowly decrease to zero, as we might infer from Figure 1f, what actually happens is that it converges to an almost constant slope, in compliance with the asymptotic value given by Equation (39), such that each stable section (between any two kicks) remains at the same height with respect to the central value $v(0,t)=0$. The average value for the height reached by each section is given by Equation (40), as is indicated by the corresponding horizontal dashed lines. This behavior was actually observed in the experiments with photons performed by Kocsis et al. [32] more than a decade ago. With the aid of weak measurements, it was possible to measure the transverse momentum, which is a quantity precisely proportional to the Bohmian instantaneous velocity field (although the latter refers to massive particles and, in the case of photons, is related to the transverse component of the wave vector under paraxial conditions [51,52]).

3.2. Bipartite-System Dynamical Behaviors

Let us now consider the dynamics generated by bipartite separable and entangled states, and, in particular, how the velocity field that rules the dynamics is influenced by these states. The latter has important implications regarding the spatial distribution of the quantum flux for each subsystem, and, hence, the corresponding (reduced) probability densities and the observation of interference traits. Unlike the analysis presented in Section 3.1, here, the joint dynamics does not facilitate an analogous space–time description, so we consider a series of illustrative snapshots instead, for both the (joint) probability density and the transverse velocity field. Actually, in the latter case, given that this field depends on both subsystem coordinates x and y, we need to consider separate plots for the velocity fields $v_x(x,y|t)$ and $v_y(x,y|t)$, specified, respectively, by Equations (55) and (56). Regarding the trajectories, we consider the trajectories projected onto the corresponding subspaces (full higher-dimensional pictures of the trajectories within the joint space are given in the Supplementary Material). Specifically, in all cases, to investigate the dynamics generated by the corresponding joint wave functions, we considered sets of trajectories with equidistant initial conditions that mapped (initially) either the state of X or the state of Y, and that played the role of flow markers. These initial conditions (21 for each wave packet) are uniformly distributed around the center of the wave packets. In the contour plots, these flow markers appear as dots with different colors along the x and y directions for better identification in the long-time regime. Thus, at $t=0$, they appear as a cross-shaped distribution (see first column panels in Figure 3), while at $t=10$ they distribute in different ways (see last column panels in Figure 3, Figure 4 and Figure 5).

In the three cases analyzed, the X subsystem was described by a two-Gaussian coherent superposition, as given by Equation (25), as a simplified version of Young’s two-slit experiment. In the case of the Y subsystem, we firstly considered a localized state, such as the Gaussian wave packet (14). Accordingly, the joint wave function read as

$$\Psi (x,y|t)=\mathcal{N}\left[{\mathcal{G}}_A(x,t)+{\mathcal{G}}_B(x,t)\right]{\mathcal{G}}_0(y,t),$$

(60)

The initial probability density is represented in Figure 3a, which effectively shows that, while Y was localized around $x_0$, X was delocalized between $x_{A,0}$ and $x_{B,0}$. In Figure 3b–d we observe the evolution of the joint probability density at $t=1$, $t=5$, and $t=10$, respectively. Due to the factorizability of the state, we could easily identify the different stages of the evolution for X and Y with the corresponding cases, shown in Figure 1a,d. The spatial dispersion of the markers was also in compliance with the corresponding subsystems. If we consider the markers distributed along the y-direction, their mutual separations increased homogeneously, as this corresponded to the Gaussian state analyzed in Section 3.1. On the other hand, because of the emergence of interference traits, if we consider the markers along the x-direction, their distribution was uneven, showing accumulations in the regions covered by maxima and voids in the regions corresponding to interference minima, separated by a distance $\Delta x=2\pi$ at $t=10$, according to Equation (36). These dynamics, though, were governed by the underlying velocity fields, as seen in the upper row panels in Figure 4 and Figure 5, for the $v^X(x,y|t)$ and $v^Y(x,y|t)$, respectively. In particular, in Figure 4a–c, we observe the snapshots for $v^X(x,y|t)$ at $t=1$, $t=5$, and $t=10$, respectively, while the analogous cases for $v^Y(x,y|t)$ are represented in Figure 5a–c. It is worth noting that, unlike the probability density, which was spatially modulated by the Gaussian wave packets, decaying to zero at some point in either direction, the velocity fields remained independent of the complementary coordinate. Thus, at time t, $v^X(x,y|t)$ kept the same value for a given x-position, regardless of the value for y. The same held for $v^Y(x,y|t)$ with respect to x. This is particularly interesting in the case of X, where its associated velocity field showed, for any value of y, the periodic alternation of nearly constant regions, in spatial intervals with a value $\Delta x=2\pi$ at $t=10$, again in compliance with Equation (36).

If we consider the case where both subsystems, X and Y, are described by identical delocalized superpositions,

$$\Psi (x,y|t)={\mathcal{N}}^2\left[{\mathcal{G}}_A(x,t)+{\mathcal{G}}_B(x,t)\right]\left[{\mathcal{G}}_A(y,t)+{\mathcal{G}}_B(y,t)\right],$$

(61)

the situation turns a bit more complicated, because now the effects of interference emerge along both directions, x and y, giving rise to the appearance of a highly organized chessboard-like pattern of local maxima modulated by a two-dimensional Gaussian envelope, as seen in Figure 3h. This can be inferred from the direct product between Equation (33) and its equivalent for the Y subsystem:

$$\rho (x,y|t)=\rho (x,t)\rho (y,t)\approx 4\left(\frac{2m^2{\sigma }_0^2}{\pi {\hslash }^2{t}^2}\right)e^{-2m^2{\sigma }_0^2(x^2+y^2)/{\hslash }^2{t}^2}{cos}^2(k_{\infty }x/2){cos}^2(k_{\infty }y/2).$$

(62)

A series of snapshots illustrating the trend, until the two-dimensional interference chessboard-like pattern has reached full visibility, is shown in Figure 3e–h. It can be noticed that, as in the previous case, sets of marker trajectories starting with the same value of $x\left(0\right)$ or $y\left(0\right)$ do not mix (cross) with trajectories launched with other initial conditions, in spite of the intricate pattern describing the asymptotic behavior of the probability density. The different colors assigned to the markers along the x and y directions for the trajectories associated with the four different distributions make evident this behavior. Nonetheless, the dynamical explanation can be found again in the underlying structure displayed by the velocity field, shown in Figure 4e–h for $v^X(x,y|t)$, and in Figure 5e–h for $v^Y(x,y|t)$. As before, the factorizability of the joint state gives rise to two totally independent velocity fields, although both have the same structure, which asymptotically displays the known alternation of stable regions separated by sudden changes in the velocity at even distances of $\Delta x=2\pi$, for $t=10$. According to this velocity pattern, the markers start acquiring an uneven spatial distribution as they become more and more dispersed.

Next, let us take a look at the topology displayed by the trajectories followed by the markers considered in the two cases analyzed so far. A graphical representation on the dynamics for the three types of bipartite states here considered is provided in Figure 6. In the first column panels, we consider projections of the x and y trajectories associated with the wave function (60). In Figure 6a, we have the x and y-components of trajectories specified with initial conditions of the type $\left[x\right(0),0]$ for the X subsystem, with $x\left(0\right)$ evenly distributed within the two regions covered by the two wave packets along the x-direction (see horizontal markers in Figure 3a). As can be seen, while the x-component (solid black lines) was in agreement with the behavior exhibited by the analogous one-dimensional system, described by Equation (25) and discussed in Section 3.1, the y-component (red solid line) remained constant [$y=y\left(0\right)=0$] at any time. For any other value of $y\left(0\right)$, we would have obtained a similar behavior, with all trajectories showing an identical behavior (for the same initial condition) along x, but with the whole set moving along the vertical in compliance with the trajectory corresponding to a single Gaussian (see Supplementary Materials). In Figure 6a′, however, the behavior is a bit different. In this case, the initial conditions were of the type $[x_{A/B,0},y\left(0\right)]$, with $y\left(0\right)$ evenly distributed within the region covered by a single wave packet along the y-direction (see markers in Figure 3a along y). We observe that the y-component of the trajectories proceeded in agreement with the behavior studied in Section 3.1 for a single Gaussian wave packet. However, the x-component of the trajectories indicated wiggling behavior, since they were not computed at $x\left(0\right)=0$, but at either $x\left(0\right)=x_{A,0}$ or $x\left(0\right)=x_{B,0}$, and, hence, there is going to be a contribution to the motion of the x-component of the velocity field. As a direct consequence of these behaviors, determined by inspecting the trajectory dynamics, either in the X subspace or in the Y one, we note that there was a strong connection between separability, a formal aspect, and diffusion, a physical behavior. Specifically, if we consider a slice of the joint probability density at a constant value of x (or y), it looks the same at any time as that any other slice taken at another constant value of x (or y). The associated (full) velocity field allowed such splittings at a dynamical level, while the trajectories showed local (instantaneous) features of the diffusive process as time proceeded. This consequence was also confirmed in the case wherein X and Y were both described by delocalized wave functions, such as in (61). The corresponding trajectories (x and y-components, as before) in the X-subspace are plotted in Figure 6b. In this case, due to the symmetry along both directions, the sets of trajectories with initial conditions $[x_{A/B,0},y\left(0\right)]$ evolved following the same dynamics as their counterparts with initial conditions $[x\left(0\right),y_{A/B,0}]$. Consequently, we did not observe any crossing, as before, as can be more clearly seen in the enlargement shown in Figure 6b′. This resulted in the chessboard-type structure displayed by the joint probability density.

Let us now focus on the entangled state described by the wave function (43). As seen in Figure 3i, the initial probability density was only important around $(d/2,-d/2)$ or around $(-d/2,d/2)$. Physically, this meant that either X went through slit A and Y through B, or X through B and Y through B, disregarding any other possibilities. In particular, both subsystems passed through the same slit, which is what the diagonal terms, $(d/2,d/2)$ or around $(-d/2,-d/2)$, describe. The distance between centroids was, thus, larger than those between any two distributions in the factorizable cases ($\sqrt{2}d$ for the entangled state, to be compared with d for the factorizable state). Thus, the effects of the overlapping started becoming apparent at later times (no effect is observed in Figure 3j, compared with the homologous separable cases shown in Figure 3b or Figure 3f). Furthermore, because of the larger distance, the number of interference maxima increased and, hence, the distance between minima decreased by a factor of $\sqrt{2}$. These interferential traits, however, could only be seen along the diagonal joining the centroids, while no interference traits were apparent either along x or y (not, at least, in a prominent manner). This translated into interference suppression in the corresponding subspaces or decoherence, discussed in Section 2.3.

A priori, the above behavior could be considered a minor issue associated with the non-separability of the joint state, since interference was still there, although along the anti-diagonal direction joining the two centroids. However, when we examined the corresponding velocity fields, $v^X(x,y|t)$ or $v^Y(x,y|t)$, displayed in Figure 4g–i and Figure 5g–i, respectively, it was not possible to distinguish unique features associated with one direction or the other since early stages of the time evolution. Note that, when we look at these fields along the anti-diagonal direction, they are dramatically interconnected at any time, without allowing us to clearly separate or distinguish two different contributions, remaining invariant along x or y, as happened with the factorizable cases. Of course, asymptotically we note that both velocity fields developed stable regions flanked by other regions with sudden changes in the velocity, in compliance with the behavior observed in the probability density. Yet, this is not enough to set an unambiguous separability criterion. Accordingly, the markers disseminated in an uneven way everywhere (see Figure 5i), which explicitly indicated the strong dependence that the dynamics of X and Y had on both coordinates, x and y.

If we now consider the trajectories associated with the markers of the entangled state, we find a significant difference with respect to the separable states. The trajectories associated with the X subsystem are plotted in Figure 6c (given the symmetry of the state, the trajectories for the Y subsystem are identical). As can readily be noticed, the renowned Bohmian non-crossing rule breaks down in the reduced subspaces as a consequence of the stretching of the velocity fields along the direction perpendicular to the anti-diagonal. The x-components of the trajectories launched with $y=0$ (solid black line) show clearly how, beyond $t\approx 2$, crossing started, i.e., they passed through the same position at the same time. The y-components, on the other hand, instead of evolving along the same path as one of the x-components, in analogy to the trajectories associated with the two-Gaussian superposition (see Figure 6b′), also started diverting and crossing the same point at the same time (see enlargement in Figure 6c′). In this way, trajectories launched from the same subspace slice (i.e., the same constant value for y, in the case of the X subsystem) could now explore different regions (slices) of the joint space. This is the reason why simply annihilating the interference term in a reduced model does not ensure that trajectories cross the symmetry axis of the setup, as was shown in [39], because the dynamics still contain contextual information on the influence of both slits. In order to find the correct answer, the reach of the information about the non-crossed slit must be rather limited, or gradually suppressed, as shown in [40], in compliance with the evolution here observed for the entangled trajectories.

4. Concluding Remarks

The above discussion poses an interesting question. Typically, in a single-particle two-slit experiment there is no possibility to distinguish which slit was crossed by the particle impinging on the detector. However, Bohmian trajectories provide us with a model, which tells us that there are immiscible streams in the direction of the interference maxima that do not mix. This is possible because there is an underlying velocity field that channels the swarms of particles (probability flux) along given (quantized) directions, leaving voids between every two neighboring ones. Such velocity fields, as we have seen, are well-defined within the standard quantum formalism (i.e., Schrödinger’s picture), as is the link between the probability density and the quantum flux, being two well-defined quantum quantities. This is the reason why the experimental outcomes reported in [32], and the trajectories inferred from them, are consistent. In other words, a trajectory model is perfectly valid to describe quantum fluxes and is, therefore, compatible with quantum mechanics (which does not mean that real particles should be identified with the markers following the trajectories; this still remains an empirically unproved matter).

Now, in the bipartite case, because of the wandering motion displayed by the tracers, also mediated by the corresponding velocity field, the detection of particles associated with the X subsystem is ambiguous, since trajectories may come either from one slit or the other, because the non-crossing rule is not valid anymore within the reduced subspace. This result explains why the experimental data reported in [25] do not disprove the feasibility of a Bohmian picture, rather showing that it is compatible with standard quantum mechanics when the correct assumptions are considered about the (full) system analyzed. Note that, following the reasoning here, this is again a direct consequence arising from the underlying velocity field that relates the probability density with the quantum flux, and not from additional postulates or interpretations. In other words, the trajectory picture provided by Bohmian mechanics is consistent with standard quantum mechanics simply because all the main elements that it makes use of are contained within the latter; the trajectories only constitute an extra refinement that allow us to determine, at a local level, the evolution of the probability flux across the available phase space.

The discussion here has turned around bipartite Young-type interference experiments, considering the particular case of Gaussian states because of analytical convenience. However, we would like to stress that the method is applicable to a wide variety of problems involving bipartite systems, because both the general expressions and the procedure are not constrained to a specific representation. In this way, it does not matter whether the two parties are spatially separated, as is the case here, or, conversely, they are bound by an interaction potential. The only difference with respect to the case here studied is the loss of analyzability, in general, which requires resorting to numerical equations of motion.