Experimental certification of millions of genuinely entangled atoms in a solid

TL;DR

This paper reports the experimental certification of 16 million genuinely entangled atoms in a solid, demonstrating the robustness and scalability of multipartite entanglement in quantum memories for future quantum networks.

Contribution

It introduces a new entanglement witness method applicable to large-scale solid-state systems, enabling certification of massive multipartite entanglement even with losses.

Findings

Certified 16 million entangled atoms in a solid-state quantum memory.

Developed a scalable entanglement witness based on collective emission.

Demonstrated the robustness of multipartite entanglement at unprecedented scales.

Abstract

Quantum theory predicts that entanglement can also persist in macroscopic physical systems, albeit difficulties to demonstrate it experimentally remain. Recently, significant progress has been achieved and genuine entanglement between up to 2900 atoms was reported. Here we demonstrate 16 million genuinely entangled atoms in a solid-state quantum memory prepared by the heralded absorption of a single photon. We develop an entanglement witness for quantifying the number of genuinely entangled particles based on the collective effect of directed emission combined with the nonclassical nature of the emitted light. The method is applicable to a wide range of physical systems and is effective even in situations with significant losses. Our results clarify the role of multipartite entanglement in ensemble-based quantum memories as a necessary prerequisite to achieve a high single-photon process fidelity crucial for future quantum networks. On a more fundamental level, our results reveal the robustness of certain classes of multipartite entangled states, contrary to, e.g., Schr\"odinger-cat states, and that the depth of entanglement can be experimentally certified at unprecedented scales.

Tables

Table 1. Table 1: Results for entanglement depth K 𝐾 K . Depending on the level of modeling the inefficiencies of the experimental setup, different values for p 1 subscript 𝑝 1 p_{1} and hence for K 𝐾 K are obtained (cf. Fig.  2 (b)). By sampling p 1 , p 2 subscript 𝑝 1 subscript 𝑝 2 p_{1},p_{2} and N 𝑁 N around the measured values within the estimated uncertainties, we calculate the expected entanglement depth K 𝐾 K . The values in the last columns are lower bounds on K 𝐾 K with confidence 3 ​ σ = 99.7 % 3 𝜎 percent 99.7 3\sigma=99.7\% .

Level of modeling	$p_1$	$K$	$K-3\sigma$
(i) raw data	0.0023(3)	$4.76\times {10}^5$	$7.54\times {10}^4$
(ii) after re-mission	0.013(2)	$1.64\times {10}^7$	$3.72\times {10}^6$
(iii) before re-mission	0.016(2)	$2.46\times {10}^7$	$5.24\times {10}^6$
(iv) after absorption	0.16(1)	$3.23\times {10}^9$	$2.09\times {10}^9$

Equations

$$H_{\mathrm{int}}=\sum _{j,\vec{k}}{e}^{-i\vec{k}\cdot {\vec{r}}_j}{a}_{\vec{k}}{\sigma }_{+}^{(j)}+e^{i\vec{k}\cdot {\vec{r}}_j}{a}_{\vec{k}}^{†}{\sigma }_{-}^{(j)},$$

$$|1⟩\to |D_1⟩\propto \sum _j{e}^{-i\vec{k}\cdot {\vec{r}}_j}|g\dots g{e}_{j}g\dots g⟩.$$

$$|\psi ⟩=|{\varphi }_1⟩\otimes \cdots \otimes |{\varphi }_M⟩,$$

$$|D_2⟩\propto \sum _{j<l}{e}^{-i\vec{k}\cdot ({\vec{r}}_j+{\vec{r}}_l)}|g\dots g{e}_{j}g\dots g{e}_{l}g\dots g⟩,$$

$$U^{†}{a}_{\vec{k}}U=\sqrt{\eta }{N}^{-1/2}{S}_{-}^{\vec{k}}+\sqrt{1-\eta }{a}_{\vec{k}},$$

$$|\psi ⟩=\underset{i=1}{\overset{M}{⨂}}\left(a_i|d_0⟩+b_i|d_1⟩+c_i|d_2⟩\right),$$

$$p_1=\frac{|A{|}^2}{M}{\left|\sum _i\frac{b_i}{a_i}\right|}^2$$

$$p_2=\frac{|A{|}^2}{M^2(1-\frac{1}{N})}{\left|\sqrt{2}\sum _{i<j}\frac{b_i{b}_j}{a_i{a}_j}+\sqrt{1-\frac{1}{K}}\sum _i\frac{c_i}{a_i}\right|}^2$$

$$p_2^{\min }(p_1,M)=\underset{\psi :p_1=\text{const}}{\min }{p}_2.$$

$$f(\{a_i,b_i{\}}_i)=f_2(\{a_i,b_i{\}}_i)+\lambda [f_1(\{a_i,b_i{\}}_i)-C]$$

$$p_1^{\mathrm{sym}}=M{a}^{2M-2}{b}^2$$

$$p_2^{\mathrm{sym}}=a^{2M}{\left(\frac{1}{\sqrt{2}}(M-1)\frac{b^2}{a^2}+\frac{c}{a}\right)}^2,$$

$$J=\frac{2d{e}^{-d}}{1-e^{-2d}}.$$

$$\text{SNR}=\frac{\eta {\left|\alpha \right|}^2}{{\left|\alpha \right|}^2/N+\delta },$$

$$\text{SNR}=\frac{\eta {\left|\alpha \right|}^2}{{\left|\alpha \right|}^2/(N(1-e^{-1/R{T}_1}))+\delta }.$$

$$\sum _{i<j}\frac{b_i{b}_j}{a_i{a}_j}=\frac{1}{2}{\left(\sum _i\frac{b_i}{a_i}\right)}^2-\frac{1}{2}\sum _i{\left(\frac{b_i}{a_i}\right)}^2.$$

