Spectral- and time-multiplexing are currently explored to generate large multipartite quantum states of light for quantum technologies. In the continuous variable approach, the deterministic generation of scalable entangled states requires the generation of a scalable number of squeezed modes. Here, we demonstrate the simultaneous generation of 21 squeezed spectral modes at the repetition rate of our laser, i.e., 156 MHz. We exploit the full repetition rate and the pulse shaping of a femtosecond light source to combine, for the first time, frequency- and time-multiplexing in multimode squeezing. This paves the way for the implementation of multipartite entangled states that are both scalable and fully reconfigurable.

## I. INTRODUCTION

Multimode quantum light fields enable continuous-variable (CV)-based quantum information technologies. In particular, light fields squeezed into different spatial, spectral, and temporal shapes are a versatile and promising resource for measurement-based quantum computing (MBQC),^{17,28} quantum-enhanced sensing,^{16,34,46} multiparty communication protocols,^{4,9} as well as simulation^{30} and quantum-enhanced machine learning protocols.^{31}

The leverage of the continuous-variable approach, which encodes quantum information in the amplitude and phase values of the light field, is the deterministic generation of multipartite entangled states using nonlinear optical processes. CV-cluster states, a necessary resource for MBQC,^{17,28,36,37} are the largest entangled structures by from spectral and temporal modes of squeezed light.^{5,9,11,25,45} Beyond the amount of generated squeezing, two crucial features of CV entangled quantum resources are the scalability, i.e., the number of entangled parties, and the reconfigurability, i.e., the capability to reshape the entanglement links at will. The CV entangled states are generated by mixing several squeezed optical modes using linear-optics operations or, more generally, via mode-basis changes. As expected, the scalability depends on the number and the production rate of the initial squeezed modes. Differently from the discrete variable case,^{20,40,43} the number of required squeezed modes also strongly depends on the shapes of the entangled structures.^{10,17} It is then crucial to easily access basis changes to reshape the quantum links.

Here we demonstrate the generation of multimode squeezed states that are multiplexed both in the temporal and spectral degrees of freedom. While temporal multiplexing already enabled the generation of the largest CV-cluster states,^{5,25} multimode squeezing in the spectral modes of a femtosecond source associated with well-controlled ultrafast pulse shaping techniques enables the full reconfigurability of the entanglement links as well as mode-selective non-Gaussian operations.^{9,35} Both elements are present in our experimental setup: time multiplexing is reached by measuring squeezing in individual laser pulses emitted at a repetition rate of 156 MHz, and 21 squeezed spectral modes are generated for each pulse.

The quantum resource is obtained via spontaneous parametric down conversion (SPDC) in a waveguide in a single-pass configuration. Nonlinear waveguides are promising platforms for the controlled generation of both single- and multimode quantum states.^{1,8,14,18} Indeed, they have large nonlinearities and can be engineered to tailor the modal decomposition of the generated quantum state.^{2} In particular, periodically poled Potassium Titanyl Phosphate (PPKTP) waveguides provide suitable phase-matching conditions for modal engineering compatible with mode-selective femtosecond homodyne detection.^{38} Periodically poled waveguides have demonstrated fully integrated^{21} and large^{22} single-mode squeezing. In particular, PPKTP crystals have shown controlled squeezing generation in amplified schemes.^{32,42} While this work focuses on the generation of a large number of squeezed spectral modes, several integrated platforms can be engineered to limit the number of spectral and temporal modes involved in order to increase the level of squeezing per mode.^{11,27,44} In the temporal domain, state-of-the-art experiments have been performed to generate cluster states^{5,25} and boson sampling experiments.^{27}

Here, we multiplex squeezed states for the first time, both in time, i.e., at a high repetition rate (156 MHz), and in reconfigurable spectral shapes (ultrafast pulse shaping). Moreover, we show entanglement correlations between optical-frequency bands. The demonstrated resource, composed of a nonlinear waveguide and an ultrafast mode-selective homodyne with large electronic bandwidth, is the pivotal building block for scalable and reconfigurable networks (see Fig. 1). As shown in the lower part of Fig. 1, by combining two waveguide sources or more^{5,25} with delay lines and linear optics, the entanglement structure in the optical frequency^{9} can be extended in the time (pulse-based) domain.

