Coulomb and exchange interaction effects on the exact two-electron dynamics in the Hong-Ou-Mandel interferometer based on Hall edge states

The electronic Hong-Ou-Mandel interferometer in the integer quantum Hall regime is an ideal system to probe the building up of quantum correlations between charge carriers and it has been proposed as a viable platform for quantum computing gates. Using a parallel implementation of the split-step Fourier method, we simulated the antibunching of two interacting fermionic wave packets impinging on a quantum point contact. Numerical results of the exact approach are compared with a simplified theoretical model based on one-dimensional scattering formalism. We show that, for a realistic geometry, the Coulomb repulsion dominates over the exchange energy, this effect being strongly dependent on the energy broadening of the particles. We define analytically the spatial entanglement between the two regions of the quantum point contact, and obtain quantitatively its entanglement-generation capabilities.


I. INTRODUCTION
In the electronic counterpart of the Hong-Ou-Mandel (HOM) experiment two indistinguishable electrons impinge on the opposite sides of an half-reflecting quantum point contact (QPC), acting as the beam splitter. Contrary to photons 1 , the antisymmetry of the two-electron wave function entails in ideal conditions a zero bunching probability 2,3 . However, in experiments 4,5 the energybroadening of electrons induces a non zero Pauli dip. This has been related to the interplay between the energy selectivity of the QPC and the energetic dispersion of the wave packets 6 , or to the decoherence induced by charge fractionalization 7,8 .
A reliable and robust HOM interferometer is essential in the flying qubit implementation of quantum computing gates based on Hall edge-state interferometers 9 . Indeed, it requires a system that is invariant to small perturbations and robust against typical scattering mechanisms in semiconductor devices, such as with phonons, impurities and background electrons. Experiments 10,11 show that coherent transport of electrons can be achieved in the Integer Quantum Hall (IQH) regime 12 , which is generated by an intense perpendicular magnetic field applied to a confined two-dimensional electron gas (2DEG). This produces edge states, chiral 1D channels following the borders of the device, where an electron can be injected and propagates without being backscattered 13,14 . Edge-channel based nanodevices have been implemented to test electron self-interference in the Mach-Zehnder interferometer [15][16][17][18][19][20][21] , to observe violation of Bell's inequality and two-qubit correlations in the Hanbury-Brown Twiss interferometer 22,23 , and to devise quantum logic gates as quantum erasers or which path detectors [24][25][26] . Moreover, the full quantum tomography of a fermionic particle has been performed by measuring the current noise in a Hanbury-Brown Twiss interferometer 27,28 . Although fractional quantum Hall effect 29 offers many ad-ditional opportunities to exploit coherent quasiparticle interference 30-32 , we limit our present study to the integer regime.
Regarding the electronic HOM effect, recent experiments highlight the presence of charge fractionalization 29,33,34 in the propagation of single excitations in edge channels, so that the coherence of the travelling qubit is not preserved 35,36 . This effect originates at bulk filling factor 2 due to inter-channel interactions that destroy the coherence of the injected Landau quasi-particles [37][38][39] .
As recently proposed, strategies can be implemented to quench this source of decoherence, e.g. the introduction of top gates to loop the second channel 40 or to increase the distance between the two copropagating states 15 . Contrary to the previous scenario, in the present work, we analyze and simulate a device operating at bulk filling factor 1, i.e. only the first Landau level is energetically available 41 . Rather than considering the injection of a steady-state current in the Hall interferometer 5,23,42,43 , we simulate our flying-qubit as enconded in a Gaussian wave packet of edge states 15,16 , with an energy well above the Fermi sea, as recently proposed using quantum dot pumps with time-dependent confining barriers 44 . The choice of a Gaussian state minimizes the computational burden, since its spatial spreading during the propagation is lower with respect to other kinds of excitations, as Levitons [45][46][47][48] . However, we expect that our results do not depend on the shape of the wave packets, at least at a qualitative level, as long as their components transmitted and reflected by the quantum point contact have a similar distribution 16 .
In our approach, the time-dependent wave function is evolved by means of a parallel implementation of the split-step Fourier method 15,16 in a two-dimensional potential landscape, that reproduces the field generated by top gates in the typical GaAs/AlGaAs heterostructure. The energy broadening of the single-particle is directly included in the propagating state and time is an intrinsic variable of our simulations, so that we can access the dynamical properties of the two-particle system. It should be noticed that traditional approaches in the literature bypass the huge computational load for such simulations by using scattering matrices in effective 1D schemes, which proved not to fully capture the interplay between two-electron correlations and the realistic geometry of the device, as for electron bunching. Here, we privilege the exact solution by developing a scalable parallel numerical solver of the time-dependent Schrödinger equation for two particles in a 2D realistic geometry. Moreover, the use of a two-particle Hamiltonian in our simulations allows us to easily include the exact Coulomb interaction of the two charges, in order to explore how electron-electron repulsion interplays with the exchange interaction of the fermionic system.
We initially describe, in Section II, our time-dependent numerical method to simulate two-particle transport in edge channels and the potential landscape we compute to reproduce the HOM interferometer. In Section III(A), we obtain the dynamical bunching probability in presence of exchange symmetry in the two-particle wave function, while in Section III(B) we study the effect of Coulomb repulsion between the two charges. Finally, in Section III(C) we present our measurement of the dynamical spatial entanglement between the two anti-bunched regions of the electronic HOM interferometer. In Section IV, we draw our conclusions.

