Massively-multiplexed generation of Bell-type entanglement using a quantum memory

High-rate generation of hybrid photon-matter entanglement remains a fundamental building block of quantum network architectures enabling protocols such as quantum secure communication or quantum distributed computing. While a tremendous effort has been made to overcome technological constraints limiting the efficiency and coherence times of current systems, an important complementary approach is to employ parallel and multiplexed architectures. Here we follow this approach experimentally demonstrating the generation of bipartite polarization-entangled photonic states across more than 500 modes, with a programmable delay for the second photon enabled by qubit storage in a wavevector multiplexed cold-atomic quantum memory. We demonstrate Clauser, Horne, Shimony, Holt inequality violation by over 3 standard deviations, lasting for at least 45 {\mu}s storage time for half of the modes. The ability to shape hybrid entanglement between the polarization and wavevector degrees of freedom provides not only multiplexing capabilities but also brings prospects for novel protocols.

E fficient and robust entanglement generation between distant parties remains one of the most fundamental steps towards practical implementations of quantum enhanced protocols [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] , enabling amongst others quantum secure communication 1,10,[16][17][18][19][20] or distributed quantum computing 21 . The optical domain remains a platform of choice for current state-of-the-art experimental demonstrations of quantum communication protocols; however, a fundamental exponential transmission loss limits the attainable distance between the parties to several tens of km 22,23 . While direct amplification of the quantum states is forbidden by the no-cloning theorem, noise-tolerant quantum repeaters have been proposed to distribute entanglement over extended distances by generating entanglement between intermediate stations (repeater nodes) and performing the entanglement swapping protocol 24,25 . Quantum memories play a significant role in such a task facilitating adaptive strategies, and enabling hierarchical architecture of the network with neighbor pairs of nodes waiting for each other, to accomplish entanglement the generation step before trying the entanglement swap.
While there are proposals for memory-less error-correctionbased schemes 26,27 , a significant technological advancement is still required, making the quantum-memory-based repeaters the most promising solution for near-term quantum networks. The performance of such repeaters heavily relies on the efficiency and lifetime of quantum memories, and a significant effort has been devoted to their improvement 28,29 . State-of-the-art quantum memories with either high retrieval efficiency [30][31][32][33][34] or long storage times [35][36][37][38] have been demonstrated. Even though, practical quantum repeater networks will most likely require a combination of those results with each other in a parallelized or multiplexed architecture 3,7,[39][40][41][42] , in particular facilitated by multimode quantum memories [42][43][44] and optical switching networks. Importantly, there has been a tremendous progress in developing either spatially 44,45 , temporally 40,[46][47][48] , spectrally 49 or hybrid 42 multimode quantum memories. Recently demonstrated wavevector multiplexing constitutes a viable alternative offering both an unprecedented number of modes and versatile in-memory processing capabilities 43,50 . In our group, we have recently proposed and analyzed theoretically a scheme to use the wavevector multiplexing paradigm with a coldatomic quantum memory as a feasible platform for a quantum repeater 51 . Roughly speaking, during entanglement generation the large number of modes allows for successful completion of this first step in almost all cases.
Here, we experimentally demonstrate a high-rate bipartite polarization entanglement generation employing ca. M ≈ 550 pairs of photonic modes of the wavevector multiplexed quantum memory (WV-MUX-QM), with the possibility of delaying a second photon in pair for up to tens of μs. The entanglement is verified by a mode-resolved Bell state measurement certifying non-locality of generated states. The wavevector multiplexing approach is compatible with modern two-photon quantum repeater protocols robust to global phase fluctuations, and particularly suitable with promising multi-mode quantum optical channels employing multicore fibers 52,53 or free-space transmission 54 . Our scheme constitutes an important step towards quasi-deterministic entanglement generation 51 , which for the achieved parameters would improve the quantum communication rate roughly M-fold.
We note that the Bell state generation has been demonstrated before in few-mode quantum memories 55,56 , and recently with a similar mode conversion setup and a long lifetime of 1 ms 57 , albeit only in three modes.

