Introduction {#Sec1}
============

The question whether quantum superposition and entanglement can exist in macroscopic systems has been at the heart of foundational debates since the beginnings of quantum theory^[@CR1]--[@CR5]^ and inspired many experiments^[@CR6]--[@CR16]^. One important type of entangled states corresponds to a single excitation that is delocalized over many systems. Such Dicke states^[@CR17]^ have been created for individual photons^[@CR18],\ [@CR19]^ and cold atoms^[@CR20]--[@CR22]^, reaching up to three thousand atoms^[@CR23]^. In solids, superradiance associated with a Dicke state was recently demonstrated for two superconducting qubits^[@CR24]^.

It is well known that, in the ideal case, a Dicke state is created whenever a single photon is stored in an atomic ensemble^[@CR25]^. However, in real experiments, neither the initial single-photon state nor the quantum storage process is perfect. It is not a priori obvious whether the multiparty entanglement will survive under these conditions. Atomic-ensemble-based quantum memories have been studied intensively both in atomic gases^[@CR26]--[@CR32]^ and in solid-state systems^[@CR33]--[@CR37]^. One widely used quantum storage method is the atomic frequency comb (AFC) quantum memory^[@CR38]--[@CR41]^. An AFC is a collection of atomic ensembles with different, equally spaced resonance frequencies. Such AFCs can be conveniently generated in rare earth ion-doped crystals through optical pumping. Each individual ensemble represents one "tooth" of the comb.

Here we create entanglement between these teeth by the absorption of a single photon. More precisely, there is a small probability of absorbing more than one photon, and it is essential to take this into account when attempting to quantify the multiparty entanglement. We derive a criterion that allows us to show that multiparty entanglement of over two hundred ensembles is present for our experimental conditions.

Results {#Sec2}
=======

For most of this work, we focus on the entanglement that is generated between the teeth, rather than within each tooth. It is then possible to treat each tooth as a single two-level system (a "qubit") with collective states $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| 0 \right\rangle$$\end{document}$ (where all atoms in the tooth are in the ground state) and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| 1 \right\rangle$$\end{document}$ (where a single atom in the tooth is excited). Ideally, the absorption of a single photon by an AFC consisting of *N* teeth creates the lowest-order Dicke state, widely known as the W state,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _N} = \frac{1}{{\sqrt N }}\left( {\left| {100...0} \right\rangle + \left| {010...0} \right\rangle ... + \left| {000...1} \right\rangle } \right),$$\end{document}$$where the first term corresponds to the case in which the first tooth has absorbed the photon, etc. Our theoretical and experimental approach to demonstrating this multiparty entanglement is based on the re-emission of the single photon from the AFC, which is due to a collective interference effect. As time passes, the different terms in the above equation acquire different phases $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${{\rm{e}}^{i{\delta _j}t}}$$\end{document}$, where *δ* ~*j*~ = *j* *Δ* stands for the detuning of the *j*th tooth relative to the lowest-frequency (*j* = 0) tooth, *j* runs from 0 to *N* − 1, and *Δ* represents the angular frequency spacing between the teeth. The photon has a high probability of being re-emitted only at the "echo" times when all the phase factors are the same, i.e., at times that are integer multiples of 2*π*/*Δ* ^[@CR38],\ [@CR39]^. This echo emission is thus an interference effect in time that is very similar to spatial multi-slit diffraction (see Fig. [1](#Fig1){ref-type="fig"}). This analogy is at the heart of our approach, which we now describe in detail.Fig. 1Principle of our approach. Our approach is based on the analogy between re-emission of a photon from an atomic frequency comb and multi-slit diffraction. **a** Light passing through a mask with a spatially periodic structure, i.e., a diffraction grating, travels different optical path lengths depending on the part of the grating through which it was transmitted. This results in sharp constructive interference and broad destructive interference in momentum space. **b** Light directed into an atomic ensemble with a spectrally periodic absorption profile (AFC) is absorbed in these atoms, causing them to oscillate at their (different) resonant frequencies. This results in an interference in time, manifested via sharp peaks in the re-emission probability at well-defined times ("echos"). The absorption "teeth" in the AFC are analogous to the slits in the diffraction grating. **c** Principle of the experiment. A single photon is created with the help of a photon pair source and a heralding detector. The single photon is stored in an AFC and detected after retrieval from the AFC. The echo contrast in combination with the single-photon character of the source can be used to find a bound on the minimum number of entangled teeth

Derivation of a lower bound for the entanglement depth {#Sec3}
------------------------------------------------------

Our method for demonstrating the entanglement between the teeth is inspired by the fact that, for multi-slit interference, the number of participating slits can be inferred from the sharpness of the interference pattern. We consider the ratio of the maximum photon re-emission probability in the first echo to the re-emission probability averaged over one period 2*π*/*Δ*,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$R = \frac{{P\left( {\frac{{2\pi }}{{\Delta }}} \right)}}{{\frac{{\Delta }}{{2\pi }}{\int}_{{\textstyle{\pi \over {\Delta }}}}^{{\textstyle{{3\pi } \over {\Delta }}}} P(t){\rm{d}}t}}.$$\end{document}$$We refer to *R* as the "echo contrast". The photon emission probability *P*(*t*) is proportional to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left\langle {{S_ + }(t){S_ - }(t)} \right\rangle$$\end{document}$ ^[@CR39],\ [@CR42]^, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${S_ - }(t) = \mathop {\sum}\nolimits_j {{e^{i{\delta _j}t}}S_- ^j}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$S_ - ^j = {\left| 0 \right\rangle ^j}{\left\langle 1 \right|^j}$$\end{document}$ is the dipole operator for tooth *j*, and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${S_ + }(t) = S_ - ^\dag (t)$$\end{document}$. Here we are using the Heisenberg picture, where observables rather than states are time dependent. At the time of the first echo, *t* ~e~  *=* 2*π*/*Δ*, all teeth are in phase, giving $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${S_ - }\left( {{t_{\rm{e}}}} \right) = \mathop {\sum}\nolimits_j {S_ - ^j \equiv {S_ - }}$$\end{document}$. On the other hand, averaging the re-emission probability over a time interval 2*π*/*Δ* centered at the echo time as in Eq. ([2](#Equ2){ref-type=""}) leads to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\frac{{\Delta }}{{2\pi }}{\int}_{{\textstyle{\pi \over {\Delta}}}}^{{\textstyle{{3\pi } \over {\Delta}}}} {\rm{d}}t\left\langle {\mathop {\sum}\nolimits_{j,l} {S_ + ^jS_ - ^l{{\rm{e}}^{i(l - j)\Delta {\kern 1pt} t}}} } \right\rangle ,$$\end{document}$ which is non-zero only for *j* = *l*. Then, the denominator of Eq. ([2](#Equ2){ref-type=""}) results in the expression $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left\langle {\mathop {\sum}\nolimits_j {{{\left| 1 \right\rangle }^j}{{\left\langle 1 \right|}^j}} } \right\rangle ,$$\end{document}$ i.e., the sum of the excitation probabilities for each tooth, yielding $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$R = \frac{{\left\langle {{S_ + }{S_ - }} \right\rangle }}{{\left\langle {\mathop {\sum}\nolimits_j {{{\left| 1 \right\rangle }^j}{{\left\langle 1 \right|}^j}} } \right\rangle }}$$\end{document}$.

We now show that the echo contrast *R* is closely related to the "entanglement depth"^[@CR43]^. A (generally mixed) quantum state of *N* qubits has entanglement depth at least equal to *M* if it cannot be decomposed into a convex sum of product states with all factors involving less than *M* entangled qubits, i.e., at least one of the terms needs to be an *M*-qubit entangled state.

First, we consider the case where exactly one photon is absorbed by the AFC. In this case, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$R = \left\langle {{S_ + }{S_ - }} \right\rangle .$$\end{document}$ Let us suppose that *R* is found experimentally to have a value of *M*. This value can be achieved by the state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M} \otimes {\left| 0 \right\rangle ^{ \otimes \left( {N - M} \right)}},$$\end{document}$ i.e., a Dicke state involving *M*-teeth and no excitation in the remaining *N* − *M*-teeth. Let us note that *S* ~+~ *S* ~−~ is permutation invariant (as is *R* in the general case), so all permutations of the teeth are equivalent for our purpose. The above state has an entanglement depth of *M*. All other states in the single-excitation subspace giving *R* = *M* involve more than *M* entangled teeth, i.e., they are of the form $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\mathop {\sum}\nolimits_{j = 1}^N {{c_j}\left| {{{0...1}_j}...0} \right\rangle }$$\end{document}$ with more than *M* non-zero coefficients *c* ~*j*~ that are not all equal and fulfilling $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left\langle {{S_ + }{S_ - }} \right\rangle = {\left| {\mathop {\sum}\nolimits_{j = 1}^N {c_j}} \right|^2} = M.$$\end{document}$ This is due to the fact that, in the single-excitation subspace of the *M*-teeth Hilbert space, the Dicke state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M}$$\end{document}$ is the only eigenstate of *S* ~+~ *S* ~−~ with a non-zero eigenvalue (namely *M*), see "Methods" section. Thus an experimental value *R* = *M* implies an entanglement depth of at least *M* in the case of perfect absorption of a single photon.

Next, we consider the more realistic case of absorption of a single photon with probability *P* ~1~ \< 1. In this case, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$R = \left\langle {{S_ + }{S_ - }} \right\rangle {\rm{/}}{P_1}.$$\end{document}$ Let us again suppose that the experiment gives *R* = *M*. This can be achieved with a mixed state or coherent superposition state that has a probability weight *P* ~1~ for $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M} \otimes {\left| 0 \right\rangle ^{ \otimes \left( {N - M} \right)}}$$\end{document}$ and a weight 1 − *P* ~1~ for $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| 0 \right\rangle ^{ \otimes N}}.$$\end{document}$ With the same line of arguments as for the case of the perfect absorption (*P* ~1~ = 1), the entanglement depth is again at least *M*.