## II. RESULTS

### A. Spectrally multimode squeezed light

We use spontaneous parametric down conversion (SPDC) in a periodically poled Potassium Titanyl Phosphate (PPKTP) waveguide to generate spectrally multimode squeezed states of light at a repetition rate of 156 MHz. The train of 22-fs pulses emitted by the light source, is converted via second harmonic generation (SHG) and is used to pump a 1-mm long PPKTP waveguide. The waveguide is designed for type 0 phase-matching, corresponding to the largest nonlinear coefficient in KTP. While most multimode squeezed sources of light to date were obtained using resonant cavities,^{5,9,11,25,39} we reach high pump intensity via light confinement into a waveguide. This allows for the generation of trains of quadrature-squeezed pulses in a single passage.

*m*can be expressed in an eigenbasis of squeezed frequency modes,

^{39}

^{,}

*S*

_{j}(

*ω*), that we call

*supermodes*, and the {

*λ*

_{j}} are related to the amount of squeezing in the spectral modes. In a collinear degenerate SPDC, where signal and idler fields are indistinguishable, the eigenmodes can be retrieved using a Bloch–Messiah reduction (equivalent to a Schmidt decomposition),

^{19,26}where their shape depends on the phase-matching and the pump spectrum. In a Gaussian approximation, i.e., an unchirped pump with a Gaussian-shaped spectrum, the supermodes can be approximated with Hermite–Gauss (HG) functions.

^{33}

The high nonlinearity of KTP, the pump intensity, and the length of the poled region set the total amount of squeezing that can be produced, while its distribution into a given number of spectral modes depends on the phase-matching configuration, the spectral width of the pump, and the waveguide dimensions.^{12,38} The estimated Schmidt parameter in our configuration is K = 98. We chose such a configuration to show the ability of the system to generate a high number of optical spectral modes that can be independently measured via homodyne detection.

We reveal the amount of squeezing in multiple spectral modes by shaping the spectrum of the homodyne reference beam, named the local oscillator (LO).^{9,39} After light detection, we obtain the variance of the spectral-mode quadrature by measuring the sideband intensity noise on a spectrum analyzer. The squeezed (anti-squeezed) quadrature displays noise below (above) the shot noise.

Table I shows the squeezing measurement results and the corresponding spectral shape of the LO. We measured squeezing in 21 Hermite-Gauss (HG) supermodes (up to HG_{20}) with 0.47 dB in the leading HG_{0} mode, which has an 18-nm bandwidth (FWHM), in good agreement with the simulations. The width of Hermite-Gauss functions is known to increase with their order *n*; therefore, the number of supermodes that we can experimentally access is limited by the 40-nm spectrum of the LO. As a result, the measured squeezing level of higher order modes decreases faster than in numerical simulations due to a poorer match between the set LO mode and the theoretical supermode. Since we chose a configuration with a large number of squeezed spectral modes, we expected the squeezing per mode to be low, despite the experiment being low-loss. Indeed, by keeping the pump power constant as well as the length of the poling region, i.e., the global strength of the nonlinearity, numerical simulations^{38} showed that by modifying the spectral width and the waveguide dimensions, we can access configurations with a lower number of squeezed modes and, consequently, a larger value of squeezing per mode. In such configurations, we then have a trade-off between the number of generated squeezed modes and the amount of squeezing per mode.^{12,38}

Mode . | Sqz (dB) . | Asqz (dB) . | |
---|---|---|---|

HG_{0} | −0.47 | 0.55 | |

HG_{1} | −0.33 | 0.42 | |

HG_{2} | −0.23 | 0.35 | |

HG_{3} | −0.20 | 0.28 | |

HG_{4} | −0.17 | 0.30 | |

HG_{5} | −0.18 | 0.30 | |

HG_{6} | −0.18 | 0.27 | |

HG_{7} | −0.14 | 0.26 | |

HG_{8} | −0.15 | 0.26 | |

HG_{15} | −0.09 | 0.17 | |

HG_{20} | −0.08 | 0.16 |

Mode . | Sqz (dB) . | Asqz (dB) . | |
---|---|---|---|

HG_{0} | −0.47 | 0.55 | |