II. PHYSICAL SYSTEM AND NUMERICAL MODEL
A. Single-particle edge states as a basis We simulate the dynamics of two charges −e with an effective mass m * , that propagate in a 2DEG in the IQH regime. Here, a perpendicular magnetic field B = (0, 0, B) T is included in the single-particle Hamiltonian H(x, y) by using the Landau gauge A = (0, Bx, 0). In presence of a confining lateral barrier V (x), translationally-invariant in theŷ-direction, the singleparticle Hamiltonian is diagonalized by Ψ n,k (x, y) = e iky ϕ n,k (x), where the index n refers to the Landau level, k to the wave vector in the direction of propagation and ϕ n,k (x) diagonalizes the effective 1D Hamiltonian in the transverse direction: with x 0 (k) = − k eB and ω c = eB m * . In proximity of the confining potential (V (x) = 0), the bending of the Landau levels generates conductive channels called edge states. Due to the chirality of the edge state ϕ n,k (x), an electron initialized in it can not be back-reflected by a potential roughness on its path, unless it is scattered by a narrow QPC to the counterpropagating state ϕ n,−k (x) at the opposite side of the confined 2DEG 13,16 . In the electronic HOM experiment, particles are initialized in counterpropagating edge states at opposite sides of a beam splitter, as shown in Fig. 1. In our geometry, only the first Landau level is filled, i.e. only ϕ n=1,k (x) is numerically computed from Eq. (1) by means of LA-PACK routines.
B. Electron as a time-dependent superposition of single-particle edge states To include the energy dispersion of the electron state, edge states with different wave vector k are combined, so that the single-particle wave packet reads where F σ (k, k 0 , y 0 ) = 4 σ 2 /2π 3 e −σ 2 (k−k0) 2 e −iky0 is a Gaussian weight function centered at k 0 = − eB x 0 and with a dispersion σ in theŷ-direction. In order to avoid numerical errors induced by unrealistic step-like potentials, the wave packets are initialized next to confining barriers with a more realistic shape: where F τ (x) = (exp(τ x) + 1) −1 is characterized by the smoothness τ , height V b , and α and β label the two initialization regions of the device (Fig. 1). The indistinguishability of the two wave packets ψ α (x, y) and ψ β (x, y) at t = 0 is ensured by the equivalences τ α = τ β and V α b = V β b , while the opposite direction of propagation is guaranteed by the the symmetry between the confining barriers of the 2DEG (V α (x) = V β (−x)). The counterpropagating states ψ α (x, y) and ψ β (x, y) are therefore characterized by an opposite wave vector of propagation k α 0 = −k β 0 and initial central position y α 0 = −y β 0 , so that they impinge on the beam splitter simultaneously.
C. Numerical solution of the two-particle Schrödinger equation in four spatial dimensions plus time.
By assuming a symmetric spin part of the wave function, the Slater determinant is finally computed from the orbital states ψ α (x, y) and ψ β (x, y), so that exchange symmetry is included in the two-electron state: Together with accounting for the realistic geometry of the device, the explicit computation of the two particle wave function in the 4D configuration space allows one to exactly reproduce the system dynamics in presence of electron-electron interaction V 12 (x 1 , y 1 , x 2 , y 2 ).
However, the memory burden to allocate the 4D wave function increases the computational cost with respect to the single-electron version of our software 15,16 . The present numerical simulations entail indeed distributedmemory parallelization techniques 49 and the support of high-performance computing facilities; to face the problem of memory allocation, we partition the (x 2 , y 2 ) domain of the configuration space between interconnected supercomputing nodes with the MPI paradigm, while a parallel version of the Split-Step Fourier method 50 evolves the distributed two-particle wave function.
The latter time-dependent algorithm has already been validated to study electron transport in single-qubit Hall interferometers in Ref. 15 and Ref. 16 . Here, we use its extension to the two-particle case. In the Landau gauge, A i = (0, Bx i , 0) with i = 1, 2, the two-particle Hamiltonian is split by using Trotter-Suzuky factorization and included in the the evolution operatorÛ (t, 0) for a total  Fig.3(a) at t = 1.2 ps ( ). This provides the geometrical parameter in Eq.(10) for our QPC, a = 60 ± 1 nm. (Bottom) Unconditional probability of particle 1 P (x1, y1, t ) = dx2dy2|Ψ(x1, y1, x2, y2, t )| 2 at different time steps t in the potential landscape of the HOM. The snapshots show that the stationary regime is achieved between 0.9 ps and 1.2 ps. Note that due to exchange correlation, the unconditional probability of particle 1 shows the evolution of both wave packets. evolution time t = N δt: F y1,y2 (F −1 y1,y2 ) and F x1,x2 (F −1 x1,x2 ) are parallel 2D Fourier transforms (antitrasforms) performed on the singleparticle real-space coordinates x i , y i (i = 1, 2). Fourier transforms exploit the locality of the modified kinetic op- , while the single-particle V 1 , V 2 and twoparticle V 12 potential operators are computed exactly in the real-space.
The two-particle evolution in the HOM experiment is then affected by the presence of the beam splitter, to partition the impinging wave packets in a reflected and transmitted component with the same probability in ideal conditions. In this geometry, the beam splitter is realized by a QPC, i.e. a narrow constriction in the potential confining barrier, that partially scatters the traveling charge to the counterpropagating edge channel. We design the QPC by using potential barriers as in Eq. (3), so that the potential profile reads where V b is the height of the confining potential, x L , x R the left and right side of the barrier along thex-direction, and y T and y B the top and bottom side of the barrier along theŷ-direction, respectively. Such parameters of the single-electron potential V (x, y) in Eq. (7) have been tuned to transmit and reflect the impinging singleelectron wave packet in Eq. (2) with 50% probability, as achieved in the Mach-Zehnder interferometer of Ref. 16 . The latter paper shows the operability of the QPC in the IQH regime and the consequence of the distinctive Fermilike energy selectivity of the present beam splitter. In this geometry, if a proper initial position x 0 (k 0 ) of the wave packet ensures an integrated transmission probability of about 50% at the QPC, the Fermi-like energy selectivity of the partitioner induces not-complete overlap between the transmitted and reflected components of the singleelectron wave packet in the k y space. This feature not only affects the visibility in early implementations of the Mach-Zehnder interferometer, but, as shown in the next section, it also has a key role in the apparent violation of the Pauli exclusion principle in the HOM experiment 6 .