The situation becomes more complex if there is a chance for the AFC to absorb more than one photon. In our experiment, the most significant higher-order probability is *P* ~2~, the probability of absorbing two photons, which is very small but non-zero. We now have to consider states with higher-order components, e.g., the fully separable state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left( {\alpha \left| 0 \right\rangle + \beta \left| 1 \right\rangle } \right)^{\mathop { \otimes }\nolimits_ N}}$$\end{document}$ with *α* ^2^ + *β* ^2^ = 1. Since this state has the Dicke state of *N* teeth, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _N},$$\end{document}$ as its component in the single-excitation subspace, it can give very large values of the echo contrast up to *R* ≃ *N*, which happens in the limit of small *β*. On the other hand, in this limit, the two-excitation probability is $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P_2} \simeq P_1^2{\rm{/}}2.$$\end{document}$ In our experiment, the values of *P* ~1~ and *P* ~2~ satisfy $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P_2} \ll P_1^2{\rm{/}}2,$$\end{document}$ see below. It is clear from this example that the values of *P* ~1~ and *P* ~2~ are important for deciding whether large values of *R* imply large values of the entanglement depth. We have numerically determined a lower bound for the entanglement depth, *M*, as a function of the echo contrast, *R*, conditioned on the experimental values of *P* ~1~ and *P* ~2~ (see "Methods" section). This bound can be very well approximated in our regime by the relation:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$M  >R - \frac{{\sqrt {2{P_2}} }}{{{P_1}}}N.$$\end{document}$$

Note that the bound gives *M* \> 0 (which is consistent with the correct value *M* = 1) for the fully separable state discussed above, for which *R* = *N*, but also $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\frac{{\sqrt {2{P_2}} }}{{{P_1}}}$$\end{document}$ = 1. Below we utilize the obtained numerical bound to find the minimum entanglement depth corresponding to the experimental value of the echo contrast.

