text:
\section{Introduction}
\label{sec:intro}

The consistent inclusion of Lorentz symmetry breaking in a curved space is certainly one of the most important open problems within studies of this phenomenological idea.  Traditionally, the Lorentz symmetry breaking has been introduced in two manners, 1) the explicit one, where a constant vector (tensor) is introduced in the theory from the very beginning, and 2) the spontaneous one, where this constant vector (tensor) arises as the vacuum expectation value of some dynamical field (see \cite{KosGra} for a detailed discussion of approaches to Lorentz symmetry breaking in gravity).

It should be { emphasized} that while the explicit approach is very convenient in the case of linearized gravity, there are essential difficulties with its application in a full-fledged scenario. Indeed, unlike in flat space, in curved space it is highly problematic to define constant vectors (or tensors) because the condition of vanishing covariant derivatives of a given tensor imposes constraints on the background metric which are difficult or even impossible to satisfy. Otherwise, quantum corrections in curved-space extensions of known Lorentz-breaking field theory models will involve an infinite tower of new terms proportional to covariant derivatives of ``constant" Lorentz-breaking tensors, which makes the calculations much more complicated (see {\it e.g.} \cite{ShaBP}). Another difficulty related with this approach is that in a curved space-time, the group of general covariant transformations plays a double role, being not only the extension of the Lorentz group but also the gauge group. As a result, introducing terms that explicitly break Lorentz symmetry will imply violation of the gauge symmetry as well, which makes the spontaneous symmetry breaking approach much more appropriate (that see discussion in \cite{KosMew}).

 In {a curved space-time}, the most convenient form of { introducing} spontaneous Lorentz symmetry breaking is the bumblebee model, first  proposed  in \cite{KosGra}, where the breaking mechanism is implemented via a new dynamical vector field with a nontrivial potential characterized by a {continuous set of minima}. The first studies of modifications of known gravitational solutions within the bumblebee gravity have been carried out in \cite{Bertolami}.
Various aspects of the bumblebee gravity have been studied in numerous papers, {\it e.g.} \cite{Seifert,Maluf,Godel,Casana:2017jkc}, all of which assume that the underlying geometry is of (pseudo)-Riemannian type. However, in the absence of any empirical evidence supporting that the space-time structure is { necessarily} Riemannian, or otherwise\footnote{Here Riemannian means that the connection is completely specified by the metric, as opposed to being an {\it a priori } independent geometrical entity.}, in the high energy regime  \cite{Will:2014kxa,Berti:2015itd,Daniel:2010ky,Krawczynski:2012ac,Damour:1991rd,Kramer:2006nb,Clifton:2005aj,Abbott:2016blz}, it is legitimate to explore other alternatives. In particular, it is well-known that the metric-affine (or Palatini) formulation of gravity theories beyond GR yields field equations which are inequivalent to those obtained in the purely metric approach, and it has been argued that quantum effects could be encoded at low energies by a non-Riemannian counterpart of the space-time geometry \cite{Olmo:2015bha,Olmo:2015wwa,Lobo:2014nwa,Hossenfelder:2017rub}. This fact motivates us to initiate the exploration of Lorentz breaking phenomenology of bumblebee gravity in its metric-affine formulation, in which the connection is treated as a geometrical object a priori independent of the metric. In this sense, it is worth mentioning that Lorentz-breaking extensions of both the Standard Model and General Relativity have already been tested experimentally, leading to stringent constraints on Lorentz-breaking parameters, see {\it e.g.} \cite{datatable,Kostelecky:2015dpa,Bailey:2014bta,Kostelecky:2016kfm,Kostelecky:2010ze,Kostelecky:2008in,Bailey:200fd}.

It has been recently shown that some metric-affine theories of gravity beyond GR have a peculiar behavior that makes them depart from their metric counterpart and become potentially testable via elementary particle interactions. In the particular case of minimally coupled Ricci-Based Gravity theories (RBG's), those in which the gravity Lagrangian is a (projective invariant) function of the metric and the Ricci tensor, one finds that the space-time metric picks up two types of contributions, one coming from the integration over the matter sources and which is mainly responsible for the space-time curvature, and another coming from the energy-momentum density of the local sources. The latter contribution is responsible for the existence of a nonzero non-metricity tensor, $Q_{\mu\alpha\beta} =\nabla_\mu g_{\alpha\beta}$, in regions where the stress-energy of the matter fields is not covariantly constant. The origin of these contributions in the metric can be traced back to the fact that these theories admit an Einstein frame (via a non-conformal transformation) in which the matter sector features new non-linear interactions and the gravity sector is described by GR. Thus, from the original frame of the theory (RBG frame), these non-linear interactions appear encoded in the space-time metric and are responsible for the non-vanishing non-metricity tensor. The existence of the Einstein frame representation for RBG theories was devised in \cite{Afonso:2017bxr,BeltranJimenez:2017doy}, and it has been explicitly proven and used for applications for different matter sectors: perfect and anisotropic fluids, scalar fields and non-linear electrodynamics \cite{Afonso:2018bpv,Afonso:2018hyj,Afonso:2018mxn,Afonso:2019fzv,Delhom:2019zrb}.
The deviation from the Riemannian condition $\nabla_\mu g_{\alpha\beta}=0$ in the RBG frame is intimately related to a departure of the space-time metric from its Minkowskian form whenever matter fields are present even if the effects of curvature ({\it i.e.} Newtonian and post-Newtonian corrections) are negligible. As a result, particle physics experiments can be used to place strong constraints on the parameters of those gravity models \cite{Latorre:2017uve,Delhom:2019wir}.

In this work we will consider a so far unexplored route, in which the bumblebee model plays a key role. In this model the gravity Lagrangian is represented by the Einstein-Palatini GR Lagrangian and an additional (projectively invariant) non-minimal coupling term between the bumblebee field and the affine Ricci tensor. This sort of non-minimal coupling has not yet been considered in detail within metric-affine theories. In addition, keeping the gravitational part exactly as in GR avoids the generation of non-metricity induced by the stress-energy tensor of the matter fields, which allows to disentangle the different contributions that the non-metricity may have. As we will see, the resulting theory admits an exact formal solution for the independent connection which leads to the emergence of a non-metricity tensor generated by the non-minimal bumblebee coupling. The bumblebee field will then become coupled to the rest of the matter fields present in the theory due to its non-minimal coupling with the gravitational sector. Thus, even though the theory is initially formulated following the postulates of metric theories of gravity so as to satisfy the Einstein equivalence principle \cite{Will}, in the end the resulting coupling of the bumblebee field to the matter sector implies a direct violation of this principle, as one would expect in a Lorentz violating theory. Besides deriving the field equations and discussing their resolution, in this paper we will also discuss the weak field, post-Minkowskian limit with a focus on the modified dispersion relations and stability of scalar and Dirac fields.

The structure of our paper is as follows. In section 2, we present a general family of Ricci-Based Gravity theories (RBGs) with non-minimal couplings between matter and the connection, introducing the  metric-affine bumblebee model as a specific implementation. After discussing the field equations of the bumblebee model, in section 3, we study the weak-field, post-Minkowskian limit and the effective dynamics when curvature effects are negligible as compared to the contribution of non-metricity. Section 4 is devoted to the derivation and discussion of the modified dispersion relations and stability corresponding to scalar and spinor fields. We conclude with a summary and discussion of the results.

\section{Metric-affine Bumblebee model as a non-minimally coupled RBG}

\subsection{General case in the RBG framework.}

The analysis of metric-affine RBG theories has focused  so far on modifications of the gravitational sector but keeping a minimally coupled matter sector, {\it i.e.} with the matter fields interacting only with the metric, not with the connection. Such theories can be described by an action of the form
\begin{eqnarray}
\mathcal{S}=\frac{1}{2\kappa^2}\int d^4 x\,\sqrt{-g}F(g^{\mu\nu},R_{(\mu\nu)}(\Gamma))+ \mathcal{S}_{m}(g_{\mu\nu},\Psi_i) \ ,
\label{MinCoupRBG0}
\end{eqnarray}
where $\kappa^2=8\pi G$, $\Psi_i$ denotes a collection of minimally coupled matter fields, and $R_{(\mu\nu)}(\Gamma))$ represents the symmetrized Ricci tensor of the independent connection $\Gamma^\alpha_{\mu\nu}$. We use units $\hbar=c=1$.  A straightforward way of {generalizing} the above action is to allow for a non-minimal coupling between some matter fields $\Phi_i$ and the connection via $R_{(\mu\nu)}(\Gamma))$, such that (\ref{MinCoupRBG0}) turns into
\begin{eqnarray}
\nonumber\mathcal{S}&=&\frac{1}{2\kappa^2}\int d^4 x\,\sqrt{-g}F(g^{\mu\nu},R_{(\mu\nu)}(\Gamma),\Phi_i)+ \mathcal{S}_{m}(g_{\mu\nu},\Psi_i),\\
\label{MinCoupRBG}
\end{eqnarray}
where the $\Phi_i$ have minimal kinetic terms in the sense of \cite{Delhom:2020hkb}. The detail that only the symmetrized part of the Ricci tensor enters the action guarantees the existence of a projective symmetry, which turns out to be relevant for the stability of the theory by ensuring the absence of ghost-like degrees of freedom in RBG theories \cite{BeltranJimenez:2019acz,Jimenez:2020dpn} as well as in metric-affine scalar-tensor theories \cite{Aoki:2019rvi}.

In the metric-affine framework, metric and connection are treated as independent fundamental fields. Accordingly,  their field equations are obtained by extremizing the action without imposing any {\it a priori} relation between $\delta g_{\mu\nu}$ and $\delta \Gamma^\alpha_{\mu\nu}$. By varying upon the above action, one finds that the metric and connection field equations of the RBG action with non-minimally coupled matter fields are formally identical to those of the minimally coupled version, namely
\begin{eqnarray}
\frac{\partial F}{\partial g^{\mu\nu}}-\frac{1}{2}F g_{\mu\nu}&=&\kappa^2 T^M_{\mu\nu};\\
\nonumber\nabla^\Gamma_\lambda \left(\sqrt{-h}h^{\mu\nu}\right)&=&\sqrt{-h}\bigg[T^{\mu}{}_{\lambda\alpha}h^{\nu\alpha}+T^{\alpha}{}_{\alpha\lambda}h^{\mu\nu}+\\
&+&\frac{1}{3}T^{\alpha}{}_{\alpha\beta}\delta^\mu{}_{\lambda}\bigg]\label{conneqRBG},
\end{eqnarray}
where $T^{\alpha}{}_{\mu\nu}=2\Gamma^{\alpha}{}_{[\mu\nu]}$  is the torsion tensor, $T^M_{\mu\nu}$ is the usual stress-energy tensor of the minimally-coupled matter sector $T_{\mu\nu}^{M}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L}_{M})}{\delta g^{\mu\nu}}$, and we have defined a new (inverse) metric $h^{\mu\nu}$ via $\sqrt{-h}h^{\mu\nu}=2\kappa^2 \sqrt{-g}\frac{\partial F}{\partial R_{(\mu\nu)}}$. Following \cite{BeltranJimenez:2017doy}, it is possible to show that the torsion terms on the right-hand side of (\ref{conneqRBG}) can be eliminated by exploiting the projective symmetry of the theory. After doing this, one finds that, up to a {non-physical} projective mode, the solution for the connection is given by the Levi-Civita connection of the metric $h_{\mu\nu}$, namely
\begin{equation}
{\Gamma}^{\mu}{}_{\alpha \beta}=\frac{1}{2}h^{\mu \lambda}\left(\partial_{\alpha} h_{\beta \lambda}+\partial_{\beta} h_{\lambda \alpha}-\partial_{\lambda} h_{\alpha \beta}\right) \ .
\end{equation}
We then see that the inclusion of matter fields that are non-minimally coupled through the Ricci tensor is still compatible with the fact that in RBG theories the connection is the Levi-Civita connection of some metric, which greatly simplifies the process of solving the field equations.

\subsection{Building the Einstein frame of generalized RBGs with matter couplings to $R_{\mu\nu}$.}

The fact that the connection can be solved as the Levi-Civitta connection of a metric $h_{\mu\nu}$ raises the question of whether these generalised RBGs also admit an Einstein frame representation. As we show next, the answer is positive and the Einstein frame can be constructed with the standard procedure (see {\it e.g.} \cite{BeltranJimenez:2017doy}) of linearising the action by introducing an auxiliary (symmetric) tensor field $\Sigma_{\mu\nu}$.
 Let us sketch how the Einstein-frame is reached for a generalised RBG of the form (\ref{MinCoupRBG}). Consider the action
\begin{eqnarray}
\nonumber\tilde{\mathcal{S}}&=&\frac{1}{2\kappa^2}\int d^4 x\,\sqrt{-g}\left[\frac{\delta f}{\delta \Sigma_{\mu\nu}}(R_{(\mu\nu)}-\Sigma_{\mu\nu})+f \right]+\\
&+& \mathcal{S}_{m}(g_{\mu\nu},\Psi_i) \ ,
\label{MinCoupRBGauxiliary}
\end{eqnarray}
where $f=F(g^{\mu\nu},\Sigma_{\mu\nu},\Phi_i)$ indicates that we have replaced all instances of $R_{(\mu\nu)}$ by $\Sigma_{\mu\nu}$. The field equations of $\Sigma_{\mu\nu}$ are algebraic, and imply that $\Sigma_{\mu\nu}=R_{(\mu\nu)}$ provided that $\frac{\delta^2 f}{\delta\Sigma^2}\neq0$. Hence, integrating out $\Sigma_{\mu\nu}$ we see that the {\it new} action (\ref{MinCoupRBGauxiliary}) is physically equivalent to the original one (\ref{MinCoupRBG}).

