In recent quantum optical continuous-variable experiments, the number of fully inseparable light modes has drastically increased by introducing a multiplexing scheme either in the time domain or in the frequency domain. Here, modifying the time-domain multiplexing experiment reported in the work of Yokoyama *et al.* [Nat. Photonics **7**, 982 (2013)], we demonstrate the successive generation of fully inseparable light modes for more than one million modes. The resulting multi-mode state is useful as a dual-rail continuous variable cluster state. We circumvent the previous problem of optical phase drifts, which has limited the number of fully inseparable light modes to around ten thousands, by continuous feedback control of the optical system.

## I. INTRODUCTION

Quantum entanglement, which is often counterintuitive and was originally questioned by Einstein, Podolsky, and Rosen (EPR),^{1} is now recognized as a versatile resource for quantum information protocols.^{2,3} In order to meet various requirements for various applications, an important technological development is to increase the number of available entangled qubits. In particular, measurement-based quantum computation (MBQC) requires the cluster type of large-scale entangled state,^{4,5} where individual components of the cluster state must be accessible by measuring devices. The number of entangled qubits with individual accessibility seems to be currently limited to a moderate scale.^{6–8}

Shifting to the continuous variable (CV) counterpart,^{9,10} recently the number of fully inseparable qumodes of light fields drastically increased by introducing multiplexing either in the time domain^{11–13} or in the frequency domain.^{14–17} The successes of the CV approach are based on the deterministic nature of field-quadrature squeezing. This is a strong advantage over the photonic qubit-based approach, where photons or photons pairs are generated with low probabilities and rare success events are postselected.^{8,18} On the flip side, the disadvantage of the CV MBQC has been considered as finite squeezing in resource cluster states, which results in accumulation of noises during computation.^{19} However, there has still been a possibility of a special encoding which can circumvent this problem, and indeed, a fault-tolerant threshold corresponding to about −20.5 dB of squeezing is recently proven by Menicucci^{20} for the Gottesman-Kitaev-Preskill (GKP) type of encoding.^{21} Although the threshold of the squeezing levels is still unreached, it is good news that CV cluster states with finite squeezing can be a sufficient resource for quantum computation.

The multiplexing allows many qumodes to be assigned to a single beam, which drastically reduces complexity and physical size of the optical system from the conventional implementation assigning a single qumode to a single optical beam.^{22–26} Although both of the time-domain multiplexing and the frequency-domain multiplexing would be important developmental directions, the time-domain multiplexing is advantageous because of time translation symmetry. With the time-domain multiplexing, a cluster state of qumodes can be an arbitrarily long chain in the longitudinal direction by extending the operating time of the cluster-state generator.

In spite of the in-principle unlimitation of the time-domain-multiplexing, the previous experimental demonstration of a dual-rail CV cluster state was limited in the number of qumodes,^{13} due to a technical reason as follows. The optical setup is a large Mach-Zehnder interferometer as depicted in Fig. 1, including a long optical delay line for a shift of qumodes on a rail. The mathematical derivation of the cluster state with this setup is contained in Sec. S1 of the supplementary material. In order to operate this cluster-state generator, the relative phases at all interference points must be properly locked. For this purpose, modulated bright beams are injected into the optical paths of the cluster state as phase probes, and their classical interference signals are exploited as error signals for the feedback control. However, the modulated bright beams are noisy, which obscures objective quantum-level correlation signals of the cluster state. The noises were circumvented in the previous demonstration by chopping the bright beams and by detecting the cluster state when the bright beams are absent. Due to the above previous situation, the optical system was uncontrollable during the cluster-state generation. Phase drifts gradually accumulated and quantum correlations gradually degraded at the rear part of the cluster state. Finally, the full-inseparability criterion^{27} became unsatisfied after the creation of around 16 000 qumodes within 1.3 ms.^{13}

