Based on the discussions in response to the answers provided, it seems to me that the issue to be resolved is a simple matter of looking at the unitary process represented in the specific case where interference is reproduced.
I am going to use a mathematical approach. Without math the discussion would just be too cumbersome and I think it is the lack of mathematics that causes the confusion.
We start just after the slit before the nonlinear crystal. There the state is (ignoring normalization)
$$ |\psi\rangle = |s_1\rangle + |s_2\rangle , $$
where $s_1$ and $s_2$ represent the two slits. Because they are in superposition, they would produce the well-known interference associated with a double slit experiment. However, when we then send it through the nonlinear crystal, we get
$$ |\psi\rangle = |H\rangle_1|V\rangle_1 + |H\rangle_2|V\rangle_2 , $$
where $H$ and $V$ represent horizontal and vertical polarization of the down-converted photons, and the subscripts indicate the slits. These two terms are now separated and then recombined in different ways for the different detections. The ones we are interested in are those that would be able to reproduce interference. Let's assume the interference observed at $D_0$ is observed in the horizontal polarization. So the vertical polarization state are manipulated to give the associated measurements.
For the interference require the sums and difference (still ignoring normalization):
$$ |+\rangle = |V\rangle_1 + |V\rangle_2 , $$
$$ |-\rangle = |V\rangle_1 - |V\rangle_2 , $$
which means that
$$ |V\rangle_1 = |+\rangle + |-\rangle , $$
$$ |V\rangle_2 = |+\rangle - |-\rangle . $$
Now I can replace the last two equations into the expression for the state. Then it becomes, with some rearrangement of terms
$$ |\psi\rangle = (|H\rangle_1 + |H\rangle_2)|+\rangle
+ (|H\rangle_1 - |H\rangle_2)|-\rangle . $$
Here the states $|+\rangle$ and $|-\rangle$ are measured at $D_1$ and $D_2$, respectively, in coincidence with the measurement at $D_0$, which would then respectively measure $|H\rangle_1 + |H\rangle_2$ and $|H\rangle_1 - |H\rangle_2$. Clearly, they produce the interference that we are looking for. They process is simply because of unitary transformation applied to the one system in an entangled state which then becomes a unitary transformation on the other system.
Discussion: The nonlinear parametric down-conversion process which converts the incoming photon after the slit into a pair of entangled photons is a unitary process. What this means is that a pure state (which we have assumed for the photon after the slit) remains a pure state even though it changed. Therefore, it remains coherent. As a result it can maintain momentum and energy conservation that are required for such a parametric process. It is in fact these conservation properties that give rise to entanglement in the wave vectors and frequencies.
So when we select a single photon after the slit, which is then converted into two photons by the BBO, it contains the full information in the bi-photon wave functions that allows us, through the appropriate manipulation of one of the two photons, to produce the interference pattern as a single-photon wave function at $D_0$. However, since it is just one photon, we can only observe one dot. To see the complete interference pattern, we need to repeat the experiment to build up the statistics from which the interference pattern emerges.
That is the only reason for the ensemble. Every photon produces its own interference. (As Dirac phrased it: every photon only interferes with itself.) But to see this pattern we need lots of photons to build up the statistics.