Discrete-event simulation of an extended Einstein-Podolsky-Rosen-Bohm experiment

We use discrete-event simulation to construct a subquantum model that can reproduce the quantum-theoretical prediction for the statistics of data produced by the Einstein-Podolsky-Rosen-Bohm experiment and an extension thereof. This model satisfies Einstein's criterion of locality and generates data in an event-by-event and cause-and-effect manner. We show that quantum theory can describe the statistics of the simulation data for a certain range of model parameters only.


I. INTRODUCTION
The Einstein-Podolsky-Rosen thought experiment was introduced to question the completeness of quantum theory [1], "completeness" being defined in Ref. 1. Bohm proposed a modified version that employs the spins instead of coordinates and momenta of a two-particle system [2], and is experimentally realizable [3][4][5][6][7][8][9][10]. A key issue in the foundations of physics is whether there exist "local realist" models that yield the statistical results of the quantum-theoretical description of the Einstein-Podolsky-Rosen-Bohm (EPRB) experiment.
In this paper, we take, as operational definition of a local realist model, any model for which • Events are the basic building blocks of any DES model, just as points are the basic building blocks of Euclidean geometry. In DES an event is a defined concept, represented by a model variable taking a particular value (e.g. a bit changing from zero into one) at a specific point in time. In contrast to quantum theory, there is no need to invoke the elusive wave function collapse to "explain" how quantum theory may eventually be reconciled with the fact that a measurement yields a definite yes/no answer, or to appeal to Born's rule.
• It may be difficult to analyze the behavior of a DES model by means of differential equations, probability theory, or other mathematical techniques of theoretical physics. Of course, we may use e.g. probability or quantum theory to model the statistics of the data produced by a DES.
• It is not practicable to perform a DES without using a digital computer. A digital computer itself is a physical device that changes its internal state (all the bits of the CPU and memory) in a discrete, step-by-step (clock cycle) manner. Therefore, a DES algorithm running on a digital computer (which we assume is error-corrected and operating flawlessly) can be viewed as a metaphor for an idealized experiment on a physical device (the digital computer) [54]. All aspects of such an experiment are under the control of a programmer. In the context of EPRB experiments, this means that any loophole [55] can be opened or closed at the discretion of the programmer. For instance, the so-called contextuality loophole, which is impossible to avoid in a laboratory EPRB experiment [49], can be trivially closed in a DES (see below).
• The outcomes of genuine laboratory experiments are subject to unknown influences but in a DES on a digital computer (operating flawlessly), there are no such influences. If there were, it would not be possible to exactly reproduce the results of a DES. Therefore, DES is "the experiment" to confront a theory with facts obtained under the same premises as those on which the theory is based.
• Although a DES algorithm changes the state of a physical device (the digital computer), the events and variables in a DES are only metaphors for the "real" detector click, etc. On the other hand, once it has been established that a DES of a subquantum model yields the correct results, one could build a macroscopic mechanical device that performs exactly the same as DES.
• DES on a digital computer complies with the notion of realism, meaning that at any time during the DES, the internal state of the digital computer is known exactly, all variables of the simulation model taking definite values. Of course, we can always "hide" an algorithm and data on purpose. For instance, we can do this to create the illusion that the "visible" data is unpredictable (a standard technique to generate pseudo-random numbers).
• In a digital computer, there are no signals traveling faster than light. Therefore, on the most basic level, the internal operation of a digital computer complies with Einstein's notion of local causality. However, there is nothing that prevents us from performing an acausal analysis of the data. For instance, if we generate and store a sequence of numbers and then wish to compute the sum, we may do this by summing the numbers in the reverse order of how they were generated. This trivial example shows that one has to distinguish between the generation of the raw data and the processing of this data. For the purpose of constructing a local realist DES model, it is essential that the process that generates the raw data complies with Einstein's notion of local causality. It is not permitted to accumulate data, perform e.g. a discrete Fourier transform or compute acausal correlations, and use the results to describe a quantum physics experiment. While both these techniques are very useful for a wide variety of data processing tasks [56], they are "forbidden" in a DES of a subquantum model.
• Consistency of the DES methodology demands that a subquantum model for, say, a beam splitter, must be re-used, without modification, for all experiments in which this beam splitter is used. Our DES approach seems to satisfy this requirement of consistency, at least for a vast number of fundamental quantum-physics experiments with photons and neutrons [51]. Our motivation for considering both the EPRB and extended EPRB (EEPRB) experiments is to scrutinize the consistency of the DES approach for this category of experiments.
Finally, to head off misunderstandings, the DES models that we construct do not, in any way, make use of the quantum-theoretical predictions for the statistics of the data. Instead, a DES builds up these statistics by an event-by-event, cause-and-effect, Einsteinlocal process. Under appropriate conditions, these statistics can be described by quantum theory, while in other cases they cannot (see below).