Results
Multimode generation of Bell states. Robust quantum repeater protocols based on two-photon interference require high-rate generation of photon-atom entanglement [58][59][60] , where a y H (a y V ) corresponds to a photonic creation operator for horizontal (vertical) polarization and S y H (S y V ) is the creation operator for an atomic excitation (spin-wave) associated with the emission of an H (a V) photon. In a basic scenario, two atomic ensembles labeled H and V generate horizontally and vertically polarized photons, respectively, to further have them combined on a polarizing beamsplitter; however, such a scheme requires four atomic ensembles per repeater node (a two-ensemble interface per each neighbor) and complicates any multiplexing or multimode platform. Alternatively, with a lightmatter interface supporting 2M single-polarization modes, M photonic modes can have H polarization (V polarization) and be superimposed onto the other orthogonally polarized M modes, provided the modes are coherent with each other. While we introduce this idea in the context of a wavevector multiplexed quantum memory (WV-MUX-QM), such an approach could be suited to other multimode systems.
The WV-MUX-QM is a high-optical-density cold-atomicensemble-based quantum memory with several hundreds of wavevector or equivalently photonic emission angle modes. The modes are interfaced via spontaneous off-resonant Raman scattering to probabilistically generate a two-mode squeezed state of atomic excitations (spin-waves), and a write-out (Stokes) photons in each mode, as detailed in the ref. 43 . The probability of generating a single pair of a spin-wave and a write-out photon per mode shall be denoted by χ ≪ 1. Importantly, the intermode coherence of WV-MUX-QM memory has been previously verified 61 .
Bell state generation across many modes. Experimentally, we employ the wavevectors (or emission angles) of photons as the degree of freedom for the M modes. The memory is interfaced with spatially large write/read beams with a well defined and the fixed angle corresponding to wavevectors k W and k R jk R j ¼ À k W jk W j for write and read beam, respectively. In the memory writing process, a spin-wave-photon pairs across many angular modes are generated. The quantum state of a single pair of this kind can be described as a coherent superposition of plane-wave contributions: where k W − k w = K represents a unique spin-wave wavevector for each (further called write-out-w) generated photon's wavevector k w and A represents the rectangular field of view in the wavevector space. Additionally, the write-out photons coming from the memory have circular polarization that is further transformed to vertical one, which we explicitly denote by V subscript.
In a single write laser shot many pairs can be generated, and the full quantum state of the memory and write-out field is: where χ M is a total (across M modes) pair generation probability, which for small χ is just χ M ≈ χM. The state Ψ j i describes many spin-wave-photon pairs distributed across all possible modes including multiphoton contributions to a single mode. However, as the probability of such an event is χ 2 (for a chosen mode), for small χ they are not likely to happen, even with relatively high average of total excitations per shot. Therefore, we can focus on a single photon-spin-wave state ψ j i from which we generate the WV-MUX Bell state.
To generate the polarization degree of freedom (DoF) Bell state, we split the field of view containing 2M modes in half and superimpose the two resulting regions. Prior to the superimposition, one half is sent through half-waveplate to change polarization of the write-out photons form vertical (V) to horizontal (H). After the superimposition we end up with a WV-MUX polarization DoF Bell-like state: where the optical phase between generated write-out H and V photons with k w wavevector is φ w (k w ). The spin waves (atomic excitations) can be read-out after a programmable delay: where X ∈ {H, V} and b y X ðk r Þ denotes a creation operator for the read-out photon with polarization X and wavevector k r , and ϕ(k r , k w ) characterizes the correlation between write-out and read-out photons. The operation brings our state to: dk w dk r ϕðk r ; k w ÞB y ðk r ; k w Þ vac j i; ð5Þ where we introduced a WV-MUX Bell state creation operator B y ðk r ; k w Þ, with φ(k r , k w ) = φ w (k w ) + φ r (k r ) and φ r (k r ) being an additional phase difference for orthogonally polarized read-out photons.
We note that with anti-collinear write and read beam configuration ( k R jk R j ¼ À k W jk W j ) the write-out-read-out wavevectors on average satisfy k w jk w j ¼ À k r jk r j . Since the read-out of a spin-wave is facilitated by light-atom interaction in a finite ensemble, the angular spread of read-out photons around specific k r corresponding to a registered k w is inversely proportional to the transverse size of the atomic cloud. Assuming Gaussian profile for the transverse atomic density, we get a Gaussian mode function (conditional on the choice of the write-out wavevector k w ) for the read-out photon with k r wavevector, described by the correlation function: with the covariance matrix C, that determines the correlation strength between write-out and read-out photons in the k-space.
The correlation size gives the number of memory modes in terms of Schmidt decomposition: where α is a geometric factor, which for a rectangular area with transversally-Gaussian correlation is α = 0.565 43 , and A = L x × L y denotes the observed rectangular area in the wavevector space (note that in total we observe 2A which is divided into H and V part) and λ 1 , λ 2 correspond to eigenvalues of the C matrix (e.g., λ = σ 2 for one dimensional Gaussian parameter σ). Finally, we can rewrite ψ B in terms of the two-photon (writeout-read-out pair) subspace and separate the wavevector and polarization photonic DoFs, obtaining: with a polarization DoF state: We note that the state given by Eq. (8) has a non-trivial interdependence between the wavevector and polarization DoF via the wavevector-dependent phase φ w (k w ) + φ r (k r ) (see "Methods" section), which is experimentally feasible to be arbitrarily shaped e.g., by placing a spatial light modulator (SLM) in the farfield of the atomic cloud.
Throughout this work we choose the coordinates so that the photons are correlated as follows: where we denote k i = (x i , y i ); i ∈ {r, w}.
Experimental Bell-state preparation. To experimentally generate a state given by Eq. (8) across M ≈ 550 modes, we employ a WV-MUX-QM with a specifically designed Mach-Zehnder interferometers (MZI) placed in the far-field of the atomic cloud (see "Methods" section). A simplified experimental setup is depicted in Fig. 1. Each MZI allows dividing the emission cone (wavevector range) of write-out (read-out) photons to H and V modes and to superimpose the two parts. Furthermore, the wavevectordependent phase φ w (k w ), φ r (k r ) can be shaped by tilting one of the MZI mirrors. The MZI consists of three mirrors, a half-wave plate (HWP) and a polarizing beam-splitter (PBS). All modes have initially H polarization which is not reflected at PBS. The first mirror reflects half of the emission cone (half of the wavevector modes), which is then routed by the second mirror into a port of PBS and transmitted through. The other half of the emission cone passes through HWP, which rotates the polarization to V, and is reflected by the third mirror into the second port of the PBS. V-polarized half is reflected at the PBS, hence the two parts are superimposed in wavevector (spatial) DoF. Importantly, both parts undergo the same number of reflections. Effectively, for each pair of modes MZI connects two wavevector modes into a product of a single wavevector mode and two polarization modes.
Each MZI works as a mode converter between the wavevector and polarization DoF. The wavevector DoF enables parallelization or multiplexing for efficient quantum repeater protocols 51 , while the polarization DoF facilitates Bell state generation and enables robust two-photon protocols for entanglement creation and connection between repeater nodes [58][59][60] .
We note that the conversion between wavevector (or spatial) and polarization DoFs has been demonstrated before, also in the context of Bell state generation 41,48,57,62 , albeit with a few-mode system.
Bell-state measurement. To quantify the entanglement of experimentally generated states we perform a wavevectorresolved Bell-state measurement (BSM) with the write-out and read-out arm corresponding to the two parties-Alice and Bobusually considered in the Bell setting. With a linear MZI phase, we select measurement bases for Alice and Bob lying on the equator of the Bloch sphere. Such a choice ensures the BSM visibility remains constant regardless of the wavevectordependent phase. The local Alice and Bob generalized measurement operators {Θ(k)} k (ξ) are given by with Π ξ ¼σ x cos ξ þσ y sin ξ, whereσ x ,σ y are the Pauli operators and ξ parametrizes the measurement (e.g., ξ = 0, ξ = π/2 correspond to a measurement in diagonal and circular polarization bases, respectively). Write-out and read-out paths undergo a projective polarization measurement on a beam-displacer (BD) with ξ adjusted by halfwave plates (HWP). The two output ports ± of the BD are observed with a single-photon sensitive I-sCMOS camera 63 , which resolve individual wavevector modes with 1 px corresponding to 2.38 rad/mm (see "Methods" section). Importantly, some of the events registered with the camera correspond to dark counts, photons from the multiexcitation component OðχÞ or misalignment of the experimental setup. To model those imperfections we assume a depolarizing channel 51 over the polarization DoF, which transforms: where the visibility V in general depends on the wavevectors of write-out and read-out photons.
Expected value and Bell parameter. In the experimental demonstration, we select a linear MZI phase which for maximally correlated write-out and read-out wavevectors k w = k r = k allows us to observe a family of Bell-like states Let us calculate the average outcome of local Alice's and Bob's measurements for a single pair of write-out and read-out wavevectors Θ(k w , ξ w ) ⊗ Θ(k r , ξ r ). Post-selecting outcomes with a registered pair of photons at k r and k w , we get the expected value to be where the cosine is given by the trace over polarization DoF:

Selecting two bases
A ¼ fξ ð1Þ w ; ξ ð2Þ w g, B ¼ fξ ð1Þ r ; ξ ð2Þ r g for Alice and Bob each, respectively, we can formulate the Bell parameter in the wavevector space: For an optimal selection of bases the Bell parameter depends only on the visibility To violate Clauser, Horne, Shimony, Holt (CHSH) inequality ? and certify non-classical correlations one needs Vðk r ; k w Þ>1= ffiffi ffi 2 p .
Visibility. BSM visibility indirectly quantifies the quality of the generated entangled state for further entanglement distillation protocols. Entanglement of distillation, i.e., the number of maximally entangled states that can be distilled per a copy of the generated state, would provide a resource-oriented characterization; however, it is difficult to calculate in a generic case. Optimistically, one may use another entanglement monotone such as entanglement of formation, concurrence or negativity to upper bound the entanglement of distillation. For a Werner state, given by RHS of Eq. (12), all those monotones can be calculated analytically 64 and depend only on the visibility. Furthermore, concurrence and negativity are in this case linear functions of V, hence we will further focus on the visibility as a figure of merit quantifying the entanglement of generated states. In addition to experimental setup imperfections, the influence of multiphoton excitations and noise either from dark counts or stray photons constitutes a fundamental factor limiting the visibility. Importantly, these factors also affect the Glauber second-order intensity cross-correlation between write-out and read-out photons g (2) (k r , k w ). As shown in "Methods" section, the visibility can be written as: Interestingly, this result facilitates estimating BSM visibility via g (2) function measurements which do not require the BSM setup. Importantly, with the BSM setup in-place, the g (2) function lets us observe wavevector-resolved interference in the correlation patterns, facilitating observation of the phase profile at MZIs (see the Supplementary Material Sec. S1 and S2 for more details about the wavevector space correlations).
Experimental BSM. To observe the mode-resolved Bell parameter SðA; B; k r ; k w Þ we shall look at maximally correlated pairs of write-out and read-out wavevectors arg max jϕðk w ; k r Þj 2 . As detailed in the Supplementary Material Sec. S3, for each point in sum coordinates [y s ≡ (y w + y r )/2, x s ≡ (x w + x r )/2], we bin the coincidences in a rectangular region 2nσ x × 2nσ y around (0, 0) point in difference (y w − y r , x w − x r ) coordinates. Importantly, this affects the second-order correlation function and thus also the BSM visibility which reads: where the binning factor is and whereB w ¼ B w =M (B r ¼ B r =M) denotes the noise probability per mode in the write-out (read-out) arm, before detection. The visibility is expected to be uniform in sum coordinates (x s , y s ). For further analysis we choose n = 1 (F(n) ≈ 0.825) as a trade-off between decreasing the visibility (as n grows) and gathering a larger statistics for each point. The choice of Alice and Bob bases A; B is optimal for a subset of wavevector modes. While, we note that a constant MZI phase would enable a choice of bases simultaneously optimal for all modes, linear phase provides experimentally feasible characterization in terms of the BSM visibility. Figure 2 depicts subsequent analysis stages leading to the Bell parameters in sum coordinates Sðx s ; y s Þ. For Fig. 2a-d we choose a single pair of bases ξ ð2Þ w ; ξ ð2Þ r and illustrate coincidence maps with different combinations of polarization measurement outcomes in write-out/read-out arms. Figure 2e-h depict expected values for the BSM measurement with each combination of Alice and Bob bases e.g., Fig. 2e corresponds to the following operation on maps depicted in Fig. 2a- Wavevector-resolved Bell parameter, depicted in Fig. 2i, is obtained according to Eq. (16) by taking a combination of BSM expected values. The linear MZI phase is clearly visible in all maps, modulating the results along y s . Hence, we take an average of the Bell parameter along x s . As depicted in Fig. 2j, the averaged Bell parameter hSðx s ; y s Þi x s is sinusoidal with a fitted amplitude of 2.60 ± 0.19 yielding CHSH violation by over three standard deviations.
BSM with memory BSM visibility. The Bell parameter under a choice of optimal bases is directly proportional to the BSM visibility, making the visibility a viable figure of merit for the generated write-out-read-out bipartite states. Furthermore, from the perspective of entanglement distillation protocols, the BSM visibility determines the entanglement monotones measuring the ebit content of a single generated state. Importantly, entanglement distribution protocols such as quantum repeater nets can greatly benefit from the delayed release of the read-out photon, which allows improving the protocol success rate via, among others, hierarchical architectures 25,59 or multiplexing of the read-out photonic mode 51 . Hence, we perform a series of measurements with increasing memory time t ∈ [0.3, 60.3] μs and with a very quickly changing linear MZI phase, which allows us to retrieve wavevector-resolved BSM visibility from the sinusoidal patterns in coincidence maps.
Spin-wave decoherence. Fundamentally, the BSM visibility at larger memory times is limited by the spin-wave decoherence. For WV-MUX-QM's light-matter interface, we select the so-called clock transitions insensitive to the first order to stray magnetic fields and fix the quantization axis by introducing external constant magnetic field along the cloud (z-axis); hence, the decoherence mechanism in our case mainly due to thethermal atomic motion. Inevitably, random displacement of atoms distorts the spatial structure of a spin-wave 43 . The decoherence quantified as an average overlap with the initial spin-wave state is Gaussian with the characteristic time τ depending on the spinwave wavevector and given by: w ; ξ ð2Þ w g and fξ ð1Þ r ; ξ ð2Þ r g bases. Each expected value is a linear combination of coincidences between ± ports. i Bell parameter Sðx s ; y s Þ retrieved by fitting a sine function (with x s and y s spatial frequencies) to the expected value for each combination of bases and taking the Bell linear combination of the results. j Bell parameter averaged over x s (red curve) along fitted sine (blue dashed curve) with amplitude 2.60 ± 0.19 indicating CHSH violation by more than three standard deviations σ. Subsequent shadings correspond to σ and 2σ. Solid lines indicate CHSH inequality violation threshold.
with γ depending on the atomic mass of 87 Rbm, cloud temperature T and Boltzmann constant k B . Spin-wave decoherence can be accounted for by plugging into Eq. (21) a time-dependent retrieval efficiency η r ðt; kÞ ¼ η r ð0Þ expðÀt 2 =τðkÞ 2 Þ: ð26Þ Visibility model. In our experimental setup the MZI used to divide the atomic emission cone into H and V polarized parts superimposes a mode with the lowest modulus wavevector min jk H j from the H part onto the highest modulus wavevector max jk V j from the V part and vice versa. As a consequence, spin-waves corresponding to different polarization parts decohere with a different rate τ(|k H |) ≠ τ(|k V |) effectively further deteriorating BSM visibility. Importantly, it could be amended by an improved MZI setup. Let us denote k min ¼ minðjk H j; jk V jÞ, k max ¼ maxðjk H j; jk V jÞ. Since longer wavevectors are associated with faster decoherence, let as also denote τ min ¼ τðk max Þ, τ max ¼ τðk min Þ and Δ ¼ τ À1 min þ τ À1 max . As demonstrated in "Methods" section, the additional visibility reduction is: With the increasing storage time the decoherence deteriorates the visibility. As detailed in Methods, the storage-time-dependent version of Eq. (21) can be written as: The total visibility is given by the product of Eq. (28) and Eq. (27).
BSM visibility maps. Experimentally, the BSM visibility can be retrieved from coincidence measurements with a quickly changing linear Mach-Zehnder inteferometer phase. The coincidence map for each combinations of measurement ports (s r , s w ) = {(+, +), (+, −), (−, +), (−, −)} is smoothed with one-dimensional Gaussian filters (σ = 1 px for y i direction and σ = 10 px for x s direction) and for each x s a row of data along y s is selected. Each row is divided into overlapping 50 px segments with an 1 px step and to each segment we fit a model given by a cosð2πf y y s Þ þ b and obtain visibility as a/b, which we assign to the y s coordinate of the visibility map given by the center position of the 50 px segment. Figure 3 depicts experimental data along a fitted visibility model given by the product of Eq. (28) and Eq. (27),   4 Persistence of high mean visibility. Visibility has been averaged over wavevector modes hVðx s ; y s ; tÞi x s ;y s for increasing memory times t. a Blue curve represents the fraction of modes violating CHSH inequality, while the red curve corresponds to the visibility averaged over this subset of modes. b Visibility averaged over all modes (purple curve) yields good agreement with the model based on Glauber cross-correlation g (2) (red dashed curve, cf. Eq. (20)). Lines between points serve as guide to the eye. Errorbars represent extreme cases in the ensemble of wavevector-modes, with top/bottom of errobar corresponding to the best/worst visibility ± 1 s. d. obtained from the fit. assuming τ(k) = γ/|k|. Hatched regions correspond to visibility above 1= ffiffi ffi 2 p , which predicts CHSH violation while additional isolines are drawn at levels from 0.6 downward with a step of 0.1. The model fit yielded with γ corresponding to a temperature of T ¼ 47 þ5 À4 μ K . For completeness, we note that an independent measurement, described in Supplementary Material Sec. S4, yielded a consistent temperature of T ¼ 48 þ6 À5 μ K , the χ-dependent noiseB r ðχÞ ¼ 0:131 ± 0:015 and a read-out efficiency of η r (0) = 0.405 ± 0.015.
Mode-averaged visibility. Obtained visibility maps let us study the properties of an average memory mode. As depicted in Fig. 4a, up to ca. 30 μs nearly all modes would violate CHSH, while for 45 μs it remains true for half of the modes. The average visibility amongst the CHSH violating modes remains fairly constant with increasing memory time. Figure 4b depicts very good agreement of the visibility averaged over all modes V h i with the modeaveraged prediction from the measured Glauber cross-correlation ðg ð2Þ ðtÞ À 1Þ=ðg ð2Þ ðtÞ þ 1Þ . Visibly, above 45 μs the performance of a significant number of modes has severely deteriorated. Nevertheless, modes with a lower modulus wavevector can have significantly longer lifetimes. Unfortunately, the number of modes with a given modulus wavevector |k| is roughly proportional to |k|, hence there are few modes with a long lifetime.
Wavevector-dependent decoherence rate. For quantum memories utilizing a dense cold ensemble of atoms, it has been demonstrated that thermal motion constitutes the main source of spinwave decoherence 38 , provided the so-called clock transitions-in the first order robust to external magnetic field fluctuations-are utilized for light-matter interface. Characteristically, decoherence due to thermal motion leads to Gaussian temporal profile of the retrieval efficiency, as given by Eq. (26), with a characteristic time of τ(|k|) = γ/|k|. A vast range of wavevector modes simultaneously utilized in the WV-MUX-QM provides a unique opportunity to precisely observe experimentally the decoherence time τ(|k|) dependence on the wavevector length |k|. As depicted in Fig. 5, experimental results are in a good agreement with the predictions of thermal-motion-induced decoherence.

