Dynamics and Hall-edge-state mixing of localized electrons in a two-channel Mach-Zehnder interferometer

We present a numerical study of a multichannel electronic Mach-Zehnder interferometer, based on magnetically-driven non-interacting edge states. The electron path is defined by a full-scale potential landscape on the two-dimensional electron gas at filling factor two, assuming initially only the first Landau level as filled. We tailor the two beam splitters with 50% interchannel mixing and measure Aharonov-Bohm oscillations in the transmission probability of the second channel. We perform time-dependent simulations by solving the electron Schroedinger equation through a parallel implementation of the split-step Fourier method and we describe the charge-carrier wave function as a Gaussian wave packet of edge states. We finally develop a simplified theoretical model to explain the features observed in the transmission probability and propose possible strategies to optimize gate performances.


I. INTRODUCTION
The concrete implementation of quantum information devices is facing a notable development, mainly based on superconducting 1 and single ion 2 qubits. Alternative approaches based on electronic states in semiconductor devices seem also to be particularly promising due to their scalability and their potential to be integrated with traditional electronic circuitry. However, decoherence represents a major problem for semiconductor devices due to the existence of several scattering sources for electrons in solids, as phonons, impurities, and electron-electron interactions. Specifically, for a flyingqubit implementation 3-5 of a quantum gate, the onset of environmental interactions would destroy the coherence of the traveling electron wave packet (WP) on very short timescales.
Topologically protected edge states (ESs) are able, in principle, to prevent the loss of coherence of the electron state by embedding it in a subspace that is invariant to small perturbations and is robust against the above scattering mechanisms 6 . For this reason, single electrons in ESs are emergent candidates for the implementation of quantum logic gates 7,8 . The most notable example of such states consists of a two-dimensional electron gas (2DEG) subject to an intense transverse magnetic field driving the system into the integer quantum Hall (IQH) regime. In this case, ESs are one-dimensional chiral conductive channels localized at the boundaries between allowed and forbidden regions for the free conduction-band electrons, where the latter can propagate for long distances (larger than 10 µm) 9,10 without being backscat-tered, due to the chirality of the channels. Thus, the IQH regime is an ideal platform for electronic interferometry aimed at quantum information processing, which however requires the realization of semiconductor nanodevices able to manipulate edge channels. Due to the analogy with the corresponding optical systems, this class of systems is often termed electron quantum optics devices. Several examples of the latter have been realized experimentally, such as Mach-Zehnder interferometers (MZI) 11,12 , Fabry-Pérot interferometers [13][14][15] , Hong-Ou-Mandel interferometers [16][17][18][19][20] and Hanbury Brown-Twiss interferometers 21,22 , and their application as quantum erasers 23 or which-path detector 24 has been tested. Numerical simulations based on stationary-state 9,25,26 or time-dependent 8,27,28 approaches have been essential to understand the experiments, but a time-dependent modeling of a whole IQH device, aiming at the proposal of a quantum gate, is lacking.
A new promising architecture for a multichannel MZI has been proposed in Ref. [7]. Whereas previous MZIs were mainly based on counterpropagating channels 29 and the typical Corbino geometry 14 , this device is characterized by a smaller loop area, which reduces the effects of phase-averaging, and, most important, strong scalability, which allows to concatenate in series a number of devices. The system under study operates at filling factor two, contrary to the previous proposal presented in Ref. [8], where the device was operating at filling factor one, with a single channel reflected/transmitted by a quantum point contact. Indeed, numerical studies show that it is possible to realize coherent superposition of edge channels with sharp potential barriers 30 , and that the interchannel mixing coefficient can be arbitrarily tuned by using arrays of top gates 9,31,32 . This mechanism has been applied experimentally to spin-resolved edge states 32,33 , with an additional in-plane magnetic field to couple the two spin channels. However, an essential drawback of this idea is the large spatial extension of this beam splitter (BS) and the difficulty in fine-tuning the device operation. Instead of using spin-resolved ESs, our research focuses on a multichannel MZI where the two non-interacting copropagating channels belong to different Landau levels (LLs), and the formation of the qubit does not require a resonant condition 9 . As a consequence, in the present device the length in which the two channels need to run on the same edge (i.e. the region at filling factor 2) is limited to a small BS region, as detailed in Appendix A. Additionally, we suppose to encode the qubit in a propagating WP of ESs 8 , with a Gaussian shape, whose injection protocol for quantum dot pumps has been recently proposed in Ref. [43]. This choice corresponds to the experimental situation where a quantum dot pump 36,37 injects the single-electron WP in the ES of an interferometer, with an energy well above the Fermi energy of the device. Though a number of implementations of ES interferometers are based on Lorentzian or exponential WPs 34,35,41,42 , we found that, in a noninteracting picture, the qualitative results of our simulations are valid also for alternative shapes of the initial wave function. To be more specific, in Appendix B we show that our approach applies also to the case of a WP with a Lorentzian distribution in energy 34 .
In order to solve the time-dependent Schrödinger equation, we use the split-step Fourier method together with the Trotter-Suzuki factorization of the evolution operator 28 , with a parallel implementation of the simulation code 38 . Specifically, we study the real-space evolution of the particle state, observing the dynamics of a carrier inside the MZI. We additionally perform support calculations with the Kwant software 39 , which solves the scattering problem in a steady-state picture for a single energy component of the WP. After optimizing the device and the performance of the quantum gate operation, we measure the transmission probability from the first to the second channel. We vary both the length mismatch of the two paths, defined by the width of the mesa W , and the orthogonal magnetic field B, to observe Aharonov-Bohm (AB) oscillations 11,40 in the transmission amplitude. We finally relate the variations of transmission probability in the two outbound channels to the device geometry and compare exact numerical results to a simplified theoretical model based on the scattering matrix formalism.