III. EXTENDED EINSTEIN-PODOLSKY-ROSEN-BOHM EXPERIMENT: THEORY
A. Thought experiment Figure 1 shows the layout of an extended Einstein-Podolsky-Rosen-Bohm (EEPRB) experiment with spin-1/2 particles [26]. A source is emitting a pair of particles in two spatially separated directions towards beam splitters BS1 and BS2. In this idealized experiment, all beam splitters are assumed to be identical, performing selective (filtering) measurements [57,58]. Selective measurements allow us to attach an attribute with definite value (e.g. the direction of the magnetic moment or of the polarization) to the particle. For instance, assuming that BS1, BS3 and BS4 perform ideal selective measurements, a particle leaving BS1 along path S 1 = +1 (S 1 = −1) will always leave BS3 (BS4) along path S 3 = +1 (S 3 = −1) if c = a. In case of selective measurements, we can attach an attribute to the particle, the value of this attribute being given by S 1 . We use the same procedure for attaching attributes to particles leaving the other beam splitters.
Particles leaving BS3,...,BS6 are registered by detectors. All detectors are assumed to be identical and to have a 100% detection efficiency (we relax this assumption later). The binary variables x i, j = 0, 1 for i = 1, . . . , 4 and j = 1, 2 (see Fig. 1) indicate which of the four detectors at the left ( j = 1) and right ( j = 2) fire. For each pair of emitted particles, exactly one of the detectors on the left and exactly one of the detectors on the right side of the source will register a particle. This implies that for j = 1, 2, only one of x 1, j , x 2, j , x 3, j , and x 4, j can be nonzero. The four new variables defined by all take values +1 or −1 (see Fig. 1). Clearly, S 1 and S 3 (S 2 and S 4 ) encode, in a unique manner, the path that the left (right) going particle took. In the following, we use the S's to formulate the DES model. For practical reasons, most laboratory EPRB experiments are carried out with photons [3][4][5][6][7][8][9][10], the polarization of the photons playing the role of the spin. From a quantum-theoretical viewpoint, there is no loss of generality in doing so because mathematically, the description of the photon polarization is in terms of Pauli-spin matrix algebra. In the following, to keep the discussion concise and concrete, we only focus on (E)EPRB experiments that employ the photon polarization as the "quantum system" of interest.
A source is emitting a pair of photons in two spatially separated directions towards beam splitters BS1 and BS2. BS1 sends one photon of the pair to either BS3 or BS4. BS2 sends the other photon of the pair to either BS5 or BS6. In an ideal model, all beam splitters are identical. Each beam splitter represents a combination of wave plates, an electro-optical modulator (EOM), and a polarizing beam splitter (PBS), see for example