Let us now define the inverse metric $h^{\mu\nu}$ by
\begin{equation}
\sqrt{-h}h^{\mu\nu}\equiv \sqrt{-g} \frac{\delta f}{\delta \Sigma_{\mu\nu}}.
\end{equation}
Generally, we can see the above equation as an implicit definition of $\Sigma_{\mu\nu}$, which should be seen as a function of $(g_{\mu\nu},h_{\mu\nu},$ $\Phi_i)$. Thus, by solving algebraically in $\Sigma_{\mu\nu}$, we would be able to write $\Sigma_{\mu\nu}(g^{\mu\nu},h^{\mu\nu},\Phi_i)$. We can then perform a field redefinition in (\ref{MinCoupRBGauxiliary})
 by substituting $\Sigma_{\mu\nu}(g^{\mu\nu},h^{\mu\nu},\Phi_i)$ to arrive at
 \begin{eqnarray}
\nonumber\tilde{\mathcal{S}}&=&\frac{1}{2\kappa^2}\int d^4 x\,\left[\sqrt{-h}h^{\mu\nu}R_{\mu\nu}+{\mathcal{U}}(g^{\mu\nu},h^{\mu\nu},\Phi_i) \right]+\\
 &+&\mathcal{S}_{m}(g_{\mu\nu},\Psi_i),
\label{MinCoupRBGalmostEF}
\end{eqnarray}
where
\begin{eqnarray}
\nonumber{\mathcal{U}}(g^{\mu\nu},h^{\mu\nu},\Phi_i)\equiv\sqrt{-g}\left(f\left.-\frac{\delta f}{\delta\Sigma_{\mu\nu}}\Sigma_{\mu\nu}\right)\right|_{\Sigma_{\mu\nu}=\Sigma_{\mu\nu}(g^{\mu\nu},h^{\mu\nu},\Phi_i)}.
\end{eqnarray}
Now, we see that the field equations for $g^{\mu\nu}$ are also algebraic and, therefore, we can formally solve them to obtain $g^{\mu\nu}(h^{\mu\nu},\Phi_i,\Psi_i)$. Thus, integrating $g^{\mu\nu}$ out by means of this solution we arrive to the Einstein-Hilbert representation of (\ref{MinCoupRBG}), given by
\begin{equation}
{\mathcal{S}}_{\rm EH}=\frac{1}{2\kappa^2}\int d^4 x\,\sqrt{-h}h^{\mu\nu}R_{\mu\nu}+ \bar{\mathcal{S}}_{m}(h_{\mu\nu},\Phi_i,\Psi_i),\label{EFaction}
\end{equation}
where
\begin{eqnarray}
\nonumber\bar{\mathcal{S}}_{m}(h_{\mu\nu},\Phi_i,\Psi_i)&=&{\mathcal{S}}_{m}\big[g^{\mu\nu}(h^{\mu\nu},\Phi_i,\Psi_i),\Psi_i\big] +\\
\nonumber &+& \int d^4 x\,{\mathcal{U}}\big[g^{\mu\nu}(h^{\mu\nu},\Phi_i,\Psi_i),h^{\mu\nu},\Phi_i\big].
\end{eqnarray}
 We can clearly see here that the above action (\ref{EFaction}), which by construction is equivalent to \ref{MinCoupRBG}, formally describes a set of Einstein-like equations for the metric $h_{\mu\nu}$ coupled to a matter sector described by $\bar{\mathcal{S}}_{m}$ which will feature {\it new} interactions between the $\Phi_i$ and $\Psi_i$ sectors.

\subsection{The metric-affine bumblebee model as a particular case.}

As a particular case of the class of theories discussed above, we find the metric-affine version of the curved space-time bumblebee model (see \cite{Casana:2017jkc}), which is defined by an action of the form
\begin{eqnarray}
\nonumber\mathcal{S}&=&\int d^4 x\,\sqrt{-g}\Big[\frac{1}{2\kappa^2}\Big(R(\Gamma)+\xi B^{\alpha} B^{\beta} R_{\alpha\beta}(\Gamma)\Big)-\\
&-&\frac{1}{4}B^{\mu\nu}B_{\mu\nu}-V(B^{\mu}B_{\mu}\pm b^2)\Big] + \mathcal{S}_{m}(g_{\mu\nu},\Psi_i) \ .
\label{bumblebee}
\end{eqnarray}
 A key feature of this model is that the bumblebee field $B_\mu$ has a non-zero  vacuum expectation value (VEV) that spontaneously breaks Lorentz symmetry by introducing a privileged space-time direction. This field is coupled non-minimally to the space-time geometry via { the} $B^\mu B^\nu R_{\mu\nu}$ term. The parameter $\xi$ characterizes the strength of this non-minimal coupling between the bumblebee field and the affine connection through $R_{\mu\nu}(\Gamma)$, and we use the minimal coupling prescription for its kinetic term, so that the field strength of the bumblebee field is  defined as the exterior derivative of $B_{\mu}$ ({\it i.e.} $B_{\mu\nu}=(\mbox{d}B)_{\mu\nu}$).

We see that by identifying $\Phi_i\equiv B_\mu$ we find
\begin{align}
& f_B=(g^{\mu\nu}+\xi B^\mu B^\nu)\Sigma_{\mu\nu},\\
& \sqrt{-h}h^{\mu\nu}=g^{\mu\nu}+\xi B^\mu B^\nu,\label{auxiliarymetric}\\
& \mathcal{U}_B=-\frac{1}{4}B^{\mu\nu}B_{\mu\nu}-V(B^\mu B_\mu\pm b^2),\label{auxiliarypotential}
\end{align}
which would allow us to write the Einstein-frame action corresponding to \ref{bumblebee} once the metric $g_{\mu\nu}$ is solved in terms of the matter fields $(B_\mu, \Psi_i)$.

 By taking variations of Eq. (\ref{bumblebee}) with respect to the metric, the connection, and the bumblebee field, we obtain the following field equations:
\begin{align}
\nonumber&R_{(\mu\nu)}(\Gamma)-\frac{1}{2}g_{\mu\nu}\bigg(R(\Gamma)+\xi B^{\alpha}B^{\beta}R_{\alpha\beta}(\Gamma)\bigg)+\\
&+2\xi\bigg(B_{(\mu}R_{\nu)\beta}(\Gamma)\bigg)B^{\beta}=\kappa^2  T_{\mu\nu};\label{Riccieq}\\
&\nabla_{\lambda}^{(\Gamma)}\bigg[\sqrt{-g}\bigg(g^{\mu\nu}+\xi B^{\mu}B^{\nu}\bigg)\bigg]=0;\label{connectioneq}\\
&\nabla_{\mu}^{(g)}B^{\mu}_{\ \nu}=-\frac{\xi}{\kappa^2}B^{\beta}R_{\beta\nu}(\Gamma)+2 V^{\prime}B_{\nu},\label{bumblebeeeq}
\end{align}
where the prime in $V'$ denotes {derivative} of $V$ with respect to its argument. Here $T_{\mu\nu}$ is given by $T_{\mu\nu}=T_{\mu\nu}^{B}+T_{\mu\nu}^{M}$, where we have defined
\begin{eqnarray}
T_{\mu\nu}^{B}&=& B_{\mu\sigma}B_{\nu}^{\ \sigma}-\frac{1}{4}g_{\mu\nu}B^{\alpha}_{\ \sigma}B^{\sigma}_{\ \alpha}-V g_{\mu\nu}+2V^{\prime}B_{\mu}B_{\nu}.
\end{eqnarray}

\subsection{Solving the connection equation}

Upon the identification
\begin{equation}\label{auxiliarymetric}
\sqrt{-h}h^{\mu\nu}=\sqrt{-g}(g^{\mu\nu}+\xi B^{\mu}B^{\nu}),
\end{equation}
one can follow \cite{BeltranJimenez:2017doy} to show that Eq.(\ref{connectioneq}) is equivalent to Eq.(\ref{conneqRBG}) up to an irrelevant projective mode. Thus, as explained above, the connection can be written as the Levi-Civita connection of $h_{\mu\nu}$ . In order to find the explicit relation between the space-time metric $g_{\mu\nu}$ and the metric $h_{\mu\nu}$,  let us rewrite (\ref{auxiliarymetric}) in matrix form as
\begin{equation}\label{auxiliarymatrix}
\sqrt{-h}\hat{h}^{-1}=\sqrt{-g}\hat{g}^{-1}\bigg(\hat{I}+\xi \hat{BB}\bigg),
\end{equation}
where the hat denotes a matrix, such that $\hat{h}^{-1}$ and $\hat{h}$ are the matrix representations of $h^{\mu\nu}$ and  $h_{\mu\nu}$, respectively. Taking the determinant of (\ref{auxiliarymatrix}) we find that $
{h}={g}\det{({I}+\xi {BB})}$ (\ref{auxiliarymetric}), and plugging this back into (\ref{auxiliarymetric}) we arrive at
\begin{equation}
h^{\mu\nu}=\frac{1}{\sqrt{\det{({I}+\xi {BB})}}}\bigg(g^{\mu\nu}+\xi B^{\mu}B^{\nu}\bigg) \ .
\label{cch}
\end{equation}
By inverting the above relation, we also have that
\begin{equation}
h_{\mu\nu}=\sqrt{\det{({I}+\xi {BB})}}\bigg(g_{\mu\nu}-\frac{\xi}{\det{({I}+\xi {BB})}}B_{\mu}B_{\nu}\bigg).
\label{ch}
\end{equation}
Using that $\det{({I}+\xi {BB})}=1+\xi X$, with $X\equiv B^\mu B_\mu$, we can finally write
\begin{eqnarray}
h^{\mu\nu}&=&\frac{1}{\sqrt{1+\xi X}}\bigg(g^{\mu\nu}+\xi B^{\mu}B^{\nu}\bigg),\label{cch1}\\
h_{\mu\nu}&=&\sqrt{1+\xi X}\bigg(g_{\mu\nu}-\frac{\xi}{1+\xi X}B_{\mu}B_{\nu}\bigg).
\label{ch1}
\end{eqnarray}
From this last result, one finds that
\begin{equation}\label{eq:gmn_0}
g_{\mu\nu}=\frac{1}{\sqrt{1+\xi X}}h_{\mu\nu}+\frac{\xi}{1+\xi X}B_{\mu}B_{\nu} \ ,
\end{equation}
which provides an algebraic relation between $g_{\mu\nu}$ with $h_{\mu\nu}$ and $B_\mu$, though the scalar $X\equiv g^{\mu\nu}B_\mu B_\nu$ still contains an explicit dependence on $g^{\mu\nu}$. This dependence can be eliminated by noting that $Y\equiv h^{\mu\nu}B_\mu B_\nu=X\sqrt{1+\xi X}$, which completely solves the problem. Hence, we can integrate $g_{\mu\nu}$  out of the action by performing a field-redefinition  in terms of $h_{\mu\nu}$ and $B_\mu$. This is useful in order to physically interpret the different elements that contribute to the space-time metric $g_{\mu\nu}$. In fact,
using \ref{eq:gmn_0}, we can write $g^{\mu\nu}$ in terms of $h^{\mu\nu}$ and $B_\mu$ in $\mathcal{U}_B$ and $\mathcal{S}_m$, thus finding the Einstein-frame representation of the action (\ref{bumblebee}), given by
\begin{eqnarray}\label{EinsteinBumblebee}
\tilde{\mathcal{S}}_{BEF}&=&\int d^{4} x \sqrt{-h}\frac{1}{2 \kappa^{2}} R(h)+\overline{\mathcal{S}}_{m}\left(h_{\mu \nu}, B_{\mu}, \Psi\right),
\end{eqnarray}
with
\begin{align}
\nonumber&\overline{\mathcal{S}}_{m}\left(h_{\mu \nu}, B_{\mu}, \Psi_i\right)\equiv{\mathcal{S}}_{m}\big(g^{\mu\nu}(h_{\mu \nu},X,B_\mu), B_{\mu}, \Psi_i\big)-\\
\nonumber&-\frac{1}{4} \bar{B}^{\mu \nu} B_{\mu \nu}-V\left(\bar{B}^{\mu} B_{\mu} \pm b^{2}\right)+\bar{V}\left(B_{\mu}, X, h_{\mu \nu}\right)
\end{align}
where any barred tensor indicates that its indices are raised with  $h^{\mu\nu}$.

In this new (Einstein frame) representation it becomes apparent that the bumblebee model can be interpreted as GR coupled to a modified matter sector in which all the matter fields couple to the bumblebee, which also presents new self-interactions encoded in the $\bar V$ term. According to this, the metric $h_{\mu\nu}$ satisfies the Einstein equations coupled to a highly non-linear matter sector. This means that $h_{\mu\nu}$ will depart from the Minkowski metric only in regions where the Newtonian and post-Newtonian effects are expected to be relevant, {\it i.e.} regions with a strong gravitational field. As a result, as it follows from (\ref{eq:gmn_0}), the metric $g_{\mu\nu}$ will not only describe the two propagating degrees of freedom of the gravitational field through $h_{\mu\nu}$, but it will also encode information on the local value of the bumblebee field via a conformal factor and a disformal term proportional to $B_{\mu}B_{\nu} $.