Here we demonstrate a new time-domain-multiplexing experiment, where a dual-rail CV cluster state with full inseparability is deterministically generated without limitation in the number of qumodes. The previous limitation explained above is removed by continuing the feedback control during the generation of a cluster state. Instead, the noises of the modulated bright beams are eliminated electrically after homodyne detections. The resulting effective squeezing levels are −4.3 dB, which are not as high as −5.0 dB for the first thousands of qumodes in the previous demonstration,^{13} but in contrast the initial squeezing levels continue without degradation. The squeezing levels are sufficiently higher than −3.0 dB of the full-inseparability criterion.^{13,27} In principle, there is no limitation in the number of fully inseparable qumodes. We test the full inseparability of about 1.2 × 10^{6} qumodes, which are generated within 100 ms. The reason why we stop at about 1.2 × 10^{6} qumodes is simply the data size for verification. In order to estimate the effective squeezing levels, we repeatedly acquire homodyne signals of the 1.2 × 10^{6} qumodes, by which the data size reaches 100 GB.

In Sec. II, we summarize theories regarding quantum correlations in the dual-rail CV cluster state, full inseparability criterion, and longitudinal mode functions. In Sec. III, we describe experimental methods. In Sec. IV, we show experimental results and make discussions. In Sec. V, we summarize our results. For self-containedness, in Sec. S1 of the supplementary material, we describe derivation of the quantum correlations in the dual-rail CV cluster state for an ideal case, and in Sec. S2 of the supplementary material, we describe derivation of inequalities for full inseparability.

## II. THEORY

### A. Quantum correlations in the cluster state

The procedure of producing the dual-rail CV cluster state is as follows. First, we produce a series of equal-time two-mode squeezed states (ordinary EPR states) in two beams [(ii) of Fig. 1]. This can be done by combining two single-mode squeezed states at a balanced beamsplitter. The time-slot width of the multiplexing is denoted by *T*. Next, we add delay to half-parts of the EPR states by the time-slot width *T*, by unbalancing the optical path lengths [(iii) of Fig. 1]. Finally, we combine the staggered EPR states by another balanced beamsplitter [(iv) of Fig. 1]. The resulting extended EPR (ExEPR) state is equivalent to the dual-rail CV cluster state up to local phase redefinitions.^{13}

In Sec. S1 of the supplementary material, the above procedure is mathematically followed for the ideal, infinitely squeezed case. A convenient way of representing the state, utilized in Sec. S1 of the supplementary material, is to specify the state by its nullifiers^{28} in the form of a linear combination of single-mode quadrature operators $ x \u02c6 \Lambda , k $ and $ p \u02c6 \Lambda , k $. Here, Λ ∈ {*A*, *B*} denotes the spatial mode index, and *k* ∈ ℕ denotes the temporal mode index, as depicted in (iv) of Fig. 1. Quadrature operators obey a commutation relation similar to that of position and momentum operators, $ [ x \u02c6 \Lambda , k , p \u02c6 \Lambda \u2032 , k \u2032 ] =i\u0127 \delta \Lambda \Lambda \u2032 \delta k k \u2032 $, where *δ*_{qq′} is the Kronecker delta. Sec. S1 of the supplementary material, is summarized in the following set of nullifiers:

The ideal ExEPR state is specified by the set of these nullifiers via the relations $ X \u02c6 k | ExEPR \u3009 =0$ and $ P \u02c6 k | ExEPR \u3009 =0$ for all *k*. These nullifiers also express the quantum correlations, which are experimentally observable. Note that all the nullifiers are mutually commutable, thus giving a simultaneous eigenstate.

In reality, the resource squeezing levels are finite. Therefore, when the nullifiers $ X \u02c6 k $ and $ P \u02c6 k $ are measured with respect to an experimental ExEPR state, there are residual noise variances,

Here, the bracket $ \u3008 O \u02c6 \u3009 $ represents the mean value of an observable $ O \u02c6 $. The squeezing parameters *r*_{x,k} and *r*_{p,k} correspond to those of the first and second single-mode squeezing sources in (i) of Fig. 1, respectively, when the subsequent transformations and measurements are ideal, which is apparent from calculations in Sec. S1 of the supplementary material.

Then, a natural question is what levels of squeezing are required for the entanglement. Here we resort to the van Loock-Furusawa criterion for full inseparability.^{13,27} A sufficient condition for full inseparability is

This condition corresponds to −3.0 dB squeezing of each nullifier. The derivation is contained in Sec. S2 of the supplementary material.

### B. Longitudinal mode functions