B. Classical Electrodynamics
It is instructive to first consider the case in which the detector signal is linearly proportional to the intensity of the impinging light. This case is covered by classical optics, described by Maxwell's theory of electrodynamics.
According to empirical evidence, the intensity of light passing through a polarizer is given by Malus' law I = I 0 cos 2 (φ − ψ), where φ is the polarization of the light beam and ψ is the rotation of the polarizer, both relative to a laboratory frame of reference. I 0 is the intensity of the incident light.
We assume that the source emits "special" randomly polarized light towards BS1 and BS2, special in the sense that the difference between the polarizations of the two beams φ 0 is fixed in time. Then, using Malus' law for BS1 and BS2, the correlated intensity for one particular, random realization of the polarization angle φ is given by where I 0 denotes the light intensity of a single beam. Integrating over all polarizations φ with a uniform density 1/2π yields Repeated use of Malus' law and exploiting the fact that the polarizations of the two beams leaving a beam splitter are orthogonal, the correlated intensity of the different beams is then given by The moments of the correlated intensity Eq. (5) are For φ 0 = π/2 (orthogonally polarized beams) and I 0 = 1, the explicit expressions for the two-S correlations are The factor 1/2 which appears in Eq. (4) and in four of the six second moments K i j is characteristic of the correlation of two light intensities. Here and in the following, we use the hat on top of the symbols to emphasize that the expressions have been obtained from a theoretical model.

C. Quantum theory
Classical electrodynamics describes the intensity of light and does not discriminate between individual events. In contrast, quantum theory can be used to describe the statistics of events, particularly in cases, such as the EPRB experiment, where detectors can discriminate between them.
For a pair of photons whose polarizations are described by the singlet state, Appendix A shows that the joint probability to observe one photon in the path labeled by S 1 and the other one in the path labeled by S 2 is given by [59] where Z denotes a valid proposition that represents all conditions under which the experiment is performed with the exception of a and b. Note that for φ 0 = π/2, Eq. (4) differs from Eq. (8) through the factor 1/2 only. In Appendix A, we show that the joint probability to observe one photon in the path labeled by (S 1 , S 3 ) and the other in the path labeled by (S 2 , S 4 ) is given by As already mentioned in the introduction, a subquantum model for the EEPRB experiment must not make use of Eq. (9) to generate the quadruples (S 1 , S 2 , S 3 , S 4 ).
Note that the only nontrivial difference between Eq. (5) and Eq. (9) is that in the former case, the absolute value of the prefactor of the S 1 S 2 term never exceeds 1/2 whereas in the latter case, it is equal to minus one. The expressions of the second moments of Eq. (9) read which are not all equal to the corresponding expressions of the K's, see Eq. (7). From Eq. (9) it follows that Clearly, in order to have a subquantum model generate data that agrees either with Maxwell's theory Eq. (5) or quantum theory Eq. (9), we only have to construct a subquantum model in which we can control the prefactor of the S 1 S 2 term in Eq. (5). Thinking of light as a collection of photons, in sections III D and IV, we explain how this control naturally results from the simple fact that we have to classify individual events as photons or something else, whereas in the "classical" case this classification is not an issue. At this point, it should be mentioned that within the context of the classical and quantum theory of light, changing the prefactor 1/2 of the S 1 S 2 term in Eq. (5) is a subtle issue, intimately related to the amount of second-order coherence one can observe by measuring either intensities or by counting clicks of a detector [60]. A discussion of this important issue is out of the scope of this paper and we refer the reader who is interested in these aspects to the in-depth analysis given in Ref. 60. In our paper, Eqs. (5) and (9) are only used to provide the classical/quantum results which any valid subquantum model for the (E)EPRB experiment has to reproduce.
An important feature of this EEPRB experiment is that all the correlations that are required to test for violations of Bell/Clauser-Horne-Shimony-Holt (CHSH) inequalities [12,61] are obtained in a single run (instead of three/four runs) of the experiment [26]. The EEPRB experiment does not suffer from the contextuality loophole [49]. As 0 ≤ P(S 1 , S 2 , S 3 , S 4 |a, b, c, d, Z) ≤ 1, it follows directly that all Bell-type inequalities, including all variants of the CHSH inequality, can never be violated [62]. This is easily seen by evaluating the sum for example. In words, the quantum-theoretical description of the EEPRB experiment predicts that all Bell/CHSH inequalities are satisfied, in stark contrast to the case in which the correlations that enter the Bell/CHSH inequalities are computed from the quantum-theoretical description of the EPRB experiment.