II. PHYSICAL SYSTEM
In our simulations, a conduction-band electron with charge −e and an effective mass m * moves in a 2DEG on the xy-plane and it is immersed in a uniform magnetic field B = (0, 0, B) along the z-direction. We describe the effect of the magnetic field on the charge carrier in the Landau gauge A = B(0, x, 0), that simplifies the definition of the initial state moving along the y-direction. The potential modulation induced by a polarized metallic gate pattern on the 2DEG is reproduced by the local potential landscape V (x, y) reported in Fig. 1. The Hamiltonian of a conduction-band electron results to bê In order to solve the time-dependent equation of motion, we initialize the electron in a region of the device where V (x, y) = V (x), i.e. where the Hamiltonian shows translational symmetry along the y-direction (region I in Fig. 1), and its eigenstate can be factorized as ϕ n,k (x)e iky . Following the standard description of the IQH effect, the exponential term e iky describes a delocalized plane wave along y, while the function ϕ n,k (x) is the eigenstate of the 1D effective Hamiltonian where ω c = − eB m * is the cyclotron frequency and is the center of a parabolic confining potential depending on the wave vector k along the y-direction. The discrete eigenvalues E n of the effective Hamiltonian (2) are the LLs. If the potential V (x, y) is not present, the LL energies do not depend on k, while the system eigenfunctions ϕ n correspond to the Landau states. On the contrary, the presence of a confining step-like potential V (x) modifies the LL band structure introducing a dispersion on the wave vector, such that E n = E n (k). In detail, the shape of the eigenstates ϕ n,k (x) changes significantly when the center of the parabolic confinement x 0 approaches the nonzero region of V (x). For a fixed k, the state ϕ n,k (x)e iky is a current-carrying ES characterized by a net probability flux in the y-direction.
Our time-dependent approach (described later in Section III) requires to go beyond the delocalized description of the wave function. In fact, we choose for the initial state a specific value of the n index, and we combine linearly the corresponding ESs on k in order to form a minimum-uncertainty WP. From a computational perspective, our choice of a minimum-uncertainty WP as the initial state avoids numerical instabilities and minimizes the real-space spreading of the wave function. Specifically, the particle is initialized in region I of Fig. 1 in the first edge channel (n = 1) and is described by where the weight function F (k) is the Fourier transform of a Gaussian function along y and it entails the wave vector localization around k 0 . The Gaussian envelope defines a finite extension around the central coordinate y 0 , while the localization around x 0 (k) is defined by the function ϕ 1,k (x). Unlike our previous work 8 , the energy dispersion of the WP includes the energies of the first two edge channels (n = 1, 2) so that, in principle, both can be occupied, even though only the first one is initially filled, as indicated by Eq. (4). The transition between low-potential and highpotential regions in the x-direction at fixed y occurs at x b , and we model it by a Fermi-like function. In particular, in the initialization region (region I in Fig. 1) the potential is assumed to depend only upon x, according to the expression where F τ (x) = (exp(τ x) + 1) −1 is the Fermi distribution with a smoothness given by the broadening parameter τ , and V b represents the energy of the forbidden region. Taking the potential of Eq. (6), the eigenstates ϕ n,k (x) of the effective Hamiltonian are computed numerically. Moving forward along the positive y-direction, the potential profile assumes also a dependence on y, such that the two edge channels, whose paths are defined by V (x, y), constitute an MZI. On the border between region I and II in Fig. 1, the WP impinges on a sharp potential dip, which acts as a BS and redistributes the wave function on the first two available channels n = 1, 2. Then, the potential mesa in the middle of region II forces the two channels to follow different paths that accumulate a relative phase. Proceeding further, between region II and III, a second BS produces the interference between the two parts of the wave function. In region III and IV, we introduce an additional mesa and an imaginary potential, respectively, as a measurement apparatus to remove the electron probability from channel n = 1 alone. As a consequence, the norm of the final wave function represents the total transmission probability of the interferometer from the first to the second channel, P tot 21 .