In our demonstration, we use continuous-wave (CW) local oscillators (LOs) for homodyne detections. Therefore, the homodyne detection system continuously outputs real values $ x \u02c6 \Lambda det ( t ) $ or $ p \u02c6 \Lambda det ( t ) $ with Λ ∈ {*A*, *B*}, which correspond to instantaneous quadrature operators $ x \u02c6 \Lambda ( t ) $ or $ p \u02c6 \Lambda ( t ) $ but filtered by the detection system with finite bandwidth. There is flexibility in how we define the actual qumodes via weight functions when integrating these continuous signals. However, we must be careful about the orthogonality of qumodes, because the longitudinal modes are deformed by the introduced electric filters, as explained below.

The quadrature operator $ x \u02c6 \Lambda , k $ of each qumode, specified by a longitudinal mode function *f _{k}*(

*t*) =

*f*

_{0}(

*t*−

*kT*) which is preferred to be contained in the

*k*-th time slot

*kT*≤

*t*≤ (

*k*+ 1)

*T*, is connected to the instantaneous quadrature operators $ x \u02c6 \Lambda ( t ) $ via the relation

However, what we can directly choose in the experiment is not the optical longitudinal mode function *f _{k}*(

*t*) but a computational weight function

*g*(

_{k}*t*), with which we perform weighted integration,

This integration is totally a postprocessing after acquisition of full-time quadrature values $ x \u02c6 \Lambda det ( t ) $. The situation is the same for the conjugate quadrature operator $ p \u02c6 \Lambda , k $. Note that the integration is actually replaced by a summation in accordance with the sampling rate of data acquisition.

If the response of the homodyne detection system is sufficiently flat, we can identify *g _{k}*(

*t*) with

*f*(

_{k}*t*). However, this is not the case in this demonstration, in contrast to the previous demonstration,

^{13}because of the electric filters applied to the homodyne signals for the reasons explained in Sec. I and in more detail in Sec. III B. They are connected as

*g*(

_{k}*t*) ∝ (

*e**

*f*)(

_{k}*t*) via the impulse response function

*e*(

*t*) of the detection system including the electric filters, where * denotes convolution. Even if we choose the weight function

*g*(

_{k}*t*) so that they are nonzero only in the time slot

*kT*≤

*t*≤ (

*k*+ 1)

*T*, by which the weight functions themselves are orthonormal ∫

*g*(

_{k}*t*)

*g*

_{k′}(

*t*)

*dt*=

*δ*

_{kk′}, the orthonormality of actual optical modes ∫

*f*(

_{k}*t*)

*f*

_{k′}(

*t*)

*dt*=

*δ*

_{kk′}is not guaranteed in general. Note that the orthonormality of the mode functions {

*f*(

_{k}*t*)}

_{k∈ℕ}is equivalent to the orthonormality of the qumodes via the relation

where *δ*(*t*) is the Dirac delta function.

The filter is a combination of both high-pass and low-pass filters, and the parameters are described in Sec. III B. A filter delays signals in general. Therefore, when the response *e*(*t*) is deconvoluted from the weight function *g _{k}*(

*t*), the optical mode function

*f*(

_{k}*t*) is stretched in the forward direction, which may ruin the mode independence to some extent. Especially, the effect of a high-pass filter is not negligible in our setup, which is expected to produce a long tail in the forward direction. In order to suppress the non-orthogonality of the actual optical modes, a possible method is to use a weight function localized to the latter part of the time slot. This is a workable method, but is not sufficient in our case, as will be discussed by using autocorrelations of shot noises in Sec. IV B.

Therefore, here we use a weight function which has both positive and negative parts. By adjusting the balance between the positive part and the negative part, the long tail produced by the high pass filter is canceled. The actual un-normalized weight functions utilized in the experimental data analysis are shown in Fig. 2. The un-normalized weight functions are the following functions:

discretized at intervals of 10 ns. The width of the weight function window *T _{w}* is 120 ns, and the envelope fall time 1/

*γ*is 40 ns. The balance between the positive part and the negative part is adjusted through an optimization parameter

*t*, which is determined to be 2 ns. The center of the weight function is

_{c}*t*=

_{k}*t*

_{0}+

*kT*, where

*k*∈ ℕ is the temporal index as usual. The mode interval

*T*is 160 ns, and

*t*

_{0}= 95 ns in Fig. 2. Note that we do not have to care about the normalization of the weight functions in the actual data analysis, because the obtained quadrature values are normalized by using corresponding shot noises acquired in the same conditions.