D. Practical realization: photon identification problem
The exposition in subsection III A assumes that each emitted pair of particles triggers exactly two detectors, namely only one of the four detectors at OS1 and one of the four detectors at OS2. In a laboratory experiment with Stern-Gerlach magnets and magnetic billiard balls, this assumption may hold true. However, it is not at all evident to have a source which only creates correlated pairs of elementary particles such as photons which upon hitting a detector, will trigger exactly one detector at OS1 and exactly one detector at OS2.
With the exception of two experiments [9,10], EPRB experiments with photons use time coincidence to identify photon pairs [3][4][5][6][7][8]. The two EPRB experiments [9,10] that do not rely on time coincidence employ local, adjustable voltage thresholds to identify photons. This procedure is mathematically equivalent to attaching a local time tag to each particle, or to using time coincidence [63]. Therefore, in the following, we only discuss the subquantum model that uses local time tags as the vehicle for identifying pairs of photons. The modifications required to deal with voltage thresholds are trivial [63].
As explained above, a minimal theoretical model of a laboratory (E)EPRB experiment with photons should include a procedure to identify (pairs of) photons. Specifically, not including the data by which the photons and/or pairs are identified opens the socalled photon identification loophole [63]. By design, the EPRB experiments that claim to be loophole free [8][9][10] all suffer from this loophole.
As shown in Fig. 1, observation stations OS1 and OS2 are equipped with local clocks C1 and C2, respectively. The time t 1 (t 2 ) at which a detector in OS1 (OS2) fires is read off from the local clock C1 (C2). The clocks C1 and C2 are synchronized before the source starts to emit pairs of particles and, being ideal clocks, remain synchronized for the duration of the experiment. Similarly, the frames of reference of OS1 and OS2 are aligned before the source starts to emit pairs of particles and do not change afterwards.
Concretizing the aim of this paper, in the next section, we describe a local realist model of the (E)EPRB experiment that reproduces the statistical predictions of quantum theory given by Eqs. (8) and (9). In formulating the DES model, we call the agents that carry the information from the source to the observation stations "photons" and use the language of optics.