Experimental results {#Sec4}
--------------------

The principle of our experimental approach is shown in Fig. [1c](#Fig1){ref-type="fig"}, and the detailed experimental set-up is shown in Fig. [2a](#Fig2){ref-type="fig"}. It consists of a source of heralded single photons, an AFC, and a time-resolved detector to register the retrieved photons. The heralded single-photon source is implemented by spontaneous parametric down conversion (SPDC), which probabilistically down converts photons at 523.5 nm into pairs of photons with wavelengths centered around 795 and 1532 nm. Provided a single pair was generated, detection of a 1532 nm photon heralds the presence of a 795 nm photon. Measuring the cross-correlation function of the source, shown in Fig. [2b](#Fig2){ref-type="fig"}, allows us to determine the mean photon number per mode *μ*, where the mode is defined by the coherence time of the pump laser. We find *μ* = (1.1 ± 0.1) × 10^−3^ (see "Methods" section). This low mean photon number indicates that the probability of generating multi pairs is very low (of the order 10^−6^ per mode), thus confirming the nearly single-photon nature of our heralded source. This can also be seen from the heralded autocorrelation function in Fig. [2b](#Fig2){ref-type="fig"}. In combination with a precise estimation of the loss at various places of the set-up and the absorption efficiency in the AFC (see "Methods" section), this allows us to determine *P* ~1~ and *P* ~2~ using a detailed theoretical model of the set-up (see "Methods" section). We find *P* ~1~ = (3.50 ± 0.03) × 10^−3^ and *P* ~2~ = (2.55 ± 0.23) × 10^−8^, where the uncertainties in *P* ~1~ and *P* ~2~ are dominated by the uncertainties in the optical depth of the AFC and in the mean photon number per mode, respectively.Fig. 2Experimental set-up. **a** *Left*: Source of heralded single photons. A 1047 nm continuous wave (CW) laser is frequency doubled and subsequently downconverted with two different periodically poled lithium niobate (PPLN) nonlinear crystals. This probabilistically generates photon pairs with central wavelengths at 795 and 1532 nm, which are separated using a dichroic mirror (DM). A Fabry--Perot (FP) cavity spectrally filters the 795 nm photons to 6 GHz in order to match the spectral acceptance bandwidth of our AFC. Wave plates allow adjusting the polarization of the photons to maximize interaction with the AFC. Similarly, the 1532 nm photons are filtered down to 10 GHz using another FP cavity before being detected by a superconducting nanowire single-photon detector (SNSPD)---this detection heralds the presence of the 795 nm photon. *Right*: AFC preparation. CW laser light at 795 nm, which is frequency and amplitude modulated by an acousto-optic modulator (AOM) and a phase modulator (PM), is used to optically pump the Tm atoms in a cryogenically cooled Tm:LiNbO~3~ bulk crystal to create an AFC. Micro electro-mechanical switches (MEMS) allow switching between optical pumping and storage of single photons. An avalanche photodiode (APD) is used to detect the 795 nm photons after interacting with the AFC. The electrical signal from both detectors (APD and SNSPD) is analyzed using a time-to-digital converter (TDC) to register the difference in arrival times of the two detection signals. **b** Characterization of the single-photon source. Cross-correlation function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{ab}}}^{(2)}$$\end{document}$ of the two downconverted photons at 1532 nm (mode a) and 795 nm (mode b) as a function of the pump power. The *error bar* indicates the standard deviation derived from the Poissonian statistics of the photon detection events. We also show the heralded autocorrelation function of the 795 nm photon $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{bb}}}^{(2)}$$\end{document}$ for one value of the pump power. **c** A sample trace of an AFC (*solid line*) with *N* = 30 teeth and bandwidth *B* = 6 GHz and the spectral-density profile of the 795 nm photon (*dashed line*)

The AFC for the 795 nm photons is implemented in a bulk Tm:LiNbO~3~ crystal. Optical pumping is used to spectrally shape the inhomogeneously broadened ^3^ *H* ~6~ → ^3^ *H* ~4~ absorption line of the Tm atoms into a series of absorption peaks (teeth) spaced by angular frequency *Δ*, see Fig. [2c](#Fig2){ref-type="fig"}. The number of teeth *N* in the AFC is given by *N* = *B*∗2*π*/*Δ*, where *B* is the bandwidth of the AFC (fixed at 6 GHz); and *N* can be changed by modifying *Δ*. After a storage time of *t* ~e~ = 2*π*/*Δ*, the 795 nm photons are retrieved and detected by a single photon detector whose signal is sent to a time-to-digital converter (TDC). This generates a time histogram of the 795 nm photon detections such as the one in Fig. [3a](#Fig3){ref-type="fig"}.Fig. 3Measured echo contrast and bounds for entanglement depth. **a** An example of a time histogram for a storage time of 68 ns. The peak value of the echo (*E*), marked by the *green x*, is extracted from fitting the echo to a Gaussian profile, which reduces the impact of TDC sampling noise, and the *blue region* around the echo represents the time interval over which counts are averaged (*A*). The *red region*, which precedes both the transmitted and echo photon, is used to extract the background noise level. **b** Experimental values of echo contrast *R* as a function of the number of teeth *N*. *Blue circles* represent the values of *R* obtained from raw data, *red open circles* are obtained after noise subtraction, and *green squares* after noise subtraction and deconvolution of the detector response (see "Methods" section). *Error bars* indicate standard deviations derived from the Gaussian fitting and Poissonian statistics of the photon detection events. **c** Numerical bounds for the entanglement depth *M* as a function of echo contrast *R* for the experimental values *P* ~1~ = 3.5 × 10^−3^ for the single-excitation probability, *N* = 564 for the number of teeth, and *P* ~2~ = 2.6 × 10^−8^ for the double-excitation probability (*solid line*). For comparison, we also show the bounds for the same *P* ~1~ and *N*, but for *P* ~2~ = 0 (*short-dashed line*), *P* ~2~ = 2.6 × 10^−9^ (*dashed line*), and *P* ~2~ = 2 × 10^−7^ (*dash-dotted line*). The experimental values of the echo contrast *R* are shown as a *blue dot* (raw data), *red dot* (after noise subtraction), and *green square* (after noise subtraction and detector deconvolution) respectively

The echo contrast *R* can be directly extracted from the histogram by taking the ratio of the peak value *E* of the echo to the average of the counts *A* in a time interval centered on the echo and equal to the storage time, i.e., *R* = *E*/*A*. Figure [3b](#Fig3){ref-type="fig"} shows a plot of the ratio *R* as a function of the number of teeth *N* in the AFC. We see that the value of *R* obtained from the raw data (*solid blue circles*) first increases with *N*, reaches a maximum of 70.6 ± 1.4 for *N* = 408, and then saturates. The main factors that limit the value of *R* are on the one hand the limited retrieval efficiency of the AFC and detector jitter, both of which limit the peak value *E* reached by the echo, and on the other hand noise coming from multi-pair emissions of the source and detector dark counts, both of which increase the denominator *A* of the echo contrast *R*. After noise subtraction (*red circles*) and additional deconvolution of the finite detector response time (*solid green squares*), we find a maximum *R* of 127.6 ± 4.3 and 256.7 ± 8.7, respectively, both for *N* = 564. The linearity of the echo contrast with *N*, after compensating for the above effects (see "Methods" section for details), is consistent with the expectation, since greater values of *N* correspond to larger Dicke states.

In Fig. [3c](#Fig3){ref-type="fig"}, we plot our numerical bounds for the entanglement depth *M* as a function of echo contrast *R* for *N* = 564 and taking into account the experimental values of *P* ~1~ = 3.5 × 10^−3^ and *P* ~2~ = 2.6 × 10^−8^, see also Eq. ([3](#Equ3){ref-type=""}). We include other values of *P* ~2~ for comparison. (Note that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P_2} \ll {P_1}$$\end{document}$ in all cases, so *P* ~1~ is kept essentially constant.) We conclude that the noise-and-jitter-corrected echo contrast of *R* = 256.7 ± 8.7 implies at least 229 ± 11 entangled teeth, where the uncertainty in the entanglement depth is dominated by the uncertainty in *R*.

Each of the entangled teeth above consists of many atoms (of the order 10^9^, see "Methods" section). As a consequence, each state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| 1 \right\rangle$$\end{document}$ in Eq. ([1](#Equ1){ref-type=""}) is a type of Dicke state itself. To illustrate this, we first create an AFC with *N* = 9 (*red trace* in Fig. [4a](#Fig4){ref-type="fig"}) and extract a minimum entanglement depth of 5 from the data. Using only the atoms that form one of the teeth in the first AFC, we then create a second---nested---AFC, again with *N* = 9 (see Fig. [4b](#Fig4){ref-type="fig"}). We experimentally determine the minimum entanglement depth for this secondary AFC to be 4, see Fig. [4c, d](#Fig4){ref-type="fig"}. This procedure of subdividing an ensemble, while restricted by the laser linewidth in our experiment, is fundamentally only limited by the homogeneous linewidth of the atoms. In principle, it would be possible to demonstrate entanglement of a very large number of entangled subsystems (and even subsubsystems, etc.) in this way.Fig. 4Entangled subsystems within each ensemble. **a** An AFC with a bandwidth of 6 GHz with *N* = 9 teeth, where each tooth is 0.33 GHz wide. **b** The plot shows an AFC with a narrower bandwidth of 0.33 GHz that is nested within a single tooth of the broadband AFC. This AFC is also created with *N* = 9 teeth. **c** Echo contrast for the broadband and narrowband AFCs. **d** Bounds on the entanglement depth as a function of echo contrast for the experimental values of *P* ~1~ = 1.1 × 10^−2^ and *P* ~2~ = 2.4 × 10^−7^ for the broadband AFC with *N* = 9, and *P* ~1~ = 8.8 × 10^−5^ and *P* ~2~ = 1.6 × 10^−11^ for the narrowband AFC with *N* = 9. It can be inferred from these two plots that the entanglement depth is at least 5 for the broadband and at least 4 for the narrowband AFC. *Error bars* in the last two plots indicate standard deviations derived from the Gaussian fitting and Poissonian statistics of the photon detection events

Discussion {#Sec5}
==========

The remaining mismatch between the obtained values of *R* and the maximum possible value *N* can be explained by other experimental limitations, such as the imperfect creation of the AFC and the Lorentzian spectral-density profile of the photons, which limit the contribution to the interference coming from the teeth that are further detuned. The state that is actually generated in the experiment is thus of the form $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\mathop {\sum}\nolimits_{j = 1}^N {{c_j}\left| {{{0...01}_j}0...0} \right\rangle }$$\end{document}$ (not all *c* ~*j*~ being equal) rather than a perfect W state. However, this only means that our bound on the entanglement depth is conservative, because a state with unequal coefficients has to involve a greater number of entangled teeth in order to achieve the same value of *R*. The most relevant type of decoherence in our experiment is irreversible dephasing due to the finite spectral width of each tooth. This effect accumulates over time and thus affects the observed value of *R*. The entangled state at the time of absorption is therefore likely to have had a higher entanglement depth than what can be inferred from observing the echo, which is emitted some time later.

Besides its fundamental interest, the multi-partite entanglement that we have demonstrated may be useful for quantum metrology^[@CR23]^, provided that the storage efficiency can be increased substantially, which is possible, e.g., using cavities^[@CR35]^. The AFC system allows one to address individual or sets of teeth in frequency space, which in principle makes it possible to characterize the created quantum state in more detail. The large ratio of inhomogeneous to homogeneous linewidths in rare-earth-doped crystals suggests that it should be possible to create entanglement between larger numbers of teeth, possibly up to one hundred million. The length and the doping concentration of rare earth materials can also be increased, allowing one to additionally increase the number of atoms in each tooth. Storing entangled photons in AFCs created in separated crystals would offer the possibility to study the nature of multi-partite entanglement at yet a higher nesting level. Furthermore the present approach in principle allows the storage of more than one photon, which would enable the creation of more complex entangled states such as higher-order Dicke states.

We note that entanglement between many individual atoms in a solid is demonstrated in a parallel submission by F. Fröwis et al., using a related, but complementary approach based on the directionality of the echo emission from an atomic frequency comb.

Methods {#Sec6}
=======

Only the Dicke state has a non-zero eigenvalue {#Sec7}
----------------------------------------------

Here we prove the statement that in the single-excitation subspace of the *M*-teeth Hilbert space the Dicke state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M}$$\end{document}$ is the only state with a non-zero eigenvalue. For this purpose, it is useful to view each qubit as a spin-1/2 system, such that the dipole operators $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| 1 \right\rangle \left\langle 0 \right|$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| 0 \right\rangle \left\langle 1 \right|$$\end{document}$ are now viewed as spin raising and lowering operators, respectively. In this spin representation, the zero-excitation state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {0...0} \right\rangle$$\end{document}$ has total spin *M*/2 and *S* ~*z*~ projection −*M*/2, and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M}$$\end{document}$ (which is also completely symmetric and thus belongs to the irreducible representation with maximum total spin) has total spin *M*/2 as well, but *S* ~*z*~ projection −*M*/2 + 1. In addition to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M},$$\end{document}$ there are *M* − 1 other basis states spanning the single-excitation subspace. They all have the same *S* ~*z*~ projection as the Dicke state (i.e., *S* ~z~ = −*M*/2 + 1), but total spin *M*/2 − 1 since they are not fully symmetric. They have the lowest possible *S* ~*z*~ projection value that is compatible with their total spin value, and for that reason, they are annihilated by the spin lowering operator *S* ~−~. Hence, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M}$$\end{document}$ is the only eigenstate of *S* ~+~ *S* ~−~ in the single-excitation subspace that has a non-zero eigenvalue.