$$\begin{array}{rl}{p}_2& =\frac{|A{|}^2}{M^2}|\frac{1}{\sqrt{2}}{\left|\sum _i\frac{b_i}{a_i}\right|}^2\\ & {-\frac{1}{\sqrt{2}}\sum _i{e}^{2i({\phi }_i-\bar{\phi })}{\left|\frac{b_i}{a_i}\right|}^2+\sum _i{e}^{i({\vartheta }_i-2\bar{\phi })}\left|\frac{c_i}{a_i}\right||}^2.\end{array}$$

$$p_1=\frac{A^2}{M}{\left(\sum _i\frac{b_i}{a_i}\right)}^2$$

$$p_2=\frac{A^2}{M^2}{\left(\sqrt{2}\sum _{i<j}\frac{b_i{b}_j}{a_i{a}_j}+\sum _i\frac{c_i}{a_i}\right)}^2.$$

$$\begin{array}{rl}\frac{\partial f}{\partial a_i}& =\frac{f+\lambda C}{a_i}\\ & -\frac{A}{a_i}\left(\frac{\sqrt{2}{b}_i}{a_i}\sum _{j\ne i}\frac{b_j}{a_j}+\frac{a_i}{c_i}+\frac{c_i}{a_i}+\lambda \frac{b_i}{a_i}\right),\end{array}$$

$$\frac{\partial f}{\partial b_i}=\frac{A}{a_i}\left(\sqrt{2}\sum _{j\ne i}\frac{b_j}{a_j}-\frac{b_i}{c_i}+\lambda \right)$$

$$\frac{\partial f}{\partial \lambda }=f_1-C.$$

$$\frac{b_i}{c_i}=\sqrt{2}\sum _{j\ne i}\frac{b_j}{a_j}+\lambda ,$$

$$\frac{1}{a_i{c}_i}=\frac{f+\lambda C}{A}.$$

$$b_i\left(\frac{1}{c_i}+\sqrt{2}\frac{1}{a_i}\right)=\sqrt{2}\sum _j\frac{b_j}{a_j}+\lambda .$$

$$a_i{c}_i=a_j{c}_j$$

$$b_i\left(\frac{1}{c_i}+\frac{\sqrt{2}}{a_i}\right)=b_j\left(\frac{1}{c_j}+\frac{\sqrt{2}}{a_j}\right)$$