IV. SUBQUANTUM MODEL
Before describing all the components of the subquantum model, we recall the basic strategy that we adopt in constructing such a model. As quantum theory describes the most ideal version of the (E)EPRB experiment and as our aim is to show that the subquantum model reproduces the results of the former, we construct a DES of the most ideal version of the (E)EPRB experiment. The DES model for the EEPRB experiment that we describe next contains the ideal implementation of EPRB laboratory experiments.
In concert with our general strategy to set up the subquantum model, we assume that the source emits pairs of photons only and this at regular time intervals ∆. The time at which the nth pair is emitted is given by T n = n∆. The time it takes for a photon to travel from the source to BS1 or BS2 is assumed to be constant and the same for all photons traveling to OS1 and OS2. We denote this time of flight by T TOF . Similarly, the time of flight from BS1 to BS3 or BS4 (BS2 to BS5 or BS6) and the time of flight from BS3, BS4, BS5 or BS6 to the corresponding detectors are denoted by T TOF and T TOF , respectively. The total time of flight is then T TOF = T TOF + T TOF + T TOF .
In the DES model, the polarization of the photon traveling to BS1 (BS2) is represented by a two-dimensional unit vector The angle φ is chosen to be uniformly random from the interval [0, 2π). As x T 1 · x 2 = 0, the polarizations of the photons of each pair are orthogonal, that is, they are maximally anticorrelated and randomly distributed over the unit circle.
In the following, we specify the DES rules for BS1 only. The rules for the other beam splitters are identical and are obtained by a simple change of symbols. In the DES model, the operation of beam splitter BS1 is defined by the rules where 0 < r < 1 denotes a uniform pseudo-random number. Here and in the following, it is implicitly understood that a new instance of the pseudo-random number r is generated with each invocation of an equation in which r appears. The unit vector x 1 denotes the polarization of the photon leaving BS1. It is not difficult to see that the model defined by Eq. (11) generates S 1 = +1 (S 1 = −1) events with a relative frequency given by cos 2 (φ − a) (sin 2 (φ − a)), that is, Eq. (11) produces data that is in concert with Malus' law if the polarization of the incident photon is constant in time.
Optical components such as wave plates and EOMs contain birefringent material which changes the polarization by retarding (or delaying) one component of the polarization with respect to its orthogonal component. In the DES, this retardation effect is accounted for by assuming that as a photon passes through a beam splitter, it may suffer from a time delay which may depend on the direction of the beam splitter relative to the polarization of the photon.
Obviously, the law of retardation in the subquantum model cannot be derived from Maxwell's theory or quantum theory. We can only find the subquantum law of retardation by trial and error. Fortunately, from earlier work we already know the subquantum law of retardation for the EPRB experiment [64,65] and we only need to extend this law slightly to have the DES reproduce the quantum-theoretical results for both the EPRB and EEPRB experiment. Specifically, for BS1, the two DES rules for the subquantum law of retardation read where x 1 (or equivalently φ ) is the polarization of the incoming photon, 0 < r < 1 is another uniform pseudo-random number, T max is an adjustable parameter specifying maximum retardation, and α > 0 is an adjustable parameter controlling the dependence of the retardation on the difference between the photon polarization x 1 and the orientation of the beam splitter a.
The new features are the last factor in Eq. (12), β > 0 being an adjustable parameter, and the rule u 1 ← x 1 which updates the two-dimensional vector u 1 . The initial value of u 1 can be any vector that has a norm less than or equal to one. This vector is attached to the beam splitter and may be thought of as representing (on a subquantum level) the electrical polarization of the material [51].
The purpose of the factor |(1 − x 1 · u 1 )/2| β in Eq. (12) is to turn off the generation of random retardation times if the polarization of the incoming photons is constant. To see how this works, first consider the case that the polarization of the incoming photons is constant, say x 1 = x 1 . Then, after the first photon has passed by, u 1 = x 1 and |(1 − x 1 · u 1 )/2| β = 0 for all photons that follow. Next, assume that the polarization of the photons entering BS1 is randomly distributed over the unit circle. Then, because x 1 and u 1 (which is equal to the polarization x 1 of the previous photon) are independent, |(1 − x 1 · u 1 )/2| β is just a random variable in [0, 1] multiplying τ EPRB . The idea to store and use the value of the polarization of the previous photon has also been used to reproduce, by DES, the quantum-theoretical results for a large variety of single-photon and single-neutron experiments [50,51]. The capability of the subquantum model to disable the generation of random retardation times if the polarization of the incoming photons is constant is essential for reproducing the quantum-theoretical results of both the EPRB and EEPRB experiments with the same subquantum model.
At this point, it may be of interest to mention that on the basis of the statistics (i.e. averages and correlations) only, it is not possible to make statements about the uniqueness of the subquantum law of retardation. As a matter of fact, in the case at hand, replacing Eq. (12) by the rules works equally well. Equation (14) defines a deterministic learning machine (DLM) [50,51], which learns, event-by-event, the time average of the polarizations x 1 carried by the photons. The speed and accuracy by which u 1 approaches the time average of the x 1 's is controlled by the parameter 0 < γ < 1 [51]. The order in which Eqs. (13) and (14) are executed is irrelevant. The DLM defined by Eq. (14) is the same as the one that has been used to reproduce, by DES, the quantum-theoretical results for a large variety of single-photon and single-neutron experiments [50,51]. If the polarizations of the incoming photons are constant, say x 1 , and a certain number (depending on γ) of photons has passed by, we have u 1 ≈ x 1 and |(1 − u 1 · u 1 )/2| β ≈ 0. If the polarizations of the photons entering BS1 are randomly distributed over the unit circle, u 1 → 0 and |(1 − u 1 · u 1 )/2| β ≈ 1. Then, just as in the case of Eq. (12), the factor |(1 − u 1 · u 1 )/2| β in Eq. (13) is used to turn off the generation of random retardation times if the polarizations of the incoming photons are constant.
In our idealized experiment, all detectors are assumed to be identical and to have a 100% detection efficiency. After a photon has passed BS3 or BS4 (BS5 or BS6), it may trigger one and only one detector in OS1 (OS2), symbolized by one of x 1,1 , x 2,1 , x 3,1 , or x 4,1 (x 1,2 , x 2,2 , x 3,2 , or x 4,2 ) being equal to one and the other ones being equal to zero. The time t 1 (t 2 ) at which a detector in OS1 (OS2) fires is read off from the local clock C1 (C2). These local clocks C1 and C2 are synchronized before the source starts to emit pairs of photons and, being ideal clocks, remain synchronized for the duration of the experiment. For the nth emitted pair, the arrival times are given by where we have attached the subscript n to keep track of which pair of the emitted pairs we are dealing with. For each pairemission event n = 1, . . . , N, Eqs. (11) -(16) generate the data (S 1,n , S 3,n ,t 1,n ) and (S 2,n , S 4,n ,t 2,n ). Note that OS1 and OS2 only share the angle of polarization φ characterizing the pair of photons, nothing else. In each triple in (S i,n , S i+2,n ,t i,n ), we replace the time variable by a (local) binary variable w i,n to indicate whether a detection event is classified as a photon (w i,n = 1) or not (w i,n = 0). Specifically, the rule to decide whether a detection event corresponds to the observation of a photon or of something else is given by where W is the time window (an adjustable parameter). We emphasize that the decision process defined by Eq. (17) collected by OS1 and OS2, respectively. From the data sets S 1 and S 2 , we compute the single-and two-particle averages ∑ N n=1 w 1,n w 2,n S i,n S j,n ∑ N n=1 w 1,n w 2,n , (i, j) = (1, 2), (1,3), (1,4), (2,3), (2,4), (3,4) , (19) without (K's) and with (E's) the photon identification process in place.  Below, we specify the DES parameters that have been used, discuss the data shown in Figs. 2 and 3, and provide additional information about simulation data that we do not show.
• The number of emitted pairs is N = 1000000 per setting (a, b, c, d).
• The maximum retardation time was chosen to be T max = 5000 (dimensionless units). Pairs of particles are emitted with a time interval ∆ > 2T max . In line with our strategy to perform an ideal experiment, this choice eliminates the possibility of misidentifying pairs and also ensures that at each instant of time, there is only one photon in transit to OS1 and only one other photon en route to OS2. For T TOF we can take any nonnegative value. In fact, from Eq. (17) it follows that the actual values of ∆ and T TOF do not enter in the DES algorithm.
• In Fig. 2(a), we show the DES results for the case without photon identification (that is if W > T max or α = β = 0). Then the DES reproduces the results of Maxwell's theory, by an event-by-event process [66]. Specifically, if the polarization of the incoming photon is constant, the DES model of the beam splitter itself generates data according to Malus' law. The DES model of the EEPRB experiment with randomly polarized light produces data that, ignoring statistical fluctuations, is in excellent agreement with Eq. (7), see Fig. 2(a).
• Using a local time window W = 1 (dimensionless units) and for α = 4 and β = 1/2, the event-by-event process yields results (see Fig. 2(b)) for the E's which are in excellent agreement with quantum theory (data for the first, third and fourth moments are not shown). The ratio of identified photon pairs to emitted pairs depends on a and b and varies between approximately 11% for a = b and 0.1% for |a − b| = π/4. For W = 8, this ratio changes to approximately 18% for a = b and 0.8% for |a − b| = π/4, while the agreement with quantum theory is still very good (data not shown).
• The data obtained by identifying photons using time coincidence instead of local time windows are almost the same and are therefore not shown. In the limit W /T max → 0, N → ∞, and α = 4, it has been proven analytically that the DES model of the EPRB experiment yields the correlation E 12 = − cos 2(a − b) exactly [65].
• Using the same pseudo-random sequence for each choice of settings renders the DES compliant with the notion of a counterfactually definite theory. In this case, the DES results (data not shown) are, for all practical purposes, the same as those obtained by using different pseudo-random sequences for each choice of settings. Operating in this mode, the subquantum model of the EPRB experiment does not suffer from the contextuality loophole [49] nor of any other known loopholes [55]. This confirms the conclusion of an earlier work [54], which adopted a different approach to realize counterfactually-definite compliant simulations. The demonstration that there exist both counterfactually-definite and non-counterfactually-definite compliant computer models for the EPRB experiments that produce results in complete agreement with those of quantum theory implies that, for the case of EPRB experiments, counterfactual definiteness is not incompatible with quantum physics [63].
• In the DES, it is trivial to account for the detection efficiency 0 ≤ η ≤ 1. For each detection event, we generate a pseudorandom number r and remove the detection event from the data set if r > η. We find that the only effect of reducing η is to increase the statistical fluctuations (data not shown). The agreement with quantum theory is not affected.
• In Fig. 3, we show the DES results for the case in which the source emits photons with the same polarization chosen randomly. In this case, the DES reproduces the results of Maxwell's theory ( Fig. 3(a)). A corresponding quantumtheoretical result does not exist, see appendix B. However, the DES data are in excellent agreement with a non-quantum probabilistic theory in which the factor 1 − S 1 S 2 cos 2(a − b) in Eq. (9) is replaced by 1 + S 1 S 2 cos 2(a − b) , see Fig. 3(b).
• Replacing the rules Eq. (12) by the rules Eqs. (13) and (14), and repeating the DES with the same value of α and β also yields data (not shown) that are in excellent agreement with the quantum-theoretical description, for γ in the range [0.1, 0.98].
• Our DES model also reproduces the theoretical results (see Appendix C) if the two photons of each pair have a fixed polarization (results not shown). In Maxwell's theory, this case is described by light beams with fixed polarizations. In quantum theory, this case is described by an (uncorrelated) product state.