From the decomposition (\ref{eq:gmn_0}) and the fact that the independent connection satisfies $\nabla^\Gamma_\alpha h_{\mu\nu}=0$, it follows that the non-metricity tensor $Q_{\alpha\mu\nu}=\nabla^\Gamma_\alpha g_{\mu\nu}$ is {non-trivial} and entirely due to the derivatives of the bumblebee field. Since this field is expected to have a non-trivial VEV that spontaneously breaks Lorentz invariance, this is an example of a gravitationally generated non-metricity tensor that can develop a VEV. In contrast, in RBGs with minimally coupled matter, the non-metricity is associated to derivatives of the stress-energy tensor of the matter fields, which vanish in vacuum. A constant background of non-metricity was assumed in \cite{Foster:2016uui}, and experimental constraints to all its possible effective couplings to fermions and photons were derived from Lorentz violation searches in Earth laboratories. Since minimally coupled matter fields do not couple explicitly to non-metricity\footnote{Indeed as we will see later, although the new couplings that arise in the matter sector have a relation with non-metricity, rather than coupling explicitly to the non-metricity tensor of the theory, the matter fields in this model couple to the source of non-metricity instead.}, these constraints do not apply to our model. However, we note that constraints on Lorentz-violating couplings such as those in the Standard Model Extension  \cite{MISSING} could translate into constraints on the bumblebee non-minimal coupling $\xi$. Further work in this direction is currently in progress.

\subsection{Metric and bumblebee equations}

Let us now continue with the exploration of the field equations (\ref{Riccieq}) and (\ref{bumblebeeeq}). By taking the trace of (\ref{Riccieq}), a relation between the scalar curvature and the trace of the stress-energy tensor $T=g^{\mu\nu}T_{\mu\nu}$ can be found { in the form}
\begin{equation}
R(\Gamma)=-\kappa^2 T,\label{Rs}
\end{equation}
which exactly matches the relation in GR. By contracting with one and two $B^{\mu}$ fields, we can also find the relations
\begin{align}
\nonumber &B^{\mu}R_{\mu\nu}(\Gamma)=\frac{1}{1+\xi X}\bigg(\kappa^2B^{\mu}T_{\mu\nu}-\\
&-\frac{1}{2}\left(\kappa^2T+\xi R_{\alpha\beta}B^\alpha B^\beta\right)B_{\nu}\bigg),\label{Rs3}\\
&B^{\mu}B^{\nu}R_{\mu\nu}(\Gamma)=\frac{\kappa^2}{2+3\xi X}\big(-TX+2B^{\mu}B^{\nu}T_{\mu\nu}\big)\label{Rs2}
\end{align}
 respectively, where $X\equiv g^{\mu\nu}B_{\mu}B_{\nu}$. { Substituting} the second of the above equations into the first {one}, we finally obtain
\begin{eqnarray}
\nonumber B^{\alpha}R_{\alpha\nu}(\Gamma)&=&\frac{\kappa^2}{1+\xi X}\bigg(B^{\alpha}T_{\alpha\nu}-\frac{1}{2+3\xi X}\bigg[(1+\xi X)T+\\
&+&\xi T_{\alpha\beta}B^\alpha B^\beta\bigg]B_{\nu}\bigg).\label{Rs4}
\end{eqnarray}

The above results can be plugged back into the metric field equations (\ref{Riccieq}) to obtain an expression for $R_{\mu\nu}(h)$ which only involves $T_{\mu\nu}$, $B_\mu$ and the metric $g_{\alpha\beta}$ which, recall, can be explicitly written in terms of $h_{\mu\nu}$ and $B_\mu$. Therefore, the field equations for $h_{\mu\nu}$ can be written in Einstein like form, though their explicit form is cumbersome and does not bring any useful new insight. This confirms that we can interpret the auxiliary metric as we did when (\ref{EinsteinBumblebee}) was introduced,  as the metric that  accounts for the cumulative effects of mass and energy.

Regarding  the bumblebee field equation, using (\ref{Rs4}) we are able to get rid of the $B^\nu R_{\mu\nu}$ term in (\ref{bumblebeeeq}), thus arriving at
\begin{eqnarray}
\nonumber\nabla_{\mu}^{(g)} B^{\mu\nu}&=&\left(2V^{\prime}+\frac{\xi T}{2+3\xi X} + \frac{\xi^2 B^{\alpha}B^{\beta}T_{\alpha\beta}}{(1+\xi X)(2+3\xi X)}\right)B^{\nu}-\\
&-&\frac{\xi}{(1+\xi X)}B_{\mu}T^{\mu\nu},
\label{Proca}
\end{eqnarray}
This equation is richer than its equivalent in the metric case, for which the $\xi$ corrections are absent. The new terms induce modifications on the effective potential that depend on the presence of other matter fields, thus implying new phenomenology. In particular, the effective mass of the field may change in regions with high energy densities, potentially leading to large effective masses and chameleon-like screening mechanisms inside and around massive objects.

\section{Weak gravitational field, post-Minkowskian limit}

Having found the solution of the connection equation and the explicit formal form of the space-time metric, we are now ready to explore the limit of negligible curvature, focusing our attention on how the sources of non-metricity modify the effective theory (post-Minkowskian limit). For this purpose, let us focus on non-gravitational experiments on Earth's surface, so that one can safely neglect all the corrections to the Minkowski metric coming from Newtonian and post-Newtonian corrections. Thus we require\footnote{Recall that (\ref{EinsteinBumblebee}) implies that $h_{\mu\nu}$ formally satisfies a set of Einstein-like equations, which yields the Minkowski metric as solution whenever one considers elementary particle interactions rather than astrophysical problems.}  $h_{\mu\nu}\approx\eta_{\mu\nu}$ . Let us also consider $\xi$ to be a small coupling and study its perturbative effects. From the above relations between the auxiliary and the space-time metrics (\ref{cch1},\ref{ch1}), we can write
\begin{eqnarray}
\nonumber g_{\mu\nu}&=&\eta_{\mu\nu}+\xi(B_{\mu}B_{\nu}-\frac{1}{2}B^{\lambda}B_{\lambda}\eta_{\mu\nu})+\mathcal{O}(\xi^2)\\
&=&\eta_{\mu\nu}+\varepsilon_{\mu\nu},\label{efmet}
\end{eqnarray}
where $B^{\lambda}B_{\lambda}=\eta^{\mu\lambda}B_{\mu}B_{\lambda}$ and $\varepsilon_{\mu\nu}=\xi(B_{\mu}B_{\nu}-\frac{1}{2}B^{\lambda}B_{\lambda}\eta_{\mu\nu})$ $+\mathcal{O}(\xi^2)$. Here we see that even when Newtonian and post-Newtonian corrections can be neglected, the space-time metric is locally departing from its Minkowskian value due to local contributions sourced by the bumblebee field $B_\mu$. Since all fields couple to the metric, as a result, in this post-Minkowskian approximation of our theory all the matter fields will couple to the bumblebee due to the unconventional way in which the connection mediates between the geometry and the matter.

\subsection{Effective dynamics of scalar and fermionic fields.}

We will next proceed to study the effective dynamics of scalar and spinor matter fields in this post-Minkowskian scenario. Regarding spinors, it is important to note that they provide a nonzero contribution to the connection equation via torsional terms. Those terms have been omitted in our presentation of the field equations for simplicity. A more careful analysis, along the lines of \cite{Afonso:2017bxr}, justifies our choice because in the case of bosonic fields, torsion can be trivialized by a simple choice of projective gauge. For fermions, however, the torsion picks up contributions that cannot be gauged away. However, such terms do not modify the equation satisfied by the symmetric part of the connection, which is still of the form (\ref{connectioneq}) and admits the Christoffel symbols of $h_{\mu\nu}$ as solution. Therefore, the resulting effective metric will be the same as (\ref{efmet}). The new torsional terms will appear as new fermionic contact interactions on the right-hand side of (\ref{Riccieq}) and will be Planck-scale suppressed (see {\it e.g.} \cite{Hehl:1976kj} or \cite{Jimenez:2020dpn} for the RBG case), so that we can neglect them for our purposes.

From (\ref{efmet}) it is clear that the inverse effective metric is given by
\begin{equation}
\label{efmetinv}
g^{\mu\nu}=\eta^{\mu\nu}-\varepsilon^{\mu\nu},
\end{equation}
where the indices are raised and lowered by the usual Minkowski metric. Similarly, one has $\sqrt{-g}\simeq 1-\frac{\xi}{2} B^{\mu}B_{\mu}$.
As a result, the scalar Lagrangian within this approximation is given by
\begin{eqnarray}
\label{approx}
\mathcal{L}_{sc}&=&\frac{1}{2}\sqrt{-g}(g^{\mu\nu}\partial_{\mu}\Phi\partial_{\nu}\Phi-m^2\Phi^2)\nonumber\\
&=&
\frac{1}{2}\left(1-\frac{\xi}{2} B^{\rho}B_{\rho}\right)\left\{\big[\eta^{\mu\nu}-\epsilon^{\mu\nu}\big]\partial_{\mu}\Phi\partial_{\nu}\Phi-m^2\Phi^2\right\}.\nonumber\\
\end{eqnarray}
Using the definition of $\varepsilon^{\mu\nu}$ and integrating by parts the former equation, we arrive at
\begin{eqnarray}
\nonumber\mathcal{L}_{sc}&=&-\frac{1}{2}\Phi(\Box+m^2)\Phi-\frac{\xi}{2}(B^{\mu}\partial_{\mu}\Phi)^2+\frac{m^2}{4}\xi\Phi^2B^{\mu}B_{\mu}+\\
&+&\nonumber\mathcal{O}(\xi^2),\\
\nonumber &=&-\frac{1}{2}\Phi(\Box+m^2)\Phi+\frac{\xi}{2}\Phi\Big[B^{\mu}B^{\nu}\partial_{\mu}\partial_{\nu}+\\
\nonumber&+& \big(B^\mu (\partial_\nu B^\nu)+B^\nu (\partial_\nu B^\mu)\big)\partial_\mu+\frac{m^2}{2} B^2\Big]\Phi+\mathcal{O}(\xi^2),\\
\label{31}
\end{eqnarray}
where $\Box\equiv\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}$. The interaction terms between the bumblebee and the scalar field naturally stem from special features of this metric-affine model, as explained previously.
As the bumblebee field develops a non-trivial VEV, we see that the last two terms in (\ref{31}) carry coefficients for Lorentz violation, which will be discussed in more detail later.

Regarding spinor fields, we start with the Hermitian Lagrangian
\begin{equation}
\label{spin}{\cal L}_{sp}=\sqrt{-g}\Big[\frac{i}{2}e_a{}^\mu\left(\bar{\Psi}\gamma^a\nabla^{(\Gamma)}_{\mu}\Psi-(\nabla^{(\Gamma)}_{\mu}\bar{\Psi})\gamma^a\Psi\right)-m\bar{\Psi}\Psi
\Big],
\end{equation}
where $e_a{}^\mu$ is the tetrad field satisfying $e_a{}^\mu e_b{}^\nu g_{\mu\nu}=\eta_{ab}$. Given that the connection $\Gamma$ is just the Levi-Civita connection of $h_{\mu\nu}$, in the weak field regime $h_{\mu\nu}\approx\eta_{\mu\nu}$ we can approximate $\nabla^\Gamma_\mu \Psi\approx\partial_\mu\Psi$ up to Planck scale suppressed and $\mathcal{O}(\xi^2)$ torsion corrections (see {\it e.g.} \cite{Delhom:2019wir}). Taking this into account and using that
\begin{equation}
e_{a}{}^\mu=\delta_a{}^\mu-\frac{1}{2} \varepsilon_{a}^\mu=\delta_{a}{}^\mu+\frac{\xi}{2}\left(\frac{1}{2} B^{2} \delta_{a}{}^\mu-B_{a} B^\mu\right),
\end{equation}
we find that the spinor Lagrangian can be written as
\begin{eqnarray}
\nonumber\mathcal{L}_{sp}&=&\bar\Psi(i \gamma^\mu \partial_\mu-m) \Psi-\\
\nonumber&-&\frac{\xi}{2} \bar\Psi\Bigg[\frac{i}{2} B^{2} \gamma^\mu \partial_{\mu}+i B^\mu B^{\nu} \gamma_{\mu} \partial_{\nu}+\frac{i}{2}\Big(B_{\alpha}\left(\partial_{\mu} B^{\alpha}\right)+\\
\nonumber&+&B^{\nu}\left(\partial_{\nu} B_{\mu}\right)+(\partial_\alpha B^\alpha) B_{\mu}\Big) \gamma^{\mu}-m B^{2}\Bigg] \Psi+\mathcal{O}(\xi^2).\\
\label{SpinorLagPert}
\end{eqnarray}
Before moving forward, it is worth calling attention to the fact that the combination of all the complex terms present in Eq. (\ref{SpinorLagPert}) result in a Hermitian Lagrangian as a consequence of Eq. (\ref{spin}).  Similarly to the scalar sector, the $\mathcal{O}(\xi)$ terms in the above Lagrangian will contribute to Lorentz violation coefficients when the bumblebee acquires a VEV. We will next proceed to analyze their physical implications.

\section{Lorentz-violating coefficients}