III. RESULTS
In the following, we present our exact numerical simulations of the HOM experiment in the IQH regime, with B = 5 T. By considering the typical parameters of GaAs (m * = 0.067m e ) for the hosting material, we tune the spatial broadening of the two indistinguishable electron wave packets (σ = 10, 12.5, 15, 17.5 and 20 nm), and we observe the interplay between the HOM geometry and two-electron correlations in different scenarios, where the exchange and/or Coulomb interaction are included.
The dynamical properties of the system are completely determined from the 4D wave function Ψ(x 1 , y 1 , x 2 , y 2 , t), which is iteratively computed at each time step. As an example, Fig. 2 displays the conditional probability of one particle by selecting the coordinates of the other at two different positions and at two different times during the evolution in the HOM interferometer. We stress that the huge computational load of the 4D wave function requires to reduce the number of degrees of freedom in order to produce viable information, as to measure the dynamical bunching (antibunching) probability P b (P ab ), which corresponds to the joint probability of finding the two particle in the same output (opposite outputs) 6 . The real-space domain (x, y) is divided in two regions, labeled top (T ) and bottom (B) as shown in Fig. 1, so that the configuration space [x 1 , y 1 , x 2 , y 2 ] is partitioned into 4 distinguishable regions labeled S i,j , where the indexes i, j = T, B refer to the first and to the second particle respectively. The joint probability of detection at the same side of the device is while the antibunching probability is similarly computed by integrating over S T B and S BT ; bunching and anti-bunching probability are related by P ab = 1 − P b .