HG_{1} | −0.33 | 0.42 | |

HG_{2} | −0.23 | 0.35 | |

HG_{3} | −0.20 | 0.28 | |

HG_{4} | −0.17 | 0.30 | |

HG_{5} | −0.18 | 0.30 | |

HG_{6} | −0.18 | 0.27 | |

HG_{7} | −0.14 | 0.26 | |

HG_{8} | −0.15 | 0.26 | |

HG_{15} | −0.09 | 0.17 | |

HG_{20} | −0.08 | 0.16 |

Importantly, while we have revealed the number of squeezed modes, this source could be reshaped via the adaptation of the measurement device into arbitrary entangled states, making it a promising resource for CV protocols.^{9}

### B. Simultaneous time and frequency multiplexing

To verify that the set of spectrally squeezed modes (Fig. 1) is produced at each pulse *m* of the train, we characterize the signal with the temporal resolution of the individual pulses. To do so, we have developed a homodyne detector with a wide electronic bandwidth, a cutoff frequency larger than the laser repetition rate, and a large common mode rejection ratio (CMRR).^{24} We acquire the homodyne signal with a fast oscilloscope; one quadrature value is associated with a single pulse by integrating the scope signal on the pulse time window. The latter has to be carefully adjusted to precisely discriminate between consecutive pulses. Details of the experimental protocols can be found in Sec. IV.

Using this procedure, we performed a pulse-resolved (pulse-by-pulse) squeezing and antisqueezing measurement in seven HG modes (up to HG_{6}). The measured values are comparable with the non-pulse-resolved sideband squeezing measurements performed with the spectrum analyzer in Table I, as shown in Fig. 2. The small discrepancies observed between the two datasets can be explained by the lower signal-to-noise ratio in the pulse-by-pulse measurement.

Therefore, we revealed the same multimode structure in the pulse-by-pulse measurement and the sideband measurement. The results confirm that the multimode spectral structure exists at a single pulse level; this opens the way to the generation of entangled structures at the repetition rate of the laser. It is also possible to use delay lines or multiplex multiple waveguided sources^{5,25,27} to engineer more elaborate entangled structures from light pulses and spectral modes, as suggested in Fig. 1.

### C. Covariance matrix and entanglement

In Secs. II A and II B, we presented the characterization of down-converted light on the Hermite-Gauss basis of optical frequencies, where the system is described as a collection of independently squeezed states. However, it is also possible to measure this state on other bases, revealing the presence of entanglement between the modes of a chosen basis. Here, we show the reconstruction of the state covariance matrix on a frequency-band basis, measuring amplitude and phase noise correlations between different frequency bands of the output light spectrum. Such a basis is the most commonly used for describing the parametric process. It is also the one where the different modes can be easily spatially separated (via dispersive elements) and distributed.

For the measurement, we divide the LO spectrum into eight individual frequency bands of 5 nm that we call *frexels*. If we call *q*_{i} and *p*_{i} the quadratures associated with each of these frexel modes, the covariance matrix is defined as $Vxi,xj=(\u27e8\delta xi\delta xj\u27e9+\u27e8\delta xj\delta xi\u27e9)/2$, where *x* is either *q* or *p* and *δx* = *x* − ⟨*x*⟩. We measure it using the pulse shaper to address the spectrum of each frexel as well as all combinations of two frexels. This leads to the measurement of 36 homodyne traces. The reconstructed covariance matrices can be seen in Fig. 3 (see Sec. IV for details). The diagonal elements show the noise in each frequency band, while the off-diagonal terms reveal the correlations between the different frequency bands. The strong correlations observed in the antidiagonal are typical of the joint spectral structure of the generated photon-pair in the chosen type-0 SPDC.

The off-diagonal correlations can be used to evaluate entanglement between any bipartition of the eight frequency bands.^{15} We have measured entanglement in 114 of the 127 possible bipartitions. The amount of entanglement, measured via a positive partial transpose (PPT) criterion, varies along the different bipartitions. It reflects the structure of the narrow type-0 joint-spectrum of the entangled photon-pairs generated by the non-linear process in the waveguide. The quantum correlations, as also shown in the covariance matrix, are mostly shared between frequency bands symmetrically spaced from the central optical frequency, i.e., 4–5, 3–6, and 2–7. For this reason, bipartitions that do not decouple these pairs of frexel are not guaranteed to be entangled. Details can be found in Sec. IV.

