\section*{Methods}
\textbf{Setup.}
The master laser is a frequency stabilized Ti:Sapphire laser (MBR-110, Coherent) with frequency $\nu_0$ of 340.69656 THz. The frequency of the master laser is stabilized to a low-drift Fabry-P\'{e}rot Interferometer (FPI) with PDH technique. The FPI is placed in vacuum housing with temperature controlled to a stability of 5 mK. The linewidth of the master laser is below 50 kHz.

Part of the laser output is fed into a cavity enhanced frequency doubler to generate blue light. The blue light is controlled with an -200-MHz AOM1. The photon pairs are filtered with the optical gratings and the etalon, combined with SMF collection. The signal photons are collected with SMF and ready for use. Another part of the master laser output is used for AFC preparation. The AOM2 with center frequency of -250 MHz is used in double-pass configuration to down shift the laser frequency. The AOM3 and AOM4 have a center frequency of 200 MHz and are also used in double-pass configuration. Thus the center frequency of the AFC preparation matches with the photon source at 340.69646 THz.

The AFC preparation light and single-photon pulses are combined with 10:90 BS and collected with SMF. The light focus to a diameter of 100 $\mu m$ with a lens (f=250 mm). The sample is placed in a cryostat (Oxford Instruments, SpectromagPT) at a temperature of 1.5 K and with a superconducting magnetic field of 0.3 T in the horizontal direction. The two parallel Nd$^{3+}$:YVO$_4$ crystals' c-axes are placed in the horizontal direction. Because each crystal only strongly absorbs $H$-polarized light, the $H$-polarized components of input photons are stored in the first crystal, and the $V$-polarized components are stored in the second crystal after polarization rotation by the HWP. The detection efficiencies for both photons are approximately 35\% with 50-Hz dark counts.

The initial absorption depth of the sample is $d\simeq8.0$. With optimized spectral tailoring, we can prepare AFC with extremely small background absorption $d_0\leq0.05$. The Zeeman population relaxation lifetime is measured by spectral-hole burning, $T_Z=42.9$ ms in the current crystals, which is 7 times longer than that of 10 ppm doped crystals \cite{SHB}. With two-pulse photon echoes, we measured an optical coherence time $T_2=28.4$ us (homogeneous linewidth $\Gamma_h=11$ kHz). Because of the better coherence properties and the more efficient optical pumping achieved in the current sample, the AFC efficiency is significantly higher than that obtained in 10 ppm doped crystals \cite{polarization}.

The preparation and measurement timing is controlled by three arbitrary function generators (AFG, Tektronix, AFG3252). It takes 12.5 ms for AFC preparation. To avoid the fluorescence noise caused by the classical pump light, the measurement cycle begins after waiting for time $T_w=1.5$ ms after the preparation cycle is completed. $N$ trials of single photon pulses are stored in the sample in the 10-ms measurement phase. The complete preparation and measurement cycles are repeated at a frequency of 40 Hz. A schematic drawing of the timing sequence of the experiment is given in the Supplemental Information.

\textbf{Details of the photon source.}
The phase matching temperature to generate degenerate photon pairs is 44.50 $^{\circ}$C.
The power of the 440 nm pump light is 20 mW for all the measurements. With continuous pump, the coincidence rate of the unfiltered photon pairs is 1.4 MHz. The heralding efficiency is 0.20. The diffraction efficiency of the gratings (1200 lines per millimeter) is 0.88. The transmission efficiency of the etalon is 0.95. The etalon has a bandwidth of 700 MHz and a free spectral range of 50 GHz. The temperature drift of the etalon is smaller than 10 mK. After spectral filtering, the coincidence rate of the photon pairs is 6.1 kHz and the heralding efficiency is 0.10. The photon source is switched into pulsed operations for interface with the AFC. The temporal width of blue pump light is 32 ns. The transmission efficiency from the photon source to the final SPD for the signal photons is approximately 0.23.

A figure of merit for the nonclassical nature of the photon correlations is the normalized cross-correlation function $g^{(2)}_{s i}=p_{s i}/p_sp_i$, where $p_s$ ($p_i$) is the probability to detect a signal (idler) photon and $p_{s i}$ is the probability to detect a coincidence in a specified time window. In practice, $p_{s i}$ and $p_sp_i$ are determined by the number of coincidences in the time window centered on and away from the coincidence peak, respectively. Assuming that the second-order auto-correlations of signal and idler $g^{(2)}_{x}$ (where x=`s' for signal or `i' for idler) satisfying $1\leq g^{(2)}_{x} \leq2$, then the non-classicality is proved by measuring $g^{(2)}_{s i}$ greater than 2 \cite{entangle,ref1,g2}.
We note that, compared with the unfiltered source, $g^{(2)}_{s i}$ rises significantly after spectral filtering. We measured $g^{(2)}_{s i} \approx 452$ between the transmitted signal photon and the idler photon. The temporal width of the coincidence window is 2 ns and the detection window for noise is 400 ns away from the transmitted signal peak. The $g^{(2)}_{s i}$ drops after storage because of low recall efficiency, the temporal broadening and the rising contributions of dark counts. The temporal width of the coincidence window for retrieved signal is 10 ns. The detection window for related noise is placed at 400 ns away from the retrieved signal peak. An example of the measurement on the $g^{(2)}_{s i}$ is given in the Supplemental Information.

We are most grateful to A. J. Leggett for his comments on the preliminary version of this
manuscript and to C. Brukner and M. Paternostro for their feedback on the revised version. This
work was supported by the National Basic Research Program (2011CB921200), the CAS, National Natural
Science Foundation of China (Grant Nos. 11274289, 11325419, 61327901.), the Fundamental Research Funds for
the Central Universities (wk2030380004 and wk2470000011), and the EU STREP project PAPETS.