Discussion
We experimentally demonstrated the generation of polarizationentangled bipartite Bell states in ca. 550 photonic modes, with an inherent programmable delay for the second photon in a pair. Our approach harnesses hybrid atom-photon entanglement generation in a single cold high-density ensemble of Rb-87 atoms, which is the central part of the wavevector-multiplexed quantum memory (WV-MUX-QM). Ca. 1100 wavevector (angular emission) modes of the WV-MUX-QM are divided into H and V polarized photonic modes, which are superimposed pairwise forming 550 modes in a polarization superposition. Importantly, a write-out photon in one of those combined modes may have been created in an H or a V part, in each case having a different initial wavevector and thus being entangled with a collective atomic excitation-a spin-wave-in a different memory mode. The wavevectors of write-out and corresponding read-out photons are correlated with each other; hence, a write-out photon created in the H(V) part of write-out emission cone is accompanied by a read-out photon also in the H(V) part, albeit of the read-out emission cone. Hence, the read-out modes need to be super-imposed in the same way as write-out modes. This way, a write-out photon and a correlated read-out photon may be either both H polarized or both V polarized, constituting a polarization entangled state. Fundamentally, this concept realizes a conversion between wavevector and polarization degrees of freedom.
We have demonstrated CHSH inequality violation by more than 3 standard deviations with the Bell parameter reaching S ¼ 2:60 ± 0:19. Via wavevector-resolved Bell state measurement (BSM) with varying storage times, we established experimentally and with a theoretical model the BSM visibility-quantity proportional to the Bell parameter under the optimal choice of bases -behavior with wavevector and temporal resolution. After 45 μs storage time, ca. 50% of modes still indicate CHSH violation, while after 60 μs it amounts to around 10%. Importantly, with single-copy distillation protocols 65 , even slightly entangled states can be probabilistically distilled into Bell states, hence the BSM visibility averaged over the modes-which reflects the average ebit content of a generated state-is a significant figure of merit. In our experiment, it has fallen merely to around 50% at 60 μs storage time, as compared to the initial ca. 80% for immediate readout. With the current experimental parameters, we estimate the probability to generate and detect at least one atom-photon Bell state across all pairs of modes to be 1 − (1 − ηχ) M ≈ 35% (89% when only the detection and not the filtering system efficiency is contained in η).
Quantum repeater networks constitute the most versatile application of entanglement amending the fundamental exponential loss of photonic channels. While the performance and feasibility of WV-MUX-QM platform for quantum repeater networks has been discussed elsewhere 51 , we note here that obtaining hybrid entanglement is essential for two-photon protocols 58-60 which solve the phase stability issues inherent to the DLCZ protocol 25 . From the perspective of feasible quantum repeaters, the 800 nm photons are suboptimal having an order of magnitude higher fibertransmission losses than telecom; however, conversion to telecom photons as well as light-atom interfaces at telecom wavelengths have been demonstrated 2,66-68 . Furthermore, quantum repeater protocols generally do not require phase stability between different modes, rendering free-space transmission a viable alternative. Importantly, WV-MUX-QM platform offers versatility beyond entanglement generation, such as intramemory processing of stored spin-wave states 69 -including operations on non-classical single-excitation states 50 , continuous variable and temporal processing 70,71 or ultra-narrow band temporal imaging 72 . Furthermore, the rich atomic structure of alkali atoms may be employed to design various light-matter interfaces suited for a particular application 73 .
Finally, we note that our experimental demonstration is far from the ultimate capabilities of wavevector-multiplexed memories. Cold atomic memories with higher coherence times 37 and read-out efficiencies 31 has been demonstrated. The wavevector multiplexing principle could be applied with atomic clouds captured in a dipole trap or an optical lattice, which can prolong the memory lifetime to the sub-second regime 36,62 . However, the simplest improvement would be to observe a larger angular range of photonic emission from the memory. As estimated in our theoretical work 51 , the number of modes attainable with standard optical components and commercial CMOS image sensors would reach over 5000 enabling feasible high-rate entanglement generation.