From the reconstructed covariance matrix, we can also recover the squeezing values associated with the uncorrelated squeezed modes of the system, i.e., the supermodes. To do so, we diagonalize the covariance matrix. The results are shown in Fig. 4, where the diagonal elements are compared with the squeezing values shown in Fig. 2. Small deviations are due to the larger resolution (larger signal-to-noise ratio) of the measurement in the supermode basis compared to the 8-pixel basis.

We have shown that the basis change, via mode-selective homodyne detection, can be used to build/modify the entanglement correlation between optical frequency modes. This is an essential precondition for building the networks of entangled modes shown in the lower right part of Fig. 1.

## III. DISCUSSION

We demonstrated a single-pass source producing 21 independently squeezed spectral modes at 156 MHz. This allows for the simultaneous multiplexing of the temporal and spectral degrees of freedom of deterministic quantum resources. We produced the largest number of spectral modes of squeezed vacuum measured via homodyne detection in a single-pass configuration. Tailoring of the phase-matching condition and the pump spectrum can be used to increase the amount of squeezing per spectral-mode by decreasing the number of squeezed modes. We can in fact access only a limited fraction of the predicted 98 modes in the actual configuration, given the large spectral span of the high-order modes. The main result of this work is the simultaneous multiplexing of both degrees of freedom; this constitutes a versatile building block for the three-dimensional structures that are required to build a scalable quantum computer.^{7} Large two-dimensional CV cluster states^{5,25} have been generated by combining time modes from multiple squeezing sources, and a three-dimensional reconfigurable structure obtained via time multiplexing of a squeezing source has been recently exploited in a boson sampling setting.^{27} Here, we demonstrate a resource that, by offering multiple squeezed modes within a given pulse, can be shaped as a three-dimensional entangled structure (as shown in the scheme of Fig. 1). The reconfigurability in our scheme is offered by multiplexing via spectrally resolved homodyne detection.^{3} Therefore, different two-dimensional spectral structures can be concatenated in the time domain as fast as the repetition rate of the pulsed laser source (156 MHz). Moreover, spectrally mode-selective non-Gaussian operations,^{35} needed for quantum computing protocols,^{7,13} can be implemented on selected nodes of the 3-D structure.

## IV. METHODS

### A. Experimental setup

The laser source is a femtosecond mode-locked Titanium-Sapphire (Ti:Sa) (Femtosource Synergy+ from Spectra-Physics), delivering trains of 22-fs pulses at a repetition rate of *f*_{rep} = 156 MHz (*τ*_{rep} ≃ 6 ns). The corresponding optical spectrum is centered at 795 nm with a FWHM of about 42 nm. At the laser output, the average power is 1 W. The IR beam is split into three beams: 30% is sent to a pulse shaper and serves as a reference for the homodyne detection; 10% is used to seed the SPDC; and the final 60% is used to pump a 1-mm long BiBO (Bismuth Triborate) crystal, critically phase-matched for the second-harmonic generation. After the SHG, the pulsed source is centered at 397.5 nm and has a spectrum of 0.7 nm (FWHM). The up-converted light pumps a periodically poled KTP (Potassium Titanyl Phosphate) waveguide and undergoes a collinear type-0 SPDC process generating a train of quadrature-squeezed pulses. The temperature of the waveguide is set to its value of optimal conversion efficiency (89.1 °C). The pump and seed beams are coupled inside the waveguide with an achromatic lens, while the exiting beams are collimated with an IR-coated lens, and the residual pump light is filtered out with a dichroic mirror. Inside the waveguide, the temporal matching of the seed and the pump pulses is achieved by optimizing the phase-sensitive parametric amplification to ∼20%. Such a value, although consistent with the presence of squeezing in the system, cannot be used to evaluate the squeezing values of the spontaneous process, as the spectrum of the amplified seed is not that of an eigenmode of the process.

The quantum state generated by the parametric process is characterized via homodyne detection. To unveil the multimode structure of the squeezed light, we use a pulse shaper to modify the LO spectrum according to the spectral mode that we want to measure. The overall transmission of the setup is 70%, and it can also be used to correct for undesired spectral dispersion.