III. NUMERICAL SIMULATIONS
As previously observed, our time-dependent numerical simulations model the evolution of a localized WP representing the propagating carrier. Our method allows to directly observe the dynamics of carrier transport in the time domain and to assess the effects of real-space localization on it. This approach does not require the diagonalization of the Hamiltonian of the whole device, which can be a very demanding task for such a large system. Indeed, we only perform the diagonalization of the 1D effective Hamiltonian in Eq. (2) with the addition of the confining potential V (x) in region I.
Once the particle is initialized, we solve the timedependent Schrödinger equation by using a parallel implementation of the split-step Fourier method, based on the recursive application of the evolution operator U (δt) = e − i Ĥ ·δt to the initial wave function Ψ I (x, y; t = 0): The kinetic and potential contributions to the Hamiltonian of Eq. (1) can be split in three parts: where the kinetic terms are defined bŷ Then, we use Trotter-Suzuki factorization to split the evolution operator, separating the kinetic and potential contributions. In order to exploit the diagonal nature ofT 1 andT 2 on the reciprocal space and of V (x, y) on the real space, we apply alternated Fourier transforms F x(y) and anti-Fourier transforms F −1 x(y) along the x(y)direction. The evolution operator assumes finally the following form: The split-step Fourier method requires a careful choice of the small time step δt. In particular δt where ∆ x(y) and v are the real-space grid spacing in the x(y)direction and the group velocity, respectively. Furthermore, to avoid aliasing effects, δt V . Consistently with the previous requirements, we select an iteration time δt = 10 −16 s. We take the initial state of Eq. (7) with n = 1, σ = 60 nm and centered at x 0 = −50.9 nm, y 0 = −800 nm. We consider GaAs parameters for the hosting material, namely m * = 0.067m e . Furthermore, we use a 2048 × 4096 simulation grid. Numerical simulations are performed on a domain including the whole device to study AB oscillations of the transmission amplitude P tot 21 , while a reduced domain is used to study each component of the MZI and optimize gate performances, as reported in the following.