On the other side of the coin, the wavy weight functions lead to increased contribution of higher frequency components, compared to the previous demonstration.^{13} However, the squeezing level degrades at higher frequencies with our squeezing sources, due to a finite bandwidth of optical cavities. The resulting correlation levels of −4.3 dB are not as high as the previous −5.0 dB for this reason.

## III. EXPERIMENTAL METHODS

### A. Optical setup

Figure 3 shows the experimental setup. The light source is a continuous-wave (CW) Ti:sapphire laser (SolsTiS, M Squared Lasers) operating at the wavelength of 860 nm and the power of about 2 W. About 1 W of the laser power is sent to a second harmonic generation (SHG) cavity. It is a bow-tie shaped cavity with a round-trip length of 500 mm, containing a KNbO_{3} (KN) crystal as a nonlinear optical medium. The generated CW second harmonic light is used as pump beams for two optical parametric oscillators (OPOs) operating below threshold. The power of the pump beam for each OPO is about 150 mW. Each OPO is a bow-tie shaped cavity with a round-trip length of 230 mm, containing a periodically poled KTiOPO_{4} (PPKTP) crystal (Raicol Crystals) as a nonlinear optical medium. It generates a squeezed light beam, which is treated as a stream of single-mode squeezed vacuum states. The bandwidth of the squeezing basically corresponds to the bandwidth of the OPO cavity, whose half-width at half-maximum is 17 MHz.

The two squeezed light beams from the OPOs are sent to an asymmetric Mach-Zehnder interferometer, composed of two balanced beamsplitters and an optical delay line. The delay line is in one of the two arms of the Mach-Zehnder interferometer, and implemented with an optical fiber. The length of the optical fiber is about 30 m, from which the time slot width *T* of the multiplexing is about 160 ns. Both ends of the fiber are anti-reflection coated, and the coupling efficiency is maximized by using special fiber aligners (FA1000S, First Mechanical Design). Total amount of losses regarding the optical fiber, including the coupling inefficiency and the propagation losses, was about 11%.

In order to characterize the resulting ExEPR state, the two output beams are each subject to a balanced homodyne detection. The bandwidth of the homodyne detectors without the electric filter is more than 30 MHz. The power of a CW local oscillator (LO) beam for each homodyne detection is 10 mW. The transverse mode of the LO beams is cleaned and stabilized by a cavity. The interference visibility between the ExEPR beams and the LO beams were about 97% on average. The quantum efficiency of the homodyne photodiodes is about 99%.

### B. Modulations, filters, and data acquisition

In order to probe the optical phases and thereby to lock the relative phases of all optical interferences by feedback control, modulated bright probe beams are utilized. At each interference point, there is a corresponding photo-detector, and demodulation of the detector signal gives an error signal, which is fed back to a piezoelectrically actuated mirror, as depicted in Fig. 3. First, a phase modulation at 32 MHz is commonly added by an electro-optic modulator (EOM) just after the Ti:sapphire laser output. This modulation is utilized for the Pound-Drever-Hall lock^{29} of the optical cavities including the OPOs. Next, two probe beams are phase modulated by piezoelectric transducers (PZTs) at 231 kHz and 326 kHz, respectively, and then injected into the individual OPOs from the back side of a highly reflective mirror composing of the OPO. The power of the probe beams is adjusted to about 2.5 *μ*W at the output of the OPOs. These lower-frequency modulations are utilized for the lock of the following three. First, the phase-sensitive parametric amplification at each OPO is locked to the point minimizing the probe beam power, by which the phase of the squeezed light is associated with that of the probe beam. Second, the two interferences of the Mach-Zehnder interferometer are locked, where the beat signal of the two modulations is exploited by demodulating at the difference frequency of 95 kHz. Third, the LO phases are locked so that the homodyne detectors measure either the objective $ x \u02c6 \Lambda det ( t ) $ or $ p \u02c6 \Lambda det ( t ) $ quadrature values. The measurement basis is selected via the choice of the demodulation frequency between 231 kHz ad 326 kHz. These modulations of bright probe beams become unwanted noises in the homodyne detections, which wipe out the objective quantum correlations and thus must be removed somehow.