The shaped LO is interfered with on a balanced beam splitter with the seed beam first for alignment, then with the signal. Their relative phase is scanned with a piezo-electric actuator placed in the LO path. The visibility of the interference fringes of the LO and the seed beam is maximized to reach the equivalent mode-matching efficiency *η*_{mm} = 90%, a value that accounts for both the temporal and spatial overlap, while the spectral overlap between the signal and the LO is optimized by the pulse-shaping. A scheme of the experiment is shown in Fig. 5.

### B. Reconstruct the covariance matrix

*q*and

*p*are uncorrelated, thus the covariance matrix

**V**takes a block-diagonal form

**= {**

*x***q**,

**p**}. We denote

*i*+

*j*as the mode combination of the frexels

*i*and

*j*and we define its quadrature operators as $x\u0302i+j=(x\u0302i+x\u0302j)/2$, provided that the optical power detected by the frexels

*i*and

*j*is the same. Given that the mean value of the quadrature is ⟨

*q*⟩ = ⟨

*p*⟩ = 0, we compute the off-diagonal elements with the formula

*q*quadrature) or the maxima (for the

*p*quadrature) of the homodyne traces. To recover the uncorrelated squeezed modes of the system, i.e., the supermodes, the covariance matrix needs to be diagonalized. To do so, we use a Gram–Schmidt process to orthogonalize the eigenmodes associated with the antisqueezed modes, and we use them to construct the matrix applied to the covariance matrix to operate the basis change. The relative sign between the eigenvalues of

**V**

_{qq}(equivalently

**V**

_{pp}) stands for the relative phase between the squeezing ellipses of the different supermodes. The uncertainty of the squeezing values stems from the uncertainty of the squeezing or anti-squeezing peaks measured with the spectrum analyzer.

### C. Entanglement between different bipartitions

^{41}We measure the PPT value for each of the possible 127 bipartitions, i.e., different ways to divide the eight frequency bands into two sets. This is done by carrying out a partial transpose operation on the density operator of the state that maps its covariance matrix

**V**into $V\u0303=\Lambda V\Lambda $, where Λ is the operation that changes the sign of the

*p*coordinate of one of the two bipartitions. A necessary condition for separability of the state reads $P=V\u0303\u2212iJ\u22650$, where

**J**is the 2

*n*× 2

*n*symplectic form

**P**(PPT value) for each possible bipartition is shown in Fig. 6, where we see that inseparability holds for 90% of the bipartitions. Moreover, we can distinguish four different bands of PPT values: a lower, a lower-intermediate, an upper-intermediate, and an upper band. The lower band is the band with the highest degree of entanglement (lower PPT value), and it corresponds to a partitioning choice where frexels 4 and 5 are separated in two different bipartitions. The lower-intermediate band contains bipartitions that separate frexel 3 and 6 (but 4 and 5 are in the same one), while the upper-intermediate band sees frexels 2 and 7 in separate bipartitions (but 4–5 and 3–6 are in the same one). Finally, the upper band is mostly composed of non-entangled bipartitions, where the above-mentioned pairs of frexels are kept coupled together and not split into two different bipartitions. The most external frequency bands of the spectrum (frexels 0 and 7) have a less significant role.

### D. Measuring squeezing pulse-by-pulse

Pulse-by-pulse squeezing measurement in the temporal domain is performed using a fast oscilloscope, which allows us to directly measure quantum fluctuations of the homodyne signal. It requires careful alignment of the homodyne detector, which is achieved by precisely distributing the optical power applied to each photodiode (2.5 mW at most), maximizing the temporal overlap between input pulses, and adjusting the bias voltage applied to the photodiodes. These steps are followed to minimize the CMRR (Common Mode Rejection Ratio) at 156 MHz.