Deriving bounds on the entanglement depth {#Sec8}
-----------------------------------------

In order to derive lower bounds on the entanglement depth *M* for given *R*, *P* ~1~ and *P* ~2~, it is convenient to---equivalently---derive upper bounds on *R* for given *M*, *P* ~1~ and *P* ~2~. The fully separable example discussed in the text shows that *R* can be much greater than the entanglement depth *M* once *P* ~2~ is different from zero. For a given *M*, we are allowed to use entangled states of (up to) *M*-teeth as the individual factors, where the factors *α* $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| 0 \right\rangle$$\end{document}$ + *β* $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| 1 \right\rangle$$\end{document}$ in the fully separable example correspond to *M* = 1. Keeping in mind that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P_1},{P_2} \ll 1$$\end{document}$ in our experiment, it is clear that there should be a significant vacuum component $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| 0 \right\rangle$$\end{document}$ ^⊗*M*^ in each factor. As for the non-vacuum part, the Dicke state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M}$$\end{document}$ is clearly a good choice because it is the single-excitation state that maximizes $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$R = \left\langle {{S_ + }{S_ - }} \right\rangle {\rm{/}}\left\langle {\mathop {\sum}\nolimits_j {{\left| 1 \right\rangle }^j}{{\left\langle 1 \right|}^j}} \right\rangle$$\end{document}$ in each *M*-teeth subspace. In fact, it is the optimum choice. Higher-order Dicke states have greater $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left\langle {{S_ + }{S_ - }} \right\rangle$$\end{document}$, but not *R*, because they also have greater values for the total number of excitations $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left\langle {\mathop {\sum}\nolimits_j {{\left| 1 \right\rangle }^j}{{\left\langle 1 \right|}^j}} \right\rangle$$\end{document}$. Moreover---in contrast to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _M}$$\end{document}$---they require a non-zero *P* ~2~ (alternatively stated, they make a non-zero contribution to *P* ~2~), whose value is one of our constraints. Higher-order non-Dicke states are clearly suboptimal. We conclude that the optimal state should contain factors of the form $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\alpha {\left| 0 \right\rangle ^{ \otimes M}} + \beta {\left| W \right\rangle _M}.$$\end{document}$ However, it is not immediately apparent how many such factors there should be. For example, consider the states $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _1}} \right\rangle = \left( {{\alpha _1}{{\left| 0 \right\rangle }^{ \otimes M}} + {\beta _1}{{\left| W \right\rangle }_M}} \right) \otimes {\left| 0 \right\rangle ^{ \otimes N - M}}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _2}} \right\rangle = {\left( {{\alpha _2}{{\left| 0 \right\rangle }^{ \otimes M}} + {\beta _2}{{\left| W \right\rangle }_M}} \right)^{ \otimes 2}} \otimes {\left| 0 \right\rangle ^{ \otimes N - 2M}}.$$\end{document}$ The latter state has a component $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _{2M}}$$\end{document}$ in the single-excitation subspace and will thus achieve a greater value of *R* compared to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _1}} \right\rangle$$\end{document}$. But it also makes a contribution to *P* ~2~ because of the cross term that is proportional to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\beta _2^2$$\end{document}$, whereas $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _1}} \right\rangle$$\end{document}$ makes no such contribution. This suggests that in general, the optimal state may be a mixed state of the form:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\rho = \mathop {\sum}\limits_i {q_i}\left| {{\psi _i}} \right\rangle \left\langle {{\psi _i}} \right|,$$\end{document}$$with $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\sum} {q_i} = 1$$\end{document}$, *q* ~*i*~ ≥ 0, and$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _i}} \right\rangle = {\left( {{\alpha _i}{{\left| 0 \right\rangle }^{ \otimes M}} + {\beta _i}{{\left| W \right\rangle }_M}} \right)^{ \otimes i}} \otimes {\left| 0 \right\rangle ^{ \otimes (N - iM)}},$$\end{document}$$where 1 ≤ *i* ≤ *k*, such that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$k = \left\lfloor {N{\rm{/}}M} \right\rfloor .$$\end{document}$ Let us first consider the case where *N* is divisible by *M*, in which case one simply has *k* = *N*/*M*. Note that it is optimal for the coefficients of the (non-vacuum) factors in each $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _i}} \right\rangle$$\end{document}$ to be identical, in order to maximize the overlap with the state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _{iM}}$$\end{document}$ in the single-excitation subspace and thus maximize *R*. The ratio *R* and the probabilities *P* ~1~ and *P* ~2~ for the state *ρ* are given by:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$R = \frac{M}{{{P_1} + 2{P_2}}}\mathop {\sum}\limits_{n = 1}^k {n^2}{q_n}\left( {\beta _n^2} \right){\left( {1 - \beta _n^2} \right)^{n - 1}},$$\end{document}$$ $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P_1} = {q_1}\beta _1^2 + \mathop {\sum}\limits_{n = 2}^k n{q_n}\left( {\beta _n^2} \right){\left( {1 - \beta _n^2} \right)^{n - 1}},$$\end{document}$$ $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P_2} = \mathop {\sum}\limits_{n = 2}^k \frac{{n(n - 1)}}{2}{q_n}\left( {\beta _n^4} \right){\left( {1 - \beta _n^2} \right)^{n - 2}}.$$\end{document}$$

To derive the upper bounds on *R*, we maximize Eq. ([6](#Equ6){ref-type=""}) with the constraints of Eqs. ([7](#Equ7){ref-type=""}) and ([8](#Equ8){ref-type=""}), where the values of *P* ~1~ and *P* ~2~ are obtained from our experiment. We use a global numerical search algorithm that runs a local nonlinear programming solver, which finds the maximum value of a constrained nonlinear multivariable function for multiple start points. It finally reports the global maximum value of *R* for a given *M*.

The numerical maximization shows that only the *q* ~1~ and *q* ~*k*~ components in Eq. ([4](#Equ4){ref-type=""}) are non-zero in the optimal case, and the latter is close to one. This can be understood physically. The state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _k}} \right\rangle$$\end{document}$ has the largest possible number of non-separable states of size *M*, giving a state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _N}$$\end{document}$ in the single-excitation subspace, which maximizes $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left\langle {{S_ + }{S_ - }} \right\rangle$$\end{document}$. Thus, the weight of this state, *q* ~*k*~, should be as large as possible, and indeed it comes out close to 1 in the maximization. However, the size of *β* ~*k*~ is constrained by the value of *P* ~2~, which is very small in our case. As a consequence, this state only makes a small contribution to *P* ~1~. The remaining contribution to *P* ~1~ is best provided by a state that does not increase *P* ~2~, i.e., $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _1}} \right\rangle$$\end{document}$.

In general, *N* is not always divisible by *M* and one has *N* = *kM* + *k*′, such that *k*′ \< *M* (while $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$k = \left\lfloor {N{\rm{/}}M} \right\rfloor$$\end{document}$). In this case, the same arguments apply with the small difference that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {{\psi _k}} \right\rangle = {\left( {{\alpha _k}{{\left| 0 \right\rangle }^{ \otimes M}} + {\beta _k}{{\left| W \right\rangle }_M}} \right)^{ \otimes k}} \otimes \left( {{\alpha _{k'}}{{\left| 0 \right\rangle }^{ \otimes k'}} + {\beta _{k'}}{{\left| W \right\rangle }_{k'}}} \right).$$\end{document}$ The optimum values for *α* ~*k*′~ and *β* ~*k*′~ are such that the state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| W \right\rangle _N}$$\end{document}$ is obtained in the single-excitation subspace.

Figure [5](#Fig5){ref-type="fig"} shows the maximum possible contrast, *R*, as a function of entanglement depth, *M*, constrained on our experimental values of *P* ~1~ and *P* ~2~. In this figure, the curve is close to a straight line. This can be explained by noting that only *q* ~1~ and *q* ~*k*~ are non-zero $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left( {{q_1} + {q_k} \simeq 1} \right),$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P_2} \simeq \frac{1}{2}{k^2}{q_k}\beta _k^4$$\end{document}$ is very small, which makes the contribution of *q* ~*k*~ to *P* ~1~ negligible ($\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$k\beta _k^2 \ll {P_1}$$\end{document}$ for any value of *q* ~*k*~) resulting in $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P_1} \simeq {q_1}\beta _1^2.$$\end{document}$ Thus, one can obtain $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$R = {\textstyle{1 \over {{P_1} + 2{P_2}}}}\left( {M{P_1} + \sqrt {2\left( {1 - {q_1}} \right){P_2}} N} \right),$$\end{document}$ which gives $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${R_{{\rm{max}}}} \simeq M + {\textstyle{{\sqrt {2{P_2}} } \over {{P_1}}}}N$$\end{document}$ and thereby leads to Eq. ([3](#Equ3){ref-type=""}).Fig. 5Bound for maximum value of echo contrast *R* as a function of entanglement depth *M*. We assume *P* ~1~ = 3.5 × 10^−3^ and *P* ~2~ = 2.6 × 10^−8^ for the single-excitation and double-excitation probability, respectively, and the total number of ensembles (teeth) *N* = 564. The experimental value of the echo contrast *R* is shown as a *green square*

Theoretical model {#Sec9}
-----------------

The atomic state of the AFC system after absorption depends on the photon statistics of the source, on transmission loss, and on inefficiencies such as inefficient single-photon heralding and absorption in the AFC. We now describe a theoretical model that allows extracting *P* ~1~ and *P* ~2~ from experimental data.

In our experiment, an intense laser beam passes through a spontaneous parametric down conversion crystal (SPDC) that converts photons into entangled pairs of photons traveling in different spatial modes a and b (see Fig. [6](#Fig6){ref-type="fig"}). A non-number-resolving photon detector D~a~, henceforth referred to as the heralding detector, is placed in mode a, and the AFC system followed by another detector D~b~ is located in mode b. A detection by the detector D~a~ heralds the presence of at least one photon in mode b.Fig. 6Visualization of our mathematical model for the set-up. This figure visualizes our mathematical model of the set-up of Fig. [1c](#Fig1){ref-type="fig"}; *ξ* ~*n*~ is the amplitude for creating *n* pairs, *η* ~a~ is the overall transmission and detection probability in mode a, D~a~ is the heralding detector, *η* ~b~ is the transmission before the memory in mode b, *η* ~w~ is the write efficiency of the atomic frequency comb memory, *η* ~t~ is the transmission and detection probability after the memory in mode b, and D~b~ is the detector in mode b. All losses and inefficiencies are modeled as beam splitters

The combined state of the two downconverted photons in modes a and b can be written as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| \psi \right\rangle = \mathop {\sum}\limits_{n = 0}^\infty \frac{{{\xi _n}}}{{n!}}{\left( {{a^\dag }} \right)^n}{\left( {{b^\dag }} \right)^n}{\left| 0 \right\rangle _{\rm{a}}}{\left| 0 \right\rangle _{\rm{b}}},$$\end{document}$$where$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| {{\xi _n}} \right\rangle ^2} = \frac{{{\mu ^n}}}{{{{(\mu + 1)}^{n + 1}}}}$$\end{document}$$is the thermal photon number distribution. Here *a* ^†^ (*b* ^†^) is the photon creation operator in mode a (b), $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| {{\xi _n}} \right|^2}$$\end{document}$ is the probability of creating *n* pairs of photons (one photon per pair in mode a and one in b), and *μ* is the mean photon pair number per temporal mode (per coherence time of the pump laser). Let us note that our source is in fact slightly temporally multi mode, which implies that the distribution should be somewhere in-between thermal and Poissonian. However, here we conservatively assume thermal statistics, which leads to slightly larger values of *P* ~2~ and thus slightly lower values of the entanglement depth extracted from the measured data (Fig. [3](#Fig3){ref-type="fig"}).

Due to loss, a photon created in mode a will reach the detector D~a~ with probability of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{a}}^*}}} < 1.$$\end{document}$ It will be detected with a probability given by the detection efficiency $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{D}}_{\rm{a}}}}}$$\end{document}$. We model the combined (limited) channel transmission and detector efficiency using a beam-splitter with transmission probability of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{\rm{a}}} = {\eta _{{{\rm{a}}^*}}} \times {\eta _{{{\rm{D}}_{\rm{a}}}}}$$\end{document}$ followed by a perfect detector. Thus, the the creation operator for mode a transforms as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${a^\dag } \to \sqrt {{\eta _{\rm{a}}}} {a^\dag } + \sqrt {1 - {\eta _{\rm{a}}}} a_{\rm{l}}^\dag ,$$\end{document}$$where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$a_{\rm{l}}^\dag$$\end{document}$ is the creation operator in the loss mode of this hypothetical beam-splitter. With this transformation we can write:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left( {{a^\dag }} \right)^n} = \mathop {\sum}\limits_{k = 0}^n \left( {\begin{array}{*{20}{c}} n \\ k \\ \end{array}} \right){\left( {\sqrt {{\eta _{\rm{a}}}} {a^\dag }} \right)^k}{\left( {\sqrt {1 - {\eta _{\rm{a}}}} a_{\rm{l}}^\dag } \right)^{n - k}}.$$\end{document}$$The state of the loss mode a~l~ and the mode b together, after a *k*-photon detection in the detector D~a~, becomes:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| {\bar \psi } \right\rangle _k} = \mathop {\sum}\limits_{n = k}^\infty \sqrt {\left( {\begin{array}{*{20}{c}} n \\ k \\ \end{array}} \right)} {\xi _n}{\sqrt {{\eta _{\rm{a}}}} k}{\sqrt {(1 - {\eta _{\rm{a}}})^{\left( {n - k} \right)}}}{\left| {n - k} \right\rangle _{{{\rm{a}}_{\rm{l}}}}}{\left| n \right\rangle _{\rm{b}}}.$$\end{document}$$Considering the fact that D~a~ is a detector that does not resolve photon numbers, we sum over all possible values of *k* ≥ 1 in mode a and trace out the loss channel from the state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\left| {\bar \psi } \right\rangle _k}$$\end{document}$. This results in the state of the mode b after heralding by D~a~ as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{*{20}{l}}\\ {{\rho _{\rm{b}}}} & = {\frac{1}{{\cal N}}\mathop {\sum}\limits_{k = 1}^\infty \mathop {\sum}\limits_{n = k}^\infty \left( {\begin{array}{*{20}{c}} n \\ k \\ \end{array}} \right){{\left| {{\xi _n}} \right|}^2}\eta _{\rm{a}}^k{{\left( {1 - {\eta _{\rm{a}}}} \right)}^{(n - k)}}{{\left| n \right\rangle }_{\rm{b}}}\left\langle n \right|} \hfill \\ \\ {} \hfill & = {\mathop {\sum}\limits_{n = 1}^\infty {p_n}{{\left| n \right\rangle }_{\rm{b}}}\left\langle n \right|,} \hfill \\ \end{array}$$\end{document}$$with$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\cal N} = \mathop {\sum}\limits_{k = 1}^\infty \mathop {\sum}\limits_{n = k}^\infty \left( {\begin{array}{*{20}{c}} n \\ k \\ \end{array}} \right){\left| {{\xi _n}} \right|^2}\eta _{\rm{a}}^k{\left( {1 - {\eta _{\rm{a}}}} \right)^{(n - k)}},$$\end{document}$$ $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${p_n} = \frac{{\mathop {\sum}\limits_{k = 1}^n \left( {\begin{array}{*{20}{c}} n \\ k \\ \end{array}} \right){{\left| {{\xi _n}} \right|}^2}\eta _{\rm{a}}^k{{\left( {1 - {\eta _{\rm{a}}}} \right)}^{(n - k)}}}}{{\cal N}}.$$\end{document}$$

In mode b, each photon experiences loss (before entering the AFC) that is modeled by a beam-splitter with transmission probability *η* ~b~. Each photon is then absorbed by the AFC with probability *η* ~w~. Otherwise, it is transmitted through the AFC and detected by the detector D~b~, similarly modeled as a beam-splitter with transmission probability of *η* ~t~ followed by a perfect detector. Thus, the creation operator for mode b transforms as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{*{20}{l}}\\ {{b^\dag }} & \to {\sqrt {{\eta _{\rm{b}}}{\kern 1pt} {\eta _{\rm{w}}}} {b^\dag } + \sqrt {{\eta _{\rm{b}}}{\kern 1pt} \left( {1 - {\eta _{\rm{w}}}} \right){\kern 1pt} {\eta _{\rm{t}}}} b_{\rm{t}}^\dag } \hfill \\ \\ {} & {} { + \sqrt {{\eta _{\rm{b}}}{\kern 1pt} \left( {1 - {\eta _{\rm{w}}}} \right)\left( {1 - {\eta _{\rm{t}}}} \right)} b_{{\rm{tl}}}^\dag + \sqrt {1 - {\eta _{\rm{b}}}} b_{\rm{l}}^\dag ,} \hfill \\ \end{array}$$\end{document}$$where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$b_{\rm{t}}^\dag ,$$\end{document}$ $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$b_{\rm{l}}^\dag ,$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$b_{{\rm{tl}}}^\dag$$\end{document}$ are the creation operators in the transmission mode, loss mode before the AFC, and the loss mode after the AFC, respectively. Hence, the associated optical state becomes:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rho _{\rm{b}}} = \mathop {\sum}\limits_n {p_n}\left| {{\Phi _n}} \right\rangle \left\langle {{\Phi _n}} \right|,$$\end{document}$$where$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{*{20}{l}}\\ {\left| {{\Phi _n}} \right\rangle } & = {\frac{1}{{\sqrt {n!} }}\left( {\sqrt {{\eta _{\rm{b}}}{\kern 1pt} {\eta _{\rm{w}}}} {b^\dag } + \sqrt {{\eta _{\rm{b}}}\left( {1 - {\eta _{\rm{w}}}} \right){\kern 1pt} {\eta _{\rm{t}}}} {\kern 1pt} b_{\rm{t}}^\dag } \right.} \hfill \\ \\ {} & {} {{{\left. { + \sqrt {{\eta _{\rm{b}}}\left( {1 - {\eta _{\rm{w}}}} \right)\left( {1 - {\eta _{\rm{t}}}} \right)} b_{{\rm{tl}}}^\dag + \sqrt {1 - {\eta _{\rm{b}}}} {\kern 1pt} b_{\rm{l}}^\dag } \right)}^n}\left| 0 \right\rangle .} \hfill \\ \end{array}$$\end{document}$$Conditioning on "no detection" in the *b* ~t~ mode (since a detection would correspond to a photon that was not absorbed in the AFC) results in the density matrix$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\tilde \rho _{\rm{b}}} = \frac{1}{{\cal M}}\mathop {\sum}\limits_n {\bar p_n}\left| {{{\tilde \Phi }_n}} \right\rangle \left\langle {{{\tilde \Phi }_n}} \right|,$$\end{document}$$where$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{l}\\ \left| {{{\tilde \Phi }_n}} \right\rangle = \frac{1}{{{{\left( {{\eta _{\rm{b}}}{\kern 1pt} {\eta _{\rm{w}}} + {Z^2}} \right)}^{n/2}}}}\left| {{{\bar \Phi }_n}} \right\rangle ,\\ \\ \left| {{{\bar \Phi }_n}} \right\rangle = \frac{1}{{\sqrt {n!} }}{\left( {\sqrt {{\eta _{\rm{b}}}{\kern 1pt} {\eta _{\rm{w}}}} {b^\dag } + Z{x^\dag }} \right)^n}\left| 0 \right\rangle ,\\ \\ {{\bar p}_n} = {p_n}{\left( {{\eta _{\rm{b}}}{\kern 1pt} {\eta _{\rm{w}}} + {Z^2}} \right)^n},\\ \\ {\cal M} = \mathop {\sum}\limits_n {p_n}{\left( {{\eta _{\rm{b}}}{\kern 1pt} {\eta _{\rm{w}}} + {Z^2}} \right)^n}.\\ \end{array}$$\end{document}$$We have used$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{l}\\ Z = \sqrt {{\eta _{\rm{b}}}\left( {1 - {\eta _{\rm{w}}}} \right)\left( {1 - {\eta _{\rm{t}}}} \right) + \left( {1 - {\eta _{\rm{b}}}} \right)} ,\\ \\ {x^\dag } = \frac{1}{Z}\left( {\sqrt {{\eta _{\rm{b}}}{\kern 1pt} \left( {1 - {\eta _{\rm{w}}}} \right)\left( {1 - {\eta _{\rm{t}}}} \right)} b_{{\rm{tl}}}^\dag + \sqrt {1 - {\eta _{\rm{b}}}} b_{\rm{l}}^\dag } \right).\\ \end{array}$$\end{document}$$We can rewrite the state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\tilde \rho _b}$$\end{document}$ in the number basis of the b mode and the x mode (defined in the previous equation) as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{c}\\ {{\tilde \rho }_{\rm{b}}} = \frac{1}{{\cal M}}\left\{ {\mathop {\sum}\limits_{n = 1}^\infty {p_n}\mathop {\sum}\limits_{r,r' = 1}^n \sqrt {\left( {\begin{array}{*{20}{c}} n \\ r \\ \end{array}} \right)\left( {\begin{array}{*{20}{c}} n \\ {r'} \\ \end{array}} \right)} {{\sqrt {{\eta _{\rm{b}}}{\eta _{\rm{w}}}} }^{r + r'}}{Z^{2n - r - r'}}} \right.\\ \\ \left. {{{\left| {r,n - r} \right\rangle }_{{\rm{b,x}}}}\left\langle {r',n - r'} \right|} \right\}.\\ \end{array}$$\end{document}$$

The probability *P* ~*r*~ that *r* photons are present in the AFC is calculated as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{*{20}{l}}\\ {{P_r}} = {_{\rm{b}}\left\langle r \right|{\kern 1pt} {\rm{Tr}}{{\kern 1pt} _x}{{\tilde \rho }_{\rm{b}}}{{\left| r \right\rangle }_{\rm{b}}}=}\ {\frac{1}{{\cal M}}\mathop {\sum}\limits_{n = r}^\infty {p_n}\left( {\begin{array}{*{20}{c}} n \\ r \\ \end{array}} \right){{\left( {{\eta _{\rm{b}}}{\eta _{\rm{w}}}} \right)}^r}{Z^{2n - 2r}}} \hfill \\ \end{array}$$\end{document}$$Therefore, we can obtain the values of *P* ~1~ and *P* ~2~ by measuring *μ*, *η* ~a~, *η* ~b~, *η* ~w~, and *η* ~t~.

Experimental estimation of *μ* {#Sec10}
------------------------------

The second-order cross-correlation function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g(\tau )_{{\rm{ab}}}^{(2)}$$\end{document}$ gives information about the photon number distribution of a two-mode field. Specifically, for the downconverted photons in mode a (1532 nm) and mode b (795 nm) at zero time delay (*τ* = 0), it can be written as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{ab}}}^{(2)}(0) = \frac{{{p_{{\rm{ab}}}}(0)}}{{{p_{\rm{a}}}{p_{\rm{b}}}}},$$\end{document}$$where *p* ~a~ (*p* ~b~) is the probability for a detection in mode a (b) and *p* ~ab~(0) the probability to detect a coincidence in a temporal window centered at *τ* = 0. Experimentally, we obtained $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{ab}}}^{(2)}$$\end{document}$ as the ratio between the coincidence rate *C* ~ab~ at *τ* = 0 (where photons in mode a and *b* exhibit maximum correlations) and the coincidence rate at a delay *τ* larger than the coherence length of the photons in mode a and b (where photon creation, and hence detection, in mode a and b is completely independent).

All experiments were performed at maximum pump power of ≈500 μW, for which we measured $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{ab}}}^{(2)}(0) = 884 \pm 50.$$\end{document}$ In the limit where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${g_{{\rm{ab}}}} \gg 1$$\end{document}$, the cross-correlation function can be written as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${g_{{\rm{ab}}}}(0) = \frac{1}{\mu }$$\end{document}$ ^[@CR44]^. This allows us to estimate a value of *μ* = (1.1 ± 0.1) × 10^−3^. Note that the uncertainty in *μ* is the dominant contribution to the uncertainty of *P* ~2~ given in the manuscript. Our coincidence measurements were performed using home-made logic electronics with a detection window of 5 ns. A characterization of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{ab}}}^{(2)}$$\end{document}$ as a function of the pump power can be seen in Fig. [2](#Fig2){ref-type="fig"}.

Experimental estimation of *η*~a~, *η*~b~, and *η*~t~ {#Sec11}
-----------------------------------------------------

To estimate *η* ~a~ and *η* ~b~, we change the set-up represented in Fig. [6](#Fig6){ref-type="fig"} and place the detector D~b~ (with detection efficiency of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{D}}_{\rm{b}}}}}$$\end{document}$) before the AFC, i.e., the photons in mode b only pass through the loss channel with probability $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{b}}^*}}}$$\end{document}$ before they reach D~b~. In this configuration, both *η* ~a~ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{b}}^*}}}$$\end{document}$ can be estimated from the coincidence detection rate *C* ~ab~ and the single detection rates *S* ~a~ and *S* ~b~. Since $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\mu \ll 1$$\end{document}$ in our case, we can write:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{l}\\ {C_{{\rm{ab}}}} = \mu {\eta _{{{\rm{a}}^*}}}{\eta _{{{\rm{D}}_{\rm{a}}}}}{\eta _{{{\rm{b}}^*}}}{\eta _{{{\rm{D}}_{\rm{b}}}}}{\rm{/}}{\tau _{\rm{p}}},\\ {S_{\rm{a}}} = \mu {\eta _{{{\rm{a}}^*}}}{\eta _{{{\rm{D}}_{\rm{a}}}}}{\rm{/}}{\tau _{\rm{p}}},\\ {S_{\rm{b}}} = \mu {\eta _{{{\rm{b}}^*}}}{\eta _{{{\rm{D}}_{\rm{b}}}}}{\rm{/}}{\tau _{\rm{p}}},\\ \end{array}$$\end{document}$$where *τ* ~p~ is the coherence time of the pump laser. From the above relations, we find:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{*{20}{l}}\\ {{\eta _{\rm{a}}} = \frac{{{C_{{\rm{ab}}}}}}{{{S_{\rm{b}}}}},} \hfill \\ \\ {{\eta _{{{\rm{b}}^*}}} = \frac{{{C_{{\rm{ab}}}}}}{{{S_{\rm{a}}}{\eta _{{{\rm{D}}_{\rm{b}}}}}}}.} \hfill \\ \end{array}$$\end{document}$$

We experimentally find *η* ~a~ = 11.0% and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{b}}^*}}} = 5.3\%$$\end{document}$, for which we used $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{D}}_{\rm{b}}}}} = 0.60$$\end{document}$. In the set-up depicted in Fig. [6](#Fig6){ref-type="fig"}, we can write $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{\rm{b}}} = {\eta _{{{\rm{b}}^*}}} \times {\eta _{{{\rm{c}}_{\rm{i}}}}}$$\end{document}$, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{c}}_{\rm{i}}}}}$$\end{document}$ is the probability with which a photon that passes through the loss will enter the AFC. Note that the background absorption of the AFC *d* ~0~, as shown in Fig. [7](#Fig7){ref-type="fig"}, can be treated as loss and is hence included in $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{c}}_{\rm{i}}}}}$$\end{document}$. Since *d* ~0~ varies for AFCs with different *N*, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{c}}_{\rm{i}}}}}$$\end{document}$ and hence *η* ~b~ also varies. For instance, the AFC with *N* = 564 and *B* = 6 GHz for which entanglement depth is calculated in Fig. [3b](#Fig3){ref-type="fig"}, results in *η* ~b~ = 1.06%. For the two AFCs shown in Fig. [4](#Fig4){ref-type="fig"} with *N* = 9 and bandwidths of *B* = 6 GHz and *B* = 0.33 GHz, we estimate *η* ~b~ = 1.99% and *η* ~b~ = 0.96%, respectively.Fig. 7A sample trace of an atomic frequency comb. This comb has *N* = 9 teeth and and a bandwidth *B* = 6 GHz. The sample is derived from an average over traces of 50 experimental runs; *d* ~1~ indicates the peak-to-peak optical depth that constitutes the AFC while *d* ~0~ is the background optical depth. The finesse *F* of the AFC is given by the ratio *Δ*/*γ*, where *Δ* is the separation between the teeth and *γ* is the linewidth of the teeth