Let us now outline the physical effects related to a spontaneous breaking of Lorentz symmetry by the bumblebee VEV $\langle B_\mu\rangle=b_\mu$ in a weak gravitational field. Generally, observables which couple to $b_{\mu}$ will be sensitive to the spontaneous breaking of Lorentz symmetry by the bumblebee field. Since the present model displays non-minimal couplings between the bumblebee field and the matter sources through the non-Riemannian part of the connection, there will arise several Lorentz-violating (LV) coefficients in the effective matter sector once the affine connection has been integrated out.

 A straightforward and simple method to examine the effects of the Lorentz violation is to consider the VEV being a fixed-norm vector or, similarly as in Minkowski space, to take all their components to be constant. Such a choice is reasonable because we are dealing with laboratory experiments. Nonetheless, for the sake of completeness, the non-constant  $b_{\mu}$ case is presented in \ref{app1}. Therefore, from now on we restrict our analysis to the standard $b_{\mu}$ constant case \cite{KosMew}. We note that the allowed values for $b_{\mu}$ have been restricted from various experiments, see {\it e.g.}\cite{datatables} and references therein. At the same time, let us note that the LV coefficients generated by the non-metric part of the connection in metric-affine bumblebee gravity are analog to some of the parameters of the LV vector sector of the Standard Model Extension. Therefore, data from table D15 in \cite{datatable} allows us to constrain the norm of $b_{\mu}$ to be less than $10^{-43}$\, GeV.

\subsection{Scalar field}
Let us first study the Lagrangian for a scalar field propagating on top of the bumblebee background, which takes the form
\begin{eqnarray}
\nonumber\mathcal{L}_{sc}&=&-\frac{1}{2}\Phi(\Box+m^2)\Phi+\frac{\xi}{2}\Phi\Big[(s^{\mu\nu}\partial_{\mu}\partial_{\nu})+\frac{1}{2}m^2 b^2\Big]\Phi+\\
&+&\mathcal{O}(\xi^2),\label{ScalarLagBackground}
\end{eqnarray}
where $s^{\mu\nu}=\xi b^{\mu}b^{\nu}$. The $\mathcal{O}(\xi)$ terms will typically induce LV coefficients through the VEV of the bumblebee field. The $s^{\mu\nu}$ term constitutes a modification of the standard kinetic term which can be encoded in an effective metric for the scalar field of the form $g_{\rm eff}^{\mu\nu}=\eta^{\mu\nu}-\xi s^{\mu\nu}$. Hence, a ``wrong'' signature of the LV coefficient $s^{\mu\nu}$ could trigger ghost-like instabilities around strong enough bumblebee backgrounds. Note however that in that case, the perturbative expansion would break down since $\xi b^2$ would be $\mathcal{O}(1)$, and a full non-perturbative analysis would be required. The correction to the mass term in (\ref{31}) can also be encoded in an effective mass of the form $m_{\rm eff}^2=m^2(1-(\xi/2) b^2)$ which could also trigger tachyonic-like instabilities for a space-like bumblebee VEV (again non-perturbative effects could play a non-negligible role).

In order to explore potential instabilities in more detail, let us analyze the particular cases of a time- and space-like constant $b_\mu$ by working out their respective dispersion relations. Starting with a time-like VEV $b_{\mu}=[b,0,0,0]$ we find the dispersion relation
\begin{equation}\label{disprelscalartimelike}
E^2=\vec{p}^2+m^2+\xi b^2 \left(\vec{p}^2+\frac{1}{2}m^2\right)+\mathcal{O}(\xi^2).
\end{equation}
This dispersion relation is healthy for positive values of $\xi$ but could potentially develop ghost- and tachyonic-like instabilities for negative values of $\xi$, although in this case the higher-order terms might be relevant to the discussion. Considering a space-like choice for $b_\mu=[0,\vec{b}]$ we end up with
\begin{equation}\label{disprelscalarspacelike}
E^2=\vec{p}^2+m^2+\xi \left(\frac{1}{2} m^2 b^2+(\vec{b}\cdot\vec{p})^2\right) +\mathcal{O}(\xi^2).
\end{equation}
Again, this dispersion relation is healthy for positive values of $\xi$. For negative values of $\xi$ a tachyonic-like instability as well as a ghost-like instability  (in directions which are non-orthogonal to $\vec{b}$) could potentially arise, and again the higher-order terms might be relevant to the discussion if they appear.

\subsection{Dirac field}

Let us now turn our attention to the spin $1/2$ fields. To explore the physics of our interest in a more convenient way, we will work with the decomposition of LV coefficients that is more commonly used in the literature \cite{Colladay:1996iz,Colladay:1998fq}. To that end, let us rewrite the weak-field spinor action \ref{SpinorLagPert} as follows
 \begin{equation}
\mathcal{L}_{sp}=\bar{\Psi}\left(i\Gamma^{\mu}\partial_{\mu}- M\right)\Psi,
\label{SpinorLagShort}
\end{equation}
where ${\Gamma}^{\mu}$ and ${M}$ are elements of the 16-dimensional Clifford algebra defined by the Dirac gamma matrices. We can thus expand them in the usual basis of this algebra as
\begin{equation}
\begin{split}
{\Gamma}^{\mu}&={e}^{\mu}I+(\delta^\mu{}_\alpha+{c}^{\mu}{}_{\alpha})\gamma^{\alpha}+{d}^{\mu}{}_{\alpha}\gamma_{5}\gamma^{\alpha}+i{f}^{\mu}\gamma_{5}+\\
&+\frac{1}{2}{g}^{\mu}{}_{\lambda\alpha}\sigma^{\lambda\alpha},\\
{M}&={m}_{\rm eff}I+{a}_{\mu}\gamma^{\mu}+{k}_{\mu}\gamma^{\mu}\gamma_{5}+\frac{1}{2}l_{\mu\nu}\sigma^{\mu\nu},
\label{LVSpinorDecomp}
\end{split}
\end{equation}
where ${c}^{\mu}{}_{\alpha}$, ${d}^{\mu}{}_{\alpha}$, ${e}^{\mu}$, ${f}^{\mu}$, ${g}^{\mu}{}_{\lambda\alpha}$, ${a}_{\mu}$, ${k}_{\mu}$ and $l_{\mu\nu}$ are LV coefficients. Comparing \ref{SpinorLagPert} to \ref{SpinorLagShort} and \ref{LVSpinorDecomp}, we find the non-zero LV coefficients
\begin{eqnarray}
{c}^{\mu}{}_{\alpha}&=&-\frac{\xi}{2}\Big(\frac{1}{2}b^2\delta^{\mu}{}_\alpha+b^\mu b_\alpha\Big),\\
{m}_{\rm eff}&=&m\left(1-\frac{\xi}{2} b^2\right).
\label{coeff}
\end{eqnarray}
We see that within the metric-affine bumblebee model, the LV coefficients that appear provide a modification of the fermionic mass through $m_{\rm eff}$ {and} a modification of the standard kinetic term through ${c}^{\mu}_{\alpha}$.

 In general, these will introduce modifications in the dispersion relation of spin $1/2$ fields. To that end, we first notice that the modified Dirac equation is
 \begin{equation}\label{SpinorEqShort}
\left(i\Gamma^{\mu}\partial_{\mu}- M\right)\Psi=0\ ,
\end{equation}
 and multiplying on the left by $\left(i\Gamma^{\mu}\partial_{\mu}+M\right)$ we arrive at
\begin{align}
\begin{split}
\big(-\Gamma^\mu\Gamma^\nu\partial_{\mu}\partial_{\nu}+i[M,\Gamma^\mu]\partial_{\mu}-M^2\big) \Psi=0 \ .
\end{split}
\end{align}
By using now the relations
\begin{align}
&\{\Gamma^\mu,\Gamma^\nu\}=2\eta^{\mu\nu}-\xi(b^2\eta^{\mu\nu}+2b^\mu b^\nu)+\mathcal{O}(\xi^2),\\
&[ M,\Gamma^\mu]=\mathcal{O}(\xi^2),\\
& M^2=m^2\left(1-\xi b^2\right)+\mathcal{O}(\xi^2),
\end{align}
we find the following dispersion relation
\begin{align}
\nonumber&\bigg[E^{2}\big(1-\xi (b^{2}+b_0^2)\big)+ 2\xi b_{0} E (\vec{b} \cdot \vec{p}) -\\
&-\left(\vec{p}^{2}+m^{2}\right)\left(1-\xi b^{2}\right)-\xi(\vec{p} \cdot \vec{b})^{2}\Bigg] \Psi=0 \ .
\end{align}
 As we did for the scalar field, let us particularize for a constant bumblebee background of both time- and space-like types. Starting with a time-like VEV $b_{\mu}=[b,0,0,0]$ we get the dispersion relation
\begin{eqnarray}
E^{2}=(1+\xi b^2)(\vec{p}^2+m^2)+\mathcal{O}(\xi^2).
\label{Spinordispreltimelike}
\end{eqnarray}
While for positive values of $\xi$ this dispersion relation is perfectly well-behaved, for negative values of $\xi$ there could be instabilities if $b^2\xi\geq1$, but higher order terms would be non-negligible in this case and \ref{Spinordispreltimelike} would not be trustable anymore. For a (constant) space-like vector $b_{\mu}=[0,\vec{b}]$ we obtain
\begin{eqnarray}
E^{2}=\vec{p}^2+m^2-\xi (\vec{p}\cdot\vec{b})^2+\mathcal{O}(\xi^2) \ .
\end{eqnarray}
For negative values of $\xi$ this is well behaved, though if it is positive then ghost-like instabilities could arise beyond the perturbative level ({see} \cite{Kostelecky:2000mm} for a discussion on the typical energy scales at which these instabilities become relevant).

\section{Summary and conclusions}

In this work we have formulated the bumblebee model within the metric-affine formalism.  By solving the equations of motion of the connection, we have been able to express  the metric as a function of this field and of an auxiliary metric which accounts for the standard effects of the gravitational interaction. To our knowledge, this is the first time that a solution for the connection in a curvature-based metric-affine gravity theory with spontaneously broken Lorentz symmetry has been found. The methods used to solve the connection are analog to those commonly employed in Ricci-based gravity theories (see {\it e.g.} \cite{Afonso:2017bxr,BeltranJimenez:2017doy,Bernal:2016lhq}), but its qualitative properties are rather different from those theories. In particular, while in RBGs the metric picks up local corrections that depend on the stress-energy tensor of the matter fields, here those corrections are entirely determined by the bumblebee field itself. This { occurs} due to the {absence} of higher-curvature terms in the action and to the non-minimal coupling of the bumblebee to gravity via the Ricci tensor.

From our analysis it follows that, unlike in RBGs, where the non-metricity is given by gradients of the energy-momentum density, in our model it is the gradient of the bumblebee field that generates a non-vanishing non-metricity. The fact that the bumblebee has a VEV which breaks Lorentz symmetry allows for a background non-metricity that could fit with the proposal in \cite{Foster:2016uui}. This is the first gravitationally-induced non-metricity model with a VEV that we are aware of.

An immediate effect of the metric dependence on the bumblebee field is that all the matter fields couple to it. Since this coupling is not a gauge one it can generate, for instance, the coupling of a vector to a neutral scalar. In the weak-field limit,  we have seen that the resulting theory looks like a bumblebee coupled to matter with non-linear interactions in Minkowski space. Therefore, this theory is naturally treated as an effective theory where the role of the energy scale is played by $\xi^{-1/2}$ (see {\it e.g.} \cite{Georgi}). Typically, for the VEV solution we have checked that the theory can present ghost and tachyonic-like instabilities depending on the range of the non-minimal parameter $\xi$ and the explicit form of the VEV.
{ A natural continuation of this study consists in calculating quantum corrections in this effective theory.} The search for astrophysically relevant solutions beyond the weak-field limit presented here is also another research avenue currently under consideration. \\

\begin{acknowledgements}
This work was partially supported by Conselho
Nacional de Desenvolvimento Cient\'{\i}fico e Tecnol\'{o}gico (CNPq) and by the Spanish projects FIS2017-84440-C2-1-P (MINECO/FEDER, EU), i-LINK1215 (CSIC), and the projects SEJI/2017/042 and PROMETEO/2020/079 (Generalitat Valenciana). The work by A. Yu. P. has been supported by the CNPq project No. 301562/2019-9. PJP would like to thank the Brazilian agency CAPES for financial support (PDE/CAPES grand, process 88881.17175/2018-01) and Department of Physics and Astronomy, University of Pennsylvania, for the hospitality. A. D. is supported by a PhD contract of the program FPU 2015 (Spanish Ministry of Economy and Competitiveness) with reference FPU15/05406, and would like to thank the Departamento de F\'{\i}sica, Universidade Federal da Para\'{\i}ba for hospitality.
\end{acknowledgements}


text:
\section{}

\section{Introduction}
Within a quantum-thermodynamical approach~\cite{Gemmer2009} a single subsystem can be shown to relax to a thermal state if
appropriately embedded in a quantum environment: The pertinent temperature is determined by the spectral density of this environment,
which, for sufficiently large modular embeddings and for a given working point in energy space, shows an approximate exponential
behavior. Beyond a design of this density there appears to be no immediate possibility of further manipulation.

