Quantum entanglement, one of the most fascinating and important features in quantum theory, is widely appreciated as an essential ingredient in quantum computations [1–5]. Simulations of quantum entanglement through optical approaches were investigated both theoretically and experimentally [6–18]. A quantum bit can be represented by a distinct path or space mode of a classical field in an interferometric setup as classical optics analogies [9–18]. However, a *n*-qubit system with 2* ^{n}* basis states must be represented by 2

*distinct paths or modes of a classical field. These simulations are usually not effective due to an exponential increase in required physical resources correlated with the addition of quantum bits [9–19]. The drawback can be traced to a lack of a rigorous tensor product structure of the system [3, 19, 20].*

^{n}In this letter, we present an effective simulation of quantum entanglement of *n* quantum bits by using an analogy of classical fields modulated with pseudorandom phase sequences (PPSs). Based on the properties of PPSs, we proved that the *n* fields modulated with *n* different PPSs constitute a *n*2* ^{n}*-dimensional Hilbert space with a tensor product structure, which differs significantly from those in classical simulations that were executed lacking a tensor product structure [9–18]. By using an optical interferometric setup, PPSs yield not only random measurement results, but also an ensemble model to define the ensemble average and correlation functions [21]. The PPSs, derived from orthogonal pseudorandom sequences, are widely applied to code division multiple access (CDMA) communication technology as a way to distinguish different users [22–24]. A set of pseudorandom sequences is generated by using a shift register guided by a Galois field

*GF*(

*p*) that satisfies orthogonal, closure, and balance properties [23]. In this letter, we utilize a m-sequence of period

*N*− 1(

*N*=

*p*) generated by a primitive polynomial of degree

^{s}*s*over

*GF*(

*p*) and apply it to 4-ary phase shift modulation, a well-known modulation format in wireless and optical communications [22, 24]. Next we generate a PPS set $\text{\Xi}=\left\{{\lambda}^{(1)},{\lambda}^{(2)},\dots {\lambda}^{(N)}\right\}$ over

*GF*(4), where each

*λ*

^{(i)}is a phase sequence with

*N*phase units and time slots: ${\lambda}^{(i)}=[{\lambda}_{1}^{(i)},\text{\hspace{0.17em}}{\lambda}_{2}^{(i)},\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}{\lambda}_{N}^{(i)}]$, while

*λ*

^{(1)}is an all-0 sequence and other sequences can be generated by using following method [23, 25]: (1) given a primitive polynomial of degree

*s*over

*GF*(4), a base sequence of a length 4

*− 1 is generated by the linear feedback shift register; (2) other sequences are obtained by cyclic shifting of the base sequence; (3) by adding zeroes to the sequences, the occurrence of any element equals to 4*

^{s}*− 1; (4) mapping the elements of the sequences to [0, 2π]: 0 mapping 0, 1 mapping π/2, 2 mapping π, and 3 mapping 3π/2.*

^{s}We first consider two orthogonal modes (polarization or transverse modes), $|0\rangle $ and $|1\rangle $, of a classical field. A simulation state can be expressed as a mode superposition: $|\psi \rangle =\alpha |0\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\beta |1\rangle $, where $|\alpha {|}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}|\beta {|}^{2}=1,\text{\hspace{0.17em}}(\alpha ,\beta \in \u2102)$. All of the mode superposition states span a Hilbert space, where we will explore properties associated with this special classical field. By introducing a map $f:\lambda \text{\hspace{0.17em}}\to \text{\hspace{0.17em}}{e}^{i\lambda}$ on the set of Ξ, we obtain a phase sequence set $\text{\Omega}=\left\{{\phi}^{(j)}\right|{\phi}^{(j)}={e}^{i{\lambda}^{(j)}},j=1\dots N\}$, and with which were written a superposed state corresponding to *n*th sequence:

*i*=

*j*, and 0 otherwise. In fact, the map

*f*corresponds to phase modulations of PPSs of Ξ onto the classical field.

A PPS map *f* constitutes a phase ensemble, wherein each phase unit represents a single simulation, and measurement of a physical quantity is a result of ensemble average. Similar to that in quantum mechanics, ensemble average and correlation measurement can be defined [22–24]. In the quadrature demodulation, each code obtained in a sequence unit of a PPS can be considered as a single measurement. The sequence number of the PPS's unit can be used to label the ensemble. Different from the ergodicity hypothesis of quantum mechanics, the ergodicity of PPS is determined and much more efficient.

Given the properties of the PPSs and the Hilbert space, the inner product of any two states and their orthogonal property can be obtained by:

*k*th units of

*λ*

^{(i)}and

*λ*

^{(j)}, respectively. Based on above properties, the classical fields modulated with different PPSs are independent and distinguishable. Figure 1 shows construction pathway of simulation states, generated by unitary transformed from initial states—the mode superposition of classical fields with PPS

*λ*

^{(j)}. Furthermore, a general form of a simulation state can be constructed from $|{\psi}_{n}\rangle $:

Following the pathway in Figure 1, a simulation state $|\Psi \rangle $ is obtained, denoting with a direct product of $|{\psi}_{n}\rangle $:

*N*2

*. Generally a simulation state takes the form:*

^{N}*N*2

*coefficients. It is obvious that the Hilbert simulation space is greater than what is required for simulation of quantum state. To obtain a space the same size as that in quantum mechanics, either restrictions or a proper measurement need to apply [21].*

^{N}PPS provides not only the tensor structure and space needed for quantum state simulation, it also yields the property that an entangled state cannot be expressed in terms of direct product of tensors by using PPS properties and phase ensemble average. In the following we use density matrix to illustrate this feature. We assume that a simple type of simulation state of *N* fields can be expressed as:

*ρ*can be calculated by:

In addition to the fact that a quantum entanglement cannot be expressed in terms of direct tensor product, quantum entanglement also makes a correlation measurement different. The correlation analysis on the simulation states is necessary because the nonlocal correlation with Bell's inequality and equality criterion is the most fundamental property of quantum entanglement. In order to perform the correlation analysis, a correlation measurement $\widehat{P}$ on $|\psi \rangle $ is given:

*α*,

*β*are set to be $1/\sqrt{2}$, yielding $\overline{P}(\theta )=\mathrm{cos}(\theta )$. Furthermore, we generalize $\widehat{P}$ to the case of

*N*fields:

*ρ*:

Key to an effective simulation of quantum entanglement is that the physical resources for the simulation do not increase exponentially with number of particle. In the following we discuss analysis of computation complexity. A simple unitary transformation, *NOT* gate, is used as an example to show computation complexity. Starting with a single field $|{\psi}_{n}\rangle ={e}^{i{\lambda}^{(n)}}({\alpha}_{n}|0\rangle +{\beta}_{n}|1\rangle )$, applying a unitary transformation switching $\widehat{U}:|0\rangle \leftrightarrow |1\rangle $ to decomposes PPS into each phase unit: $\widehat{U}|{\psi}_{n}\rangle \to \left[{e}^{i{\lambda}_{k}^{(n)}}\widehat{U}({\alpha}_{n}|0\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\beta}_{n}|1\rangle )|k=1\dots N\right]={e}^{i{\lambda}^{(n)}}({\alpha}_{n}|1\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\beta}_{n}|0\rangle )$. For each phase unit, its computation is the same as that in quantum computation, therefore computation for *N* phase units equals *N* times of quantum computation of each phase unit. We can extend unitary transformations to simulation states with *N* fields:

*N*phase units ${e}^{i{\lambda}_{k}^{(j)}}$ and

*N*time slots, therefore the required computation is

*N*times that of quantum computation, but 2

*times is unnecessary [3].*

^{N}*Two-particles Bell states*: Consider the case that the modes $|1\rangle $ in the states $|{\psi}_{a}\rangle $ and $|{\psi}_{b}\rangle $ similar to Equation (1) are exchanged by a mode exchanger constituted by mode splitters and combiners [21, 26]. The exchange yields the following states:

The simulation state $|\Psi \rangle $ is obtained:

Then we obtain the results of the fields in the correlation measurement $\overline{P}({\theta}_{a},k)=\mathrm{cos}({\theta}_{a}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\gamma}_{k}^{(a)})$ and $\overline{P}({\theta}_{b},k)=\mathrm{cos}({\theta}_{b}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\gamma}_{k}^{(b)})$, where ${\gamma}_{k}^{(a)},\text{\hspace{0.17em}}{\gamma}_{k}^{(b)}$ are the *k*th units of the RPSs ${\gamma}^{(a)}$ and ${\gamma}^{(b)}$, respectively. Then the correlation function is:

*π*/2, respectively, when Bell's inequality is maximally violated.

