\section{Collapse-induced dephasing and diffusion of a two-mode BEC}\label{app:CSL_collSpin}

For $N$ indistinguishable bosonic particles in two modes, represented by annihilation operators $\mathsf{c}_{a,b}$, we can make use of the collective spin representation \cite{Arecchi1972}, with $J=N/2$ and spin operators $\mathsf{J}_x=(\mathsf{c}_a^\dagger\mathsf{c}_b+\mathsf{c}_b^\dagger\mathsf{c}_a)/2$, $\mathsf{J}_y=(\mathsf{c}_a^\dagger\mathsf{c}_b-\mathsf{c}_b^\dagger\mathsf{c}_a)/2i$, and $\mathsf{J}_z=(\mathsf{c}_a^\dagger\mathsf{c}_a-\mathsf{c}_b^\dagger\mathsf{c}_b)/2$.
In the presence of CSL, the time evolution of the collective spin state representing the condensate can be described by the master equation \cite{Javanainen1997PhaseDispersion,schrinski2019macroscopicity}
\begin{align}\label{eq:LioSup}
\partial_t\rho=\frac{1}{i\hbar}[\epsilon\mathsf{J}_z+\hbar\zeta\mathsf{J}_z^2,\rho]+\mathcal{L}\rho
\end{align}
The first term in the Hamiltonian describes the free rotation around the $z$-axis of the generalized Bloch sphere at an angular frequency given by the energy difference $\epsilon$ between the two involved modes. The second term accounts for atom-atom interactions in the condensate to lowest order close to the equator, $\hbar \zeta=2 ({\rm d} \mu/ {\rm d} n )_{n=0}$, where $n$ is the population difference (i.e.~two times the $\mathsf{J}_z$-eigenvalue) and $\mu$ is the chemical potential \cite{Javanainen1997PhaseDispersion}. CSL in second quantization contributes the Lindblad generator
\begin{align}\label{eq:TwoModeLindblad}
\mathcal{L}\rho=&\lambda\frac{ m^2}{\rm u^2}\,\frac{r_{C}^3 }{\pi^{3/2}} \int {\rm d}^3 {q}\, e^{-q^2r_C^2}\nonumber\\
&\times\left[\mathsf{A}(\mathbf{q})\rho\mathsf{A}^\dagger(\mathbf{q})-\frac{1}{2}\left\{\mathsf{A}^\dagger(\mathbf{q})\mathsf{A}(\mathbf{q}),\rho\right\}\right],
\end{align}
with $\mathsf{A}(\mathbf{q})=\int d^3{p}\,\mathsf{a}^\dagger(\mathbf{p})\mathsf{a}(\mathbf{p-\hbar q})$ and $\mathsf{a} (\mathbf{p})$ the particle annihilation operator in momentum representation. Restricting the dynamics to the two relevant modes of the interferometer, we can expand the Lindblad operator as
\begin{align} \label{eq:TwoModeLindbladOp}
\mathsf{A}(\mathbf{q}) =& \sum_{j,k \in \{a,b\}} \underbrace{\langle \psi_j|e^{i \mathbf{q}\cdot\boldsymbol{\mathsf{r}} }|\psi_k \rangle}_{=:W_{jk} (\mathbf{q})} \mathsf{c}_j^\dagger \mathsf{c}_k \\
=&\left[ W_{aa}(\mathbf{q}) + W_{bb}(\mathbf{q}) \right] \frac{N}{2} + \left[ W_{aa}(\mathbf{q}) - W_{bb}(\mathbf{q}) \right] \mathsf{J}_z \nonumber \\
& + \left[W_{ab} (\mathbf{q}) + W_{ba} (\mathbf{q}) \right] \mathsf{J}_x + i \left[W_{ab} (\mathbf{q}) - W_{ba} (\mathbf{q}) \right] \mathsf{J}_y. \nonumber 
\end{align}
Here, $|\psi_{a,b}\rangle$ denote the single-particle wave functions of the two modes, $\boldsymbol{\mathsf{r}}$ the respective position operator, and $W_{jk} (\mathbf{q}) = W_{kj}^{*} (-\mathbf{q})$ the overlap matrix elements of the momentum displacement operator. Note that the $\mathsf{J}_y$-term vanishes for bound states with real-valued wave functions.
The terms that we have omitted here would describe incoherent hopping of atoms between the condensate and other undetected modes, which causes atom loss. Formally, we can account for this CSL-induced heating effect in the data analysis of an experiment by introducing conditional outcome probabilities, given that the detected $N$ atoms have remained in the condensate at all times. This way, the CSL bounds would not depend on the additional classical observation of depletion \cite{schrinski2019macroscopicity}. Moreover, since particles lost from the condensate (or, less likely, regained ones) will always increase the phase uncertainty $\sigma_\varphi$, our omission of the depletion effect only underestimates the influence of CSL. Additionally, we will ignore the $z$-coordinate as the condensate is elongated in this direction  and  integrated out upon readout \cite{Schmiedmayer2013,Schmiedmayer2014}.