Putting aside explicitly time-dependent Hamiltonians, an interesting way to circumvent this limitation consists in the application of quantum measurements.
Measurements play an outstanding role in quantum mechanics, since they provide a connection between the abstract theory and
experimentally accessible quantities. However, in contrast to classical measurements,
quantum measurements usually influence the measured system. A pertinent
example for the influence of periodic measurements is the well-known Zeno effect~\cite{Misra1977,Itano1990}, i.e., the suppression
of decay by fast repeated measurements. Recently, also
the possibility of cooling a TLS coupled to an oscillator bath by applying periodic quantum non demolition measurements
has been discussed~\cite{Erez2008,Gordon2009,Bensky2010,Gordon2010}.

In a recent Letter~\cite{Jahnke2010} we discussed the influence of periodic measurements in a quantum-thermodynamical
setting and pointed out how different concepts of statistical mechanics may arise from this ``observed quantum thermodynamics''.
Here, we present a detailed discussion of the effect of periodic measurements in a generalized environment model, allowing for a
detuning between system and environment. As we will see, this permits us to prepare the embedded system in almost arbitrary thermal
states --- even the creation of inversion becomes possible.

\section{Model \label{sec:model}}
A typical quantum-thermodynamical model consists of the system of interest and its environment, with which it is allowed to exchange
energy via a weak coupling.
The model we investigate here belongs to the class of systems with a modular environment.
Following~\cite{Jahnke2010} we consider a system $S$ (a TLS with energy splitting $\dS$) described by the Pauli operator $\op{\sigma}_z$ and coupled via $\op{V}$ (strength $\lambda$) to a quantum environment $B$ with Hamiltonian $\op{H}_B$,
\begin{equation}
    \op{H}_\text{tot}=\frac{\delta}{2}\op{\sigma}_z\otimes \op{1}_B+\op{1}_S\otimes\op{H}_B+\lambda\op{V}\,.
\end{equation}

The environment consists of $n$ spins (e.g., a paramagnetic salt) with energy splitting $\dB\equiv\dS+\Dd$ each, i.e., we allow for
a detuning $\Dd$ between system and environment.
We will show below that this detuning $\Dd$ is key for preparing the TLS in almost arbitrary thermal states.

Without spin-spin interaction the degeneracy of the energy level $E_k=\dB k$ with $k$ spins up is given by $N_k={{n}\choose{k}}$.
Around a given $k_0$, the working point in energy space, the degeneracy structure can be approximated by an exponential of the form
\begin{equation}
\label{eq:degeneracy}
 N_k\approx\mathcal{N}e^{\beta E_k}\,.
\end{equation}
For $n>k_0\gg1$ one finds (see Appendix \ref{sec:thermalizingenvironment})
\begin{equation}
\label{eq:beta}
 \beta(k_0)\approx\frac{1}{\dB}\ln\left( \frac{n}{k_0}-1\right)\,.
\end{equation}
Including a weak interaction within the environment, each energy level broadens into an energy band of width $\Delta \varepsilon_k\ll\dB$.

For the interaction between system and environment we choose
\begin{equation}
\label{eq:V}
 \op{V}=\left( \op{\sigma}^{+}+\op{\sigma}^{-}\right) \otimes \op{B}\,,
\end{equation}
where $\op{\sigma}^{+}=\ketbra{1}{0}$ is the creation operator and $\op{\sigma}^{-}=\ketbra{0}{1}$ the annihilation operator for the TLS
and the environmental part reads
\begin{eqnarray}
\label{eq:B}
 \op{B}=\sum_{k}\sum_{n_k,m_{k+1}}&&C_{k+1,k}(m_{k+1},n_{k})\ketbra{n_k}{m_{k+1}}\nonumber\\
&\quad+&C_{k,k+1}(n_k,m_{k+1})\ketbra{m_{k+1}}{n_k}\,.
\end{eqnarray}
To keep the model as unbiased as possible, we do not further specify the interaction but consider the $C_{i,j}$ to be a set of
random matrices whose entries are taken from a Gaussian distribution and are normalized to
$\overline{|C_{i,j}(a,b)|^2}=(N_i N_j)^{-1/2}$. This choice of the normalization together with a sufficiently small $\lambda$
guarantees the weak coupling for arbitrary states.
Due to Hermiticity of $\op{V}$, we have $C_{i,j}=C_{j,i}^\dagger$.

Initially we assume a product state for the total system $\op{\rho}_\text{tot}(0)=\op{\rho}(0)\otimes\op{\rho}_{B}^0(k_0)$,
where the state for the environment $B$ is given by
\begin{equation}
 \label{eq:rho0B}
   \rB^0(k_0)=\frac{1}{N_{k_0}}\sum_{n_{k_0}}\ketbra{n_{k_0}}{n_{k_0}}\,,
\end{equation}
i.e., only band $k_0$ is occupied.

Under the above conditions and for $\Dd=0$ the subsystem $B$ acts as a thermalizing environment on the TLS:
Subject to undisturbed Schr\"odinger evolution of the total system (i.e., without measurements), the TLS would relax
into a state with inverse temperature $T^{-1}=\beta(k_0)$ given by~(\ref{eq:beta}) independently of its initial state
 $\op{\rho}(0)$~\cite{Gemmer2009, Breuer2006, Gemmer2006}.

In the following, we want to study the effect of external observation (by means of periodic measurements) for this
quantum-thermodynamical resonant case as well as for the case of a detuned environment (i.e., $\dB\neq\dS$).

\section{Effective environment}
Without measurements, the environmental spins have to be in resonance with the TLS, i.e., the detuning has to vanish ($\Dd=0$)
in order to allow for a thermalization of the TLS. However, this no longer needs to be the case if we include
periodic measurements.
We thus analyze the dynamics of the system with arbitrary detuning $\Dd$ disturbed by periodic measurements of the environmental
magnetization, i.e., the system first evolves under pure Schr\"odinger dynamics for some time $\Dt$. Then a measurement of the
environmental magnetization (or rather its energy) is executed, which means we perform a projection on one of the energy bands $k_1$.
The pertinent projection operator reads
\begin{equation}
 \Ptot{k_1}=\op{1}_S\otimes\P{k_1}\,
\end{equation}
with
\begin{equation}
 \label{eq:projectionoperator}
 \P{k_1}=\sum_{n_{k_1}}\ketbra{n_{k_1}}{n_{k_1}}\,,
\end{equation}
where we sum over all levels of the measured energy band $k_1$.
This leads to a state of the TLS given by
\begin{equation}
 \label{eq:rhocojump}
 \rS(1)=\frac{\Tr[B]{\Ptot{k_1}\rtotbef(1)}}{\Tr[B]{\P{k_1}\rBbef(1)}}\,,
\end{equation}
which can be rewritten as
\begin{equation}
 \label{eq:rhocojump2}
 \rS(1)=\rSbef(1)+\frac{\Tr[B]{\Ptot{k_1} \op{C}_{SB}(1)}}{\Tr[B]{\P{k_1}\rBbef(1)}}\,,
\end{equation}
where $\rSbef$ denotes the state of the TLS before the measurement and the second term is the so-called
co-jump~\cite{Granzow1998} caused by the correlations  $\op{C}_{SB}:=\rtotbef-\rSbef\otimes\rBbef$
between system and environment straight before the measurement. Since the measurement is incomplete (measurement of an energy band
and not of a single level), there may still be some correlations left after the measurement. However, the effect of these correlations
is negligible in the case of weak interaction and small band width.
That is, the state of the total system after the measurement can be approximated by a product state. For simplicity, we introduce a
further approximation: Since we are not interested in details like the occupation probabilities of the single levels within the
energy bands of the environment, we apply some kind of coarse graining by replacing the typically rather complicated state of
the environment after the measurement by a simple one of the form~(\ref{eq:rho0B}). This can be done, provided the environment is
sufficiently large.
Thus, the new state of the total system after the measurement is approximated by
\begin{equation}
 \rtot(1)\approx\rS(1)\otimes\rB^0(k_1)\,,
\end{equation}
i.e., a form similar to the initial total state but with a new state $\rS(1)$ for the system and possibly a different energy band being
occupied for the environment (if a different energy band has been measured).
This procedure can now be repeated periodically. In order to analyze what happens with the TLS due to these measurements, we first
have to calculate the state of the TLS after a single measurement~(\ref{eq:rhocojump}).
This can be done by means of perturbation theory~(see Appendix~\ref{App:derivation}), which leads to three possible states of the TLS
corresponding to the three possible measurement results:

In the case of measuring the same energy as in the previous measurement, we obtain the occupation probability of the ground state up
to second order in the interaction strength $\lambda$:
\begin{equation}
\label{eq:rho000expansion}
 \rc{00,0}(j)\approx\rc{00}(j-1)\left( 1-4\lambda^2\rc{11}(j-1)(\textup{e}^{-\beta\dB/2}-\textup{e}^{\beta\dB/2})\left( \sinA-\sinB\right)\right)
\end{equation}
with
\begin{eqnarray}
 \sinA&:=&\frac{\sin^2\left(\frac{\Dd}{2}\Dt\right)}{\Dd^2}\,,\\
 \sinB&:=&\frac{\sin^2\left(\left(\dS+\frac{\Dd}{2}\right)\Dt\right)}{(2\dS+\Dd)^2}\,.
\end{eqnarray}
In case of measuring one energy band higher than before, we get
\begin{equation}
\label{eq:rho00p}
 \rc{00,+}(j)\approx\frac{\rc{11}(j-1)\sinA}{\rc{11}(j-1)\sinA+\rc{00}(j-1)\sinB } \,,
\end{equation}
whereas measuring one energy band lower leads to
\begin{equation}
\label{eq:rho00m}
 \rc{00,-}(j)\approx\frac{\rc{11}(j-1)\sinB}{\rc{11}(j-1)\sinB+\rc{00}(j-1)\sinA } \,.
\end{equation}
For the off-diagonal elements we get accordingly
\begin{eqnarray}
 \label{eq:rho100expansion}
  \rc{10,0}(j)&\approx&\rc{10}(j-1)\Bigg[1+\lambda^2\Bigg(\left(\textup{e}^{-\beta\dB/2}\rc{00}(j-1)+\textup{e}^{\beta\dB/2}\rc{11}(j-1)\right)4\sinA\nonumber\\
&&+\left(\textup{e}^{\beta\dB/2}\rc{00}(j-1)+\textup{e}^{-\beta\dB/2}\rc{11}(j-1)\right)4\sinB\nonumber\\
&&-\left(\textup{e}^{-\beta\dB/2}+\textup{e}^{\beta\dB/2}\right)\nonumber\\
&&\times\left(\frac{1-\textup{e}^{\textup{i}\Dd\Dt}+\textup{i} \Dd\Dt}{\Dd^2}+\frac{1-\textup{e}^{\textup{i}(2\dS+\Dd)\Dt}+\textup{i} (2\dS+\Dd)\Dt}{(2\dS+\Dd)^2}\right)\Bigg)\Bigg]\,.\\
\label{eq:rho10p}
 \rc{10,+}(j)&\approx&\rho_{01}(j-1)\frac{1+\textup{e}^{2\textup{i} \dS \Dt}-2\textup{e}^{\textup{i}\dS \Dt}\cos((\dS+\Dd)\Dt)}{4(2\dS\Dd+\Dd^2)\left(\rc{11}(j-1)\sinA+\rc{00}(j-1)\sinB\right)}\,,\\
\label{eq:rho10m}
 \rc{10,-}(j)&\approx&\rho_{01}(j-1)\frac{1+\textup{e}^{2\textup{i} \dS \Dt}-2\textup{e}^{\textup{i}\dS \Dt}\cos((\dS+\Dd)\Dt)}{4(2\dS\Dd+\Dd^2)\left(\rc{11}(j-1)\sinB+\rc{00}(j-1)\sinA\right)}\,.
\end{eqnarray}

As one can see, measuring a different band causes significant changes in the occupation probabilities of the TLS: For example, if $\rho_{00}$ was
close to $1$ at the previous measurement it will be almost $0$ after such a measurement, and vice versa. Thus, one cannot expect
that a single TLS will relax to a stable thermal attractor state as is the case without measurements in the quantum thermodynamical model.
Indeed, after sufficiently many measurements the TLS will always be in its ground or its excited state, i.e., passing through
some kind of quasi-classical trajectory.

Instead, one may ask what happens to the TLS in ensemble average.
Therefore, we first have to weight each possible outcome with the probability to obtain the corresponding measurement result.
These probabilities are given by~(\ref{eq:TrPkrhoB}):
\begin{equation}
  p_{k_{j}}(j)=\Tr[B]{\P{k_j}\rBbef(j)}\,,
\end{equation}
which equals the occupation probabilities of the bands $k_{j}=\{k_{j-1},k_{j-1}+1,k_{j-1}-1\}$ before the measurement.

Thus, for a given energy band $k_{j-1}$ determined at the preceding measurement $j-1$, the probabilities for measuring one band higher, one band
lower, or the same band, respectively, read
\begin{eqnarray}\label{eq:probabilitybandup}
 p_{+}(j)&\approx&4\lambda^2\textup{e}^{\beta\dB/2}\left(\rho_{11}(j-1)\sinA+\rho_{00}(j-1)\sinB\right)\,,\\
 \label{eq:probabilitybanddown}
 p_{-}(j)&\approx&4\lambda^2\textup{e}^{-\beta\dB/2}\left(\rho_{00}(j-1)\sinA+\rho_{11}(j-1)\sinB\right)\,,\\
 \label{eq:probabilitybandequal}
 p_{0}(j)&\approx&1- p_{+}(j)- p_{-}(j)\,.
