Observation of gauge in artificial magnetic fields
Disclaimer: I do particle physics / cosmology, so this is definitely outside my field, apply grains of salt to this answer appropriately.
I think Reference [29] (Lin et al, arxiv reference: 1008.4864) honestly does a better job of explaining what is going on (which makes sense, the impression I get is that 1008.4864 is a foundational paper in this subfield).
The gist, as I understand it from Lin et al, is that the hamiltonian for the system they are interested in--neutral atoms in a BEC state coupled to two intersecting lasers--can be written like:
\begin{equation} H_{BEC} = (\vec{p} - \vec{p}_{min})^2, \end{equation} where $\vec{p}_{min}$ is a quantity that can be controlled by the experimenters (adjusting the lasers).
Lin et al note that this is formally similar to the hamiltonian for a charged particle moving in a background electromagnetic field \begin{equation} H_{charged} = (\vec{p} - q \vec{A})^2 + q \phi, \end{equation} where $\vec{A}$ and $\phi$ are the vector and scalar potentials respectively. Of course, $H_{charged}$ is gauge invariant.
To implement this analogy, they introduce an analogue vector potential $\vec{A^*}$ and scalar potential $\phi^*$ (they don't really introduce $\phi*$ but let's run with it for now). Then Lin et al identify \begin{equation} \vec{p}_{min} = q^* \vec{A^*}, \end{equation} and also work in a gauge where \begin{equation} \phi^* = 0. \end{equation} Then in that gauge, the analogue electric field is given by $\vec{E^*} = -\dot{\vec{A^*}}$.
Thus, Lin et al point out that there are effects from setting up a time dependent $\vec{A}$ (ie a time dependent $\vec{p}_{min}$). From the perspective of the analogue gauge theory, this is due to the fact that the analogue electric field $\vec{E}^*$ is non-zero. The analogue electric field is a gauge invariant, observable quantity.
The passage you cite (from 1503.08243, Kennedy et al) suggests that the effect they measure comes from a time dependent $\vec{A^*}$. Again this would lead to a non-zero $\vec{E^*}$.
Of course, from the perspective of the analogue gauge theory, they are free to perform a gauge transformation, and they must get the same answer because physical observables must be gauge invariant (this point is ironclad--if the rest of this answer is wrong, this one point can't be wrong unless the analogy to gauge fields completely breaks down).
However, a gauge transformation will necessarily turn on $\phi^*$. In other words, $\phi^*$ will be nonzero in any other gauge. This requires one to change the original Hamiltonian to take this into account. What will still be true is that \begin{equation} H_{BEC}=(\vec{p}-\vec{p}_{min})^2 = (\vec{p} - q\vec{A^*})^2 + q \phi^* \end{equation} so obviously in this new gauge we can no longer identify $q\vec{A^*} = \vec{p}_{min}$. I think this is really what Kennedy et al are getting at, this relationship between $\vec{A^*}$ and $\vec{p}_{min}$ is not gauge invariant.
When the new hamiltonian is used correctly, the final answer will be the same in any gauge.
However I think that actually showing this works would be overkill--the bottom line I think is that $\vec{E^*}$ is non-zero, so everyone in the end is making measurements of a gauge invariant quantity.
Update 7/4
So I had a chance to look at this more. In the end I just have to disagree with the quote you cited--the observables do not depend on the gauge. However the gauge they choose is particularly nice, and finding a manifestly gauge invariant formulation of what they are doing might not be worth it. The bottom line is that once you gauge fix, every combination of operators you write down is gauge invariant (since there is no gauge freedom left), and therefore there is guaranteed to be a gauge invariant combination of operators that reduces to the combination you wrote down in the gauge you picked. In other words, one completely valid way to describe a gauge invariant quantity is to say what it looks like in a well defined gauge. What I think is going on is that the observable Kennedy et al are measuring (column density in momentum space) is very natural in one gauge fixed version of the problem, but finding the manifestly gauge invariant version would be unnecessarily complicated.
More details:
Powell et al (1009.1389) is a really good paper that discusses the theoretical aspects of what is going on in setups like the one used in Kennedy et al. The underlying formalism you need is gauge theory on a lattice. The basic idea is that there are fermion fields living on a given lattice site that create particles at that lattice site. In Kennedy et al these are referred to as $a_{m,n}$ (where $m,n$ are indices on a 2D lattice). There are also link fields, which are given by Wilson lines that connect the lattice site $(m,n)$ to the lattice site $(m',n')$ \begin{equation} W_{(m,n),(m',n')} = \exp \left(i\int_{(x_m,y_n)}^{(x_{m'},y_{n'})} d\vec{x}\cdot\vec{A}\right) = e^{i\phi_{(m,n),(m',n')}} \end{equation} The last line works because this is a $U(1)$ gauge field, so the integrals are numbers.
Both the $a_{m,n}$ operators, and the phases $\phi_{m,n}$, transform under gauge transformations. The gauge transformations occur on each site independently, so we can write the parameters of the gauge transformation as $\lambda_{m,n}$. The $a$ operators transform as \begin{equation} a_{m,n} \rightarrow e^{i \lambda_{m,n}} a_{m,n} \end{equation} and the Wilson lines transform as \begin{equation} W_{(m,n),(m',n')} \rightarrow e^{i \lambda_{m,n}} W_{(m,n),(m',n')} e^{-i \lambda_{m',n'}} \end{equation} or equivalently \begin{equation} \phi_{(m,n),(m',n')} \rightarrow \phi_{(m,n),(m',n')} + \lambda_{m,n} - \lambda_{m',n'} \end{equation} As a particle theorist / cosmologist, I am more familiar with the above formulas in their continuum form, where I would call $a$ by the name $\psi$, so $\psi(x)\rightarrow e^{i\lambda(x)}\psi(x)$ and $W(x,y) \rightarrow e^{i \lambda(x)} W(x,y) e^{-i\lambda(y)}$.
One key observation is that the "hopping Hamiltonian" from Kennedy et al is gauge invariant \begin{equation} H = -t \sum_{m,n} a_{m+1,n}^\dagger e^{i\phi_{(m+1,n),(m,n)}} a_{m,n} \end{equation} which can be seen using the transformation rules above. Incidentally, my particle-y instincts are to think of the above as a discretized version of the fermionic part of the QED lagrangian, $\bar{\psi} i \gamma^\mu D_\mu \psi$, where $D$ is the gauge covariant derivative.
The fact that the hopping Hamiltonian is gauge invariant really means that nothing physical is going to end up depending on the choice of gauge. To the extent that this Hamiltonian describes the system Kennedy et al are measuring, nothing can end up depending on a choice of gauge because the underlying hamlitonian describing all of the dynamics does not. (This could be broken, for example, if (1) the approximation of the full dynamics of the system by this gauge theory breaks down in a way that breaks gauge invariance, or (2) if the way that the experimental apparatus is coupled to the BEC breaks the gauge symmetry. I am assuming both of those don't happen--if they do that is more the fault on the experimental side than the theoretical side, ie it is a boring breaking of gauge invariance).
For example, the number of particles on each site \begin{equation} n_{m,n} = a_{m,n}^\dagger a_{m,n} \end{equation} The number of particles is gauge invariant (as you can check from the rules and physically has to be the case since the number of particles is observable).
Both Powell et al and Kennedy et all find it convenient to work in a gauge where the phases only depend on one lattice direction, so \begin{equation} \phi_{(m,n),(m',n')} \rightarrow \phi_{m,m'} \end{equation} This is a very nice gauge for what they want to do. In particular, the translations in the $y$ direction commute with the Wilson line operators, but translations in $x$ do not. Their basic point is that all gauge fixings will force translation invariance to be broken somehow, so full translation invariance is not a real symmetry of the system.
Now the measurements in Kennedy et al, as far as I can tell, are really done in momentum space (from a theoretical point of view momentum space is nice for this problem because this diagonalizes the Hamiltonian). The momentum space operators are \begin{equation} \tilde{a}_{p,q} = \sum_{m,n} e^{i 2\pi (p m+qn) / N} a_{m,n} \end{equation} where $N$ is the number of lattice sites.
Things now get complicated because the momentum space operators don't have obviously nice transformation properties under gauge transformations (the gauge transformation of the $\tilde{a}_{p,q}$ will end up being some convolution of the gauge parameters with the real space operators $a_{m,n}$). This is related to the fact that the commutator of the Hamiltonian with the translation operator will be complicated in a general gauge.
So, what I think is going on is that Kennedy et al construct a column density in momentum space, which I am guessing amounts to the probability which you can compute in a given state by $\langle \tilde{a}^\dagger_{p,q}\tilde{a}_{p,q} \rangle$, where the $\tilde{a}_{p,q}$ are defined in the gauge that they describe. One frustrating thing is that I am not 100% sure what specific combination corresponds to the plots they make, so I can't be more explicit about what I'm saying, but conceptually it doesn't matter what precise combination of $\tilde{a}$'s they are plotting.
This does not make the observable gauge dependent, however it does mean that showing the gauge invariance is tricky. There is guaranteed to be some gauge invariant combination of Wilson lines and fermion operators that reduces to the combination that Kennedy et al plot, in the gauge that they pick. One avenue to discover the precise gauge invaraint combination is by guessing--if you find one gauge invariant combination that reduces to their observable, that is the correct one. Another more systematic approach is to take their observable, written in the gauge that they chose, and perform an arbitrary gauge transformation. The result will likely be messy (since the momentum space operators don't have nice transformation properties), but you are guaranteed to be able to write the result in terms of manifestly gauge invariant objects if you do everything correctly (you will probably have to add in gauge transformed combinations of operators that were zero in the original gauge, and the net goal is to cancel out all of the dependence on the gauge parameter).
In other words, once you gauge fix, you can write down any arbitrarily complicated combination you like of the operators you have and you are guaranteed to be talking about gauge invariant quantities, since there's no gauge freedom left. However, finding the manifestly gauge invariant form can be hard.
In the problem that Kennedy et al are considering, there is such a natural choice of gauge that I think they basically want to argue that there's no point in finding the gauge invariant form of what they are measuring--the main pragmatic reason to find a gauge invariant form would be if different groups were using different gauges and needed to compare their answers. The gauge invariant form could be interesting theoretically to get more insight into the system. Based on what Powell says in section II, I think the gauge invariant formulation involves studying the properties of the projective symmetry group of the system. But that would definitely be beyond the scope of an experimental paper.