Methods
Wavevector-dependent phase. The ability to shape the wave-vector dependent phase φ w (k w ), φ r (k r ) opens vast possibilities. For instance, consider the moststrongly correlated write-out and read-out wavevectors k w = k r = k. By selecting φ (k) ≡ φ w (k) + φ r (k) periodically equal to either π or 0 on a rectangular grid, we generate either a Φ − or Φ + Bell-state, respectively, depending on the wavevector k. Interestingly, the wavevector of the scattered write-out photon k is inherently random and selected by the process of spontaneous Raman scattering. Such a quantum-random sampling of generated states is equivalent to a quantum-random choice of measurement bases for Alice and Bob in the Bell setting (since either the states or the bases may be equivalently altered). After the measurement, the central WV-MUX-QM can reveal the exact phase profile φ(k) that was employed, allowing Alice and Bob to interpret their results.
The WV-MUX-QM platform. The WV-MUX-QM is based on an elongated (0.5 mm × 0.5 mm × 7 mm), high-optical-depth ensemble of rubidium-87 atoms prepared in magnetooptical trap (MOT). The optical depth (OD) amounts to 150 at the F = 1, m F = 1 → F = 2, m F = 2 transition of D2 line, which corresponds to OD = 25 at the Write laser transition (in resonance). To define the quantization axis, the cloud is kept in a constant (1 Gauss) magnetic field oriented along the propagation axis (z-axis). The MOT loading time is set to 2 ms which includes compression stage (700 μs) and polarization gradient cooling, with magnetic gradient switched off (300 μs). After MOT loading stage the atoms are optically pumped to the F = 1, m F = −1 state using 70 μs-long hyperfine pump pulse (resonant with the F = 2 → F = 2 transition of D1 line), that illuminates the cloud from four sides (along cooling beams) and has 15 mW power in total. The Zeeman sublevel pumping is achieved by illuminating the atoms along the z-axis by a circularly polarized laser beam (resonant with F = 1 → F = 1 transition of D2 line) for 55 μs from the beginning of the hyperfine pump pulse. We estimate the efficiency of the Zeeman pumping to be about 70% 43 . Longer hyperfine pump pulse duration combined with additional (5 mW, 75 μs) "clearing" pulse of the Read laser guarantees, that the storage level (F = 2, m F = 1) is emptied before the quantum memory protocol. To generate the two-mode squeezed state of spin-wave and a write-out photon, we use Raman interaction in the Λ scheme. The write-out photon and spin-wave pairs are generated using a 30 MHz red-detuned σ +polarized Write laser pulse (300 ns) on F ¼ 1 ! F 0 ¼ 2 transition of D2 line. The Write laser power is chosen to provide desired pair generation probability χ ≈ 0.01, which in our case corresponds to about 10 μW. Since we filter the write-out photons to have orthogonal polarization to the laser photons, the spin waves of our interest are created between the F = 1, m F = − 1 and F = 2, m F = 1 states. To convert (read-out) the spin waves to read-out photons we use σ − -polarized Read laser pulse (300 ns) tuned to the F ¼ 2 ! F 0 ¼ 2 transition of the D1 line. The Read laser power is chosen to give the readout pulse duration of approx. 200 ns, which corresponds to about 100 μW. The double-Λ configuration (consisting of D1 and D2 lines of rubidium-87) allows us to efficiently filter write-out and read-out photons from back-emissions, that may lead to uncorrelated detection events and spoil the measurements. Finally, to provide the best signal-to-noise ratio, the writeout and read-out photons pass through narrowband atomic filters. The filters are glass cells containing optically pumped 87 Rb that absorb residual laser light (which is separated from photons by 6.8 GHz), while being transparent to the signal (write-out or read-out) photons (see Parniak et al. for details of the filtering system 43 ).
Single-photon sensitive I-sCMOS camera. Generation of Bell states across many modes requires efficient detection of single-photons with a high spatial resolution. High number of modes-several hundreds-render arrays of single-mode singlephoton detectors practically unfeasible; however, recent development in singlephoton cameras shows the detection technology is sufficient for near-term applications 63,74 . Among single-photon camera solution, the most commonly employed are intensified CCD (ICCD) and electron-multiplying CCD (EMCCD). Neither of those solutions achieves acquisition times on the order of a few μs, which would allow real-time response to detected photons required in many protocols such as quantum repeater operation. Additionally, EMCCD cameras have a high read-out noise resulting in high dark-count rates and consequently deteriorating the fidelity of the generated entanglement. Recently commercialized, intensified scientific CMOS (I-sCMOS) cameras solve those issues to a great extent. Additionally, so-called quanta image sensors 74 offer further improvements, especially in terms of quantum efficiency.
In our experiment, we employ a custom I-sCMOS camera which involves a 2stage image intensifier (luminous gain of 5 × 10 6 ) imaged by a relay lens (f-number of 1.1, magnification of −0.44) onto a 5.5 Mpx (mega pixel), 2560 × 2160 px CMOS sensor (effective pixel pitch ca. 15 μm accounting for relay lens magnification). Image intensifier involves a GaAs photocathode which around 800 nm offers the overall detection efficiency of ca. 20%. Active multiplexing, required in many applications 51 would inherently need real-time feedback from the camera, which has not been implemented in the current experiment; however, it would be possible with one of the modern fast CMOS sensors.
Mode detection cross-talk. With the high resolution of single-photon cameras, the probability for two or more generated write-out photons to share the same (up to the camera resolution) central wavevector is given by ≈1 − n!/(n − ⌈χM⌉)! × 1/ n ⌈χM⌉ , where n = 130 × 160 is the number of observed pixels (px). In our case (χ ≈ 0.01) the probability is 4.8 × 10 −4 and we will neglect such cases. Similarly, if we focus on a single jth mode with a detected write-out photon, the probability of registering a read-out photon coming from another mode in the 3σ radius around the most-likely position of jth read-out can be approximated as 1 À ½1 À χη r η πð3σÞ 2 =n MÀ1 , where σ is the mode size in pixels, η r ≈ 40% corresponds to the read-out efficiency and η ≈ 8% to the filtering and detection system efficiency. In our case the probability amounts to around 0.17%. Therefore, in most cases, we can consider a single isolated jth mode with the results being valid for any j. The used formulae have been derived in the Supplementary Material Sec. S6.
Visibility model. Assuming low average number of detected photons per experiment n ( 1, we can approximate the photon number cross-correlation function as: g ð2Þ ðk r ; k w Þ % p w;r =ðp w p r Þ; ð32Þ where p w,r ≡ p w,r (k r , k w ) is the single-experiment probability of observing a coincidence between a write-out and read-out photon and p w ≡ p w (k w ) (p r ≡ p r (k r )) denotes the marginal probability of observing a write-out (read-out) photon. The coincidence probability can be written as p w,r = g (2) p w p r , with g (2) ≡ g (2) (k r , k w ). For a given point (k r , k w ) in the wavevector space, we measure the BSM visibility by comparing the number of coincidences during measurement set up for maximally constructive (+) or destructive (−) interference. The numbers of coincidences serve as probability estimates and hence give the visibility as Vðk r ; k w Þ ¼ ðp ðþÞ w;r À p ðÀÞ w;r Þ=ðp ðþÞ w;r þ p ðÀÞ w;r Þ: ð33Þ In (−) settings only noise coincidences are registered i.e. p ðÀÞ w;r ¼ p w p r while p ðþÞ w;r ¼ p w;r ¼ g ð2Þ p w p r . This directly gives us Eq. (20).
Storage-time dependence. Let us first consider the visibility model as given by Eq. (21) and under experimentally verified assumptions of negligible write-out noisẽ B w ( χ : We shall further employ a model for the read-out noise: Different wavevectors for superimposed H, V modes. Here we quantify the effects of different decoherence rates on the generated states. Let us denote jk H j 2 ¼ y 2 þ x 2 ; ð39Þ and jk V j 2 ¼ ðy max À yÞ 2 þ x 2 ¼ jk H j 2 þ y 2 max À 2y max y; where y max corresponds to the maximal y wavevector component entering the MZI. We shall focus on a single mode and modify the state given by Eq. (9) so that H and V parts have different (modulus) coefficients: with c 2 P ðtÞ ¼ expðÀt 2 jk P j 2 =γ 2 Þ; P 2 fH; Vg: Assuming |k H | > |k V | the temporal evolution of the state in Eq. (41) brings it closer to H j i r H j i w . If we consider a reduced state e.g., of write-out only, the evolution moves the state from the equator of the Bloch sphere to the pole. Intuitively, as the polarization measurement projects the state on some axis in the equator plane, the visibility will be reduced. Direct calculation yields withṼ ðtÞ ¼ 1= cosh½t 2 =ð2γ 2 Þ ðy 2 max À 2y max yÞ: Importantly, no assumption on the actual form of τ(k) need to be made to give an expression forṼðtÞ, as long as Eq. (26) holds. Let us denote k min ¼ minðjk H j; jk V jÞ, k max ¼ maxðjk H j; jk V jÞ. Since longer wavevectors are associated with faster decoherence, let as also denote τ min ¼ τðk max Þ, τ max ¼ τðk min Þ and Δ ¼ τ À1 min þ τ À1 max . This way, we arrive at Eq. (27). Interestingly, MZI configured to give different c H (t) and c V (t) could be used to demonstrate a violation of the so-called tilted Bell inequalities 75 by generating a family of states with a varying degree of entanglement, across the memory modes. The only required modification in the experimental setup would be to remove quarter waveplates, depicted in Fig. 1, which would amount to performing BSM with polarization operators of the formσ z cos ξ þσ x sin ξ corresponding to the projection on the Bloch sphere's meridian.

Data availability
Data presented in Figs. 2-5 have been deposited in RepOD Repository for Open Data at https://doi.org/10.18150/RJBMOD. Any other data that support the findings of this study are available from the corresponding author upon reasonable request.