\end{eqnarray}
This yields the ensemble average (indicated by an overbar) for the diagonal elements,
\begin{eqnarray}
\label{eq:rho00short}
 \avrc{00}(j)&=&p_0(j) \rc{00,0}(j)+p_+(j) \rc{00,+}(j)+ p_-(j) \rc{00,-}(j)\nonumber\\
 &\approx&\rc{00}(j-1)+4\lambda^2\Big[ \left( \textup{e}^{\beta\dB/2} \rc{11}(j-1)-\textup{e}^{-\beta\dB/2}\rc{00}(j-1)\right) \sinA\nonumber\\
&&+\left( \textup{e}^{-\beta\dB/2} \rc{11}(j-1)-\textup{e}^{\beta\dB/2}\rc{00}(j-1)\right)\sinB\Big]\, .
\end{eqnarray}
We can now iterate this
result to get the ensemble-averaged density matrix as a function of the number of measurements, i.e., the ensemble average after
measurement $j$ is used as initial state for measurement $j+1$ and so forth. The calculation for the off-diagonal elements can be found
in Appendix~\ref{app:offdiagonalelements}. It turns out that the off-diagonal elements vanish for almost any choice of the parameters
after sufficiently many measurements (an interesting exception is discussed below), i.e., the attractor state of the TLS will be a thermal state.
We note in passing that the ensemble average coincides with the infinite-time average (ergodicity).

\section{Effective relaxation controlled by $\Dt$ and $\Dd$}
To derive the analytical expression for $\avrc{00}(j)$, we first rewrite~(\ref{eq:rho00short}) as
\begin{equation}
\label{eq:averagederivation1}
\avrc{00}(j+1)-\avrc{00}(j)= -R \avrc{00}(j)+d
\end{equation}
with
\begin{eqnarray}
 R&:=&8\lambda^2\cosh(\beta\dB/2)\left(\sinA + \sinB\right)\,,\\
 d&:=&4\lambda^2\left(\textup{e}^{\beta\dB/2}\sinA+ \textup{e}^{-\beta\dB/2}\sinB\right)\,.
\end{eqnarray}
The left-hand side can be approximated by a derivative with respect to $j$, as long as $\avrc{00}$ changes only slowly
with $j$, which is guaranteed by the weak coupling ($\lambda\ll1$).
Thus, we can approximate~(\ref{eq:averagederivation1}) by a differential equation of the form
\begin{equation}
 \dnd{\avrc{00}(j)}{j}=-R \avrc{00}(j)+d\,,
\end{equation}
with the solution
\begin{equation}
\label{eq:exponentialrelaxrho00}
 \avrc{00}(j)=\left(\rc{00}(0)-\frac{d}{R}\right)\textup{e}^{-R j}+\frac{d}{R}\,,
\end{equation}
This means that the ensemble of TLSs exponentially approaches an attractor state with increasing number $j$ of measurements.

This relaxation has some interesting properties: First, it is noteworthy that due to the measurements a relaxation becomes possible
not only in the resonant case (as without measurements) but also for $\dB\neq\dS$. However, the attractor state reached
as well as the relaxation constant depend on the time $\Dt$ between the measurements and on the detuning $\Dd$.

Considering the relaxation constant
\begin{equation}
\label{eq:relaxationconstant}
 R=8\lambda^2\cosh(\beta\dB/2)\left(\frac{\sin^2\left(\frac{\Dd}{2}\Dt\right)}{\Dd^2} + \frac{\sin^2\left(\left(\dS+\frac{\Dd}{2}\right)\Dt\right)}{(2\dS+\Dd)^2}\right)\,,
\end{equation}
one finds
\begin{equation}
 \lim_{\Dt\rightarrow0}R=0\,,
\end{equation}
which means that the relaxation is slowed down for very rapidly repeated measurements. Such a suppression of decay due to fast,
periodic measurements is well known as the so-called quantum Zeno effect~\cite{Misra1977,Itano1990}.
Furthermore, $R$ also vanishes for $\Dt=n\frac{\pi}{\dS}$ and $\Dd=\frac{2m\pi}{\Dt}$, $n=1, 2,\ldots$, $m=1, 2,\ldots$.
For this special choice of the parameters (characterizing the exceptions we mentioned before), the absolute value of the off-diagonal elements also stays
constant, as discussed in Appendix \ref{app:offdiagonalelements}.
Thus, it is possible to freeze the initial state of the TLS and suppress decoherence by periodic measurements carried out even
with finite frequency.

For any other choice of the parameters, however, a thermal attractor state is obtained, which also shows remarkable features.
In general, it is given by (cf.~(\ref{eq:exponentialrelaxrho00}))
\begin{eqnarray}
 \label{eq:rho00Attr}
 \attr{00}&=&\frac{d}{R}\nonumber\\
 &=&\frac{\textup{e}^{\beta\dB/2} \sinA+\textup{e}^{-\beta\dB/2}\sinB }{2\cosh(\beta\dB/2)\left(\sinA + \sinB\right) }\,.
\end{eqnarray}
Let us first study the resonant case $\Dd=0$. Because of
\begin{eqnarray}
\lim_{\Dd\rightarrow0} \sinA&=&\lim_{\Dd\rightarrow0}\frac{\sin^2\left(\frac{\Dd}{2}\Dt\right)}{\Dd^2}=\frac{\Dt^2}{4}\,,\\
\lim_{\Dd\rightarrow0} \sinB&=&\lim_{\Dd\rightarrow0}\frac{\sin^2\left(\left(\dS+\frac{\Dd}{2}\right)\Dt\right)}{(2\dS+\Dd)^2}=\frac{\sin^2\left(\dS\Dt\right)}{4\dS^2}\,,
\end{eqnarray}
we get
\begin{equation}
 \label{eq:rho00AttrRes}
 \attr{00}(\Dt)=\frac{\textup{e}^{-\beta\dS/2}\sin^2(\dS \Dt)+ \textup{e}^{\beta\dS/2} \dS^2 \Dt^2}{2\cosh(\beta\dS/2)\left(\sin^2(\dS \Dt) + \dS^2 \Dt^2\right) }\,.
\end{equation}
The lowest possible temperature for the TLS in this resonant case is obtained for $\Dt=\frac{n\pi}{\dS}$, $n=1,2,3,\ldots$ and is given by
$T_{\text{min,res}}=1/\beta$. This is just the temperature the TLS would also get according to
quantum thermodynamics without measurements (cf.~\cite{Jahnke2010}). However, any other choice of $\Dt$ will lead to a higher temperature. In particular, for
rapidly repeated measurements we obtain
\begin{equation}
  \label{eq:rho00AttrResLimit1}
  \lim_{\Dt\rightarrow0}\attr{00}=\frac{1}{2}\,,
\end{equation}
i.e., the TLS is heated up to very high temperatures ($T_\text{eff}\rightarrow\infty$).

Finally we want to focus on the effect of a detuning $\Dd\neq0$ between system and environment: Now, the
attractor state~(\ref{eq:rho00Attr}) reads
\begin{eqnarray}
 \label{eq:rho00AttrOffres}
 \attr{00}&=&\frac{\textup{e}^{\beta\dB/2} \frac{\sin^2\left(\frac{\Dd}{2}\Dt\right)}{\Dd^2}+\textup{e}^{-\beta\dB/2}\frac{\sin^2\left(\left(\dS+\frac{\Dd}{2}\right)\Dt\right)}{(2\dS+\Dd)^2} }{2\cosh(\beta\dB/2)\left(\frac{\sin^2\left(\frac{\Dd}{2}\Dt\right)}{\Dd^2} + \frac{\sin^2\left(\left(\dS+\frac{\Dd}{2}\right)\Dt\right)}{(2\dS+\Dd)^2}\right) }\,,
\end{eqnarray}
and has two new interesting properties compared to the resonant case discussed before.
The lowest possible temperature here is reached for $\Dt=n\pi/(\dS+\Dd/2)$, $n=1, 2, 3,\ldots$, and is given by
\begin{equation}
 T_\text{min}=\frac{\dS}{\dB}\frac{1}{\beta}\,.
\end{equation}
This temperature depends on the ratio of the energy splittings of the TLS and environmental spins,
that is, the off-resonance may be used for cooling the TLS down to very low temperatures by choosing proper values
for $\Dt$ and $\Dd$.

The other remarkable feature is the limit of the lowest possible occupation probability of the ground state, i.e., the highest temperature.
As in the resonant case, rapidly repeated measurements lead to $T_\text{eff}\rightarrow\infty$. However, in the off-resonant case, it is even
possible to further increase the energy of the TLS, i.e., to create inversion. The maximum inversion is reached for $\Dt=2 n\pi/|\Omega|$, $n=1, 2, 3,\ldots$, and is given by
\begin{equation}
  \label{eq:rho00AttrOffresLimit2}
  \text{min}\left(\attr{00}(\Dt)\right)=\frac{\textup{e}^{-\beta\dB/2}}{\textup{e}^{\beta\dB/2}+\textup{e}^{-\beta\dB/2}}\,.
\end{equation}
which corresponds to the negative temperature~\cite{Landsberg1959}
\begin{equation}
 T_\text{max}=-\frac{\dS}{\dB}\frac{1}{\beta}\,.
\end{equation}
Note, that this inversion originates from the influence of the periodic measurements and not from a negative temperature state
of the spin environment (the environmental $\beta$ is positive).

By varying the time $\Dt$ between the measurements, thermal states of arbitrary temperature within these limits can be prepared.
Possible values of inverse temperatures thus cover the range between $-(\dB/\dS) \beta$ and $+(\dB/\dS) \beta$, symmetrically around
zero. (In the resonant case these values are between zero and $\beta$ only.)
Fig.~\ref{fig:attractorOffres} shows the attractor for the occupation probability of the ground state as a function of time $\Dt$ and
detuning $\Dd$.
\begin{figure}
\begin{center}
\includegraphics[scale=0.8]{attractoroffres.eps}
\end{center}
\caption{Attractor for the occupation probability of the ground state $\attr{00}$ as a function of the parameters $\Dt$ and $\Dd$. States with low
temperature are white and those with high (negative) temperature are black. Parameters used: $\delta=1\eu$, $\beta=0.75\eu^{-1}$. \label{fig:attractorOffres}}
\end{figure}

\section{Numerical verification for small environments}
The discussed behavior of the TLS under periodic measurements of the environment can be obtained not only for large environments but even
for quite small ones consisting of a few spins only. For such small environments the attractor state for the TLS
depends not only on the initial state of the environment (working point) but also on the initial state of the system, since the
environmental degeneracy structure can no longer be approximated by an exponential behavior. This means that the system no
longer shows thermodynamic behavior, which requires independence of its initial state. Nevertheless, given the initial state, we
can predict the attractor state using our analytical results.

To demonstrate this, we consider a TLS with splitting $\dS$ coupled via a random $(\sigma_x\otimes\sigma_x)$-interaction
to a seven-spin environment (splitting of the spins $\dB$), with initially two spins up.
Let us first consider the resonant case $\dS=\dB=1\eu$ with $\Dt=\pi/\dS$ and the TLS initially being in its ground state.
In this case only the initially occupied environmental ``energy band'' $E_2=2\dB$ and the band beneath ($E_1=\dB$) are involved
in the dynamics (the ground state energy is set to $E_0=0$). This means that only the ratio of the degeneracies of these bands determines the final attractor state.
Thus, because $N_1={{7}\choose{1}}=7$ and $N_2={{7}\choose{2}}=21$, we get an effective inverse temperature
$\beta_\text{eff}=\ln(N_2/N_1)\dB^{-1}=\ln(3)\dB^{-1}$ and therefore expect $\attr{00}=3/4$. Indeed, this final state is obtained in
our numerical exact simulation shown in Fig.~\ref{fig:rho00measdtpiquasiattr}.
\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{rho00measdtpiquasiattr.eps}
\end{center}
\caption{(Color online) Relaxation of the TLS given by the ensemble-averaged $\avrc{00}$ under periodic measurements of its environment
(seven spins with initially two spins up).
We obtain a relaxation to the thermal attractor state with $\attr{00}=3/4$ (dashed line).
Chosen parameters: $\dS=1\eu$, $\Dd=0$, $\Dt=\pi/\delta$. \label{fig:rho00measdtpiquasiattr}}
\end{figure}

In the off-resonant case, it is possible to obtain inversion as discussed in the last section. To illustrate this, we choose
$\Dd=0.7$ and $\Dt=2\pi/\Dd$. Again, initially we start with two spins up. Since now only the second and third bands are involved, we
obtain $\beta_\text{eff}=\ln(35/21)\dB^{-1}$, which yields according to~(\ref{eq:rho00AttrOffresLimit2})
$\attr{00}=3/8$. Again, this value is verified by the simulation, as shown in Fig.~\ref{fig:rho00measdt2pioverOmega}.
\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{rho00measdt2pioverOmega.eps}
\end{center}
\caption{(Color online) Relaxation of the TLS to a state of negative temperature due to the periodic measurements of its environment (seven spins
with initially two spins up). The final state is determined by $\attr{00}=3/8$ (dashed green line). Chosen parameters: $\dS=1\eu$, $\Dd=0.7\eu$,
$\Dt=2\pi/\Omega$. \label{fig:rho00measdt2pioverOmega}}
\end{figure}