A. Beam splitter
The BS must scatter coherently a particle WP initialized in one of the available channels to fill both LLs and leave the electron in a coherent superposition of the two outgoing channels. In order to produce the highest visibility of the interferometer, we tune the BS functionality to obtain a 50% mix. Numerical simulations based on delocalized plane-waves 30 show that a coherent edge mixing can be achieved by introducing spatial inhomogeneities on a scale smaller than the magnetic length l m , on the path of the ES. Indeed, an abrupt potential profile scatters elastically an impinging plane wave and redistributes the incoming wave function on the available states (the first two LLs in the present case), with a transmission coefficient t BS f,i from the initial i = 1, 2 to the final f = 1, 2 channels. t BS f,i depends on the energy of the incoming state, on the value of the magnetic field B and on the shape of the local potential.
Regarding our system, the above mechanism, which is represented in Fig. 2(a), is valid for each wave vector component of the particle WP, whose energy distribution is conserved along the whole device. Note that, however, the weight function F (k) depends on the local dispersion of the LLs. In particular, we aim at realizing an edgechannel superposition with equal probabilities, thus requiring a potential profile which ensures a constant transmission probability |t BS 12 (E)| 2 =|t BS 21 (E)| 2 =0.5 for each energy component E of the initial WP. We achieve such result by using the potential profile shown in Fig. 2(b). Differently from proposals based on spin-resolved ESs 31 , no resonant condition is required. Specifically, our BS consists of a square with the corners smoothed by Fermi profiles: where τ BS is the smoothing parameter. In order to evaluate the energy dependence of t BS if (E), we use the wave   Figure 3(a) also reports the energy broadening of the initial WP for a particle initialized in the first (red solid line) and second (blue dashed line) LL. We finally measure the total transmission probabilities simulating the scattering process at the BS with our time-dependent approach. Results are reported in Fig. 3(b) for the WP initialized in the first LL and in Fig. 3(c) for the WP initialized in the second one, at B= 5 T. A small scattering to the third LL, whose energy is slightly reached by the energy broadening of our initial WP, explains the discrepancy between the sum of the two scattered intensities and unity.
Finally, we perform support calculations with Kwant software 39 , simulating delocalized ESs impinging on the BS. The scattering matrix method is used to calculate the maps of Fig. 4, where the probabilities |t BS f,i | 2 are reported also as a function of the magnetic field B. The latter results confirm the transmission probabilities of Fig. 3(a) obtained with the time-dependent method and shows how B tailors the transmission coefficients.