The probability with which a photon that is transmitted through the AFC is detected by the detector D~b~, can be written as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _t} = {\eta _{{t^*}}} \times {\eta _{{{\rm{D}}_{\rm{b}}}}}.$$\end{document}$ Here, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{t^*}}}$$\end{document}$ is the probability that the photon passes through the loss channel after the AFC, reaches into the detector D~b~ which has a detection efficiency of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{{{\rm{D}}_{\rm{b}}}}}.$$\end{document}$ We estimate *η* ~t~ = 36.0%.

Experimental estimation of AFC write efficiency *η*~w~ {#Sec12}
------------------------------------------------------

The write efficiency of the AFC is the probability for the AFC to absorb a photon. It can be calculated from the effective optical depth of the AFC as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\eta _{\rm{w}}} = 1 - {{\rm{e}}^{\frac{{ - {d_1}}}{F}}},$$\end{document}$$where *d* ~1~ is the peak absorption, *F* is the finesse of the AFC, and *d* ~1~/*F* is the effective optical depth of the comb (see Fig. [7](#Fig7){ref-type="fig"}). For the different AFCs that we create, *d* ~1~ varies, and hence the write efficiency *η* ~*w*~ also varies. For instance, the write efficiencies for the broadband AFCs with number of teeth *N* = 564 and *N* = 9 are 33% and 54%, respectively. For the narrowband AFC with *N* = 9, which is shown in Fig. [4](#Fig4){ref-type="fig"}, we estimate *η* ~w~ = 0.9%.

The dominant contribution to the uncertainty in *P* ~1~ given in the manuscript comes from the uncertainty in the optical depth of around 10%. The uncertainties associated with the measurements of the other efficiencies are negligible in comparison.

Heralded autocorrelation $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${g}_{{\bf{bb}}}^{{\bf{(2)}}}{\bf{(0)}}$$\end{document}$ {#Sec13}
-----------------------------------------------------------------------------------------

The second-order zero-time autocorrelation function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{bb}}}^{(2)}(0)$$\end{document}$ is a witness of non-classicality for single-mode fields. To measure $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{bb}}}^{(2)}(0)$$\end{document}$, we add a balanced (50/50) beam-splitter (BS) before the detection of the photon in mode b (795 nm photon), and measure the probability to detect a coincidence between the two BS outputs $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left( {{p_{{{\rm{b}}_{\rm{1}}}{{\rm{b}}_{\rm{2}}}}}} \right)$$\end{document}$, as well as the single output probabilities ($\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${p_{{{\rm{b}}_{\rm{1}}}}}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${p_{{{\rm{b}}_{\rm{2}}}}}$$\end{document}$), all conditioned on the detection of a photon in mode a (1532 nm photon). We can now write the autocorrelation function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{bb}}}^{(2)}$$\end{document}$ as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{bb}}}^{(2)}(0) = \frac{{{p_{{{\rm{b}}_{\rm{1}}}{{\rm{b}}_{\rm{2}}}}}(0)}}{{{p_{{{\rm{b}}_{\rm{1}}}}}{p_{{{\rm{b}}_{\rm{2}}}}}}}.$$\end{document}$$

For a perfect single-photon state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{bb}}}^{(2)}(0) = 0$$\end{document}$, whereas for a coherent state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{bb}}}^{(2)}(0) = 1$$\end{document}$. We find $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$g_{{\rm{bb}}}^{(2)}(0) = 0.0024 \pm 0.0006$$\end{document}$, see Fig. [2b](#Fig2){ref-type="fig"}, proving the nearly ideal nature of our heralded single-photon source. This value is consistent with the value of *μ* determined from measuring the cross-correlation function above^[@CR45]^.

Background noise estimation and subtraction {#Sec14}
-------------------------------------------

The collection of our experimental data suffers from imperfections that limit the the values of *R*, but are unrelated to the detection of the heralded single photon. In our analysis, we account for two such sources of imperfection, namely background noise and limited temporal resolution of the detector, where the latter is treated in the next section. The background noise causes a constant level of detection events, which we can assess by averaging over the detection events before a heralding photon detection as indicated by the *red region* in the *insert* of Fig. [3](#Fig3){ref-type="fig"}. The background noise is constant at about 0.9 counts per 80 ps data bin for a 15 min measurement. Since this background noise is completely uncorrelated with the heralded photon, we can subtract it in order to obtain a more accurate echo contrast. Clearly, the impact of such noise subtraction will increase when the echo signal is weakened, e.g., due to increased AFC storage times or reduction of the AFC bandwidth, as evidenced by the data shown in Figs. [3](#Fig3){ref-type="fig"}, [4](#Fig4){ref-type="fig"}, and [8](#Fig8){ref-type="fig"}.Fig. 8Echo contrast depends on atomic frequency comb bandwidth. This figure shows the echo contrast *R* as a percentage of the total number of teeth as a function of the bandwidth of the atomic frequency comb. *Error bars* indicate standard deviations derived from the Gaussian fitting and Poissonian statistics of the photon detection events

We can determine the origin of the background noise from a series of independent measurements. Detector dark counts contribute about 12%, while noise due to leaked optical pumping light and spontaneous emission from atoms excited during optical pumping amounts to \~11%. Finally, the photon pair source may generate additional pairs other than that causing the heralding event, as discussed above in the context of the second-order cross-correlation function. If the 795 nm photon member of one of these additional pairs reaches the detector after the AFC, it will appear as a background noise count. This is the main contribution to the background noise at around 77%. The relative contributions to the background noise depend slightly on the AFC efficiency and loss.

Deconvolution of the detector response {#Sec15}
--------------------------------------

The temporal shape of the heralded photon is given by the Fourier transform of its spectral distribution. At the AFC input, the spectrum is defined by the 6 GHz FP filter, whereas after re-emission from the AFC the photon spectrum will depend on the overall bandwidth as well as exact shape of the AFC. However, when the heralded photon is detected, the inherent detector jitter can smear out its temporal shape thus reducing the signal-to-noise ratio and thereby the echo contrast.

The joint detector response is measured by removing all spectral filtering elements and recording the distribution of correlated detection immediately after the, now extremely broadband, pair source. We find that the joint detector response is well approximated by a Gaussian function of full-width-at-half-maximum (FWHM) Δ*t* ~D~ = 354 ps. Fitting the echoes with a Gaussian function with amplitude *E* and a FWHM of Δ*t* ~p~ allows us to deconvolute the detector response and compute a corrected signal amplitude of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$E' = E\sqrt {\Delta {\kern 1pt} t_{\rm{p}}^2{\rm{/}}\left( {\Delta {\kern 1pt} t_{\rm{p}}^2 - \Delta {\kern 1pt} t_{\rm{D}}^2} \right)}$$\end{document}$. Note that for the broadband AFCs, the output photon duration is on the order of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\Delta {\kern 1pt} t{\prime _{\rm{p}}}\sim100 - 200$$\end{document}$ ps, which is significantly smeared by the detector jitter. As expected, the deconvolution of the detector response leads to a large increase of the echo contrast for the broadband AFCs and less of an increase for narrowband AFCs, for which the temporal duration of the echoes exceeds the detector jitter---this interplay is evident in Figs. [4](#Fig4){ref-type="fig"} and [8](#Fig8){ref-type="fig"}.

Number of atoms corresponding to a single AFC tooth {#Sec16}
---------------------------------------------------

We calculate the number of atoms that correspond to a single AFC tooth using two complementary approaches. The first method utilizes the experimentally determined absorption spectrum and the Tm-atom density of the Tm:LiNbO~3~ crystal, while the second method relies on single-ion spectroscopic properties.

For the first method, we note that the integrated absorption spectrum Θ of an inhomogeneously broadened transition of an atomic ensemble is $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\Theta _{ {\rm{i}}}} = {\int} {\alpha (\nu )L{\rm{d}}\nu }$$\end{document}$, where *α*(*ν*) is the absorption coefficient resulting from all transitions that feature a resonance frequency *ν*, and *L* is the length of the medium. Similarly, the integrated absorption spectrum of a single AFC tooth is *Θ* ~t~, where the integration is taken over the tooth. If a laser beam of cross-sectional area *A* (defined by FWHM of the intensity distribution of a Gaussian beam) is sent through a medium that features atom density *n* ~d~, then the number of atoms within the beam is *n* ~d~ *LA*, and the number of atoms corresponding to a single tooth is $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$N_{\rm{t}}^{(1)} = {n_{\rm{d}}}LA\left( {{\Theta _{ {\rm{t}}}}{\rm{/}}{\Theta _{ {\rm{i}}}}} \right)$$\end{document}$.

For the second method, we calculate the optical depth *d* ~atom~ that corresponds to a single atom in a crystal using *d* ~atom~ = \[(*n* ^2^ + 2)^2^/(72*πnAσ* ^2^)\](*γ* ~s~/*Γ* ~h~), where *γ* ~s~ and *Γ* ~h~ are the spontaneous emission rate and homogeneous broadening of the transition, respectively, *n* is the index of refraction, and *σ* = *ν*/*c* ^[@CR46]^. Since *d* ~atom~ refers to an atom that features a linewidth of *Γ* ~h~, we estimate the number of atoms that correspond to a single tooth using $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$N_{\rm{t}}^{(2)} = {\Theta _{ {\rm{t}}}}{\rm{/}}\left( {{{\it{\Gamma }}_{\rm{h}}}{d_{{\rm{atom}}}}} \right)$$\end{document}$.