A. Effect of energy symmetry
In a steady-state framework, where the injection of plane waves is assumed, P b is expected to be zero for an half-reflecting beam splitter, due to the exchange symmetry of a fermionic two-particle state Ψ, even neglecting Coulomb interaction 2 . Differently, by using a time-dependent model in a 1D effective framework, it is proved 6 that the energy broadening of the single-electron state affects the bunching probability. This behavior is confirmed by our full-scale time-dependent simulations of the two-electron bunching in 2D for different wave packet sizes, reported in Fig. 3. Here, in presence of exchange symmetry only, the bunching probability of the two-fermion wave function initially varies with time, reaching its maximum when the two particles interact at the QPC. Then, the two charges leave the QPC only partially from different outputs, so that the bunching probability decreases without vanishing. As shown by comparing P b (t) to the evolution of the two-particle probability (bottom panel of Fig. 3), the stationary regime is achieved at t 1.0 ps. By increasing σ, the final bunching probability linearly decreases, so that we expect that the full antibunching is restored in the plane wave limit, i.e. σ → ∞.
The stationary trend of P b (σ) in the exclusive presence of exchange symmetry can be analytically computed by means of the one-dimensional effective model we already validated to study single-qubit interferometers 15,16 . The model relies on the chirality of edge states to map our two-dimensional system in the IQH regime with wave packets of plane waves impinging on a 1D barrier, as schematically depicted in the top panel of Fig. 4. The edge state ϕ n=1,k (x)e iky is replaced by a plane-wave |k , with an effective parabolic dispersion E(k) characterized by a magnetic mass m * B . In the reciprocal space, the single-particle wave functions are defined on the pseudospin basis {|k α , |k β }, where |k α = −|k β = |k . Neglecting Coulomb interaction, the effect of the potential barrier is described by a single-particle scattering matrix: where r(k) and t(k) are empirical equations we provide 16 to compute the energy selectivity of a realistic QPC: The a parameter depends on the smoothness τ and height V b of the QPC in Eq. (7), while γ = 4ln2 in our model. The scattering at the QPC splits the singleparticle wave packet |ψ α(β) in two contributions, that we label |α(β), R for the reflected and |α(β), T for the transmitted component: The scattered single-particle wave packets |ψ α(β) = |α(β), R + |α(β), T are inserted in Eq. (4) to compute the bunched state: where the first ket refers to particle 1, while the second to particle 2. Taking into account the orthogonality of |α(β)R and |α(β)T due to their spatial separation, the bunching probability P b = ψ bun |ψ bun is calculated: We assume the ideal condition α(β)R|α(β)R = α(β)T |α(β)T = 1/2, while the overlap integrals α(β)R|β(α)T are computed exploiting Eq. (10). In case of asymmetry between the initial positions of the two wave packets (y α 0 = y β 0 + ∆y), the bunching probability for a QPC with Σ 2 = σ 2 + a 2 /8γ reads: In the inset of Fig. 3 (top panel), numerical bunching probabilities at t = 1.2 ps are fit by means of Eq. (15) with ∆y = 0. This 1D effective model not only confirms quantitatively the results of our 2D simulations, but it also explains how the energy broadening affects the Pauli dip. According to Eq. (10), the transmitted and reflected components of the single-particle scattered wave packets are centered at different wave vectors, thus decreasing the overlap αT |βR in Eq. (14). Two strategies are therefore possible to reduce P b : the manipulation of the spatial broadening σ, as shown by Fig. 3, or a variation in the smoothness of the QPC. The former case is realized when σ → ∞ (plane-wave limit), while the latter if a → 0 (flat energy selectivity) in Eq. (8). Regarding the smoothing of the energy selectivity, a different geometry, as the beam splitter at bulk filling factor 2 in Ref. 15 , represents a possible solution.
We additionally remark that in Eq. (14) the overlap −2| αT |βR | 2 decreases the bunching probability due to the antisymmetry of the two-electron state. In case of a symmetric wave function the same term is expected to be positive (+2| αT |βR | 2 ), as shown in the bottom panel of Fig. 4. Here, we compare the dynamical bunching and antibunching probabilities of a symmetric/antisymmetric two-particle wave function that propagates in our geometry of the HOM interferometer: in the stationary regime, the two configurations are characterized by exchanged values of P b and P ab . Differently, as indicated by Eq. (15), the symmetry of the two-particle wave function is not expected to affect P b at ∆y → ∞. The stationary bunching probability for different initial displacements ∆y and σ = 15 nm is reported in Fig.5(a) (blue dots for V 12 = 0). In our time-dependent scenario, a spatial mismatch in the initial position of the two electron states (y α 0 − y β 0 = ∆y = 0) modifies the time at which each wave packet is expected to impinge on the QPC, so that the overlap between the two wave packets decreases, as well as the effect of their exchange symmetry. The classical probability of joint detection for two distinguishable particles is gradually restored with ∆y → ±∞. Eq. (15) shows that the characteristic length for this process is λ = 2 σ 2 + a 2 8γ ≈ 30 nm for σ = 15 nm, as confirmed by the fitting curve P b (∆y) (blue solid line) in Fig.5(a), which vanishes at 3λ 90 nm.