\section{Conclusion}
In this paper, we have studied the influence of an effective environment on a coupled two-level
system. We have shown analytically that periodic measurements of the environmental magnetization (energy) may lead to a freezing of the initial state, i.e., suppress decoherence,
or may (in ensemble average) yield a thermal state for the TLS, depending on the choice of the detuning and the time between the measurements. The
temperature of the attractor state reached is controlled by these two parameters, which allows for the preparation of thermal
states with almost arbitrary temperature. It is even possible to create an inversion for the TLS. The key mechanism, which underlies
this cooling or creation of inversion, is the periodic destruction of the correlations between the system and its detuned
environment. The resulting periodic reset of the interaction energy allows for a successive change of energy of the system and
environment, even if out of resonance. It is to be expected that the deviation $T_\text{eff}$ from $T$ would give way to non-thermal
attractor states in the case of larger subsystems $S$, i.e., beyond the two-level case.~\footnote{The potential generalization to
larger systems poses a challenging question, which should be addressed in future work. We expect that, in general, periodic
measurements of a thermalizing environment will only lead to the undisturbed thermal state of the system if
the repetition of measurements is sufficiently slow, such that the effect of the periodic reset of the interaction energy
is negligible.}

As demonstrated in our numerical simulations, this intriguing behavior can be observed even for quite small environments consisting
of a few spins only. In this case, the attractor exhibits dependence on the initial state of the system $S$.

In a recent paper~\cite{Alvarez2010} a nuclear magnetic resonance setup has been used for studying the influence of periodic
quantum nondemolition measurements in a related model. We think that similar experiments may also be suitable to experimentally
test the results we have derived here.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

It is intriguing to study problems of classical fluid mechanics on a relativistic level, since the character of the underlying equations may change \cite{FreistuehlerTemple}. However, as illustrated in this note, care must be taken when addressing incompressible fluids since infinite speeds of propagation may be introduced, violating the axioms of Relativity.

The study of incompressible fluids in General Relativity began with  Karl Schwarzschild's work \cite{Schwarzschild}, addressing the special case of a static fluid. A general model analogous to incompressible fluids has been introduced by Lichnerowicz \cite{Lichnerowicz}, (see also \cite{Disconzi}), by requiring that wave speeds equal the speed of light. However, this model does not capture the rigid motion of non-relativistic incompressible fluids in the relativistic setting.

Our starting point here is to assume such rigid motion by requiring the particle number density $N$ to be constant along the fluid flow,
\begin{equation} \label{incomp}
u^\sigma \partial_\sigma N = 0 ,
\end{equation}
where $u$ denotes the fluid $4$-velocity and $\partial_\sigma$ partial differentiation in some fixed coordinate system, and we always sum over repeated upper and lower indices. The relativistic Euler equations are
\begin{eqnarray}
\nabla_\mu (Nu^\mu) &=&0 \label{intro_Euler_1} \\
\nabla_\nu T^{\mu\nu} &=& 0, \label{intro_Euler_2}
\end{eqnarray}
where $T$ is the energy-momentum tensor of a perfect fluid,
\begin{equation} \label{intro_T}
T^{\mu\nu} \equiv (\rho + p) u^\mu u^\nu + p g^{\mu\nu},
\end{equation}
for $\rho$ the energy density and $p$ the pressure of the fluid. We assume a fixed background spacetime with metric tensor $g$ of signature $(-+++)$ and we raise and lower indices by contraction with $g$. We work in natural units ($c=1$) and assume the fluid $4$-velocity $u$ to be normalized to
\begin{equation} \label{normalization}
u^\sigma u_\sigma \equiv g_{\mu\nu} u^\mu u^\nu=-1.
\end{equation}
We assume that $N$, $\rho$, $u$ and $p$ are smooth (at least $C^2$) and that $N$ and $\rho$ are both strictly positive.  In our main result, Theorem \ref{Thm1}, we show that \eqref{incomp} - \eqref{intro_Euler_2} imply an elliptic equation on families of hypersurfaces, from which we expect infinite wave speeds. In light of such infinite speeds, we refer to \eqref{incomp} - \eqref{intro_Euler_2} as the incompressible \emph{covariant} Euler equations, (instead of incompressible ``relativistic'' Euler equations).

\begin{Thm} \label{Thm1}
Assume $(N,\rho,u,p)$ is a smooth solution of the incompressible {covariant} Euler equations, \eqref{incomp} - \eqref{intro_Euler_2}, such that $\rho\equiv \rho_0$ is constant and $p+\rho>0$. Then, on each of the spacelike hypersurfaces orthogonal to $u^\mu$, the function $\Phi \equiv \ln (\rho_0 + p)$ satisfies the elliptic equation
\begin{equation} \label{intro_pressure_eqn}
\Delta_u \Phi + \dot{u}^\nu \partial_\nu \Phi  + \nabla_\nu u^\mu\: \nabla_\mu u^\nu + R_{\mu\nu} u^\mu u^\nu = 0,
\end{equation}
where $R_{\mu\nu}$ is the Ricci tensor of the metric, $\dot{u}^\nu \equiv  u^\mu \nabla_{\mu} u^\nu$ and $$\Delta_u \Phi \equiv \Pi^{\mu\nu} \nabla_\mu \partial_\nu \Phi$$ for $\Pi^{\mu\nu} \equiv g^{\mu\nu} + u^\mu u^\nu$. Moreover, for $\rho >0$ not constant, $p$ satisfies the elliptic equation recorded in \eqref{pressure_eqn_general} on each hypersurface orthogonal to $u^\mu$.
\end{Thm}

A key observation for the proof of Theorem \ref{Thm1} is that $\Pi^{\mu\nu}$ projects onto the spacelike hypersurfaces orthogonal to $u^\mu$ and that $\Pi^{\mu\nu}$ is a Riemannian metric on these hypersurfaces, so that $\Delta_u$ is the Laplace operator of $\Pi^{\mu\nu}$ on each of these spacelike surfaces. From the ellipticity of \eqref{intro_pressure_eqn}, we thus expect \eqref{intro_pressure_eqn} to induce infinite speed of propagation into the incompressible covariant Euler equations, violating the axioms of Relativity.\footnote{Let us clarify our expectation: On any of the spacelike hypersurfaces orthogonal to $u$, boundary values determine solutions of \eqref{intro_pressure_eqn} by standard elliptic theory \cite{Taylor}. So a signal (specifying boundary values) send from one such hypersurface to a later such hypersurface, would determine $\Phi$ instantaneously throughout that surface. This indicates that the ellipticity of \eqref{pressure_eqn} could induce infinite speed of propagation and hence instantaneity. However, infinite wave speeds could be avoided, if the elliptic equation \eqref{intro_pressure_eqn} is suitably propagated by the evolution \eqref{incomp} -  \eqref{intro_Euler_2} imposes, (analogous to the propagation of the constraint equations in General Relativity), but it is not clear to us whether such a mechanism for propagating \eqref{intro_pressure_eqn} exists.}      However, taking the classical limit of \eqref{incomp} - \eqref{intro_Euler_2}, the equations reduce correctly to the incompressible Euler equations, c.f. Theorem \ref{Thm2}. So the incompressible covariant Euler equations could be of mathematical interest, e.g., to study problems of classical fluid dynamics.

\begin{Thm} \label{Thm2}
Assume $g_{\mu\nu}$ is the Minkowski metric. Then, replacing $\rho$ in \eqref{intro_T} by $\rho c^2$ and taking the classical limit ($c\rightarrow \infty$) of \eqref{incomp} - \eqref{intro_Euler_2}, we obtain the classical incompressible Euler equations
\begin{eqnarray} \label{incomp_Euler_classical}
\nabla \cdot {v} = 0, \ \ \ \  \
\partial_t \rho + v\cdot \nabla \rho = 0,  \ \ \ \   \
 \rho\left( \partial_t v + v\cdot\nabla v\right) + \nabla p = 0,
\end{eqnarray}
where $v\in \mathbb{R}^3$ is the fluid velocity, $\rho$ denotes the mass density of the fluid, $\nabla$ is the gradient on $\mathbb{R}^3$ and $\nabla \cdot {v}$ is the divergence of $v$. Moreover, taking the classical limit of \eqref{intro_pressure_eqn} yields
\begin{equation} \label{pressure_eqn_limit}
\Delta p + \rho_0 \: \partial_\beta {v^\alpha} \,\partial_\alpha{v^\beta} =0,
\end{equation}
where summation is over $\alpha,\beta =1,2,3$, and $\Delta$ is the Laplace operator in $\mathbb{R}^3$.
\end{Thm}

Equation \eqref{pressure_eqn_limit} is the classical equation for the pressure of an incompressible fluid, c.f. \cite[(1.84)]{MajdaBertozzi}. We prove Theorem \ref{Thm1} in section \ref{Sec_instantaneity} and leave the proof of Theorem \ref{Thm2} to section \ref{Sec_Newtonian_limit}.

\section{The elliptic equation for the pressure} \label{Sec_instantaneity}

We first establish a lemma which expresses the incompressible covariant Euler equation in a form closer to \eqref{incomp_Euler_classical}.

\begin{Lemma}
Assume $(N,\rho,u,p)$ is smooth and satisfies \eqref{incomp}. Then $(N,\rho,u,p)$ solves \eqref{intro_Euler_1} - \eqref{intro_Euler_2} if and only if $(N,\rho,u,p)$ solves
\begin{eqnarray}
\nabla_\sigma u^\sigma &=& 0, \label{incomp_rel_Euler_divu} \\
u^\sigma \partial_{\sigma}\rho  &=& 0, \label{incomp_rel_Euler_mass}\\
(\rho + p) u^\nu \nabla_\nu {u^\mu}  +\Pi^{\mu\nu} \partial_\nu p &=&0. \label{incomp_rel_Euler_momentum}
\end{eqnarray}
\end{Lemma}

\begin{proof}
Our condition of incompressibility \eqref{incomp}, keeping in mind that $N$ is assumed positive, immediately implies that \eqref{intro_Euler_1} is equivalent to \eqref{incomp_rel_Euler_divu}.

To obtain \eqref{incomp_rel_Euler_mass},  observe that the normalization of the fluid velocity \eqref{normalization} implies that
\begin{equation}  \label{1}
u^\sigma \nabla_\nu u_{\sigma}=0.
\end{equation}
From the expression for the energy-momentum tensor of a perfect fluid \eqref{intro_T} and from the identity \eqref{1}, we find that
\begin{eqnarray} \nonumber
u_\sigma\nabla_{\nu} {T^{\sigma\nu}}
&=&    u_\sigma\nabla_{\nu} \big( (\rho + p) u^\sigma u^\nu + p g^{\sigma\nu} \big) \cr
&=& (\rho + p) ( u^\nu u_\sigma \nabla_{\nu} u^\sigma  - \nabla_{\nu}  u^\nu ) - u^\nu \partial_{\nu}(\rho + p)  +  u^\sigma \partial_{\sigma} p \cr
&=& - u^\nu \partial_{\nu}\rho,
\end{eqnarray}
where we used \eqref{1} and \eqref{incomp_rel_Euler_divu} to obtain the last equality. Thus, the equation resulting from contracting \eqref{intro_Euler_2} with $u_\sigma$ is equivalent to \eqref{incomp_rel_Euler_mass}.

To obtain \eqref{incomp_rel_Euler_momentum}, observe that the normalization of the fluid $4$-velocity \eqref{normalization} implies that the projection $\Pi_{\mu\nu} \equiv g_{\mu\nu} + u_\mu u_\nu$ satisfies
\begin{eqnarray} \label{prejector-property}
\Pi_{\sigma\mu} u^\sigma =0 \ \ \ \ \ \  \ \text{and} \ \ \ \ \ \
\Pi_{\sigma\mu} v^\sigma = v_\mu
\end{eqnarray}
for any tangent vector $v^\mu$ with $v^\sigma u_\sigma=0$. Moreover, from \eqref{1} we obtain that
\begin{equation} \label{2}
\Pi_{\mu\sigma}  \nabla_\nu {u^\sigma} = g_{\mu\sigma}  \nabla_\nu {u^\sigma}  + u_\mu  u_\sigma\nabla_\nu  {u^\sigma}  =  \nabla_\nu u_{\mu}.
\end{equation}
For the energy momentum tensor of a perfect fluid \eqref{intro_T}, we find using \eqref{incomp_rel_Euler_divu}, \eqref{prejector-property} and \label{2} that
\begin{eqnarray}  \label{3}
\Pi_{\mu\sigma} \nabla_\nu {T^{\sigma\nu}}
&=& \Pi_{\mu\sigma} \nabla_\nu \big( (\rho + p) u^\sigma u^\nu + p g^{\sigma\nu} \big) \cr
&=& (\rho+p)\Pi_{\mu\sigma} u^\nu\nabla_\nu u^\sigma +\Pi_{\mu\sigma} g^{\sigma\nu} \partial_\nu p  \cr
&=& (\rho+p)  u^\nu\nabla_\nu  u_\sigma +\Pi_{\mu\sigma} g^{\sigma\nu} \partial_\nu p.
\end{eqnarray}
Raising the index with the metric, we conclude from \eqref{3} that the equations resulting from contraction of \eqref{intro_Euler_2} with $\Pi_{\mu\sigma}$ are equivalent to \eqref{incomp_rel_Euler_momentum}.
\end{proof} \noindent

Let us remark that \eqref{incomp_rel_Euler_divu} - \eqref{incomp_rel_Euler_mass} also follow by the standard replacement of $3$-vectors in \eqref{incomp_Euler_classical} by $4$-vectors. However, deriving \eqref{incomp_rel_Euler_divu} - \eqref{incomp_rel_Euler_mass} from \eqref{incomp}, it is clear that we captured the condition of rigid particle motion.  \vspace{.2cm}