VI. DISCUSSION AND SUMMARY
Laboratory EPRB experiments unavoidably require a procedure to classify a detection event as corresponding to a photon or as something else. Independent of the precise nature of this procedure (voltage threshold, local time window, time coincidence, etc.), any model that aims at describing an EPRB experiment should, from the start, account for this procedure by introducing additional variables into the description. In contrast, Bell's model, while charmingly simple, does not account for an essential aspect of laboratory EPRB experiments, namely the classification of detection events in terms of photons or something else. Consequently, any subquantum model that aims at reproducing the results of quantum theory for the (E)EPRB experiment should have features that are not included in Bell's model. As a matter of fact, a quick glance at how the data of laboratory EPRB experiments are being processed reveals that it is the photon identification process which is lacking in Bell's model. Including this process implies that correlations between events are calculated only from subsets of the data, in which case Bell's theorem does not apply. Using only subsets of the data, there is only the constraint that the correlation should, in absolute value, be less than or equal to one. Apart from that "almost everything" is possible [20,[67][68][69], including a subquantum model that, in the appropriate limit, yields the correlation of the singlet state [64,65].
Clearly, on the basis of the statistical data alone, it is not possible to reject subquantum models of the EPRB and EEPRB experiments presented in this paper. The relevant question is how a laboratory experiment can rule out or confirm that (i) a subquantum level description is possible and/or (ii) the rules by which the DES model of the beam splitter operates provide a reasonable description.
Regarding (i): If we consider it as irrelevant to ask what kind of process gives rise to the statistics of events, it seems very difficult to beat quantum theory in terms of descriptive power [70]. Therefore, it is clear that addressing (i) requires the analysis of the data on the level of individual events, without being biased by what quantum theory predicts for the statistics.
Regarding (ii): The DES model defined by Eqs. (11) and (12) or Eqs. (11), (13) and (14) produces data in concert with Malus' law, i.e. with the experiment, and therefore seems solid. The additional feature of the DES model (which allows us to reproduce the statistics of the EPRB and EEPRB experiments as given by quantum theory) is the subquantum law of retardation, defined by Eq. (12) or Eqs. (13) and (14). At first sight, there is no experimental support for such a law. However, let us look at the EPRB experiment from a slightly different perspective, namely, as a setup to characterize the response of the observation stations to very feeble, randomly anticorrelated light. Then, our DES data and also the analysis of experimental data [65] support the hypothesis that the EPRB experiment demonstrates that the statistics of the detection events that have been classified as photons depend on the settings.
What if we try to measure the retardation by an experiment that uses feeble light with fixed polarization? As explained in section IV, in our subquantum model the retardation time does not depend on the setting of the beam splitter if the polarization is constant in time. Only if the polarization is not fixed in time, our subquantum model yields retardation times that depend on the setting of the beam splitter. This is precisely what the EPRB experiment does: through the correlation, it provides information about the retardation as a function of the setting of the beam splitter. Therefore, to rule out the subquantum law of retardation used in our DES, it is necessary to perform both the experiments with feeble, randomly polarized light and with feeble light of fixed polarization.
Summarizing, we have proposed a subquantum model which satisfies Einstein's criterion of locality and which generates, event-by-event, data that agrees with the quantum-theoretical description of the Einstein-Podolsky-Rosen-Bohm and the extended Einstein-Podolsky-Rosen-Bohm experiments. This demonstration does not build on the traditional methods of theoretical physics but instead uses a digital computer and a discrete-event simulation as a metaphor for idealized, realizable laboratory experiments. For reference, we briefly review the quantum-theoretical description of the EPRB and EEPRB experiment. In this appendix, for the sake of generality, we first consider magnetic spin-1/2 particles passing through Stern-Gerlach magnets.
In the context of (E)EPRB experiments, the case of interest is a system of two spins in the singlet state, described by the density matrix A selective measurement on the spin-1/2 particle is described by the operator [58] M(S, σ σ σ , projecting a state of the spin-1/2 system onto the eigenstate of σ σ σ · x with eigenvalue S = ±1. The probabilities to observe the outcomes S 1 and S 2 in an EPRB experiment are given by [58] P respectively. For two spin-1/2 particles in the singlet state Eq. (A1), we have P(S 1 |a) = σ σ σ 1 · a = 1/2 and P(S 2 |b) = σ σ σ 2 · b = 1/2, The probability to observe the joint event (S 1 , S 2 ) is given by [58] P Performing the matrix multiplications and calculating the trace we obtain P(S 1 , S 2 , S 3 , S 4 |a, b, c, d, Z) = 1 16 1 − S 1 S 2 a · b 1 + S 1 S 3 a · c 1 + S 2 S 4 b · d .
The derivation of the quantum-theoretical description of an experiment with photon polarization instead of magnetic spin-1/2 particles is not much different, for details see Ref. 71. The upshot is that we only have to replace a · b by cos 2(a − b) etc. Thus, in the case of an EEPRB experiment that uses the polarization of the photons, the probability to observe the joint event (S 1 , S 2 , S 3 , S 4 ) is given by P(S 1 , S 2 , S 3 , S 4 |a, b, c, d, Z) = 1 16 1 − S 1 S 2 cos 2(a − b) 1 + S 1 S 3 cos 2(a − c) 1 + S 2 S 4 cos 2(b − d) .

(A8)
Appendix B: A limitation of quantum theory for two spin-1/2 particles If the random polarizations of particles that enter BS1 and BS2 are the same instead of orthogonal, the DES generates data which, within the usual statistical fluctuations, is characterized by E 1 = E 2 = 0, and E 12 = + cos 2(a − b). In this appendix, we prove that a two-particle system with single-particle averages E 1 = E 2 = 0, and pair correlation E 12 = + cos 2(a − b) cannot be described by the quantum theory of two spin-1/2 particles. As in Appendix A, we consider the general case of two spin-1/2 particles and deal with the case of photon polarization at the end.