For each squeezing measurement, we collect 2 × 10^{−6} data samples every 50 *µ*s. This corresponds to an average of 7800 pulses per squeezing value. We use the signal from one of the photodiodes as a reference to determine the time window corresponding to a single pulse. Additionally, we verify that the cross correlation function is minimized to confirm that the measured value corresponds unequivocally to the average squeezing of a single pulse. The phase of the LO is varied with a piezoelectric actuator driven much slower than the measurement rate of the oscilloscope. The spectrum analyzer is synchronized with the oscilloscope, and we measure the quadrature variance for different positions of the piezoelectric actuator. We keep the squeezed and antisqueezed quadratures, discard intermediate values, and measure the variance of the shot noise (not squeezed) to set a reference for each measurement.

### E. Design and characterization of the nonlinear waveguide

Periodic poling is used to increase the gain of the SPDC by compensating for the phase mismatch (quasi-phase-matching). The poling period depends on the crystal properties and the wavelength used in the nonlinear process. For this experiment, we used numerical simulations to determine the poling period (Λ = 3.19 *μ*m) and the interaction length (*ℓ* = 1 mm) of the nonlinear waveguide. The company AdvR, Inc. designed and manufactured the waveguide chip. It holds 30 graded-index waveguides embedded in an x-cut PPKTP crystal. The chip is 5.7-mm long to guarantee the quality of the input and output facets of the crystal, and each waveguide is made of a periodically-poled portion *ℓ*, where the parametric process occurs, followed by a non-poled (not phase-matched) portion.

The chip manufacturing process does not guarantee the production of identical waveguides.^{6} Therefore, we characterize them all by measuring the coupling efficiency and the output spectrum, and we select the best ones for our experiment. The waveguides are designed to be single mode at 795 nm but, in practice, we observed the propagation of higher-order spatial modes. We believe this is caused by the gradient of refractive index in the vertical direction of the waveguide. Adjusting the waveguide height ensures that only the fundamental mode at 795 nm can couple inside the waveguide. Additionally, we select the waveguides whose output spatial profile provides good visibility with the LO. Ultimately, we chose a 3-*μ*m-wide waveguide with a coupling efficiency of 57%.

At lower wavelengths, the waveguides are spatially multimode, but only the fundamental mode is optimal for the SPDC. Therefore, spatial mode-matching is required to shape the pump and optimize the coupling of the fundamental mode compared to higher-order modes.

### F. Design of the ultrafast pulse shaper

The pulse shaper setup is based on a folded zero-dispersion line. The spectral shaping is performed with a computer-controlled reflective spatial light modulation (SLM) from Hamamatsu (model X10468). The SLM screen has a 792 × 600 pixel matrix on which the phase masks are displayed. In this configuration, each pixel of the spatial light modulator can control the amplitude and phase of a small frequency band. The setup uses a blazed grating and a cylindrical mirror to spread out the spectrum of the incoming pulse horizontally prior to shaping, and a beam expander to optimize the area of the screen illuminated by the beam and thus increase the resolution of the setup. In our setup, all the phase masks displayed on the SLM include a blazed grating to spatially separate the shaped pulse from the incoming one. The SLM itself has 70% reflectivity. Given the chosen optical elements, the optical complexity of the pulse shaper is 110, which serves as a good indication of the number of modes that can be shaped.^{29}

### G. Wideband balanced homodyne detector

In this section, we summarize the features of a homemade wideband homodyne detector. The electrical circuit uses two Si PIN photodiodes from Hamamatsu with a 300 MHz bandwidth and 91% quantum efficiency at 795 nm. After detection, the photocurrents generated from both photodiodes are subtracted and amplified to generate the homodyne signal that we measured and analyzed with an Agilent 86142B spectrum analyzer and a Teledyne Waverunner oscilloscope. We measured a CMRR of 64 dB at 300 *μ*W. More details on the design of the detector can be found in Ref. 23.

## ACKNOWLEDGMENTS

We thank Bérengère Argence for her support in the design and realization of the homodyne detector. This work was supported by the European Research Council under the Consolidator Grant COQCOoN (Grant No. 820079) and by the Region Ile-de-France in the framework of DIM SIRTEQ.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

T.K. and F.S. contributed equally to this work.

**Tiphaine Kouadou**: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). **F. Sansavini**: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). **M. Ansquer**: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal). **J. Henaff**: Data curation (equal); Investigation (equal); Methodology (equal). **N. Treps**: Conceptualization (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). **V. Parigi**: Conceptualization (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

_{4}