$$a_i=\{\begin{array}{rl}& a\\ & \sqrt{\frac{x_0}{3}\left(1+(-1{)}^{2/3}\frac{x_1}{x_2}+(-1{)}^{-2/3}{x}_2\right)}\\ & \sqrt{\frac{x_0}{3}\left(1+\frac{x_1}{x_2}+x_2\right)}\\ & \sqrt{\frac{x_0}{3}\left(1+(-1{)}^{-2/3}\frac{x_1}{x_2}+(-1{)}^{2/3}{x}_2\right)},\end{array}$$

$$x_0=1-a^2-2\sqrt{2}ac,$$

$$x_1=1-6c^2\frac{1-c^2-\sqrt{2}ac}{x_0^2}$$

Full text

A clear picture of large-scale entanglement with its complex structure is so far not developed. It is however important to understand the role of different facets of multipartite entanglement in nature and in technical applications Ladd et al. (2010); Pezzè et al. (2016). For example, the so-called Schrödinger cat states Schrödinger (1935) are fundamentally different from a single photon coherently absorbed by a large atomic ensemble; even though both are instances of multipartite entanglement (Mandel and Wolf, 1995, chapter 16.5). The theoretical study of large-scale entanglement has to be followed by an experimental demonstration, which consists of two basic steps: the preparation of an entangled system and a subsequent appropriate measurement verifying the presence of entanglement. In the context of entanglement in large systems, the preparation of entanglement is generally much simpler than its verification. For example, single-particle measurements are often not possible and collective measurements are typically restricted to certain types and are of finite resolution. These limitations call for new witnesses that allow one to certify entanglement based on accessible measurement data.

The concept of entanglement depth Sørensen and Mølmer (2001) was shown to be meaningful for and applicable to large quantum systems. It is defined as the smallest number of genuinely entangled particles that is compatible with the measured data. This allows one to witness at least one subgroup of genuinely entangled particles in a state-independent and scalable way. Large entanglement depth was successfully demonstrated with so-called spin-squeezed and oversqueezed states by measuring first and second moments of collective spin operators Riedel et al. (2010); Gross et al. (2010); Vitagliano et al. (2014); Lücke et al. (2014); lately up of 680 atoms Hosten et al. (2016). Similar ideas were realized for photonic systems Beduini et al. (2015); Mitchell and Beduini (2014). Recently, a witness was proposed that is designed for the W state, which is a coherent superposition of a single excitation shared by many atoms Haas et al. (2014). Based on this witness, an entanglement depth of around 2900 was measured McConnell et al. (2015). However, these witnesses do not detect entanglement when the vacuum component of the state is dominant Haas et al. (2014), even though the W state is known to be quite robust against various sources of noise, in particular, against loss of particles and excitation Stockton et al. (2003). Hence, much larger values for the entanglement depth could be expected.

In this paper, we present theoretical methods and experimental data that verify a large entanglement depth in a solid-state quantum memory. A rare-earth-ion-doped crystal spectrally shaped to an atomic frequency comb (AFC) is used to absorb and re-emit light at the single-photon level Afzelius et al. (2009); Clausen et al. (2011); Saglamyurek et al. (2011); Clausen et al. (2012), where at least 40 billion atoms collectively interact with the optical field. Using the measured photon number statistics of the re-emitted light we collect partial information about the quantum state of the atomic ensemble before emission. Then, we show that certain combinations of re-emission probabilities for one and two photons imply entanglement between a large number of atoms. With the measured data from our solid-state quantum memory we demonstrate inseparable groups of entangled particles containing at least 16 million atoms.

I Results

Before discussing the experiment, we give an intuitive explanation for the appearance of large entanglement depth when a large atomic ensemble coherently interacts with a single photon (see Fig. 1(a)). Suppose that $N$ two-level atoms ($|g\rangle$ and $|e\rangle$ denote ground and excited state, respectively), couple to a light field. The quantised interaction in the dipole approximation is described by  Loudon (2000)

[TABLE]

that is, a single photon with wave vector $\vec{k}$ is annihilated by exciting atom $j$ via $\sigma_{+}\left|g\right\rangle=\left|e\right\rangle$ and vice versa. The phase is given by the scalar product between $\vec{k}$ and the position $\vec{r}_{j}$ of the atom. When an incoming light field is absorbed via interaction (1), the imprinted phase relation between the atoms serves as a memory for the direction and the energy of the absorbed photons. Without this information, a spontaneous, directed re-emission is not possible. In other words, phase coherence between the atoms is necessary in order to a have well-controlled re-emission direction Duan et al. (2001); Scully et al. (2006). Now, depending on the nature of the absorbed light, this coherence implies entanglement between the atoms or not. On the one hand, the absorption of a coherent state leads to a coherent atomic state, which is unentangled (Mandel and Wolf, 1995, chapter 16.7). On the other hand, if a single photon $|1\rangle$ is absorbed, the quantisation of the field leads to a W state (or Dicke state with a single excitation) of the atomic state (Mandel and Wolf, 1995, chapter 16.5)

[TABLE]

Then, the ensemble is genuinely multipartite entangled Stockton et al. (2003). These examples suggest a generic relation between directed emission, single-photon character of the emitted light and large entangled groups.

In our experiment, we use a neodymium-based solid-state quantum memory operating at a total read-write efficiency of 7% (see Fig. 1(b)). This memory was demonstrated to be capable of storing different types of photonic states and preserving state properties such as the single-photon character Clausen et al. (2011, 2012); Bussières et al. (2014); Tiranov et al. (2015, 2016). A heralded single photon is produced via spontaneous parametric down conversion Clausen et al. (2014) and coupled to the atomic ensemble, which was prepared in the ground state $\left|D_{0}\right\rangle=|g\rangle^{\otimes N}$. After a 50 ns delay time, the coherent excitation is spontaneously re-emitted in forward direction and detected. In practice, this optical state is not exactly a single photon. Due to losses at different levels, the state contains a large vacuum component. Also higher photon components are present. However, since directed emission and non-classical photon number statistics are largely preserved, entanglement between large groups of atoms is expected.

In order to certify this entanglement, we develop the following entanglement witness (see Methods for details). Suppose a pure state that is subdivided into a product of $M$ groups

[TABLE]

where the $|\phi_{i}\rangle$ are arbitrary. Phase coherence between the groups imply that each group has to carry some excitation. This necessarily amounts to an emission spectrum that also contains multi-photon components.

To be more specific, we consider the probabilities of the atoms emitting one and two photons, $p_{1}$ and $p_{2}$, respectively. In the low-excitation limit (see Methods), these probabilities correspond to $p_{1}=|\left\langle D_{1}\right|\left.\psi\right\rangle|^{2}$ and $p_{2}=|\left\langle D_{2}\right|\left.\psi\right\rangle|^{2}$, where

[TABLE]

that is, the phase-coherent superposition of two excitiatons. As shown in the Methods, it is possible to find the minimal $p_{2}$ for a given $p_{1}$ within the class (3) with fixed $M$. By varying $p_{1}$ and $M$ one finds a lower bound on $p_{2}$ as a function of $p_{1}$ and $M$. Given the linearity of $p_{1}$ and $p_{2}$ when mixing states like in Eq. (3) (with arbitrary grouping but lower-bounded $M$), the extension of the bound to mixed states is straightforward. Examples of such lower bounds are shown in Fig. 2. Note that the bounds are independent of $N$ if $N\gg 1$. Comparing the lower bounds with experimental data in turn gives an upper bound on $M$ and, by additionally measuring $N$, a lower bound on the entanglement depth, which simply reads $K=N/M$.

Experimentally, $p_{1}$ is obtained from the probability to measure a single re-emitted photon in the forward mode ($\vec{k}_{f}$ in Fig. 1 (b)) at a predetermined time, which we herald by the detection of the idler photon at the source. The value of $p_{2}$ corresponds to the two-photon statistics of the re-emitted light in the $\vec{k}_{f}$ mode. It is inferred from the measured autocorrelation function $g^{(2)}_{ss|i}=2p_{2}/p_{1}^{2}$ and $p_{1}$. The identification of photonic Fock states with atomic Dicke states is possible because the memory is initially prepared in the ground state and because of the photon number statistics of the source (see Methods). From the raw data, we find $p_{1}=2.3(3)\times 10^{-3}$ and $p_{2}=5(2)\times 10^{-8}$. The relatively small value of $p_{1}$ is a product of the efficiencies of the source, the memory and detectors. The partial subtraction of these losses leads to different values of the effective $p_{1}$ (see Fig. 2, table 1 and Methods): (i) raw data; (ii) subtraction of detector noise, that is, the actual value for the emitted light; (iii) subtraction of the re-emission inefficiency, that is, the single-excitation component of the atomic ensemble just before emission; (iv) subtraction of phase noise during storage, that is, overlap with the $|D_{1}\rangle$ state right after absorption. The values of $p_{2}$ directly follow using the measured autocorrelation function $g^{(2)}_{ss|i}=0.020(3)$, which is –under conservative assumptions– not affected by this modeling.

A key element in the experiment is the high-precision measurement of $N$. For this, the relation between ensemble size and directionality of the re-emitted light is exploited. The ratio of the coherent emission in the forward direction and the incoherent emission in the backward direction is a lower bound on the number of resonant atoms Scully et al. (2006). Since incoherent emission from single photons is much lower than detector dark counts, the single photon source is replaced by a bright coherent state for this measurement (see Fig. 1(b) and Methods) and we find $N\geq 4.0(1)\times 10^{10}$. The resulting $K$ depends on the level of modeling (i) to (iv) as mentioned before (see table 1). Our data analysis illustrates how decoherence and noise reduces the certifiable entanglement depth. Immediately after absorption of the single photon (iv) we have an entanglement depth of about $10^{9}$; this reduces to about $10^{7}$ just before and immediately after re-emission (iii) and (ii), respectively, while when taking into account all losses and detector inefficiencies the certifiable entanglement depth drops to about $10^{5}$ (i). In our opinion, a conservative but reasonable number of the certified entanglement depth is $10^{7}$. Indeed, on the one side the entanglement depth of $10^{9}$ in (iv) relies more on our theoretical model than on our data, while on the other side the $10^{5}$ value in (i) takes into account well-understood losses that are not part of the physical phenomenon we like to certify.

II Discussion

This work demonstrates that large entanglement depth is experimentally certifiable even with atomic ensembles beyond $10^{10}$ atoms and low detection and re-emission efficiencies. We prove that entanglement between many atoms is necessary for the functioning of quantum memories that are based on collective emission, because the combination of directed emission (i.e., high memory efficiency) and preservation of the single-photon character imply large entanglement depth.

Our results further illustrate the fundamental difference between various manifestations of large entanglement. The scales at which we observe entanglement depth seem to be completely out of reach for other types of large entanglement, such as Schrödinger-cat states Pezzè et al. (2016); Fröwis (2017).

As detailed in the Methods, our reasoning is based on two steps. First, a model-independent witness for entanglement depth is derived, which only depends on the overlap of the atomic state with $|D_{1}\rangle$ and $|D_{2}\rangle$ as well as the total number of atoms, $N$. In the experiment, we measure the probabilities $p_{1}$ and $p_{2}$ for one and two photons, respectively, emitted from the atomic ensemble. The second step consists in identifying $p_{1}$ and $p_{2}$ with the probabilities of the atomic ensemble being in the $|D_{1}\rangle$ and $|D_{2}\rangle$ state before the re-emission, respectively. This step as well as the measurement of $N$ are based on some assumptions regarding the atomic ensemble, the single-photon source and the light-matter interaction, Eq. (1). In addition, our claimed entanglement depth in the order of $10^{7}$ takes finite detector efficiencies into account. We emphasis, however, that these assumptions have been thoroughly tested in the classical and quantum regime in many previous experiments. Further note that the entanglement depth is generated by a probabilistic but heralded source. Hence no post-selection has been made in our experiment.

We report lower bounds on the minimal number of genuinely entangled atoms, which should not be confused with quantifying entanglement with an entanglement measure. Indeed, the nature of the target state, the W state $|D_{1}\rangle$, and the experimental challenges suggest that only a small amount of entanglement is present in the crystal during the storage.

We note that entanglement between many large groups of atoms in a solid is demonstrated in a parallel submission by P. Zarkeshian et al., where the coherence between these groups is revealed by analyzing the temporal profile of the re-emission in the forward direction.

**Acknowledgments ** We thank Philipp Treutlein, Otfried Gühne, Christoph Simon, Parisa Zarkeshian, Khabat Heshami and Wolfgang Tittel for useful discussions. This work was supported by the European Research Council (ERC-AG MEC) and the Swiss National Science Foundation (No. 200021_149109 and Starting grant DIAQ).

**Author Contributions ** F.F., N.B., F.B. and M.A. conceived the basic concept. F.F., C.G., N.B. and N.G. worked out the theory, supported by M.A. P.C.S., A.T. and J.L. performed the experiment and analyzed the data. F.F., A.T., N.B. and F.B. wrote the paper. All authors discussed the results and implications and commented on the manuscript at all stages.

III Methods

III.1 Heralded single photon source

The single photon used to prepare the entangled state of the atomic ensemble is generated using spontaneous parametric down conversion (SPDC). A 2 mW monochromatic continuous-wave 532 nm laser pumps a periodically poled potassium titanyl phosphate (ppKTP) waveguide to generate the signal and idler photons at 883 nm and 1338 nm, respectively. The two down-converted photons are energy-time entangled. The narrow spectral filtering of the signal (idler) photon is performed using a Fabry-Perot cavity with a linewidth of 600 MHz (240 MHz) Clausen et al. (2014). The detection of the idler photon heralds the presence of a signal photon in a well defined spectral, temporal and polarization mode. The heralded single photon is then absorbed in the quantum memory. The zero-time second-order autocorrelation of the heralded single photon before storage in the quantum memory was measured to be $g_{ss|i}=0.0055(2)$ when using a 1.2 ns coincidence window. The idler mode is detected by an InGaAs/InP single-photon detector ID220 from ID Quantique (20% detection efficiency), while the signal mode is analysed using silicon avalanche photodiodes from Perkin Elmer (30% detection efficiency).

III.2 Solid-state quantum memory

The single-photon storage is performed using a broadband and polarization-preserving quantum memory Clausen et al. (2012) realized by placing two 5.8 mm-long 75 ppm Nd3+:Y2SiO5 crystals around a 2 mm-thick half-wave plate. An atomic frequency comb is prepared on the center of the Nd3+ ions transition at 883 nm (absorption line ${}^{4}I_{9/2}\longrightarrow{}^{4}F_{3/2}$) by optical pumping. Using optical path consisting of acousto-optic and phase modulators we create a 600 MHz comb with a spacing of 20 MHz between the absorption peaks Tiranov et al. (2015). The resulting optical depth of the absorption peaks is $d=2.0\pm 0.1$. This value was doubled using a double pass propagation through the crystals (Fig. 1(b)). The overall single photon efficiency of the AFC quantum memory is 7(1)% with a 50 ns storage time. The absorption probability of one photon by the crystal was estimated to be 82(1)%, which was obtained from the probability for the single photon to be transmitted.

III.3 One- and two-photon probabilities

The one-photon and two-photon probabilities in the forward mode are obtained as follows. First, the transmission probability along the path from the photon pair source to the single-photon detectors was carefully estimated using heralded single photons. The overall transmission consists of (i) the heralding probability of the single photon before the quantum memory 19(1)%; (ii) the overall quantum memory efficiency 7(1)%; and (iii) the detection efficiency including the transmission of the system 17(1)%. From this we found that the total probability to detect a single photon reemitted from the crystal is $p_{1}=2.3(3)\times 10^{-3}$. Then, the probability $p_{2}$ can be estimated from a measurement of the zero-time second-order autocorrelation function $g^{(2)}_{ss|i}=2p_{2}/p_{1}^{2}$. We measured $g^{(2)}_{ss|i}=0.020(3)$. This value is higher than the one before the storage, which is due to spurious noise coming from the photon pair source Usmani et al. (2012). The probability to detect two photons is estimated to be $p_{2}=5(2)\times 10^{-8}$.

To connect the experimental observation of one and two photons with $p_{1}=\left\langle D_{1}\right|\rho\left|D_{1}\right\rangle$ and $p_{2}=\left\langle D_{2}\right|\rho\left|D_{2}\right\rangle$, respectively, some details have to be clarified. First, we note that the actual coupling between light and atoms is not uniform as in Eq. (1) due to position dependent field intensities and the inhomogeneous broadening of the atomic ensemble. However, it is possible to mathematically replace this by an ideal, uniform coupling with a reduced ensemble size Hu et al. (2015). We expect that replacing an ensemble by a smaller one only lowers the bounds on entanglement depth. Since $N\gg 1$ and the weak coupling between the field and a single atom, the dynamics from Eq. (1) are well approximated by a first-order expansion of the Holstein-Primakoff transform Holstein and Primakoff (1940), that is, the linear regime

[TABLE]

where $\eta$ is the transfer efficiency and $S_{\pm}^{\vec{k}}=\sum_{j=1}^{N}e^{\mp i\vec{r}_{j}\cdot\vec{k}}\sigma_{\pm}^{(j)}$ are the creation and annihilation operators for a collective atomic excitation (up to the normalization factor $N^{-1/2}$; see Hammerer et al. (2010)). All formulas in this paper are based on this approximation and the next-order correction $O(1/N)$ is omitted.

Second, the states $|D_{1}\rangle$, Eq. (2), and $|D_{2}\rangle$, Eq. (4), are not the only ones that give rise to the emission of one and two photons, respectively. Let us introduce the canonical basis $\{\left|j,m,\alpha\right\rangle_{\vec{k}}\}_{j,m,\alpha}$ for the angular momentum operators $S_{z}=\frac{1}{2}[S_{+}^{\vec{k}},S_{-}^{\vec{k}}]$ and $S^{2}_{\vec{k}}=S_{z}^{2}+\frac{1}{2}\{S_{+}^{\vec{k}},S_{-}^{\vec{k}}\}$, where $S^{2}_{\vec{k}}\left|j,m,\alpha\right\rangle_{\vec{k}}=j(j+1)\left|j,m,\alpha\right\rangle_{\vec{k}}$ and $S_{z}\left|j,m,\alpha\right\rangle_{\vec{k}}=m\left|j,m,\alpha\right\rangle_{\vec{k}}$. The third quantum number $\alpha$ labels the degeneracies for asymmetric states. For $\eta=1$, any atomic state $\left|j,-j+l,\alpha\right\rangle_{\vec{k}}$ with $N/2-j=O(1)$ (low-excitation limit) and $\vec{k}$ the forward mode transforms via Eq. (5) to the photonic Fock state $|l\rangle$. Now, assuming that the single photon source is the only source of coherent excitation, the population of the subspaces $N/2+m$ is strongly decaying with $m$, such that the overlap of the atomic state before re-emission with $|D_{1}\rangle\equiv\left|N/2,-N/2+1,1\right\rangle_{\vec{k}}$ is much larger than with the entire subspace spanned by $\{\left|j,-j+1,\alpha\right\rangle_{\vec{k}}\}_{j<N/2}$ (i.e., the nonsymmetric subspace emitting a single photon). Using the $g^{(2)}_{ss|i}$ directly measured at the source, it can be estimated that corrections taking the nonsymmetric subspace into account are much smaller than the uncertainty of the measured $p_{1}$. A similar argument applies to the two-photon emission.

Finally, the memory preparation ideally sets the atomic ensemble to the ground state. In practice, we estimate that roughly $10^{-5}\times N$ atoms are at the end of the preparation phase in the excited state without any phase coherence between them. In the linear regime, these excitations simply drop out from all calculations and can hence be safely ignored.

III.4 Number of atoms in the atomic ensemble

Another parameter is the number of atoms $N$ participating to the collective atomic mode, which was estimated with a separate measurement. The ratio between coherent (signal) and incoherent (noise) emission from the atomic ensemble was used to estimate the number of atoms coupled to the optical mode. A simple model for this signal-to-noise ratio (SNR) was developed, where the number of atoms $N$ is a free parameter. By independent characterization of the remaining other parameters and by measuring the SNR, we obtain and estimate of $N$. Intuitively, this is based on the fact that the re-emission from the atomic ensemble is enhanced in the spatial mode of the incident single photon, due to the constructive interference between all the atoms which have collectively absorbed the single photon. This probability ideally equals to $Np_{s}$, where $p_{s}$ is the probability for spontaneous emission of a single atom Scully et al. (2006). In any other mode, including the backward mode, there is no collective enhancement, and the probability of an incoherent re-emission is just $p_{s}$. Hence, the SNR is given by $N$ in the ideal case where no other source of noise is present.

In principle, one could measure the SNR from the probability to detect the heralded single photon in the backward mode. In practice, this cannot be done because $N\sim 10^{10}$. The incoherent re-emission probability is extremely small and is therefore lost in the noise due to detector dark counts and spurious light. To overcome this limitation, strong coherent state pulses with mean photon number $\left|\alpha\right|^{2}$ up to $10^{6}$ were used instead to estimate the SNR. This value is still much lower that the total number of atoms which keeps the interaction in linear regime. In this case the noise becomes less important and the true incoherent re-emission can be measured. To detect it with a low noise level we used a Picoquant silicon avalanche photodiode detector with 35% efficiency and 4 Hz dark count rate.

To perform the SNR measurement the forward $\vec{k}_{f}$ and the backward $\vec{k}_{b}$ spatial modes were used to measure signal and noise, respectively (see Fig. 1(b)). For each mode spatial filtering was performed using single-mode fibers. This allowed us to confirm that modes in both directions (forward $\vec{k}_{f}$ and backward $\vec{k}_{b}$) are probing the same volume of the QM. For this the light was sent in both directions and the coupling was aligned simultaneously for both couplers after the PBS, as shown in Fig. 1(b). Furthermore, the incoherent re-emission from the strong coherent state pulses was detected simultaneously in both modes. By applying corrections for diverse optical losses the ratio between the intensities in both modes was found to be very close to 1 as expected (see the Supplementary Information). This confirms that the forward $\vec{k}_{f}$ and backward $\vec{k}_{b}$ modes correspond to the coherent and incoherent re-emission modes defined by our model. Note that any systematic error that leads to an underestimated signal or to overestimated noise only reduces the inferred $N$ and hence leads to an underestimated entanglement depth. For example, scattered photons from the $\vec{k}_{f}$ mode that could have been mistakenly collected in the $\vec{k}_{b}$ mode would increase the noise and hence would lead to an underestimated $N$. We are not aware of any systematic error in our experiment that would let us overestimate $N$.

Using a coherent state pulse with a mean number of photons equal to $|\alpha|^{2}$, the signal is proportional to $\eta|\alpha|^{2}$, where $\eta$ is the re-phasing efficiency of the quantum memory. The incoherent re-emission in the backward mode is proportional to $|\alpha|^{2}/N+\delta$, where $\delta$ is a noise probability. Hence, the SNR is given by $\eta|\alpha|^{2}/(|\alpha|^{2}/N+\delta)$, from which $N$ can be obtained (see Supplementary Information). From this, the number of atoms was found to be $N=4.0(1)\times 10^{10}$. An estimate for the number of atoms obtained by considering the doping concentration, the length of the crytals and the size of the optical mode roughly gives $3\times 10^{11}$, which confirms at least the order of magnitude. Note, however, that the latter method comes with much larger uncertainties and therefore we rely only on the first number.

III.5 Ansatz for $M$-separability

We now give some details for the derivation of a lower bound of $p_{2}$ given $p_{1}$ and the ansatz state (3). We start with the ansatz that for every pure state decomposition of an atomic state the pure states are separable between at least $M$ groups (where each group is genuinely multipartite) and there exists at least one pure state in every decomposition that consists of exactly $M$ separable groups. Such a state is called $M$-separable. In principle, the sizes of the groups are independent from each other as long as the total number of atoms is conserved. However, we fix the group size $K$ to be constant, that is, $MK=N$ for the following reason. Our final goal are bounds on numbers of entangled atoms. This is a “min-max” problem. For every possible state the entanglement depth is the size of the largest entangled group in the state. By varying the state, our task is to find a state such that this largest group is minimized. From this, it follows that it is best to have an equal size for all groups in order to avoid few very large groups. Clearly, if we fix $N$ and $M$, $K$ does not have to be an integer. So, generally, one has to reduce the size of one group such that $(M-1)K+K^{\prime}=N$. However, we will consider many groups such that the size of a single group is in the order of or smaller than the uncertainty of $N$. Hence a detailed analysis with $K^{\prime}<K$ is not necessary.

Since we are concerned with at most two excitations in total, it is sufficient to work with pure states of the form

[TABLE]

where $|d_{k}\rangle$ here refers to Dicke states within one group, that is, symmetric superposition of $k$ excitations. With this ansatz, the probabilities read

[TABLE]

and

[TABLE]

where $A=\Pi_{i=1}^{M}a_{i}$. In the following, we ignore the corrections $1/N$ and $1/K$. While the factor $(1-1/N)^{-1}$ is arguably negligible, dropping $\sqrt{1-1/K}$ only lowers $p_{2}$ in the relevant regime (i.e., the interval $I$ discussed later).

Formally, the task is now to minimize $p_{2}$ for given $p_{1},M$ over the parameters of the ansatz state (6), that is,

[TABLE]

When this is done for all $p_{1}$, one has to find the minimum of all convex combinations for a bound on mixed states. In our case, it will turn out that the lower bounds for pure states are already convex implying that they are also valid for mixed states. Then, one can invert the results and determine the maximal $M$ that is compatible with a given pair $(p_{1},p_{2})$. This bounds the $M$-separability of the atomic state generated in the experiment.

Since the state (6) could contain further elements not contributing to $p_{1}$ and $p_{2}$, one has in general that $|a_{i}|^{2}+|b_{i}|^{2}+|c_{i}|^{2}\leq 1$ for all $i$. Thus the minimization (9) has to be done over $3M$ complex parameters. As discussed in the Supplementary Information, one easily shows that the complexity reduces to $2M$ real parameters as the optimal state has $a_{i}\geq 0,b_{i}\in\mathbbm{R}$ and $c_{i}=-\sqrt{1-a_{i}^{2}-b_{i}^{2}}$.

III.6 Lagrange multiplier

To solve Eq. (9), we use the Lagrange multiplier method, which is suitable for constrained minimization problems. In our case, we consider

[TABLE]

with $f_{1}=\sqrt{Mp_{1}}$ and $f_{2}=\sqrt{M^{2}p_{2}}$. The formulas for the partial derivatives are given in the Supplementary Information. The results are quartic equations with four solutions for every group $i>1$ that depend on the parameters of the first group $a_{1}\equiv a$ and $b_{1}\equiv b$ for all $i$.

Note that the first solution (called the symmetric solution in the following) implies that $a_{i}=a$ and $b_{i}=b$. The probabilities in this case read

[TABLE]

and

[TABLE]

with $c=-\sqrt{1-a^{2}-b^{2}}$.

Fixing $a,b$ of the first group determines the four possible solutions for every other group. A priori, the unknowns $(a,b,\lambda)$ can be found by solving the remaining equations $\partial f/\partial a=0$, $\partial f/\partial b=0$ and $\partial f/\partial\lambda=0$. Given the complexity of the equations, this is analytically not possible. Alternatively, one chooses for each group $i=2,\dots,M$ one out of four solutions, and solve Eq. (9) numerically for $(a,b)$. As a result, we find many local extrema, from which one has to choose the global minimum. Hence, we reduced the minimization problem to a finite set of possibilities. The problem is the large number of solutions. Due to symmetry, it is sufficient to determine the number of groups, $m_{j}$, that are chosen to take solution $j=1,\dots,4$ with the constraints $m_{1}>0$ and $\sum_{j}m_{j}=M$ (see Supplementary Information). In the following, we call a certain choice $C=(m_{1},m_{2},m_{3},m_{4})$ a configuration. The number of configurations scales as $O(M^{3})$.

It turns out that most of the configurations are not relevant for the following reason. For all $M$, there exist states such that $p_{1}>0$ and $p_{2}=0$. Two analytic examples are (i) $(a_{1},b_{1})=(0,1)$ and $(a_{i>1},b_{i>1})=(1,0)$, resulting in $p_{1}^{\lim 1}=1/M$ and (ii) a symmetric state where one maximizes $p_{1}$ given that Eq. (12) vanishes. The solution $p_{1}^{\lim 2}$ is a long and analytic expression which scales as $p_{1}^{\lim 2}=(eM)^{-1/2}+O(M^{-1})$. For $M>4$, $p_{1}^{\lim 2}>p_{1}^{\lim 1}$. Furthermore, the maximal $p_{1}$ is given by $p_{1}^{\max}=(1-1/M)^{M-1}$. We conclude that there is a nontrivial interval $I=[\max(p_{1}^{\lim 1},p_{1}^{\lim 2}),p_{1}^{\max}]$ where we look for the minimal $p_{2}$. In the Supplementary Information, it is shown that for large $M$, only the symmetric solution lies in $I$. With a numerical study (an unconstrained maximization over $(a,b)$), we find that this happens when $M>53$, but already for $M\gtrsim 30$, we observe that $p_{1}^{\max}-p_{1}^{\lim 2}\ll 1$ for all nonsymmetric configurations.

Numerically, we find that for all groups $M\geq 5$, the symmetric solution gives the global minimal $p_{2}$. From Eq. (11), we obtain $b^{2}=p_{1}/Ma^{2-2M}$ and insert this into Eq. (12), which has to be minimized. With the parameters $p_{1}$ and $M$, this is a single-parameter polynomial and hence a numerically stable minimization is possible. The example $M=3$ is discussed in the Supplementary Information to demonstrate a case where a nonsymmetric configuration realizes the global minimum.

Supplementary Information

.7 Fluorescence measurement

The measurement of the incoherent reemission (the fluorescence) in the backward $\vec{k}_{b}$ and the forward $\vec{k}_{f}$ modes shown on Fig 1(b) in the main text required several corrections. First, the coupling efficiencies were measured independently to be equal to more than 70% inside the optical fibre. A slight imbalance in the couplings was taken into account. Another correction was necessary due to the non-perfect transmission of the narrow bandpass filter (FWHM of 10 nm at 883.2 nm wavelength) that was used in the backward mode. The transmission was measured to be 65% at the wavelength of the heralded single photon. Finally, another correction is attributed to the PBS used in the backward mode. To apply it the polarization state of the incoherent re-emission has to be characterized. Due to the use of polarization preserving quantum memory (consisting of two crystals separated by the half-wave plate Clausen et al. (2012)) and a bow-tie configuration, the polarization state of the fluorescence emitted in backward and forward is close to a mixed state. To confirm it we performed a polarization state tomography of the emitted fluorescence. We obtained the purity of 51% and the fidelity of 95% with respect to the completely mixed polarization state. Based on this result a correction of 0.5 was used in all measurements.

In a single crystal, the absorption probability for every spatial position inside the crystal is not the same. Hence, the total probabilities to have incoherent reemission (the fluorescence) in the forward or backward directions are not equal. Assuming the total optical depth of the crystal is $d$ one can write the ratio $J$ between emission rate between two modes as

[TABLE]

However, due to the use of the double-pass configuration in our experiment, this asymmetry disappears. More precisely, the amount of the emitted fluorescence in forward and backward modes should be equal. This was verified by the direct and simultaneous measurement of the fluorescence in both the forward and backward modes. We found a ratio of at most $J=1.03(10)$ over a time interval of more than 800 µs after the absorption of the optical pulse, as shown on Fig. 3(a).

To measure the SNR the strong coherent state pulses were coupled to the QM prepared in the crystal. Generally, the SNR is fundamentally limited by the number of atoms $N$ involved to the collective re-emission process and can be expressed as

[TABLE]

where $\left|\alpha\right|^{2}$ is the absorbed mean photon number of the coherent state pulse, $\eta$ is the rephasing efficiency of the QM and $\delta$ is the contribution from the intrinsic noise of the detection system. This expression is correct for a single coherent state pulse or in the limit of the low repetition rate $R$ of the experiment ($R\ll 1/T_{1}$, $T_{1}$ being the spin relaxation time). In the general case, taking into account accumulation of the population in the excited state which comes from many coherent state pulses one can write SNR as

[TABLE]

This expression can be used to verify the source of the fluorescence and finally to estimate the number of atoms that are contributing to the collective emission. First, for the fixed repetition rate of $R=1$ MHz, the input intensity of the coherent state pulses was varied to see the influence of the noise from the detection system $\delta$ (Fig. 3(b)). The maximum measured SNR value was 72 dB for the high mean photon number at the input where the detector’s noise contribution is negligible. While decreasing $\left|\alpha\right|^{2}$ its contribution starts to be more dominant which decreases SNR. The model based on Eq. (15) and the independently measured parameters explains well the obtained data.

Next, the repetition rate $R$ was varied for the fixed $\left|\alpha\right|^{2}$ which is high enough to make the detector’s noise contribution negligible (Fig. 3(c)). In this case higher values of the SNR can be obtained due to the lower accumulation of the population in the excited state (such that $R\ll 1/T_{1}$ which corresponds to the highest SNR value). The fluorescence lifetime $T_{1}$ was measured to be 250 µs from separate measurement. The only free parameter to fit the data is the number of atoms $N$, which was estimated to be $106.0(1)\leavevmode\nobreak\ \text{dB}\approx 4.0(1)\times 10^{10}$ atoms. The maximal SNR for the measurement with the heralded single photon is expected to be higher since the spectral bandwidth of the AFC structure in this case is five times larger. This means that the number of atoms that contribute to the absorption of a single photon is larger than the measured $N$. Nevertheless, we do not correct the measured SNR but directly work with $N$.

.8 Simplification of the minimization problem

Equation (9) is a constrained minimization problem over $3M$ complex numbers. Here, we show how to reduce the complexity with a few simple arguments.

First, note that state Eq. (6) could be subnormalized, because we neglect populations in other subspaces. However, it is straightforward to see that subnormalized states give the same $(p_{1},p_{2})$ values as the renormalized state mixed with the ground state. Hence, it is sufficient to consider normalized states only.

Second, we choose $a_{i}\geq 0$ without loss of generality. For the phases of the $b_{i}$ and $c_{i}$, note that the $b_{i}$ dependent term in Eq. (8) can be written as

[TABLE]

Let us write $b_{i}=e^{i\varphi_{i}}|b_{i}|$ and $\sum_{i}b_{i}/a_{i}=e^{i\bar{\varphi}}|\sum_{i}b_{i}/a_{i}|$. Then, with $c_{i}=e^{i\vartheta_{i}}|c_{i}|$, Eq. (8) reads

[TABLE]

Clearly, the phases have to be set to $\varphi_{i}=\bar{\varphi}$ and $\vartheta_{i}=\pi+2\bar{\varphi}$ in order to minimize $p_{2}$. Without loss of generality, we set $\bar{\varphi}=0$ implying that $b_{i}\in\mathbbm{R}$ and $c_{i}\leq 0$. Note that the case where $\sum_{i<j}\frac{b_{i}b_{j}}{a_{i}a_{j}}<0$ is not interesting here because it is only possible for so-called subradiant states, that is, states with lower intensity in forward direction than the incoherent emission.

To summarize, we have $a_{i}\geq 0$, $b_{i}\in\mathbbm{R}$ and $c_{i}=-\sqrt{1-a_{i}^{2}-b_{i}^{2}}$, thus reducing the problem to $2M$ real parameters. The simplified formulas read

[TABLE]

and

[TABLE]

.9 Formulas from the Lagrange multiplier

The partial derivatives of Eq. (10) are

[TABLE]

[TABLE]

and

[TABLE]

From $\partial f/\partial b_{i}=0$, we find

[TABLE]

which we insert into $\partial f/\partial a_{i}=0$ and find

[TABLE]

Equation (23) can also be written as

[TABLE]

We notice that from $\partial f/\partial a_{i}=0$ and $\partial f/\partial b_{i}=0$ we find Eqs. (24) and (25), where the right hand side for both is independent of $i$. Hence, it follows that

[TABLE]

and

[TABLE]

for all pairs $(i,j)$. Let us fix $j=1$ and write $a\equiv a_{1}$, $b\equiv b_{1}$ and $c\equiv c_{1}$. Equations (26) and (27) thus give two equations for two unknowns $(a_{i},b_{i})$ (recall that $c_{i}$ is just a function of $a_{i},b_{i}$).

Inserting Eq. (26) into Eq. (27) allows to eliminate $b_{i}$. The remaining equation is a polynomial of degree four in $a_{i}^{2}$. One finds the four solutions

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The solutions for $b_{i}$ and $c_{i}$ follow accordingly.

.10 Asymptotic formula for $p_{1}$

To argue that for large $M$ the only relevant configuration is the symmetric one, we calculate $p_{1}$ for arbitrary configurations and show that all but the symmetric configuration are asymptotically outside the relevant interval $I$ (see Methods). To this end, note the factor $A^{2}=\Pi_{i=1}^{M}a_{i}^{2}$ in the Eqs. (18) and (19). For large $M$, this implies that almost all $a_{i}$ have to be very close to one. However, solutions two to four in Eq. (28) are close to zero for $a$ close to one. Asymptotically, we expect that $a_{i}=1-O(1/M)$ to have a finite A. To see this, consider $b=\beta/\sqrt{M}$ and $c=-\gamma/\sqrt{M}$ with $\beta,\gamma=O(1)$ and do a Taylor series of the solutions Eq. (28) around $1/M=0$. One finds

[TABLE]

and

[TABLE]

We therefore see that only configurations with $\delta m=\sum_{j=2}^{4}m_{j}=O(1)$ asymptotically give finite values of $A$. More explicitly, inserting Eqs. (32) and (33) into $p_{1}$ for a fixed configuration $C$ with $\delta m=O(1)$ gives

[TABLE]

A simple calculation shows that the maximum of Eq. (34) over all $\beta,\gamma$ and $\delta m>0$ is $p_{1}^{\max}=(e^{3}M)^{-1/2}+O(1/M)$. This means that for large enough $M$, all but the symmetric configuration exhibit a $p_{1}^{\max}$ that do not enter the nontrivial interval $I$ (which has a lower bound $(eM)^{-1/2}+O(1/M)$).

.11 Example $M=3$

We discuss an example where a nonsymmetric solution partially constitutes the global minimum. This is because for $M=2,3,4$, $p_{1}^{\lim 1}>p_{1}^{\lim 2}$ (see Methods). For $M=3$, one has $p_{1}^{\lim 1}=1/3$ and $p_{1}^{\lim 2}\approx 0.29$.

In Fig. 4, we compare the “full” numerical constrained minimization of Eq. (19) without the Lagrange multiplier method with a constrained minimization over $(a,b)$ for two relevant configurations $C$. We see that the kink for the full minimization is nicely explained by the crossing of the minimization for two different configurations.