As stated in Sec. I, these noises were eliminated at the level of the optical system in the previous work,^{13} employing periodic switching. The probe beams were alternately turned on and off by using a pair of acousto-optic modulators (AOMs), and the feedback control of the optical system was active when the probe beams were on, while the optical system was held and the ExEPR state was tested when the probe beams were off. In that demonstration, uncontrollable drifts of the optical system were negligible for first several thousands of qumodes, but gradually degraded the correlation in the ExEPR state, finally spoiling the full inseparability after creation of about 16 000 qumodes. In contrast, here we employ the method of eliminating the noises electrically after the homodyne detections, by using frequency filters. In this demonstration, the switching AOMs are always on, and the feedback control of the optical system continues during creation of the ExEPR state.

The electric filter applied to each homodyne signal in this demonstration is a serial combination of a 5th-order high pass filter with the cutoff frequency of 1.5 MHz, a 3rd-order inverse Chebyshev low-pass filter which has a notch characteristic at 32 MHz, and a 3rd-order low-pass filter with the cutoff frequency of 40 MHz for the purpose of anti-aliasing. After the electric filtering and appropriate low-noise amplification, the homodyne signals are digitized by an 8-bit oscilloscope (DPO7054, Tektronix) with the sampling rate of 100 MHz. The data acquisition frame width is 100 ms, and we acquire 1500 frames of data for each measurement quadrature in order to determine the residual noise variance of each nullifier with acceptable precision. In the same condition, we also acquire optical shot noises for references, which correspond to vacuum fluctuations. The shot noises are obtained by simply blocking the optical beam paths of the ExEPR dual rail.

Figure 4 shows the power spectra from the detection system with and without the electric filters, obtained from the fast Fourier transform (FFT) of oscilloscope data. We can see that the signals around the modulation frequencies are largely attenuated by the filters.

## IV. RESULTS AND DISCUSSION

### A. Power spectra of homodyne signals

Since now the ExEPR state is stabilized by the continuous feedback, first we discuss the power spectra of the homodyne signals. Figure 5 shows the power spectra of the quadrature signals $ x \u02c6 \Lambda det ( t ) $ and $ p \u02c6 \Lambda det ( t ) $ of the output ExEPR beams, together with that of the optical shot noise $ x \u02c6 \Lambda vac,det ( t ) $ for references, for each of the two output ports Λ ∈ {*A*, *B*}. We can see that the output of the ExEPR state generator is, in the frequency domain, oscillating between $ x \u02c6 $-squeezed states and $ p \u02c6 $-squeezed states, and that the squeezed quadratures are opposite between the two output ports. This interesting point was not discussed in the previous work.^{13}

By using the schematic in Fig. 1, we discuss how this periodic behavior in the frequency domain is understood. The situation is that two squeezed beams, one is $ x \u02c6 $-squeezed and the other is $ p \u02c6 $-squeezed, enter a Mach-Zehnder interferometer with asymmetric lengths of arms. When we focus on a single sideband frequency component, the unbalance of the arm lengths acts as just a phase shift. The amount of the phase shift is *θ* = 2*πfT*, proportional to the sideband frequency *f* and the difference between the propagation times of the two arms *T*. Then, the two-input two-output Mach-Zehnder interferometer as a whole works as a single beamsplitter with a reflectivity *R* = cos*θ* dependent on the phase shift *θ*. At the phases of *R* = 1, the $ x \u02c6 $-squeezed beam fully goes to one output port and the $ p \u02c6 $-squeezed beam the other output port, and at the phases of *R* = 0, the opposite occurs. Between the above two cases, the beamsplitter combines the two squeezed beams and makes entanglement between the output ports. These cases periodically appear in the frequency domain because the phase shift *θ* is proportional to the sideband frequency *f*. The period of the oscillation corresponds to the propagation time difference *T*.

### B. Orthogonality of qumodes

We apply the weighted integration to the above homodyne signals by using weight functions described in Sec. II B. Before that, we check the orthogonality of the qumodes by applying the weighted integration to the optical shot noises,

An unfiltered shot noise $ x \u02c6 \Lambda vac ( t ) $ is a white noise whose autocovariance is represented by a Dirac delta function. Therefore, if the integrated shot noise values of a temporal mode $ x \u02c6 \Lambda , k vac $ have some correlation with those of neighboring temporal modes $ x \u02c6 \Lambda , k + m vac $, *m* ≠ 0, the correlation is due to the overlap between the optical mode functions *f _{k}*(

*t*) and