B. The MZI
Once the coherent superposition is realized, the MZI requires that the two channels accumulate a relative phase. This can be induced by a mismatch of the path lengths or by a net flux of the magnetic field through the loop area, which is the area enclosed by the paths of the two channels. To separate the channels, we introduce an area where the potential V (x, y), which mimics the landscape of polarized top gates, has an energy value V s in between the first and the second LL. In order to avoid an unwanted mix of the two channels and to better model a real device, we create a smooth transition between the two regions by means of the following function: The smoothness of the local curvature τ s must ensure an adiabatic separation of the two edge channels 9 , with a negligible mixing among them. This creates a region (lighter blue in Fig. 1) where the filling factor is one, in contrast to the bulk filling factor of two. From a different perspective, in order to split the two channels, we exploit the relation between the real coordinate x 0 , defining the center of the WP along x, and the momentum k of the traveling particle along y, as given by Eq. (3). Indeed, the band structures of the LLs are strictly related to the shape of the potential profile, as shown in Fig. 5(a) for the region I and Fig. 5(b) for a section of the mesa structure in region II. In detail, in Fig. 5(b) it is clear that the potential step pushes upwards the local band structure, and the two LLs are then filled at different k. The elasticity of the scattering process at the mesa ensures that the first LL is filled on top of the step potential, while the second LL intersects the energy window at its bottom. The channels are therefore forced to follow a different path, whose length can be tuned by changing the width of the mesa W . The simultaneous recollection of the WPs at the second BS, which is needed to observe the interference, could be prevented by the different group velocity of the WPs in the two edge channels. Indeed the group velocity of the first channel is larger than the group velocity of the second one due to the different band structures of the two LLs. Therefore, we introduce a sort of indentation in the forbidden region on the mesa (region II in Fig. 1), in order to increase the length of the channel n = 1 and compensate this effect. Additionally, we smooth the local confining potential inside the indentation, in order to reduce the group velocity of the WP in n = 1.
Finally, the regions III and IV of the device correspond to the measurement apparatus. After the interference, the two channels are separated by an additional mesa in region III. In order to remove from the device the part of the wave function occupying the first LL after the MZI, we introduce the absorbing imaginary potential where y c defines its center, d its length, V 0 abs its maximum, and x a , x b define its spatial extension in the xdirection. This potential is represented by the gold shape in region IV of Fig. 1 at y = 500 nm and models a metallic absorbing lead on the path of the first LL. Consequently, the surviving part of the final wave function gives the probability for the electron to be transmitted in the second LL by the interference process taking place inside the device. Using the split-step Fourier method, we finally simulate the interference for different values of the orthogonal magnetic field B at W = 200 nm, and for different widths of mesa W at B = 5 T, modifying the magnetic and the dynamic phase respectively. Numerical simulations have been performed considering V b = 0.031 eV, τ b = 0.25 nm −1 for the confining potential, |x 1 − x 2 | = |y 1 − y 2 | = 20 nm and τ BS = 0.5 nm −1 at the BS, V s = 0.011 eV, τ s = 0.2 nm −1 for the mesa structure and V 0 abs = −100 eV, d = 30 nm for the absorbing potential. The numerical results are reported in Fig. 6. We observe AB oscillations in the transmission amplitude with an high visibility, defined as thanks to the optimization of the scattering process at the BS. Before discussing the results, in the following section we propose a simplified theoretical model whose predictions will help in understanding the outcomes of the exact time-dependent approach.

IV. THEORETICAL MODEL
Here we present a theoretical model based on the description of edge channels as strictly one-dimensional systems, using the scattering matrix formalism. An ES of the n th LL is represented by a plane wave along y, |k, n , with the energy dispersion of that LL, k(E, n). In order to introduce particle localization on the y-direction, our initial wave function is computed by combining different ESs of the n = 1 level, with the Gaussian weight F (k) of Eq. (5): where |E, n denotes |k(E, n), n for brevity, and |Ψ I (|Ψ III ) is the one-dimensional wave function in region I (III). We assume a bulk filling factor of two, so that n can be either 1 or 2 and represents a pseudo-spin degree of freedom. The WP in region III can be related to the initial one by describing the scattering process through the application of three operators: whereB describes the effect of a BS, andΦ the relative phase accumulated by the two channels in the mesa region. Here, differently from the full numerical simulation of the previous sections, the energy-dependence ofB and Φ is neglected for simplicity. Finally, since the absorbing potential in region IV collects the contribution of the first LL, only n = 2 survives, and the total transmission probability at the end of the device is defined by the following equation: In order to solve Eq. (17), we consider the general 2x2 matrix form of operatorsB andΦ on the pseudo-spin basis:B The phase ϕ i (i = 1, 2) includes the contributions of the magnetic (φ i ) and the dynamical (ξ i ) phases where the integration is performed along the path of the edge channel i on the mesa. The transmission coefficients b ij are related by the probability flux conservation: |b ii | 2 = |b jj | 2 .
As in the previous sections, we tune the BS to 50% transmission, so that all the coefficients |b ij | 2 = 0.5. Therefore, b 11 and b 22 only differ by a phase factor ϕ, such that In order to define a gauge-independent dynamical phase, we consider a quasi-linear dispersion of the two LLs around the central energy of the WP E 0 , and we rewrite p in terms of the constant energy E and of the group velocity v II i in region II: where ∆S i is the length of the path of channel i in the mesa region. Note that, while considering a linear dispersion is appropriate for the second LL, it represents an approximation for the first one. Such assumption is the main source of discrepancy between our exact numerical results and the present theoretical model. Using the Stokes theorem for the magnetic contribution φ i of Eq. (19), we can rewrite the total phase as with A the area enclosed by the paths of the two channels, which is tuned by changing the width W of the mesa along the x-direction. Performing the integration over the energy in Eq. (17), the total transmission probability from channel 1 to channel 2 is where the argument of the cosine Φ = eBA + ϕ exposes the dependence of P tot 21 on the magnetic field B and on the width W of the mesa. Indeed, according to the geometry of the step potential in Fig. 1, the mesa has an area A = W · L + (2δ y + L) · δ x , such that the two following definitions of Φ hold: where Φ 0 = ϕ + eB (2δ y + L)δ x and Φ 1 = ϕ. Besides, according to Fig. 1, the paths of the two channels are equivalent to ∆S 1 = 2δ y +2W +L and ∆S 2 = 2δ x +2δ y + L, and using an effective standard deviation Σ = σv II 1 /v I 1 , the total transmission probability is with W 0 containing the geometrical correction to the paths of the two edge channels.