B. Effect of Coulomb interaction
If our single-particle approach provides enough information to understand electron bunching probability in presence of exchange-symmetry only, the introduction of Coulomb interaction requires a full-scale two-particle approach. We introduce the electron-electron repulsion between the two carriers by adding to the two-particle Hamiltonian the potential where is the medium permittivity and d accounts for the finite thickness of the 2D system in the divergence at To evaluate the interplay between exchange symmetry and Coulomb interaction, we initially compute the bunching probability for a product state Ψ C = ψ α (x 1 , y 1 )ψ β (x 2 , y 2 ), where the antibunching is exclusively generated by electron-electron repulsion.
In Fig. 5(b) we report the bunching probability at final time t = 1.2 ps with σ = 20 nm for different values of the d parameter. Note that, for a Coulomb repulsion small enough (d = 100 nm), P b reaches the limit 1/2. Indeed, in absence of the exchange symmetry, the two wave packets evolve independently and the bunching probability in Eq. (14) does not contain the overlap αT |βR . P b equals 1/2 only if α(β)T |α(β)T = α(β)R|α(β)R = 1/2, namely the two initial wave pakcets are properly initialized, so that they are transmitted with 50% probability by the QPC, as in the present case. Additionally, Fig. 5(b) shows that the Coulomb-driven scattering at the QPC does not produce perfect antibunching when d → 0, but rather P b saturates by decreasing the d parameter. This differs from a one-dimensional system, where the two electron are confined on the same rail, e.g. in thex-direction. In the latter case, they are forced to get across x 1 = x 2 , so that they feel an effective infinite barrier during the scattering for d = 0. In the two-dimensional geometry, alternative paths with a finite barrier are present, and each charge is only partially reflected by the Coulomb potential for d = 0, producing partial bunching. Fig. 5(c) shows the evolution of P b (t) for the separable case with different values of σ at d = 1 nm. Finally, we add the exchange symmetry to the interacting system ( Fig. 5(d)). In our operating regime, the Coulomb repulsion dominates on the exchange interaction, which only shifts the bunching probability to lower values. The effect of the Coulomb interaction turns out to depend on the spatial broadening of the wave packet, as visible in both separable (Fig. 5(c)) and non separable (Fig. 5(d)) interacting scenarios. With contrast to the stationary P b (σ) in the inset of Fig. 3(a), in presence of electron-electron interaction the decrease of P b (σ) is not linear, but it rather saturates with increasing σ. We can interpret this as a consequence of the broadening in the charge distribu-tion, which reduces the effect of the Coulomb interaction for larger wave packets. Additionally, Fig. 5(a) provides the simulated stationary probability for different initial displacements ∆y in presence of Coulomb repulsion and σ = 15 nm: numerical data (green squares) are fit by ef f +a 2 ef f /8γ) )] (green dashed line), which corresponds to Eq. (15) with an effective σ ef f and a ef f = 2 √ 2γ Σ 2 ef f − σ 2 ef f used as fitting parameters. The numerical fit provides σ ef f 21 nm and a ef f 86 nm, that are larger than the actual ones (Fig.3 inset) without Coulomb interaction.