For the ^3^H~6~ → ^3^H~4~ transition of Tm:LiNbO~3~, the Tm-atom number density is *n* ~d~ = 1.89 × 10^19^ cm^−3^, *n* = 2.256, *Γ* ~h~ = 10 kHz, *γ* ~s~ = 2.6 kHz, and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\int} {\alpha (\sigma ){\rm{d}}\sigma }$$\end{document}$ = 497 cm^−2^ ^[@CR47]^, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\Theta _{ {\rm{i}}}} = Lc{\int} {\alpha (\sigma ){\rm{d}}\sigma }$$\end{document}$. The laser beam that is used to determine the absorption spectra is collimated and has a cross-section of *A* = *π*(80 μm)^2^, and the crystal length is *L* = 6.8 mm. For our AFC with *N* = 564 teeth, we measure *Θ* ~t~ = 4.3 MHz, giving $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$N_{\rm{t}}^{(1)} = 1.1 \times {10^9}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$N_{\rm{t}}^{(2)} = 1.7 \times {10^9}$$\end{document}$ atoms per tooth. The difference between the two estimates may be due to the measurement uncertainty or to the fact that not all transitions that contribute to the absorption line have identical properties.

Echo contrast *R* as a function of AFC bandwidth {#Sec17}
------------------------------------------------

In order to create a perfect W state, the AFC structure should be uniformly illuminated by the incoming photon. This guarantees that each AFC tooth has the same probability to contribute in the photon absorption. In the experiment, we use a 6 GHz FWHM bandwidth photon with a Lorentzian profile given by the transmission profile of the used Fabry--Perot filtering cavity (FP). The uniformity of the absorption over the different teeth should thus depend on the bandwidth of the AFC. In order to analyze this effect, we measure, for a fixed storage time (5 ns), the echo contrast *R* as a function of the prepared AFC bandwidth (see Fig. [8](#Fig8){ref-type="fig"}). Here *R* is expressed in terms of the percentage with respect to the maximum attainable *R* value given the known number of teeth created in each AFC. We observe that, as we narrow the AFC bandwidth, the relative *R* value increases. This is consistent with the expectation that we are approaching the creation of an ideal W state. Note that our derived bounds on the entanglement depth are still correct, albeit too conservative, for the case where the created state is not a perfect W state, since a state with unequal coefficients has to involve more entangled teeth in order to generate the same value of *R*.

Data availability {#Sec18}
-----------------

Experimental data that support the findings of this study have been deposited at Open Science Framework (OSF) with the accession codes "JZMCA (<https://osf.io/jzmca/>), DOI 10.17605/OSF.IO/JZMCA".

Electronic supplementary material
=================================

 {#Sec19}

Peer Review File

**Electronic supplementary material**

**Supplementary Information** accompanies this paper at doi:10.1038/s41467-017-00897-7.

**Publisher\'s note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank F. Fröwis and F. Bussières for very useful discussions and V. Kiselyov for technical support. This work was funded through NSERC, AITF and the DARPA Quiness program subaward contract number SP0020412-PROJ0005188, under prime contract number W31P4Q-13-1-0004. V.B.V. and S.W.N. acknowledge partial funding for detector development from the Defense Advanced Research Projects Agency (DARPA) Information in a Photon (InPho) program. Part of the detector research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. W.T. acknowledges funding as a Senior Fellow of the Canadian Institute for Advanced Research.

N.S. and C.S. conceived the project with some input from W.T. The theoretical approach was developed by P.Z., S.K.G., K.H., and C.S., and the calculations performed by P.Z. and S.K.G. with guidance from C.S. The experiments were developed and performed by C.D., N.S., G.H.A., P.L., M.G., and D.O. with guidance from W.T. The detectors were designed and fabricated by V.B.V., F.M., M.D.S., and S.W.N. The paper was written by P.Z., C.D., P.L., G.H.A., S.K.G., M.G.P., D.O., N.S., W.T., and C.S.

Competing interests {#FPar1}
===================

The authors declare no competing financial interests.