\noindent \emph{Proof of Theorem \ref{Thm1}:}
Assume that $\rho\equiv \rho_0$ is constant and $p+\rho>0$.  Then, introducing $\Phi \equiv \ln (\rho_0 + p)$, we write \eqref{incomp_rel_Euler_momentum} in its equivalent form
\begin{equation} \label{incomp_rel_Euler_momentum2}
u^\nu \nabla_\nu {u^\mu} +\Pi^{\mu\nu} \partial_\nu \Phi = 0.
\end{equation}
Observe that $\Pi_{\mu\nu} \equiv g_{\mu\nu} + u_\mu u_\nu$ together with \eqref{incomp_rel_Euler_divu} imply
\begin{equation} \nonumber
\nabla_\mu \big( \Pi^{\mu\nu} \partial_\nu \Phi \big) = \Delta_u \Phi + \dot{u}^\nu \partial_\nu \Phi,
\end{equation}
where we set $\dot{u}^\nu \equiv  u^\mu \nabla_\mu u^\nu$ and $\Delta_u \Phi \equiv \Pi^{\mu\nu} \nabla_\mu \partial_\nu \Phi$. Moreover, \eqref{incomp_rel_Euler_divu} together with the definition of the Riemann curvature tensor of $g$, imply
\begin{equation} \nonumber
u^\nu \nabla_\mu  \nabla_\nu {u^\mu}  =  R^\mu_{\ \sigma \mu\nu} u^\sigma u^\nu = R_{\mu\nu} u^\mu u^\nu ,
\end{equation}
where $R_{\mu\nu}$ denotes the Ricci tensor.
Thus taking the divergence of \eqref{incomp_rel_Euler_momentum2} gives us \eqref{intro_pressure_eqn}, that is,
\begin{equation} \label{pressure_eqn}
\Delta_u \Phi +  \dot{u}^\nu \partial_\nu \Phi + \nabla_\nu {u^\mu} \nabla_\mu {u^\nu}  + R_{\mu\nu} u^\mu u^\nu = 0.
\end{equation}

We now show that \eqref{pressure_eqn} is an elliptic equation on each of the spacelike hypersurfaces orthogonal to $u$. By the second identity in \eqref{prejector-property} it follows that $\Pi_{\mu\nu}$ defines a Riemannian metric on each of these hypersurfaces. Thus $\Delta_u  \equiv \Pi^{\mu\nu} \nabla_\mu \partial_\nu$ is the Laplace operator on the surfaces orthogonal to $u$. Moreover, since $\nabla_\nu u^\mu = \Pi^{\mu\rho} \nabla_\nu u_{\rho}$, we have $\dot{u}^\nu =\Pi^\nu_{\ \rho} \dot{u}^\rho$ and it follows that the differentiation $\dot{u}^\nu \partial_\nu \Phi$ is taken tangential to the hypersurface orthogonal to the fluid flow $u$. The remaining terms in \eqref{pressure_eqn} are given source terms depending only on the metric tensor and $u$. In summary, we proved that \eqref{pressure_eqn} is an elliptic equation in $\Phi$, defined on each of the spacelike hypersurfaces orthogonal to $u$.

To prove the supplement, observe that a computation similar to that leading to \eqref{pressure_eqn}, yields
\begin{equation} \label{pressure_eqn_general}
\frac{1}{\rho + p} \left(\Delta_u p +  \dot{u}^\nu \partial_\nu p\right) - \frac{\Pi^{\mu\nu}}{(\rho + p)^2}\partial_{\nu}p \; \partial_\mu(\rho+p) + \nabla_\nu {u^\mu} \nabla_\mu {u^\nu}  + R_{\mu\nu} u^\mu u^\nu = 0,
\end{equation}
which is again an elliptic equation on each of the hypersurfaces orthogonal to the fluid $4$-velocity. This completes the proof.
\hfill $\Box$  \vspace{.2cm}

\section{The Newtonian limit} \label{Sec_Newtonian_limit}

In this section we compute the classical limit ($c\rightarrow \infty$) of the incompressible covariant Euler equations, \eqref{incomp} - \eqref{intro_Euler_2}, and of the pressure equation \eqref{intro_pressure_eqn}. We thereby prove Theorem \ref{Thm2}.  \vspace{.2cm}

\noindent \emph{Proof of Theorem \ref{Thm1}:}
Assume $g_{\mu\nu}=\eta_{\mu\nu}$ is the Minkowski metric and assume coordinates $(x^0,...,x^3)=(ct,x,y,z)$, where $(x,y,z)$ are Cartesian coordinates on $\mathbb{R}^3$. To work in SI-units, we replace $\rho$ in \eqref{intro_T} by $\rho c^2$, (for $c>1$ the speed of light in SI-units), and we relate the fluid $4$-velocity $u^\mu$ to the fluid velocity $v$ by setting\footnote{By our convention $u^\mu$ is dimensionless, in contrast to the convention in \cite{Choquet,Weinberg}.}
\begin{equation} \label{techeqn1}
u^\mu = \frac{\gamma(v)}{c} \left(\begin{array}{c} c \cr v \end{array}  \right)
\hspace{.6cm} \text{for}  \hspace{.6cm}
\gamma(v) \equiv \frac{1}{\sqrt{1-\frac{v^2}{c^2}} },
\end{equation}
which is consistent with the normalization \eqref{normalization}. From \eqref{techeqn1} we find that
\begin{eqnarray} \nonumber
\nabla_\mu u^\mu
&=& \frac{\partial u^0}{\partial (ct)} + \sum_{\alpha =1}^{3} \frac{\partial u^\alpha}{\partial x^\alpha} \ = \ \frac{1}{c} \left( \frac{\partial \gamma(v)}{\partial t} + \sum_{\alpha =1}^{3} v^\alpha \frac{\partial \gamma(v)}{\partial x^\alpha} \right)  + \frac{\gamma(v)}{c} \nabla \cdot v .
\end{eqnarray}
Now, observing that
\begin{equation} \label{techeqn2}
\frac{\partial \gamma(v)}{\partial t}= - \frac{\gamma(v)^3}{c^2}\, v\cdot \frac{\partial v}{\partial t},
\ \ \  \ \ \  \ \  \ \
\frac{\partial \gamma(v)}{\partial x^\alpha} = -\frac{\gamma(v)^3}{c^2}\, v\cdot \frac{\partial v}{\partial x^\alpha}
\end{equation}
and that $\gamma(v)\rightarrow 1$ as $c\rightarrow \infty$, it follows that
\begin{equation} \nonumber
\lim\limits_{c\rightarrow \infty} \nabla_\mu u^\mu  = \nabla \cdot v.
\end{equation}
We conclude that \eqref{incomp_rel_Euler_divu} reduces to the first equation in \eqref{incomp_Euler_classical} as $c\rightarrow \infty$.

From \eqref{techeqn1}, we immediately find that
\begin{equation} \nonumber
u^\sigma \partial_\sigma \rho
= \frac{\gamma(v)}{c} \left( \partial_t\rho + v\cdot \nabla \rho \right),
\end{equation}
from which we conclude that \eqref{incomp_rel_Euler_mass} is equivalent to the second incompressible Euler equation in \eqref{incomp_Euler_classical}.

We now compute the classical limit of \eqref{incomp_rel_Euler_momentum}. For this, introducing the notation $v^\mu\equiv \frac{c}{\gamma(v)} u^\mu$, observe that \eqref{techeqn1} implies
 \begin{equation} \label{techeqn1b}
 u^\nu  \nabla_\nu u^\mu
= \frac{1}{c^2} \gamma(v) v^\mu v^\nu \partial_\nu \gamma(v) + \frac{1}{c^2}  \gamma(v)^2 v^\nu \partial_\nu v^\mu ,
 \end{equation}
from which we find that
 \begin{equation} \nonumber
 (\rho c^2 + p) u^\nu  \nabla_\nu u^\mu
 = \left(\rho + \frac{p}{c^2} \right)  \left( \gamma(v)  v^\nu  v^\mu \partial_\nu \gamma(v) +  \gamma(v)^2 v^\nu \partial_\nu v^\mu \right).
 \end{equation}
Now, since $\lim\limits_{c\rightarrow \infty} \gamma(v)=1 $ and since by \eqref{techeqn2}
 \begin{eqnarray} \label{techeqn2b}
  v^\nu   \partial_\nu\gamma(v)
 \ = \ \frac{\partial \gamma(v)}{\partial t}  + v\cdot \nabla \gamma(v) \ = \ - \frac{\gamma(v)^3}{c^2}\, v\cdot \left(  \frac{\partial v}{\partial t} + v\cdot \nabla v  \right) ,
 \end{eqnarray}
we obtain that
  \begin{equation} \nonumber
\lim_{c\rightarrow \infty} (\rho c^2 + p)  u^\nu \nabla_\nu {u^\mu}
\ = \ \rho v^\nu \partial_\nu  {v^\mu}
\ = \ \rho  \left(\begin{array}{c} 0 \cr \partial_t v + v\cdot \nabla v \end{array} \right),
 \end{equation}
where the $\mu=0$ component vanishes.
 Moreover, a straightforward computation yields
 \begin{equation} \label{techeqn4}
 \Pi^{\mu\nu} = \eta^{\mu\nu} + u^\mu u^\nu \longrightarrow  \left( \begin{array}{cc} 0 & 0 \cr 0 & \text{id}_3  \end{array}  \right),
 \end{equation}
as $c\rightarrow \infty$, where $\text{id}_3$ denotes the identity mapping on $\mathbb{R}^3$. To summarize, we showed that \eqref{incomp_rel_Euler_momentum} reduces to the third equation in \eqref{incomp_Euler_classical} as $c\rightarrow \infty$.

It remains to confirm that \eqref{intro_pressure_eqn} converges to the classical pressure equation \eqref{pressure_eqn_limit} as $c\rightarrow \infty$. To begin, observe that for $g$ being the Minkowski metric we have $R_{\mu\nu}=0$ and \eqref{intro_pressure_eqn} reduces to
\begin{equation} \label{pressure_eqn_flat}
\Pi^{\mu\nu} \partial_\mu \partial_\nu \Phi + \dot{u}^\nu \partial_\nu \Phi  + \partial_\nu u^\mu\: \partial_\mu u^\nu = 0.
\end{equation}
Substituting $\Phi \equiv \ln(\rho_0 c^2 + p)$, we compute  for the first two terms that
\begin{eqnarray} \nonumber
(\rho_0 c^2+p) \left(\Pi^{\mu\nu} \partial_\mu\partial_\nu \Phi +  \dot{u}^\nu \partial_\nu \Phi \right)
&=& \Pi^{\mu\nu} \left( \partial_\mu\partial_\nu p - \frac{ \partial_\mu p \partial_\nu p}{\rho_0 c^2 +p} + \dot{u}^\nu \partial_\nu p \right).
\end{eqnarray}
Recalling that $\dot{u}^\nu \equiv u^\mu \nabla_{\mu} u^\nu$, \eqref{techeqn1b} together with \eqref{techeqn2b} imply $\lim\limits_{c\rightarrow \infty} \dot{u}^\nu \partial_\nu p =0$, so that in light of  \eqref{techeqn4} we conclude
\begin{eqnarray} \nonumber
\lim\limits_{c\rightarrow \infty}\; (\rho_0 c^2+p) \left(\Pi^{\mu\nu} \partial_\mu\partial_\nu \Phi +  \dot{u}^\nu \partial_\nu \Phi \right) \: =\:   \Delta p.
\end{eqnarray}
To handle the last term in \eqref{pressure_eqn_flat}, use  \eqref{techeqn2}, \eqref{techeqn1b} and  \eqref{techeqn2b}, to compute
\begin{equation} \nonumber
\lim\limits_{c\rightarrow \infty}\; (\rho_0 c^2+p) \partial_\nu u^\mu \partial_\mu {u^\nu}
\; = \; \lim\limits_{c\rightarrow \infty}\; (\rho_0 c^2+p) \frac{\gamma(v)^2}{c^2} \partial_\nu v^\mu \partial_\mu {v^\nu}
\; = \; \rho_0  \lim\limits_{c\rightarrow \infty} \partial_\nu v^\mu \partial_\mu {v^\nu} ,
\end{equation}
and since $\partial_\mu v^0 = 0$,  we obtain that
\begin{equation} \nonumber
\lim\limits_{c\rightarrow \infty}\; (\rho_0 c^2+p) \partial_\nu u^\mu \partial_\mu {u^\nu}
\; = \; \rho_0   \partial_\beta {v^\alpha} \partial_\alpha {v^\beta},
\end{equation}
where summation is over $\alpha,\beta =1,2,3$. In summary, we proved that multiplying \eqref{pressure_eqn_flat} by $(\rho_0 c^2+p)$, the resulting equation converges to the sought after equation \eqref{pressure_eqn_limit} as $c\rightarrow \infty$. This proves Theorem \ref{Thm2}.    \hfill $\Box$

\section*{Acknowledgments}
I thank Blake Temple for helpful discussions.

\providecommand{\bysame}{\leavevmode\hbox to3em{\hrulefill}\thinspace}
\providecommand{\MR}{\relax\ifhmode\unskip\space\fi MR }
\providecommand{\MRhref}[2]{%
  \href{http://www.ams.org/mathscinet-getitem?mr=#1}{#2} }
\providecommand{\href}[2]{#2}