C. Spatial entanglement
We finally provide a dynamical estimation of the spatial entanglement for the antibunched configuration in presence of exchange and/or Coulomb interaction. The Hilbert space is divided into two separate subsystems, the top (T ) and the bottom (B) region in Fig. 6, so that H = H T ⊗H B . Due to the indistinguishability of the two electrons, entanglement is measured by assuming particle 1 in the top region ((x 1 , y 1 ) ∈ T ), and particle 2 in the bottom one ((x 2 , y 2 ) ∈ B), while the opposite configuration is equally entangled due to the symmetry of the two wave packets. To reduce the computational burden, we exploit the chiral nature of edge channels by projecting the two-particle wave function on 1D slices of the subsystem: y 1 and y 2 coordinates are fixed at the expected maxima of the scattered single-particle wave packet in T and B regions, in the present case y * 1 = 11.5 nm and y * 2 = −11.5 nm (black dashed lines in the inset of Fig. 6). The projected two-particle wave function, φ(x 1 , x 2 ) = Ψ(x 1 , y * 1 , x 2 , y * 2 ) is renormalized, and the density matrix reads: whose Von Neumann entropy does not depend on the subspace chosen to be traced out. Therefore, we compute the reduced density matrix on the lattice points of our domain B. We calculate the spatial entanglement by means of the Von Neumann entropy Fig. 6(a) shows the dynamics of the Von Neumann entropy of the wave-function spatial distribution in the T and B regions with exchange symmetry only in the antibunching configuration. The entanglement is quenched by an increase of the spatial distribution of the wave packet, so that we expect it to vanish in the plane-wave limit. In Fig. 6(b) we compare the stationary Von Neumann entropy for the cases of: (i) exchange symmetry alone (red solid line), (ii) a separable wave function with Coulomb interaction (green dashed line), (iii) Coulomb interaction and symmetric wave function (blue dotted line). In the latter, the Coulomb repulsion, which acts as an additional barrier, further prevents the two particles to reach the opposite regions thus damping entanglement. The Von Neumann entropy in the interacting scenario with exchange interaction does not differ from the one in the distinguishable case: the two electrons are already prevented to occupy the same coordinates due to the infinite barrier represented by V 12 at (x 1 , y 1 ) = (x 2 , y 2 ).

IV. CONCLUSIONS
We realized a full-scale numerical simulation of an electronic HOM interferometer in the IQH regime that contributes to shed light on the apparent violation of Pauli exclusion principle in two-electron bunching, by including exactly the interplay between a realistic geometry of the QPC in 2D and the exchange and correlation of the two-particle wavefunction. A full understanding of such interplay is important for the implementation of a HOM device for quantum computing protocols, e.g. to measure the degree of indistinguishability of electrons generated from different sources 11 .
In contrast to traditional approaches in literature, that are based on 1D effective models implemented in stationary frameworks and with effective scattering matrices for the QPC, we privilege the exact solution of the 4D timedependent Schrödinger equation for two particles in a 2D real space. The exact solution requires a parallel implementation of our numerical solver by means of the MPI library. Thanks to this effort, we provide dynamical measurements of the bunching probability for a two-particle system, where the electrons are initialized in Gaussian wave packets of the exact edge states confined at the barrier, coherently to the most recent single-electron injection protocols in Hall nanodevices. The time evolution of the two interacting and correlated wave packets successfully reproduce the two-electron dynamics in the present interferometer. The decrease in the bunching probability by increasing the spatial localization confirms, for a full scale 2D spatial geometry, the findings of Ref. 6, where this effect is explained with a 1D time-dependent model as a signature of the non-orthogonality between the states scattered by the potential barrier. In addition to the fullscale numerical simulation, we provide a simplified analytical model to relate the stationary bunching probability to the non-perfect overlap between the transmitted and reflected wave-packets from the QPC. This model clarifies the interplay between the spatial dispersion of the wave packet σ and the geometry of the QPC, which is encoded in the single-particle parameter Σ of Eq. (15). We show how the perfect antibunching is recovered in the plane-wave limit and point out the role of exchange symmetry by simulating the HOM experiment both for a symmetric and an antisymmetric wave function. As an additional advantage in treating exactly the two-particle scattering, we include electron-electron repulsion in the Hamiltonian to evaluate the interplay with the fermionic statistics. We observe that in a 2D real-space geometry, differently from the typical 1D scenario adopted in litera-ture, the bunching probability does not vanish for an infinite repulsive Coulomb interaction, together with a nonlinear dependence on the spatial broadening σ. Finally, a dynamical measurement of the spatial Von Neumann entropy between the top and bottom regions of the device allows us to assess the spatial entanglement between the two antibunched carriers; we found that long-range Coulomb interaction quenches the entanglement by enhancing the Pauli dip with respect to the antisymmetry alone. This study represents the starting point for the simulation of more sophisticated 2D geometries, as the conditional phase shifter 51,52 , where the Coulomb interaction entangles electrons in different edge states.