Plugging the Lindblad operator \eqref{eq:TwoModeLindbladOp} into the master equation \eqref{eq:TwoModeLindblad} reduced to the $xy$-plane, we identify two $r_{C}$-dependent geometry factors that will be relevant for the two configurations discussed here,
\begin{align}
    f_\mathrm{P} (r_{C}) &:= \frac{r_{C}^2 }{2\pi} \int {\rm d}^2 {q}\,  e^{-q^2 r_C^2} |W_{aa}(\mathbf{q}) - W_{bb}(\mathbf{q})|^2, \\
    f_\mathrm{S} (r_{C}) &:=\frac{r_{C}^2 }{2\pi} \int {\rm d}^2 {q}\, e^{-q^2 r_C^2} |W_{ab}(\mathbf{q}) + W_{ba}(\mathbf{q})|^2 .
\end{align}
For the MZI scenario in the main text, we consider two spatially distinct, but otherwise identical modes with no overlap, $\langle \mathbf{r} | \psi_b \rangle = \langle \mathbf{r} - \Delta \mathbf{e}_x | \psi_a \rangle$ such that $W_{ab} (\mathbf{q}) \approx 0$. The CSL generator \eqref{eq:TwoModeLindblad} then reduces to pure dephasing between the modes, i.e.~a phase noise channel in the collective spin representation \cite{Ma2011},
\begin{align}\label{eq:PhaseDiffusion}
\mathcal{L}\rho &= \Gamma_\mathrm{P}\left[\mathsf{J}_z\rho\mathsf{J}_z-\frac{1}{2}\left\{\mathsf{J}_z^2,\rho\right\}\right], \quad \Gamma_{\rm P} = 2\lambda \frac{m^2}{\rm u^2} f_\mathrm{P} (r_{C}).
\end{align}
The typical scaling behavior of the geometry factor $f(r_{C})$ for the dephasing rate becomes evident for Gaussian modes of width $w_x=w_y \ll \Delta_x$,
\begin{equation}
    f_\mathrm{P} (r_{C}) = \frac{1-\exp\left[-\Delta_x^2/4(w_x^2+r_C^2)\right]}{1+w_x^2/r_C^2}.
\end{equation}

If the two interfering modes overlap spatially, then $W_{ab} \neq 0$ and CSL will not only induce dephasing, but also diffusion. For the SWI configuration in the main text, we consider a condensate split between the ground and first excited state of a harmonic potential in $x$-direction, with trapping frequency $\omega$ and ground-state width $x_0 = \sqrt{\hbar/2m\omega}$. Moreover, we assume both modes to be Gaussian along the $y$-direction, with a width $w_y$, and we omit the elongated $z$-profile as argued above. This leads to
\begin{align}
    W_{aa} (\mathbf{q}) &= e^{-q_y^2w_y^2/2 - q_x^2 x_0^2/2}, \nonumber \\
    W_{bb} (\mathbf{q}) &= \left( 1 - q_x^2 x_0^2 \right) W_{aa} (\mathbf{q}), \nonumber \\
    W_{ab} (\mathbf{q}) &= i q_x x_0 W_{aa} (\mathbf{q}) = W_{ba} (\mathbf{q}), \nonumber \\
    f_\mathrm{P} (r_{C}) &= \frac{3 r_{C}^2 x_0^4 }{8\sqrt{r_{C}^2+w_y^2}(r_{C}^2 + x_0^2)^{5/2}}, \nonumber \\
    f_\mathrm{S} (r_{C}) &= \frac{r_{C}^2 x_0^2 }{\sqrt{r_{C}^2+w_y^2}(r_{C}^2 + x_0^2)^{3/2}}.
\end{align}
We have set $w_y=x_0/\sqrt{6}$ for the proposed setups in the main text. 

\section{Phase dispersion together with diffusion in phase and number difference} \label{app:DiffWigner}