*f*

_{k+m}(

*t*), via the relation

Figure 6 shows the correlation of the integrated shot noise values among neighboring qumodes,

for |*m*| ≤ 5. As we can see with trace (i), the correlation is almost canceled, |*C*_{Λ}(*m*)| < 0.01 for *m* ≠ 0, when the weight functions *g _{k}*(

*t*) in Eq. (7) and Fig. 2 is used. For comparison, trace (ii) shows the correlation when only the positive parts of the weight functions $ g k + ( t ) \u221dmax ( g k ( t ) , 0 ) $ are used. In spite of the spacing between neighboring weight functions, a negative correlation

*C*

_{Λ}(±1) ≈ − 0.13 is observed, which is mainly due to the high-pass filter. This correlation is canceled by the negative parts of the weight functions. An insight regarding this cancellation is that the wavy weight functions make little use of the low-frequency components, which are affected by the high-pass filter and cause overlap between qumodes.

Note that the fractions of a qumode invading neighboring qumodes correspond to square of the correlation coefficients |*C*_{Λ}(*m*)|^{2}, which are very small and thus experimentally negligible.

### C. Full inseparability

Finally we check the inseparability condition. Figure 7 shows the residual variances of the nullifiers, $ \u3008 X \u02c6 k 2 \u3009 $ and $ \u3008 P \u02c6 k 2 \u3009 $ in Eq. (2), together with the variances of the same observables calculated from the shot noises for references. We can clearly see the variances well squeezed below the −3 dB bound for full inseparability, for up to *k* = 6 × 10^{5}. Taking into account the two spatial mode indices *A* and *B*, the number of qumodes tested here is 1.2 × 10^{6}. As mentioned in Sec. III B, the variance of each nullifier is averaged over 1500 times. The average and the standard deviation of the squeezing levels over the 6 × 10^{5} temporal indices are − 4.3 ± 0.2 dB for $ X \u02c6 k $ and − 4.3 ± 0.2 dB for $ P \u02c6 k $. The worst squeezing level among the 6 × 10^{5} temporal indices is −3.6 dB for $ X \u02c6 k $ and −3.5 dB for $ P \u02c6 k $. They will approach the average −4.3 dB if we could increase the number of data frames for the averaging.

In contrast to the previous work,^{13} the squeezing levels do not degrade during the generation of over-one-million qumodes. This is owing to the stabilized optical system by introducing the continuous feedback.

## V. CONCLUSION

We have demonstrated the deterministic generation and verification of a fully inseparable dual-rail CV cluster state consisting of more than one million qumodes of light, by employing a time-domain multiplexing scheme. Compared to the previous work,^{13} the cluster state generator does not degrade during operation, owing to the continuous feedback of the optical system. We can in principle further increase the number of qumodes, but we stopped at around one million qumodes because of the data size for verification. Time domain multiplexing is one of the key technologies of CV information processing. We note that our scheme is potentially upgraded to two-dimensional CV cluster states, where one dimension is finite and the other dimension is unlimited, by combining two optical delay lines.^{12} Open questions from technological aspects are how to combine the cluster states with photon counters aiming at universal processing,^{28} as well as how to increase the squeezing levels toward fault tolerance. We mention recent experiments of optical filters for quantum states,^{30} as well as real-time integration of homodyne signals,^{31} which in the future may be combined with this entanglement multiplexing technologies.

## SUPPLEMENTARY MATERIAL

See supplementary material (including Refs. 32 and 33) for derivation of the quantum correlations in the dual-rail CV cluster state for an ideal case (Sec. S1 of the supplementary material), as well as for derivation of the inequalities for full inseparability (Sec. S2 of the supplementary material).

## ACKNOWLEDGMENT

This work was partly supported by PDIS and APSA commissioned by the MEXT, JSPS KAKENHI, CREST of the JST, and the SCOPE program of the MIC of Japan, and ARC CQC2T (Grant No. CE1101027) of Australia. S.Y. acknowledges support from ALPS and JSPS Postdoctoral Fellowships for Research Abroad of Japan.