Are there cases where given an ensemble, we cannot blindly replace the ensemble by another that has the same density matrix? It is known that if two ensembles result in the same density matrix, then they give rise to the same observable statistics and the two ensembles are related by a transformation described in this post.
So this leads me to ask a really naive question. Given two ensembles that give rise to the same density matrix, can we always just swap one for the other?
More specifically, suppose $\mathcal{E} = \{(p_{i}, |\psi_{i}\rangle)\}$ and $\mathcal{F} = \{(q_{i}, |\phi_{i}\rangle)\}$ result in the same density matrix.
Is there a situation where a system is modeled with ensemble $\mathcal{E}$, and blindly swapping $\mathcal{E}$ for $\mathcal{F}$ leads to wrong and fallacious results?
Let me explain why I don't think it's obvious to me that equivalent ensembles are interchangeable. For example, suppose we have system $A$ modeled by ensemble $\mathcal{E}$, and then we let it interact with system $B$ to get an entangled (possibly mixed) state. Wouldn't it matter if $A$ was modeled by $\mathcal{E}$ or $\mathcal{F}$ in this case? Or does it end up not mattering?
If it doesn't matter, is there a general proof showing that ensembles that give identical statistics continue to give identical statistics even if you let the system start interacting with other systems?
Also, what if systems $A$ and $B$ are already entangled? They have joint density matrix $\rho^{AB}$, and by taking the partial trace I can get $\rho^{A}$ for system $A$ only. Am I allowed to start attributing ensembles to system $A$ however I wish (even if $A$ and $B$ are entangled)?
 A: Let
$$
| + \rangle = \frac{1}{\sqrt{2}}( |0\rangle + |1\rangle ),
\;\;\;\;\;\;\;\;\;\;\;
| - \rangle = \frac{1}{\sqrt{2}}( |0\rangle  - |1\rangle ).
$$
Let system A consist of two particles prepared in
$$
| \psi_A \rangle = 
\frac{1}{\sqrt{2}}(|0\rangle |0 \rangle + |1 \rangle |1\rangle).
$$
while system B consists of two particles prepared in
$$
| \psi_B \rangle = 
\frac{1}{\sqrt{2}}(|+\rangle |+ \rangle + |- \rangle |-\rangle).
$$
Now suppose that in each case the second particle is passed to some third party, or sent a long way away, and we only have the first particle. The density matrix describing the first particle on its own is constructed by 'tracing out' the second particle. One can imagine that the second particle was measured, but we do not know the outcome of the measurement. The density matrix which results in case A is
$$
\rho_A =  \frac{1}{2} \left( |0\rangle\langle 0| + |1\rangle\langle 1|  \right).
$$
In case B we get
$$
\rho_B =  \frac{1}{2} \left( |+\rangle\langle +| + |-\rangle\langle -|  \right)
= \frac{1}{2} \left( |0\rangle\langle 0| + |1\rangle\langle 1|  \right) 
= \rho_A.
$$
In the absence of access to the information present at the third party, we cannot tell these two versions (A and B) apart. Same density matrix means all observations will give the same outcomes. And yet if we could get access to the third party, then we can tell case A from case B. We could carry out a two-particle measurement, for example, or apply an operation to disentangle the particles. The two cases will yield different states.
The above is readily generalized.
The overall implication is that when a density matrix is used to express information about a system in the absence of further information about other systems with which it may be entangled, then the density matrix is not capturing the whole of the information which may be expressed in the physical world, contained in the whole entangled group of systems. But if we think that there are processes, such as thermodynamically irreversible ones (e.g. the death of a cat, or an entangled partner sent into a black hole, or something like that), which make it impossible to bring the parts back together and undo the entanglement then we may argue that the density matrix of the remaining subsystem is a faithful representation of all the information expressed about it in the physical world.
A: The short answer is, if you consider the whole density matrix, the whole quantum state of the entire world, then no it really doesn't matter which way you got to that density matrix.
However, very often we are looking at a subset of the density matrix, obtained by tracing out certain degrees of freedom of some “environment” or “supersystem” which the subsystem does not have access to. Then it depends:

*

*No experiment “local” to the subsystem (meaning: no observable of the form $\hat A=\hat A_{\text{sub}}\otimes I$) can tell the difference. Indeed this can be viewed as basis for tracing out those degrees of freedom in the first place, start from $\langle A\rangle =\operatorname {Tr} \rho \hat A$ and use this to find an expression $\langle A\rangle=\operatorname {Tr} \rho_{\text{sub}} \hat A_{\text{sub}}.$


*But you can witness the difference if you can correlate with the rest of the supersystem.
A quick example
A good example to keep in mind is the delayed choice quantum eraser experiment. Recall that this consists of

*

*photons pass through a dual slit experiment

*those slits are outfitted with plates that rotate their polarizations 45°, so that the polarization difference is 90° between the two slits and we see a sum of two overlapping bell-shaped curves for the intensity on the screen, and

*we generate the photons as entangled pairs, the initial polarization is entangled with the polarization of another photon which spins round and round in a long coil of fiber-optic cable.

The photons that you see on the screen are displaying a classic density matrix $$\rho_\text{sub}= \frac12|0\rangle\langle0|+\frac12|1\rangle\langle1|,$$ by not showing these interference fringes. If you only have that data and do not correlate with the polarization in the cable, that is the best you can get.
If you correlate with the polarization in the cable, well, first you have to measure the polarisation in the cable. You will put a polaroid filter in some orientation and see if the photon is detected and sort your photon detections on the screen into two groups—was the other photon detected or not?
Now, one orientation of your polaroid filter will reveal zero information about which polarization the original photon had. That will split the photon screen into two identical graphs, both looking like two overlapping bell curves, exactly half of the original intensity profile. We say that you chose to “measure which way” the photon went, but you chose to measure this some long time after the event was logged on the screen.
Or, you rotate this filter 45° and now you get two interference fringe patterns, one looks exactly like the classic two slit interference pattern, the other is kind of staggered... Again, they both sum to the overlapping bell curves! But ignoring the staggered one with a coincidence counter we can say that we “erased the measurement”.
