How is the phase of a BEC changing when some particles are being removed? I have some trouble understanding what happens to the phase of a BEC when some particles are removed. The motivation of the question is the experiment of observing interference between two BECs. In my current picture, a BEC is a state $\Phi(r)e^{i\varphi}$, where $\Phi(r)$ is the N-fold tensorproduct of the single particle ground state. This state "lives" in the corresponding Hilbert space that is the N-fold tensor product of the single particle Hilbert space. As the tensor product is linear, the phase $\varphi$ would then just be the sum of all single-particle phases. Now this means that if we remove a single particle from the BEC, its phase is expected to jump between $[0,2\pi)$, as each particle can make a contribution $\in[0,2\pi)$. But it seems that this is not happening. Because when probing the density of two interfering BECs with absorption imaging, there is a certain exposure time. And if only very few particles escape the BEC during this exposure time, and the phase jumps significantly each time, the interference pattern should be completely blurred. So somehow the phase of a BEC seems to be stable when only a few particles are removed. My question is then "why".
 A: I had posted an answer about 10 days ago, but comments from user @Thomas made me realise that my understanding was lacking key aspects of the field. I have read quite a lot in this time and feel like I can provide another, hopefully better, answer.
Let me first start with a quote from Anderson that said:

Do two superfluids that have never ‘‘seen’’ one another possess a
definite relative phase?

to which Leggett had pointed out that the question is meaningless as long as no measurement is performed on the system.
Now, there are two ways of answering the question. The second one is correct, while the first one is kinda correct for some situations, but it is conceptually and mathematically easier and still predicts the correct experimental outcomes, and is therefore ubiquitous in the literature.
Coherent state formalism:
The typical ansatz is that the BEC wavefunction can be written as $\psi(x) = \sqrt{N_0}\mathrm{e}^{\mathrm{i}\theta(x)}$. I.e. the BEC already has "chosen" a phase, compatibly with the spontaneous symmetry breaking formalism. The broken symmetry is the $U(1)$ symmetry associated with particle number conservation. It follows that there exists a commutation relation between number and phase operators $[\hat N, \hat \theta] \propto \mathrm{i}$. The number of atoms in a BEC is hence $N_0 \pm \Delta N$, where $\Delta N$ are fluctuations. The phase $\theta$ is decently well-defined for large atom numbers as its error $\Delta \theta \sim 1/\sqrt{N} \sim 0$, so people are usually ok with this idea of the BEC "choosing" a phase.
In a BEC, all atoms occupy the same (ground) state. They are indistinguishable, so you don't need a tensor product. Actually you first start with one state $\mathrm{e}^{\mathrm{i}\theta(x)}$ which you then literally multiply $\sqrt{N_0}$ times so that $\psi$ is normalised to $N_0$.
In this formalism, removing one particle just means that the scaling factor changes from $\sqrt{N_0}$ to $\sqrt{N_0-1}$, without affecting the phase.
So, in this formalism, a BEC "chooses" a phase $\theta$. Two independent BECs choose two random phases $\theta_1$ and $\theta_2$. These two phases themselves are not observable, but, upon interference of the two BECs, an interference pattern appears with fringe separation related to the phase difference $\theta_1 - \theta_2$. This phase difference will be different for each repetition of the experiment, as each phase is random.
While convenient, easy, and intuitive, this formalism falls short on multiple aspects. One of these, as @Thomas had correctly pointed out in the comments to my previous answer, is baryon number conservation. I.e. in an experimentally trapped BEC, the total number of atoms is actually conserved. The coherent state comes with a fluctuation in atom number. This is fine as long as the thermal fraction can provide a reservoir of particles for the BEC, such that they can account for these fluctuations (with opposite sign, $\Delta N_{\mathrm{cond}}=-\Delta N_{\mathrm{exc}}$) ensuring that the total number of atoms is conserved. However, the atom number fluctuation tends to zero experimentally as $T\rightarrow 0$ meaning the coherent state cannot accurately describe the BEC at $T=0$ (100% condensed, so no reservoir to account for the fluctuations in atom number).
Fock state formalism:
The key paper for this is Relative phase of two Bose-Einstein condensates, where they explicitly show (through a thought experiment$^\dagger$) that you can still have an interference pattern (i.e. a phase difference) between two BECs even though each BEC is a Fock state, i.e. particle conserving. To re-iterate, in an actual experiment you can measure the total number of trapped atoms, so a coherent state cannot be the only explanation for the observed physical phenomena.
To quote directly from them:

The notion of phase broken symmetry [and hence coherent state stuff] is therefore not indispensable in
order to understand the beating of two condensates. On the other
hand, it provides a simple way of analyzing such an experiment, while
[...] Fock states are more difficult to handle in such a situation.
The predictions for an initial Fock state and for an initial coherent
state with random phase are therefore equivalent, but the result for
the coherent state is obtained in a much more straightforward and
intuitive manner than for the Fock state

They start with a number-specific Fock state $|N\rangle$, which therefore has a flat phase distribution. So if you remove particles, you just start with a different Fock state  $|N-1\rangle$ which still has a flat (completely random) phase distribution.
They find that both formalisms give you an interference pattern. And, actually, just from the experimental beating pattern you cannot ascertain whether each BEC was a Fock state or a coherent state to start with.
The two formalisms are associated with two different interpretations:

To summarize, we have two different points of view on the system: for
an initial coherent state, the measurement ‘‘reveals’’ the pre-existing
phase [...]; for an initial Fock state, the detection
sequence ‘‘builds up’’ the phase.


EDIT: Let me specify that all the above is applicable to a true Bose-Einstein condensate, that is, a non-interacting system. With repulsive interactions, a specific phase for the condensed fraction is chosen both because it minimises the many-body ground state energy (that can be shown to go as $\propto \cos \theta_0$, so that $\theta_0 = \pi$ is favoured), and because an interacting theory now falls within the spontaneous symmetry breaking formalism. With an order parameter, $U(1)$ breaking etc.  The (weakly) interacting Bose gas is now in the many-body ground state, which does not correspond to the macroscopic occupation of a single-particle state that would be a true BEC. But if you write the Bogoliubov vacuum in terms of real single-paticle operators, however, you see that the original (non-interacting) condensed fraction has "only" lost some atoms to excitations through quantum depletion. Now these excitations act like a reservoir of particles that justify the grand-canonical ensemble formalism whereby the phase of the condensate is fixed and its number of particles uncertain because of fluctuations with the reservoir. The total number of particles is still conserved.

$^\dagger$ Through tough thorough thought.