V. DISCUSSION
The AB oscillations simulated numerically are compared to the transmission probability P tot 21 of Eq. (29) predicted by our theoretical model. In detail, the numerical data are fit by the function for a variation of the magnetic field B [ Fig. 6(a)], and by the function  Table I. Comparison between fitting parameters for the results of exact numerical simulations and the corresponding parameters of the theoretical model of Sec. IV. The two cases of Fig. 6 are considered, namely with a variable magnetic field (left column) or mesa width (right column).
for a variation of the width W of the mesa region [ Fig. 6(b)]. The comparison between numerical and theoretical parameters is presented in Tab. I. Regarding the magnetically-driven AB oscillations in Fig. 6(a), we observe that the shape of the interference curve does not describe a perfect sinusoid, but the amplitude slightly increases with the magnetic field. Indeed, an increase of B enhances the spacing between the two LLs, reducing the unwanted interchannel mixing at the step potential, therefore increasing the oscillation visibility. Additionally, our theoretical model neglects the dependence of the transmission coefficients b if on B. Fig. 4 shows indeed that an increase of the magnetic field increases the scattering from the first to the second channel at the BS, affecting the values of A B and A * B in Eq. (30). The underestimation of the pseudo periodicity k B is induced by the approximation of the loop area A, which doesn't take into account the small difference in the x position of the two channels also in the regions with filling factor two.
The amplitude of P tot 21 (W ) in Fig. 6(b) has a damping induced by the relative dynamical phase together with the finite dimension of the wave function. Indeed, when the width W of the mesa is large enough, the two WPs do not overlap anymore and the interference is quenched 44 . Such damping was observed also in the single-channel MZI 8 , but in the present device Σ is reduced with respect to the standard deviation of the initial WP, σ. In fact, in this two-channel MZI, the smoother slope of the indentation in region II reduces the group velocity v II 1 above the mesa with respect to v I 1 . This can be interpreted as an effective dilatation of the width W in Eq. (29), determining a larger phase difference. Moreover, as shown in the inset of Fig. 6(b), the Gaussian fit describes properly the oscillation amplitude of P tot 21 only in the central region, while on the two sides the AB oscillations are larger than the predicted ones.
We measure a visibility ν = 0.87 at W = W 0 in place of 1, as a consequence of the energy dependence of the phase factors and of the scattering processes inside the device. In particular, in addition to neglecting the energy dependence of the transmission probability at the BS, our theoretical model does not take into account the unwanted interchannel mixing induced by the step encountered by the second LL when entering the mesa region. The mix is actually non zero, and it depends on the energy of the impinging WP. We expect that the high-energy components of the wave function in the second LL [top of the orange striped zone in Fig. 5(b)] are transferred more easily to the states of the first LL with the same energy and a higher group velocity, leaving sooner the scattering region.
In summary, in this paper we have investigated the transport properties of a Gaussian electronic WP in a two-channel MZI in the IQH regime. Our numerical modeling of the device required the definition of a proper potential landscape V (x, y) to ensure a high visibility of the transmission amplitude. A specific design of the BS has been used to separate the impinging state into a 50% coherent superposition of the two available channels. However, we found that the proper function of the BS is preserved when different shapes of the mixing potential are used, as we show in the Appendix B. We observed AB oscillations, relating the features of transmissionprobability amplitude to particle localization, which is inherent in our time-dependent solution. Finally, our numerical results are clarified by a simplified theoretical model based on the scattering matrix formalism and a one-dimensional model for chiral transport in edge states. We emphasize that this implementation of an MZI solves the scalability problem 7 of the single-channel MZI we studied in Ref. [8], thus potentially enabling its concatenation in series and its integration into sophisticated quantum computing architectures. The possibility to concatenate two or more MZI in series, exploiting as an input the two possible outputs of a previous interferometer, is essential for the implementation of two-qubit interferometers, as the Hanbury-Brown-Twiss one 21 , where interfering identical Gaussian WPs could be, in principle, generated from nonidentical sources 43 . In addition, the present device shows a larger visibility with respect to our previous single-channel interferometer 8 , mainly due to the weak energy selectivity of the present BS compared to the quantum point contact. Moreover, our BS does not require the resonant condition of the spin-resolved multichannel MZI proposed in Ref. [32], thus reducing the interchannel interaction induced by the spatial extension of the top gate array (Appendix A).
The functioning of our MZI does not depend on the specific shape of the WP. The choice of a Gaussian weight function is motivated by the higher control of its time evolution with respect to alternative shapes. For the sake of completeness, here we show the evolution of a WP with a Lorentzian distribution in energy, in order to mimic the emission of electrons by a mesoscopic capacitor 34 . Fig. 8 shows the initial broadening of the WP in energy and real space: the two long tails of the Lorentzian distribution produce a small filling of the states with no velocity and collect a very large number of wave vectors, thus inducing a larger spread of the WP during its evolution. Additionally, due to its very small energy peak, the wave function in real space has a long tail, that required to double the dimensionality of the initialization region. Fig. 9 shows its evolution at different time steps: the initial beam in the first edge channel (t = 0 ps) is split in a coherent superposition of the two channels by the first BS (t = 2 ps), than the mesa structure separates the component with different n (t = 10 ps), and finally the WPs are recollected at the second BS (t = 20 ps) to realize the interference. Numerical results confirm that our device is still fully operational in this case.
Finally, we present some support time-independent simulations performed with alternative shapes of the BS. We modeled a triangular and rectangular potential dip, whose profiles are reported in Fig. 10(a) and Fig. 10(b), respectively. In Fig. 10(c)-(d), the two BSs are inserted in a simple device with the leads of injection and absorption, in order to show that they produce a coherent superposition of the first and the second channel. As in the previous case, the central energy of the WP can be chosen to obtain a 50% scattering probability between the two channels. We found that for both rectangular and triangular potential dips the scattering probability from the first to the second channel computed with Kwant software shows a small variation, around 5%, for an energy dispersion of 0.2 meV, which is comparable to the energy uncertainty usually obtained in experiments 19 .     Fig. 7) to show that a coherent superposition of the first and the second channel is formed. In both cases, the beam is initialized in the first edge channel at E = 20.4 meV and it is scattered to the second channel with about (50±4)% probability by the (c) triangular potential dip and (d) rectangular potential dip.