Here we recall the phase-space method of Ref.~\cite{schrinski2019macroscopicity} to solve the time evolution of two-mode BEC states according to the master equation \eqref{eq:LioSup} in the presence of dephasing and diffusion,
\begin{align}
\mathcal{L}\rho &= \Gamma_\mathrm{P}\left[\mathsf{J}_z\rho\mathsf{J}_z-\frac{1}{2}\left\{\mathsf{J}_z^2,\rho\right\}\right]+
\Gamma_\mathrm{S}\left[\mathsf{J}_x\rho\mathsf{J}_x-\frac{1}{2}\left\{\mathsf{J}_x^2,\rho\right\}\right]
\end{align}
We consider high atom numbers $N \gg 1$ and coherent superpositions with a well-defined interferometric phase, i.e.~collective spin states that are sharply localized on the equator of the generalized Bloch sphere, $\langle \mathsf{J}_z \rangle \approx 0$ and $\Delta \mathsf{J}_{x,y,z}^2 \ll N^2$. To a very good approximation, the state is then represented by a continuous Wigner function $w_t(\varphi,n)$ in the flat phase space of the conjugate variables $n$ (population difference) and $\varphi$ (interferometric phase angle), provided the support of $w_t$ is limited to $|n| \ll N$ and stretches over $\varphi$-windows $\ll 2\pi$. The latter constraint can be alleviated with help of a trigonometric mapping of the Wigner function onto a periodic function over the equator, $\varphi \in [0,2\pi)$, but this will not be necessary for the cases studied here. In a rotating frame that absorbs the free linear evolution of the interferometric phase, $\varphi \mapsto \varphi + \epsilon t/\hbar$, the master equation \eqref{eq:LioSup} translates into the Fokker-Planck equation
\begin{align}\label{eq:WignerFunction}
\partial_t w_t(\varphi,n) \approx&-\zeta n\partial_\varphi w_t(\varphi,n)+\frac{N^2\Gamma_\mathrm{S}}{4}\partial_n^2w_t(\varphi,n)\nonumber\\
&+\frac{\Gamma_\mathrm{P}}{2}\partial_\varphi^2w_t(\varphi,n).
\end{align}
It further approximates the spin diffusion term, which would normally have an additional $\varphi$-dependence, by its angular average, resulting in an additional factor $1/2$ for the $\Gamma_{\rm S}$-dependent term. The approximation is valid as long as the diffusion rate is small, $\Gamma_{\rm S} \ll \epsilon/\hbar$, and the interference time extends over at least one free oscillation period.

Equation \eqref{eq:WignerFunction} constitutes a Gaussian channel, and every initially Gaussian Wigner function will thus remain Gaussian. In particular, the $\varphi$-marginal of the Wigner function, which represents the phase distribution of the interfering state after time $t$, will be of the form 
\begin{align}
\int dn\,w_t(\varphi,n)=\frac{1}{\sqrt{2\pi}\sigma_\varphi(t)}\exp\left[-\frac{\varphi^2}{2\sigma_\varphi^2(t)}\right].
\end{align}
More generally, we can obtain an exact solution for arbitrary initial states with help of the characteristic function in Fourier space,
\begin{equation}
    \chi_t (s,q) := \int {\rm d} \varphi {\rm d} n\, w_t (\varphi,n) e^{is \varphi + iq n}.
\end{equation}
It evolves according to
\begin{equation}
    \partial_t \chi_t (s,q) = \zeta s \partial_q \chi_t (s,q) - \left( \frac{\Gamma_{\rm P}}{2} s^2 + \frac{N^2 \Gamma_{\rm S}}{4} q^2 \right) \chi_t (s,q),
\end{equation}
which yields the solution
\begin{align} \label{eq:chit_sol}
    \chi_t (s,q) &= \chi_0 (s,q +\zeta t s) \\
    &\times \exp \left[ - \frac{\Gamma_{\rm P}t}{2} s^2 - \frac{N^2 \Gamma_{\rm S}t}{4} \left( q^2 + \zeta t q s + \frac{\zeta^2 t^2}{3} s^2 \right) \right]. \nonumber
\end{align}
The characteristic function generates the moments of the state's phase distribution as $\langle \varphi^k \rangle_t = (-i\partial_s)^k \chi_t(0,0)$. Plugging the above solution into $\sigma_\varphi^2 (t) = \langle \varphi^2 \rangle_t - \langle \varphi \rangle_t^2$ leads to the formula \eqref{eq:phaseVar_withDiff} in the main text, if we further assume that initially $\langle n\rangle_0, \langle n\varphi \rangle_0 = 0$.

For the two-step echo protocol, the solution can be constructed by first computing $\chi_{t/2} (s,q)$ according to \eqref{eq:chit_sol} and then inserting the result as the new initial condition for \eqref{eq:chit_sol} with a sign-flipped $\zeta$. The result simplifies to
\begin{equation}\label{eq:chit_sol_echo}
    \chi_t (s,q) = \chi_0 (s,q) \exp \left[- \frac{\Gamma_{\rm P}t}{2} s^2 - \frac{N^2 \Gamma_{\rm S}t}{4} \left( q^2 + \frac{\zeta^2 t^2}{12} s^2 \right) \right],
\end{equation}
which produces the results discussed in the main text.



 

\end{document}