Bell state $|{\Psi}^{+}\rangle $ differs from $|{\Psi}^{-}\rangle $ by *π* phase. Similarly, simulation of the Bell state $|{\Psi}^{-}\rangle $ is expressed as $|{\psi}_{a}^{\prime}\rangle ={e}^{i{\lambda}^{(a)}}(|0\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{e}^{i{\gamma}^{(a)}}|1\rangle )/\sqrt{2},|{\psi}_{b}^{\prime}\rangle ={e}^{i{\lambda}^{(b)}}(|0\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{e}^{i({\gamma}^{(b)}+\pi )}|1\rangle )/\sqrt{2}$. By performing the transformation ${\widehat{\sigma}}_{x}:|0\rangle \leftrightarrow |1\rangle $ on $|{\psi}_{b}^{\prime}\rangle $ of the state $|{\Psi}^{\pm}\rangle $, we obtain the simulation of the Bell state $|{\Phi}^{+}\rangle $ expressed as $|{\psi}_{a}^{\prime}\rangle ={e}^{i{\lambda}^{(a)}}(|0\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{e}^{i{\gamma}^{(a)}}|1\rangle )/\sqrt{2},|{\psi}_{b}^{\prime}\rangle ={e}^{i{\lambda}^{(b)}}(|1\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{e}^{i{\gamma}^{(b)}}|0\rangle )/\sqrt{2}$, and of $|{\Phi}^{-}\rangle $ expressed as $|{\psi}_{a}^{\prime}\rangle ={e}^{i{\lambda}^{(a)}}(|0\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{e}^{i{\gamma}^{(a)}}|1\rangle )/\sqrt{2},|{\psi}_{b}^{\prime}\rangle ={e}^{i{\lambda}^{(b)}}(|1\rangle \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{e}^{i({\gamma}^{(b)}+\pi )}|0\rangle )/\sqrt{2}$. Then their correlation functions ${E}_{{\Psi}^{-}}({\theta}_{a},\text{\hspace{0.17em}}{\theta}_{b})=-\mathrm{cos}({\theta}_{a}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\theta}_{b}),\text{\hspace{0.17em}}{E}_{{\Phi}^{\pm}}({\theta}_{a},{\theta}_{b})=\pm \mathrm{cos}({\theta}_{a}-{\theta}_{b})$ are obtained. To substitute the correlation functions into Equation (18), we also obtain the maximal violation of Bell's inequality. The violation of Bell's criterion demonstrates the nonlocal correlation of the two classical fields in our simulation, which results from shared randomness of the PPSs.

*GHZ states*: The nonlocality of the multipartite entangled GHZ states can in principle be manifest in a new criterion and need not be statistical as the violation of Bell inequality [28]. Preparing three states $|{\psi}_{a}\rangle ,\text{\hspace{0.17em}}|{\psi}_{b}\rangle $ and $|{\psi}_{c}\rangle $ similar to Equation (1), and by cyclically exchanging the modes $|1\rangle $ of the states, we obtain the states as following:

*C*= 1/4 is the normalized coefficient. If ${\theta}_{a}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\theta}_{b}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\theta}_{c}=0,\text{\hspace{0.17em}}E({\theta}_{a},\text{\hspace{0.17em}}{\theta}_{b},\text{\hspace{0.17em}}{\theta}_{c})=1$. If ${\theta}_{a}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\theta}_{b}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\theta}_{c}=\pi ,\text{\hspace{0.17em}}E({\theta}_{a},\text{\hspace{0.17em}}{\theta}_{b},\text{\hspace{0.17em}}{\theta}_{c})=-1$. By using GHZ state, the family of simple proofs of Bell's theorem without inequalities can be obtained [29], which is different from the criterion of CHSH inequality. The sign of the correlation function can be also treated as the criterion, such as the negative correlation for nonlocal and the positive correlation for local when ${\theta}_{a}=\pi /3,{\theta}_{b}=\pi /3,{\theta}_{c}=\pi /3$. We obtain that the simulation state in Equation (19) shows the negative correlation. The results are similar to the quantum case of GHZ states.

Furthermore, the simulation of GHZ state could be generalized to the case of *N* particles. By preparing *N* states similar to Equation (1) and cyclically exchanging the modes $|1\rangle $ of the states, the RPSs satisfy ${\gamma}^{(1)}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\gamma}^{(N)}=0$. We obtain the correlation function:

*i*th RPSs at the

*k*th sequence units in the correlation measurement, and $C=1/{2}^{N-1}$ is the normalized coefficient.

Using the same notion, we can obtain simulation results of other quantum entanglement states. It should be pointed out that the phase randomness provided by PPSs is different from the case of quantum mixed states. Quantum mixed states result from decoherence and all coherent superposition items disappear. In contract to the decoherence, some coherent superposition items remain in the simulation state due to the constraints of the RPSs, such as ${\gamma}^{(a)}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\gamma}^{(b)}=0,\text{\hspace{0.17em}}{\gamma}^{(a)}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\gamma}^{(b)}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\gamma}^{(c)}=0$ for the simulation of Bell states and GHZ state, respectively. These remaining items make it possible to simulate quantum entangled pure states.

As shown in the above examples, we utilized the properties of PPSs to label classical fields that are even overlapped in the same space and time. In simulation of entangled states, the resources required are the PPSs instead of classical field modes. It means that the amount of PPSs grows linearly with the number of quantum particles. According to the m-sequence theory, the number of PPSs in the set equals to the length of sequences, which means that the time resource (the length of sequence) required also grows linearly with the number of the particles.

In conclusion, a novel simulation method for quantum entanglement is presented, with its mathematical expressions and physical meanings identical to those in quantum mechanics. In the framework of quantum mechanics, the overall phase of a wave function can be ignored, as it has no contribution to the probability distribution. However, quantum entanglement must be related to two or more spatially separable quantum particles. By introducing a phase factor to superposed states with PPS properties, we conclude that quantum entanglement can be efficiently simulated by using a classical field modulated with PPSs. The research on this simulation not only provides useful insights into fundamental features of quantum entanglement but also yields new insights into quantum computation.