1 Introduction↩︎

The late–time acceleration of the Universe has been well described, at the background level, by general relativity (GR) with a cosmological constant, which acts as a uniform source with an effective dark energy equation-of-state parameter \(w=-1\) [1]–[3]. However, recent observations have begun to place strain on this minimal picture, with DESI DR2 reporting mild tensions in \(\Lambda\)CDM and an improved fit for time-dependent dark energy when BAO data are combined with CMB and supernova measurements [4]. More specifically, the preference is for an evolution in which the dark energy equation of state crosses the divide from a phantom to a quintessence behaviour at late times. Independently, the large discrepancy between quantum field theory estimates of the vacuum energy and the observed value of \(\Lambda\) [5]–[7] keeps open the possibility that cosmic acceleration is emergent or dynamical rather than a rigid constant, motivating theories that reconsider the provenance of \(\Lambda\) and seek dynamical origins for the acceleration within an enlarged gravitational sector.

One such framework is that of unimodular gravity, which reframes the cosmological constant as an integration constant fixed by boundary conditions rather than a coupling in the action [8], [9]. Whilst this does not, by itself, solve the cosmological constant problem [10], it nevertheless provides a geometric origin for \(\Lambda\), arising as an integration constant associated with enforcing a constraint on the metric determinant.

Another alternative is scalar–tensor gravity [11]–[13], in which a single spin–0 degree of freedom can modify the large-scale expansion whilst remaining compatible with local tests of gravity through screening mechanisms [14]. The additional scalar may be further used to generalise matter couplings: matter fields interact with a matter metric \(\hat{g}_{\mu\nu}\) that is related to \(g_{\mu\nu}\) by a combination of conformal and disformal transformations mediated by a scalar field \(\phi\) [15], [16]. Such couplings are particularly appealing in cosmology, as the disformal component is typically suppressed in weak-field, quasi-static regimes but can become relevant on homogeneous backgrounds [17]–[19], thereby providing a natural mechanism for dark energy.

Since unimodular gravity provides a geometric rationale for the cosmological constant, whilst scalar–tensor theories introduce additional degrees of freedom and generalised matter couplings, it is only natural to ask what becomes of this unimodular constant when matter sees a second metric. This work addresses precisely that question: given two metrics, the Einstein metric \(g_{\mu\nu}\) governing the gravitational field equations and the matter metric \(\hat{g}_{\mu\nu}\) to which non-gravitational fields couple, what is the analogue of the unimodular cosmological constant, and how does it affect the cosmological dynamics?

We show that the answer to this question depends on how unimodularity is implemented. In particular, imposing dual unimodular constraints on both \(g_{\mu\nu}\) and \(\hat{g}_{\mu\nu}\) breaks the classical equivalence between the different formalisms of unimodular gravity. Specifically, the nature of the resulting vacuum contribution, and the ability to achieve a phantom crossing, depend on how unimodularity is implemented.

The structure of this paper is as follows. In Sec. 2 we review the relevant formulations of unimodular gravity and introduce the disformal bimetric framework used throughout. In Sec. 3 we apply different unimodularisation procedures to this bimetric theory, comparing the resulting dynamics and clarifying the role played by the associated vacuum contributions. In Sec. 4 we specialise these constructions to homogeneous cosmology and examine their implications for the dark energy equation of state, including the possibility of phantom-divide crossing. In Sec. 5 we present a diffeomorphism-invariant reformulation which recovers the fixed-determinant construction in an appropriate gauge. Finally, we conclude in Sec. 6.

Throughout this paper, we adopt the mostly-plus metric signature \((-+++)\), and define the canonical kinetic term of the scalar field as \(X \coloneq -\tfrac12 g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi\). We work in natural units with \(c=1\) and define the reduced Planck mass as \(M_\text{Pl}^2 \equiv(8\pi G)^{-1 }\).

2 Preliminaries↩︎

We begin by outlining the theoretical background required for the construction of bimodular gravity. Since the framework is centred on deriving the cosmological constant as a Lagrange multiplier, it is helpful to first recall how this arises in the different formulations of unimodular gravity. For brevity, we focus on the three versions reviewed extensively in [10]. We then turn to the scalar–tensor bimetric gravity theory of [20]–[23], for it is the framework subject to unimodularisation in this work. Since this work is interested in purely the combined gravitational and biscalar effects, we neglect additional matter fields.

2.1 Unimodular Gravity↩︎

If one imposes a constraint on the metric via a Lagrange multiplier, that multiplier appears in the gravitational field equations and, for a suitable choice of constraint, can be identified with the cosmological constant. This modified formulation of GR is known as unimodular gravity, whose various formulations all yield field equations classically equivalent to GR\(+\Lambda\). The differences amongst them lie only in how the determinant constraint (and thus the constancy of the associated Lagrange multiplier) is enforced, and hence in their off-shell symmetry content.

In the original fixed-determinant formulation of unimodular gravity (UG), one constrains \(\sqrt{-g}\) to a fixed scalar density \(\varepsilon_0(x)\) by means of a Lagrange multiplier \(\lambda\), \[S_\text{UG} = \int d^4x \left[ \frac{1}{2} \sqrt{-g}R[g] - \lambda(x) \left(\sqrt{-g}-\varepsilon_0(x) \right) \right].\] It is the presence of the last term that breaks full diffeomorphism invariance. Variation with respect to \(\lambda\) imposes the constraint \(\sqrt{-g}=\varepsilon_0\), whilst variation with respect to \(g_{\mu\nu}\) yields the Einstein equations \[G_{\mu\nu} +\lambda g_{\mu\nu} = 0. \label{eq:UG95EinsteinEquation}\tag{1}\] Taking the divergence of 1 , and using the Bianchi identity \(\nabla_\mu G^{\mu\nu}=0\), we find that the multiplier must satisfy \(\nabla_\mu\lambda=0\). Thus, it is a derived result that \(\lambda\) is constant on shell and can be identified with the cosmological constant \(\Lambda \coloneq \lambda\) [24]–[27]. Hence the gravitational field equations of UG coincide with those of GR\(+\Lambda\), with \(\Lambda\) appearing as an integration constant rather than as a parameter in the action.

The covariant formulation of Henneaux and Teitelboim (HT) [28] restores full diffeomorphism invariance by replacing the fixed-determinant constraint on the metric with a constraint on the divergence of an auxiliary vector density \(\tau^\mu\), \[S_\text{HT} = \int d^4x \left[ \frac{1}{2} \sqrt{-g} R[g] - \lambda(x) \left( \sqrt{-g} - \partial_\mu\tau^\mu(x) \right) \right].\] Here, variation with respect to \(\lambda\) enforces \(\sqrt{-g}=\partial_\mu \tau^\mu\), whilst variation with respect to \(\tau^\mu\) yields \(\partial_\mu\lambda=0\), so the constancy of \(\lambda\) is imposed already at the level of the action [29]. Variation with respect to \(g_{\mu\nu}\) then gives the Einstein field equations with a more standard cosmological constant \(\Lambda \coloneq \lambda\), so the HT formulation reproduces GR\(+\Lambda\) whilst maintaining full diffeomorphism invariance.

A closely related fully diffeomorphism-invariant formalism (DUG) trades the metric constraint for a divergence condition on an auxiliary vector \(V^\mu\), \[S_\text{DUG} = \int d^4x \left[ \frac{1}{2}\sqrt{-g} R[g] - \sqrt{-g}\left(\lambda(x)+V^\mu(x)\partial_\mu \lambda(x) \right) \right],\] which is equivalent to HT up to a boundary term. Variation of \(V^\mu\) directly yields \(\partial_\mu \lambda=0\), and variation with respect to \(\lambda\) imposes \(\nabla_\mu V^\mu =1\) [10].

In all three cases, the classical field equations are on-shell equivalent to those of GR with a cosmological constant. Their differences are purely off shell, specifically in how the unimodular condition is implemented and how the constancy of \(\lambda\) is imposed—derived via the Bianchi identity (UG) or assumed at the level of the action (HT/DUG). These distinctions are classically inconsequential in the single-metric setting but will become decisive for the failure of the classical equivalence once unimodularity is applied in a bimetric framework. To that end, the next section introduces the disformal bimetric theory that will subsequently be unimodularised.

2.2 Bimetric Gravity↩︎

In the bimetric gravity theory considered here, gravity is described by an Einstein metric \(g_{\mu\nu}\) and a matter metric \(\hat{g}_{\mu\nu}\), disformally related through a scalar field \(\phi\), often referred to as the biscalar, via \[\hat{g}_{\mu\nu} = g_{\mu\nu} +B(\phi) \partial_\mu\phi\partial_\nu\phi, \label{eq:disformal95relation}\tag{2}\] with the disformal coupling \(B(\phi)\). We restrict to the branch \(B(\phi)>0\) and impose the condition \(1-2BX>0\) which ensures invertibility of the disformal map and preservation of Lorentzian signature for timelike scalar gradients. In that case, from the Sherman–Morrison formula and matrix determinant lemma, the inverse and determinant of 2 follow as \[\begin{align} \hat{g}^{\mu\nu} &= g^{\mu\nu} - \frac{B}{1-2BX}\partial^\mu\phi\partial^\nu\phi \tag{3} , \\ \sqrt{-\hat{g}} &= \sqrt{-g}\sqrt{1-2BX}, \tag{4} \end{align}\] so that the relative volume element between the two metrics is \[\nu \mathrel{\vcenter{:}}= \frac{\sqrt{-\hat{g}}}{\sqrt{- g}} = \sqrt{1-2BX}. \label{eq:mathcalV}\tag{5}\] These are precisely the objects that will be subjected to unimodular constraints in Sec. 3. It may already be obvious how the different unimodular formalisms will affect \(\nu\), and subsequently the biscalar dynamics.

For the biscalar to be dynamical, a kinetic term must be chosen. One may place \(\phi\) in the Einstein frame, where it is canonical relative to \(g_{\mu\nu}\) and luminal with respect to that frame, or in the matter frame, where conservation laws take their standard form [30]. For present purposes, it is more insightful to adopt the latter so that the gravitational sector alone sees two metrics. We thus consider the bimetric action \[S = \int d^4x \left[ \frac{1}{2}\sqrt{-g}R + \sqrt{-\hat{g}} \hat{X} \right], \label{eq:BGT}\tag{6}\] with \(\hat{X} \coloneq -\frac{1}{2} \hat{g}^{\mu\nu}\partial_\mu\phi\partial_\nu\phi =X/\nu^2\) where the latter equality is a non-trivial result. Here, the matter stress tensor is conserved in the matter frame, \(\hat{\nabla}_\mu \hat{T}^{\mu\nu}=0\), i.e. on the hatted scalar equations of motion. Indeed, by rewriting the theory in a single-metric presentation, obtained by direct substitution of 3 –5 into 6 , one finds a particular k-essence theory, \[S = \int d^4x \sqrt{-g} \left[ \frac{1}{2} R + K(\phi,X) \right], \qquad K(\phi,X) = \frac{X}{\sqrt{1-2BX}} . \label{eq:BGT95kessence}\tag{7}\] This is a change of presentation, not of physical content. In the single-metric presentation, diffeomorphism invariance implies the usual on-shell conservation law \(\nabla_\mu T^{\mu\nu}=0\). The two conservation laws are therefore not independent assumptions: they are equivalent on shell under the disformal field redefinition, although the stress tensors and covariant derivatives appearing in the two descriptions are different.¹

If one were to include matter, its action should be coupled universally to a single physical metric, \[S_{\rm m} = S_{\rm m}[\tilde{g}_{\mu\nu}, \psi_{\rm m}], \qquad \tilde{g}_{\mu\nu} = g_{\mu\nu} \quad\text{or}\quad \tilde{g}_{\mu\nu}= \hat{g}_{\mu\nu},\] in order to preserve the weak equivalence principle. The choice \(\tilde{g}_{\mu\nu}=\hat{g}_{\mu\nu}\) is natural from the original bimetric perspective, since \(\hat{g}_{\mu\nu}\) is then the matter metric. In this case, observables would then be defined in the matter frame.

3 Unimodularising Bimetric Gravity↩︎

In this section, we subject the bimetric gravity action 6 to the unimodularisation procedures of Sec. 2.1. Whilst there is no classical difference between the two formalisms in the single metric case, we find here that that classical equivalence is broken when a disformal metric is present. In order to unimodularise the disformal bimetric setup in full generality, we impose one unimodular-type condition per metric using two independent Lagrange multipliers, \(\lambda_1\) and \(\lambda_2\), and two constraint functionals \(\mathcal{C}_1[g]\) and \(\mathcal{C}_2[\hat{g}]\). The general action is then \[\begin{align} S =& \int d^4x \left[ \frac{1}{2} \sqrt{-g} R[g] + \sqrt{-\hat{g}} \hat{X} - \lambda_1 \mathcal{C}_1[g] - \lambda_2 \mathcal{C}_2[\hat{g}] \right]. \label{eq:bimodular95gravity95action95general} \end{align}\tag{8}\] The choice of constraint functionals specifies the appropriate unimodular formalism, and thus the bimetric unimodular formalism, which we now present for the three approaches discussed in Sec. 2.1.

3.1 Bimetric Unimodular Gravity Theory↩︎

We begin by subjecting 8 to the fixed-determinant constraints. Since \(\nu\) is tied to the biscalar kinetic term, the constraints have direct consequences for the biscalar dynamics. Furthermore, in single-metric UG the constancy of the Lagrange multiplier, and hence its identification with \(\Lambda\), was not assumed but derived from the Bianchi identity. We therefore close by computing the analogous relation to determine the nature of the bimetric multipliers.

3.1.1 Fixed-Determinant Bimetric Constraints↩︎

We thus define Bimetric Unimodular Gravity (BUG) theory as \[S_{BUG} \supset -\int d^4x \left[ \lambda_1\left( \sqrt{-g}-\varepsilon_1\right) +\lambda_2 \left( \sqrt{-\hat{g}} - \varepsilon_2 \right) \right], \label{eq:BUG95action}\tag{9}\] where \(\varepsilon_{1,2}(x)\) are fixed scalar densities. Variation with respect to \(\lambda_{1,2}\) enforces \[\sqrt{-g}=\varepsilon_1,\qquad \sqrt{-\hat{g}}=\varepsilon_2, \label{eq:BUG95constaints}\tag{10}\] Hence the relative volume element 5 is not dynamical, but fixed to the ratio of the fixed scalar densities \[\nu(x) = \frac{\varepsilon_2(x)}{\varepsilon_1(x)} = \sqrt{1-2B X}. \label{eq:BX-const}\tag{11}\] It follows that the biscalar kinetic density is kinematically locked to the fixed density ratio \[X = \frac{1-\nu^2(x)}{2B(\phi)}. \label{eq:kin95lock}\tag{12}\] Differentiating 12 gives \[\nabla_\mu X = -\frac{B_\phi X}{B}\nabla_\mu\phi - \frac{\nu}{B}\nabla_\mu\nu . \label{eq:app95gradX95BUG}\tag{13}\]

Now, suppose that the scalar gradient is timelike, \(X>0\). One may then define the unit four-velocity \[u_\mu=\frac{\nabla_\mu\phi}{\sqrt{2X}}, \qquad u^\mu u_\mu=-1,\] together with the orthogonal projector \[h_{\mu\nu}=g_{\mu\nu}+u_\mu u_\nu .\] The four-acceleration of the scalar flow is \[a_\mu \coloneq u^\nu\nabla_\nu u_\mu = -\frac{1}{2X}h_\mu{}^\nu\nabla_\nu X . \label{eq:app95acceleration95X}\tag{14}\] Thus the flow is geodesic if and only if the spatially projected gradient of \(X\) vanishes, \(h_\mu{}^\nu\nabla_\nu X=0\). Projecting 13 orthogonally to the scalar flow removes the first term, since it is parallel to \(\nabla_\mu\phi\). Hence the acceleration becomes \[a_\mu = \frac{\nu}{2BX} h_\mu{}^\nu\nabla_\nu\nu \neq0. \label{eq:app95BUG95acceleration}\tag{15}\] The same conclusion holds in the matter frame, since the corresponding acceleration is related by \(a_\mu=\nu^2\hat{a}_\mu\). Generic BUG therefore does not classify as dusty dark energy [31].

This also clarifies the status of the biscalar sound speed. Since the BUG scalar sector is kinematically constrained by the fixed relative volume element, it is not an ordinary freely propagating k-essence degree of freedom, and the standard k-essence expression for \(c_s^2\) is not the relevant physical propagation speed. At the same time, the Dust of Dark Energy argument for an identically vanishing physical sound speed relies on geodesic flow [31], or equivalently on the constraint reducing to \(X=X(\phi)\). In this sense BUG is neither standard k-essence nor generically identical to Dust of Dark Energy, although it contains the latter as a special case if \(\nu = \text{const}\) with \(B(\phi)\). For the remainder of this work we take \(B(\phi)=B\) to be constant and treat \(\nu(x)\) as an arbitrary fixed function, since it is this combination that gives us the simplest novel behaviour. As such, the BUG constraint reduces to \[\nabla_\mu X = -\frac{\nu}{B}\nabla_\mu\nu . \label{eq:kin95lock95deriv}\tag{16}\] Interestingly, in the special case \(\nu(x)=\nu\) is constant in addition to constant \(B\), we mimic the constant kinetic constraint of Mimetic Gravity [32], [33].

But of course, no condition has yet been imposed upon \(\lambda_{1,2}\). In single-metric UG, constancy of the multiplier is derived via the Bianchi identities, which then allows it to be identified with the cosmological constant. Therefore, the analogous relation must be derived through explicit calculation of the BUG Bianchi identity. For that, the metric and bimetric field equations must be determined.

3.1.2 The Dynamics of BUG↩︎

For a constant disformal coupling \(B\), extremising the general BMG action 8 with respect to the biscalar gives the field equation² \[\begin{align} \left[B (\hat{X} + \lambda_2) + 1 \right] \hat{\Box} \phi = - B\hat{\nabla}_\mu \hat{X} \hat{\nabla}^\mu \phi - B \hat{\nabla}_\mu \lambda_2 \hat{\nabla}^\mu \phi . \label{eq:BUG95ScalarEoM} \end{align}\tag{17}\] In BUG, however, \(\hat{X}=X/\nu^2\) and its derivative are fixed by the kinematic constraints 12 and 16 . We keep 17 in matter-frame notation for compactness.

Varying the BUG action 9 with respect to the Einstein-frame metric, expressing all quantities in the Einstein frame, and dividing through by \(\sqrt{-g}\) gives \[\begin{align} G_{\mu\nu}+\lambda_1 g_{\mu\nu} = T_{\mu\nu}, \label{eq:BUG95Einstein} \end{align}\tag{18}\] where³ \[\begin{align} T_{\mu\nu} = \left( \frac{X}{\nu} - \nu \lambda_2\right) g_{\mu\nu} + \left( \frac{1-BX}{\nu^3} + \frac{B\lambda_2 }{\nu} \right) \partial_\mu\phi\partial_\nu\phi . \label{eq:BUG95SET} \end{align}\tag{19}\] Interestingly, this effective stress tensor has the same algebraic scalar-fluid form encountered in related mimetic gravity constructions [33]. Though the resemblance is only structural, where the physical distinction lies in the constraint determining \(X(x,\phi)\) and in the behaviour of the Lagrange multipliers.

3.1.3 BUG Vacuum Contributions↩︎

It remains to determine the nature of the Lagrange multipliers \(\lambda_{1,2}\). In ordinary UG, one identifies the Lagrange multiplier with a cosmological constant by taking the divergence of the Einstein equations (see Sec. 2.1). The same calculation in the present bimetric case, however, yields a different result. Using the biscalar equation of motion 17 , re-expressed in the Einstein frame, one finds that the BUG biscalar stress tensor 19 is not separately conserved, \[\nabla_\mu T^{\mu\nu} = -\nu \nabla^\nu \lambda_2 . \label{eq:BUG95T95noncons}\tag{20}\] This comes with no surprise, since the fixed-density construction is not fully diffeomorphism invariant [35]. Taking the divergence of 18 and using the Bianchi identity therefore yields \[\nabla_\mu\lambda_1 = -\nu \nabla_\mu\lambda_2 . \label{eq:BUG95balancelaw}\tag{21}\] Thus neither BUG multiplier need be constant. Their spacetime dependence is instead tied together by the balance law 21 . Consequently, the final term in the scalar equation 17 survives.

The balance law 21 also carries an integrability condition. Since \(\lambda_1\) is a scalar, the one-form \(-\nu d\lambda_2\) must be exact [36]. Taking an exterior derivative of \(d\lambda_1=-\nu d\lambda_2\) gives \(d\nu\wedge d\lambda_2=0\), or equivalently in components, \[\nabla_{[\mu}\nu\,\nabla_{\nu]}\lambda_2=0 .\] Thus, for a prescribed relative volume element \(\nu(x)\), the second multiplier cannot be chosen as an arbitrary independent spacetime function. Rather, locally its gradient must be aligned with that of \(\nu\). This condition is of course trivially satisfied on homogeneous FLRW backgrounds, where all quantities depend only on time.

The key novel aspect of generic BUG, therefore, is that it does not possess a conserved bimodular cosmological constant. In the special case where the relative volume element is constant, \(\nu(x)=\nu\), the balance law becomes \[\nabla_\mu\left(\lambda_1+\nu\lambda_2\right)=0. \label{eq:Bi-UG95lambda95constraint}\tag{22}\] Only in this special constant-\(\nu\) case may the bracketed combination be identified with a genuine bimodular cosmological constant, albeit one built from dynamical constituents. For a general fixed scalar density ratio, however, \[\Pi(x) \mathrel{\vcenter{:}}= \lambda_1(x)+\nu(x)\lambda_2(x)\] obeys \[\nabla_\mu \Pi = \lambda_2\nabla_\mu\nu . \label{eq:varying95vacuum}\tag{23}\] It is therefore better interpreted as an effective vacuum contribution whose spacetime variation is sourced by the prescribed relative volume element.

Motivated by these properties, we shall henceforth refer to any classical realisation of the construction whose dynamics prescribe the relative volume element, constrain the kinetic term to the form \(X=X(x;\phi)\), and generate an effective vacuum contribution which need not be conserved as Bimodular Gravity.⁴

3.2 Bimetric HT Gravity Theory↩︎

We now turn to the HT implementation of the same bimetric theory. In the single-metric case, this formulation is classically equivalent to fixed-determinant unimodular gravity. In the present disformal bimetric setting, however, the two implementations do not lead to the same classical dynamics. We show that the resulting theory is therefore closer to a standard k-essence model, supplemented by a cosmological constant and a constant scalar-sector offset, whose perturbative stability can be classified in the usual way.

3.2.1 HT Bimetric Constraints↩︎

To preserve full diffeomorphism invariance in bimodular gravity, we instead implement HT constraints for both metrics, thereby defining the Bimetric Henneaux and Teitelboim (BHT) theory, whose constraint sector in the bimodular action 8 reads \[\begin{align} S_{BHT} \supset - \int d^4x \Bigg[& \lambda_1 \left( \sqrt{-g} - \partial_\mu\tau_1^\mu \right) + \lambda_2 \left( \sqrt{-\hat{g}} - \partial_\mu\tau_2^\mu \right) \Bigg]. \label{eq:eq:BHT95action} \end{align}\tag{24}\] Variation with respect to the auxiliary densities enforces the spacetime constancy of the multipliers, \[\begin{align} \partial_\mu\lambda_1 = 0, \qquad \partial_\mu\lambda_2 = 0, \label{eq:BHT95constant95lambdas} \end{align}\tag{25}\] so that \(\lambda_{1,2}\) are individually integration constants, in contrast to BUG in which the multipliers obey a balance law 21 . This is the first distinction from BUG; the second is kinematic. Variation with respect to \(\lambda_{1,2}\) imposes the corresponding bimodular constraints \[\sqrt{-g}=\partial_\mu\tau_1^\mu, \qquad \sqrt{-\hat{g}}=\partial_\mu\tau_2^\mu. \label{eq:BHT95constaints}\tag{26}\] Consequently, the relative volume element becomes \[\nu = \frac{\partial_\mu\tau_2^\mu}{\partial_\mu\tau_1^\mu} = \sqrt{1-2B X}, \label{eq:BHT95V}\tag{27}\] which, unlike BUG, remains fully dynamical. Consequently, the biscalar may keep its full dynamical freedom.

3.2.2 The Dynamics of BHT↩︎

Variation of the BHT action 24 with respect to the biscalar yields \[\begin{align} \left[B (\hat{X} + \lambda_2) + 1 \right] \hat{\Box} \phi = - B\hat{\nabla}_\mu \hat{X} \hat{\nabla}^\mu \phi . \label{eq:BiHT95ScalarEoM} \end{align}\tag{28}\] Here we have used the defining BHT property that the Lagrange multipliers are integration constants 25 . Consequently, the term proportional to \(\hat{\nabla}_\mu\lambda_2\) vanishes identically.

Variation with respect to the Einstein-frame metric gives \[\begin{align} G_{\mu\nu}+\lambda_1 g_{\mu\nu} = T_{\mu\nu}, \label{eq:BHT95Einstein} \end{align}\tag{29}\] with \[\begin{align} T_{\mu\nu} = \left( \frac{X}{\nu} - \nu \lambda_2\right) g_{\mu\nu} + \left( \frac{1-BX}{\nu^3} + \frac{B\lambda_2 }{\nu} \right) \partial_\mu\phi\partial_\nu\phi . \label{eq:BHT95SET} \end{align}\tag{30}\] This is algebraically identical to the BUG stress tensor, but its interpretation is different. In BHT, \(\lambda_1\) and \(\lambda_2\) are separately constant, and the action is fully diffeomorphism invariant. Taking the divergence of 29 therefore gives \[\nabla_\mu T^{\mu\nu}=0 , \label{eq:BHT95T95cons}\tag{31}\] so the Bianchi identity is satisfied in the standard way [34].

3.2.3 BHT Vacuum Contributions↩︎

The two constant vacuum contributions are thus separated in BHT. The first multiplier appears directly in the Einstein equations and may be identified with the cosmological constant, \[\Lambda \mathrel{\vcenter{:}}= \lambda_1 . \label{eq:BHT95lambda1}\tag{32}\] The second multiplier acts as a constant scalar potential in the matter frame, \[V_0\coloneq\lambda_2. \label{eq:BHT95lambda2}\tag{33}\] Hence BHT behaves as a unimodular gravity plus a disformally induced k-essence scalar, with the cosmological constant and scalar potential offset arising as distinct integration constants. This structure is close in spirit to mimetic k-essence [33] and generalized unimodular k-essence [37] constructions.

Evidently, BHT is a classically distinct theory to BUG. In BUG, the scalar field is non-propagating, and the two multipliers are tied by the balance law 21 , such that for a general fixed relative volume element no conserved combination \(\Pi(x)\) exists. In BHT, by contrast, the two constants remain independently conserved, whilst the scalar retains its genuine dynamical degree of freedom. Thus, BHT does not fall within the class of Bimodular Gravity. It is a more standard k-essence with a standard k-essence speed of sound which we now classify.

3.2.4 BHT Speed of Sound↩︎

For the propagating scalar field of BHT, the scalar sound speed is obtained via \[c_s^2 =\frac{K_X}{K_X + 2X K_{XX}}, \label{eq:sound95speed95general}\tag{34}\] where subscripts denote partial derivatives with respect to \(X\) [38]. For BHT, the scalar kinetic function is read off as \[\begin{align} K(\phi, X;\lambda_2) = \frac{X}{\sqrt{1-2BX}} -\sqrt{1-2BX} \, V_0. \label{eq:BMG95kessence} \end{align}\tag{35}\] Evaluating the partial derivatives with respect to \(X\) at fixed \(\phi\) and \(\lambda_2\), one finds \[\begin{align} K_X &= \frac{1-BX}{\nu^3} + \frac{B}{\nu}V_0, \tag{36} \\ K_{XX} &= \frac{B (2-BX)}{\nu^5} + \frac{B^2}{\nu^3}V_0. \tag{37} \end{align}\] Substituting 36 37 into 34 yields \[c_s^2= \nu^2\frac{ N}{D}, \label{eq:cs95explicit}\tag{38}\] where we have defined \[N \coloneq 1-BX +\nu^2 BV_0, \qquad D \coloneq 1 + B X + \nu^2BV_0.\]

The absence of gradient instabilities \(c_s^2>0\) can be satisfied on two algebraic branches. On the positive-denominator branch, \(D>0\), the condition requires \[V_0> \frac{BX-1}{B\nu^2}. \label{eq:BHT95grad95cond}\tag{39}\] There is also a negative-denominator branch, which requires \[V_0<-\frac{1+BX}{B\nu^2}, \label{eq:BHT95neg95grad95cond}\tag{40}\] whereby gradient stability would also require the numerator to be negative. For \(BX>0\), this follows immediately from 40 , since then \[(1-BX)+B\nu^2 V_0<-2BX<0 .\] For \(BX<0\), however, the negativity of the numerator is not automatic and must be imposed separately. In any case, this branch has negative kinetic coefficient and is therefore discarded as ghost-like. We therefore restrict to the branch satisfying both \(D>0\) and \(N>0\), with \(BX<1/2\), such that the biscalar is free of ghosts and gradient instabilities.

On this healthy branch, the sound speed can then be classified according to the value of \(BX\) as \[\begin{align} &\text{(i) } BX\to 1/2^-: && \nu^2\to0, && c_s^2\to0, && \text{dust-like limit}, \nonumber\\[0.3em] &\text{(ii) } 0<BX<1/2: && 0<\nu^2<1, && 0<c_s^2<1, && \text{subluminal}, \nonumber\\[0.3em] &\text{(iii) } BX=0: && \nu^2=1, && c_s^2=1, && \text{luminal}, \nonumber \\[0.3em] &\text{(iv) } BX<0: && \nu^2>1, && c_s^2>1, && \text{superluminal}. \nonumber \label{eq:BHT95sound95speed95classes} \end{align}\tag{41}\] Note that superluminality [39], [40] occurs when \(BX<0\), which for \(B>0\), requires \(X<0\), corresponding to a spacelike kinetic term.

3.3 Bimetric Fully Diffeomorphism-Invariant Unimodular Gravity Theory↩︎

Alternatively, still under the wish to preserve full diffeomorphism invariance in the bimodular framework, one may instead employ the DUG formalism as to define Bimetric Diffeomorphism-invariant Unimodular Gravity (BDUG), whose constraint sector becomes \[S_{BDUG} \supset - \int d^4x \Bigg[ \sqrt{-g} \left( \lambda_1 + V_1^\mu\nabla_\mu \lambda_1 \right) + \sqrt{-\hat{g}}\left( \lambda_2 + V_2^\mu \hat{\nabla}_\mu \lambda_2 \right) \Bigg], \label{eq:BDUG95action}\tag{42}\] where \(\nabla_{\mu}\) and \(\hat{\nabla}_{\mu}\) denote the covariant derivatives compatible with \(g_{\mu\nu}\) and \(\hat{g}_{\mu\nu}\), respectively. Variation with respect to the auxiliary vectors enforces spacetime constancy of the multipliers, \[\nabla_\mu \lambda_1 = 0, \qquad \hat{\nabla}_\mu\lambda_2 = 0, \label{eq:BHT95int95consts}\tag{43}\] so that, just as in BHT, \(\lambda_{1,2}\) play the role of strict integration constants. Varying with respect to \(\lambda_{1,2}\) yields the divergence conditions \[\nabla_\mu V^\mu_1 = 1, \qquad \hat{\nabla}_\mu V^\mu_2 = 1,\] which entirely replace the fixed–determinant constraints by normalisation conditions on the auxiliary vectors. As a consequence, the relative volume element remains fully dynamical.

The shared imposition of constant \(\lambda_{1,2}\) and fully dynamical \(BX\) means that the BHT and BDUG theories are equivalent on shell. Hence, for all intents and purposes of this work, we treat them as identical and focus on the BHT formulation for notational convenience, with its classical dynamics taken to represent those of BDUG—neither of which classifies as a Bimodular Gravity theory.

4 Cosmology↩︎

Having established the distinct dynamics of BUG and BHT, we now investigate their cosmological implications. We impose the spatially flat FLRW ansatz in the bimetric setting, derive the corresponding background equations, and examine how each implementation shapes the effective expansion history. In particular, we reconstruct the dark energy equation of state and determine whether a crossing of the phantom divide may occur in either case.

4.1 Bimetric FLRW Ansatz↩︎

We specialise to a spatially flat, homogeneous, and isotropic cosmology for both geometries, displaying explicit FLRW forms for the Einstein metric and its disformally related counterpart [22]. In the Einstein frame, \[ds^2 = -N(t)^2dt^2 + a(t)^2 (dx^2 + dy^2 + dz^2),\] and the biscalar is homogeneous, \(\phi = \phi(t)\), so \(\partial_{\mu}\phi = (\dot{\phi}, 0, 0, 0)\) (where an overdot \(\dot{}\) denotes \(d/dt\)). Using the disformal relation 2 , the components of the matter-frame metric are \[\hat{g}_{00} = -N(t)^2 + B(\phi)\dot{\phi}^2, \qquad \hat{g}_{ij} = a(t)^2\delta_{ij}.\] Hence, only the time-time component is altered between the frames, whilst the scale factor remains the same. The matter-frame line element thus takes the form \[\begin{align} d\hat{s}^2 &= -\hat{N}(t)^2dt^2 + \hat{a}(t)^2(dx^2 + dy^2 + dz^2), \end{align}\] where \[\hat{N}(t)^2 \mathrel{\vcenter{:}}= N(t)^2 - B(\phi)\dot{\phi}^2, \quad \hat{a}(t) \equiv a(t). \label{eq:disformal95lapse}\tag{44}\] Reality of the lapse requires \(\hat{N}^2>0\), i.e. \(N^2>B\dot{\phi}^2\), corresponding to the healthy domain \(1-2BX>0\). For completeness, the relative volume element is then \(\nu = \hat{N}/N\).

The proper times in each frame correspond to \[d\tau \mathrel{\vcenter{:}}= Ndt \qquad d\hat{\tau} \mathrel{\vcenter{:}}= \hat{N}dt \label{eq:propertimes}\tag{45}\] We thus define the Einstein- and matter-frame Hubble rates as \[H \mathrel{\vcenter{:}}= \frac{1}{a}\frac{da}{d\tau} = \frac{1}{a}\frac{da}{Ndt}, \qquad \hat{H} \mathrel{\vcenter{:}}= \frac{1}{a}\frac{da}{d\hat{\tau}} = \frac{1}{a}\frac{da}{\hat{N} dt},\] where, since \(d\hat{\tau} = \nu d\tau\), the two are related by \[H = \nu\hat{H}. \label{eq:H95Hhat95relation}\tag{46}\]

4.2 Cosmological Field Equations↩︎

Reducing the bimetric action with general constraints 8 to the spatially flat FLRW ansatz yields \[\begin{align} S = \int dt \Bigg[& -3 \frac{a \dot{a}^2}{N} + a^3 \frac{1}{2} \frac{\dot{\phi}^2}{\hat{N}} - \lambda_1 \left( N a^3 - \mathcal{C}_1 \right) - \lambda_2 \left( \hat{N} a^3 - \mathcal{C}_2 \right) \Bigg]. \label{eq:BiM95action95MSS} \end{align}\tag{47}\] Variation with respect to \(\lambda_{1,2}\) reproduces the constraints \[N a^3= \mathcal{C}_1, \qquad \hat{N} a^3 = \mathcal{C}_2,\] whose specific forms depend on the chosen constraints. Variation with respect to \(N\) gives the modified Friedmann equation \[3H^2 = \lambda_1 + \dfrac{1}{\nu} \left( \hat{X} + \lambda_2 \right). \label{eq:ModFriedmann}\tag{48}\] Varying with respect to \(a\) and eliminating \(\lambda_1\) using 48 gives the matter-frame Raychaudhuri equation \[\frac{d\hat{H}}{d\hat{\tau} } = \frac{B \hat{H}}{\nu^2} \frac{dX}{d\hat{\tau}} +\frac{\hat{X}}{\nu} \left[ B(\hat{X} - \lambda_2) - \frac{1}{\nu^2} \right], \label{eq:ray}\tag{49}\] whilst variation with respect to \(\phi\) gives the background biscalar equation, \[\frac{d}{dt}\left[ a^3 \frac{\dot{\phi}}{\hat{N}} \left(1 + B(\hat{X} + \lambda_2)\right) \right] = 0 , \label{eq:homo95biscalar95eom}\tag{50}\] which integrates once to give \[a^3\phi' \left(1 + B(\hat{X} + \lambda_2)\right) = C, \label{eq:homo95biscalar95eom95int}\tag{51}\] where \(C\) is a scalar-charge integration constant. Equations 49 51 together determine the background cosmology once the specific constraints are specified. We now look at the cosmologies of each of the theories in turn.

4.3 The Cosmological Dynamics of BUG↩︎

Specialising BUG to the spatially flat FLRW background, we find that the background dynamics reduce to a single first-order equation which we then use to investigate the dark energy equation of state and the possibility of phantom-divide crossing in the Einstein and matter frames. Hereafter we denote differentiation with respect to proper time by a prime \('\).

4.3.1 Homogeneous BUG Constraints↩︎

The fixed-density constraints 10 read \[N a^3 = \varepsilon_1, \qquad \hat{N} a^3 = \varepsilon_2,\] so that the relative volume element is simply the ratio of the two lapse functions, \[\nu(\hat{\tau}) = \frac{\hat{N}}{N} = \frac{\varepsilon_2(\hat{\tau})}{\varepsilon_1(\hat{\tau})}, \label{eq:BUG95nu95FLRW}\tag{52}\] which themselves are prescribed via the fixed densities \(\varepsilon_{1,2}(\hat{\tau})\). This is the defining feature of the BUG background 11 : \(\nu\) is not a dynamical variable but an externally prescribed function, and the entire scalar sector is determined by it. The kinematic constraints 12 16 fix the biscalar velocity and acceleration algebraically, \[\phi' = \sigma\sqrt{\frac{1-\nu^2}{B\nu^2}}, \qquad \hat{X} = \frac{1}{2}\phi'^2 = \frac{1-\nu^2}{2B\nu^2}, \qquad X' = -\frac{\nu}{B}\nu', \label{eq:BUG95phiprime95hatX}\tag{53}\] where \(\sigma = \pm 1\) labels the scalar branch. The evolution of the biscalar is inherited entirely from the prescribed \(\nu(\hat{\tau})\).

4.3.2 BUG Background Expansion↩︎

With \(\phi'\) and \(\hat{X}\) already fixed, the homogeneous biscalar equation 51 no longer propagates the field but instead determines the second multiplier, \[\lambda_2 = \frac{C}{B a^3\phi'} - \frac{1}{B} - \hat{X}, \label{eq:BUG95lambda295reconstruction95raw}\tag{54}\] which, using 53 , becomes \[\lambda_2 = \mu \frac{\nu}{a^3\sqrt{1-\nu^2}} - \frac{1+\nu^2}{2B\nu^2}, \qquad \mu \mathrel{\vcenter{:}}= \frac{\sigma C}{\sqrt{B}}. \label{eq:BUG95lambda295reconstruction}\tag{55}\] The first multiplier then follows from the Friedmann equation, using 46 , \[\lambda_1 = 3H^2 - \frac{\hat{X} + \lambda_2}{\nu} = 3H^2 - \frac{\mu}{a^3\sqrt{1-\nu^2}} + \frac{1}{B\nu}. \label{eq:BUG95lambda195reconstruction}\tag{56}\] Neither multiplier is therefore an independent—or even a combined—cosmological constant in the general BUG background. Both run with time, at rates set by \(\mu\), the expansion rate, and the prescribed relative volume element.

Substituting 53 and 55 into the Raychaudhuri equation 49 , every purely algebraic term cancels and the dynamics collapses to a single first-order equation, \[\hat{H}' = -\frac{\nu'}{\nu}\hat{H} - \frac{\mu}{2a^3}\frac{\sqrt{1-\nu^2}}{\nu^2}, \label{eq:BUG95Hhat95prime}\tag{57}\] which, together with \[a' = a\hat{H}, \label{eq:BUG95a95prime}\tag{58}\] closes the BUG background system.

In the Einstein frame it takes the simpler form \[H' = -\frac{\mu}{2a^3}\frac{\sqrt{1-\nu^2}}{\nu}, \label{eq:BUG95H95prime}\tag{59}\] and, using \(a' = aH/\nu\), may be recast directly as an equation in the scale factor, \[\frac{dH^2}{da} = -\mu \frac{\sqrt{1-\nu^2(a)}}{a^4}. \label{eq:BUG95H295a}\tag{60}\] The background has thus been reduced to a single quadrature. Imposing \(H(a=1) = H_0\) gives \[H^2(a) = H_0^2 + \mu\int_a^1 \frac{\sqrt{1-\nu^2(\bar a)}}{\bar a^4} d\bar a, \label{eq:BUG95H95solution}\tag{61}\] or, in dimensionless form with \(E(a) \mathrel{\vcenter{:}}= H(a)/H_0\), \[E^2(a) = 1 + \Gamma\int_a^1 \frac{\sqrt{1-\nu^2(\bar a)}}{\bar a^4} d\bar a, \qquad \Gamma \mathrel{\vcenter{:}}= \frac{\mu}{H_0^2}. \label{eq:BUG95E95solution}\tag{62}\] Once \(\nu(a)\) is given the expansion history is fully determined, with the single parameter \(\Gamma\) setting the amplitude and sign of the departure from de Sitter.

4.3.3 Frame Dependence of Phantom Crossing in BUG↩︎

The frame in which observables are defined depends on which metric the matter fields couple to. If ordinary matter is coupled to the matter metric \(\widehat g_{\mu\nu}\), \(S_{\rm m}[\widehat g_{\mu\nu}, \psi]\), then the expansion history inferred by physical observers is the matter-frame one. Since the two FLRW metrics share the same scale factor but differ in their proper times, up to an irrelevant constant normalisation of the dimensionless Hubble rates, we have \[\widehat E^2(a)=\frac{E^2(a)}{\nu^2(a)} .\] The effective matter-frame equation of state reconstructed from the background expansion is therefore \[\widehat w_{\rm BUG}(a) \mathrel{\vcenter{:}}= -1-\frac{1}{3}\frac{d\ln \widehat E^2}{d\ln a}. \label{eq:wBUG95matter95def}\tag{63}\] Equivalently, \[\widehat w_{\rm BUG}(a) = -1-\frac{1}{3}\frac{d\ln E^2}{d\ln a} +\frac{1}{3}\frac{d\ln \nu^2}{d\ln a} = w_{\rm BUG}(a) +\frac{1}{3}\frac{d\ln \nu^2}{d\ln a}. \label{eq:wBUG95frame95relation}\tag{64}\]

For the BUG solution 62 , one finds \[\begin{align} \widehat w_{\rm BUG}(a)+1 &= w_{\rm BUG}(a)+1 +\frac{1}{3}\frac{d\ln \nu^2}{d\ln a}, \tag{65} \\[0.5em] w_{\rm BUG}(a)+1 &= \frac{\Gamma}{3E^2(a)} \frac{\sqrt{1-\nu^2(a)}}{a^3}. \tag{66} \end{align}\] For any physical BUG background, \(E^2>0\) and \(0<\nu^2\leq 1\). Hence \(\sqrt{1-\nu^2}/a^3\) is non-negative, and the sign of \(w_{\rm BUG}+1\) is fixed entirely by the sign of \(\Gamma\). The Einstein-frame BUG expansion may therefore be quintessence-like for \(\Gamma>0\) or phantom-like for \(\Gamma<0\), but it cannot dynamically cross the phantom divide.

The matter-frame equation of state, however, is different. Its distance from the phantom divide is \[\widehat w_{\rm BUG}(a)+1 = \frac{1}{3} \left[ \frac{\Gamma}{E^2(a)} \frac{\sqrt{1-\nu^2(a)}}{a^3} + \frac{d\ln \nu^2}{d\ln a} \right]. \label{eq:wBUG95hat95crossing95condition}\tag{67}\] Consequently, a matter-frame phantom crossing occurs whenever \[\frac{d\ln \nu^2}{d\ln a} = -\frac{\Gamma}{E^2(a)} \frac{\sqrt{1-\nu^2(a)}}{a^3}. \label{eq:wBUG95hat95crossing}\tag{68}\] Thus the prescribed evolution of the relative volume element can drive \(\widehat w_{\rm BUG}+1\) through zero even though \(w_{\rm BUG}+1\) has a fixed sign in the Einstein frame.

This crossing is therefore not an intrinsic Einstein-frame phantom crossing. Rather, it is a frame-dependent effect induced by the changing relative volume element between the Einstein and matter metrics. If matter is universally coupled to \(\widehat g_{\mu\nu}\), then observers would reconstruct the matter-frame quantity \(\widehat w_{\rm BUG}\) and could infer a genuine background-level crossing of the phantom divide. Its origin, however, would lie in the time-dependent relative volume element between the two frames, not in a sign change of the Einstein-frame BUG source.

4.3.4 Power-Law Relative Volume Element↩︎

The residual freedom in BUG is carried entirely by the prescribed relative volume element \(\nu(a)\). To illustrate the frame dependence of the reconstructed equation of state, we consider the simple power-law profile \[\nu^2(a)=1-q a^n, \qquad 0<q a^n<1, \label{eq:BUG95powerlaw95nu}\tag{69}\] where the inequality specifies the physical domain in which the matter-frame lapse remains real. For \(0<q<1\) and \(n>0\), this condition is automatically satisfied over the interval \(0<a\leq 1\). For \(n<0\), the ansatz should instead be regarded as a late-time profile over the range in which \(q a^n<1\).

With 69 , the solution 62 reads⁵ \[E^2(a) = 1 + \Gamma\sqrt q \frac{1-a^{n/2-3}}{n/2-3}, \qquad n\neq 6. \label{eq:BUG95powerlaw95solution}\tag{71}\] The Einstein-frame equation of state therefore obeys \[w_{\rm BUG}(a)+1 = \frac{\Gamma\sqrt q}{3E^2(a)} a^{n/2-3}. \label{eq:wBUG95powerlaw95Einstein}\tag{72}\] The matter-frame correction is \[\frac{1}{3}\frac{d\ln\nu^2}{d\ln a} = -\frac{nq a^n}{3(1-q a^n)}. \label{eq:BUG95powerlaw95frame95term}\tag{73}\] Hence the matter-frame equation of state is \[\widehat w_{\rm BUG}(a)+1 = \frac{1}{3} \left[ \frac{\Gamma\sqrt q}{E^2(a)} a^{n/2-3} - \frac{nq a^n}{1-q a^n} \right]. \label{eq:wBUG95powerlaw95matter}\tag{74}\] A matter-frame phantom crossing occurs at any scale factor \(a=a_\star\) satisfying \[\frac{\Gamma\sqrt q}{E^2(a_\star)} a_\star^{n/2-3} = \frac{nq a_\star^n}{1-q a_\star^n}. \label{eq:BUG95powerlaw95crossing95condition}\tag{75}\] Since \(E^2>0\) and \(1-q a^n>0\) in the physical domain, this condition can be satisfied only when \(\Gamma\) and \(n\) have the same sign, \(\Gamma n>0\).

The direction of the crossing is then controlled by the sign of \(n\). For \(n>0\) and \(\Gamma>0\), the Einstein-frame solution is quintessence-like, whilst the matter-frame correction is negative. The matter frame may therefore cross from quintessence-like to phantom-like behaviour if the evolution of \(\nu\) becomes sufficiently important. Conversely, for \(n<0\) and \(\Gamma<0\), the Einstein-frame solution is phantom-like, whilst the matter-frame correction is positive. This is the branch on which a matter-frame observer can reconstruct a transition from phantom-like to quintessence-like expansion.

To illustrate this frame-dependent crossing explicitly, we plot in Fig. 1 the Einstein- and matter-frame equations of state for the same BUG background solution. We choose the inverse-power branch \(n<0\) with \(\Gamma<0\), for which the Einstein-frame source is everywhere phantom-like. As expected from 74 , the Einstein-frame equation of state remains below \(-1\) throughout the interval shown and therefore does not cross the phantom divide. By contrast, the matter-frame reconstruction does cross from phantom-like to quintessence-like behaviour. This crossing is not caused by a sign change of the Einstein-frame BUG source, but by the additional matter-frame contribution induced by the evolving relative volume element. Thus the same physical background can appear crossing or non-crossing depending on the frame in which the equation of state is reconstructed.

4.4 The Cosmological Dynamics of BHT↩︎

We now turn to the BHT realisation of bimodular gravity. Unlike BUG, the HT constraints do not prescribe the relative volume element, but instead enforce the two bimodular multipliers to be independent constants. The biscalar therefore remains fully dynamical, and the resulting theory is a more standard shift-symmetric k-essence model supplemented by a cosmological constant and a constant scalar-sector offset. We shall see that, whilst this branch can give non-trivial k-essence cosmology, it cannot dynamically cross the phantom divide.

4.4.1 Homogeneous BHT Constraints↩︎

Under the FLRW ansatz, with \(\tau_i^\mu=\tau_i^\mu(t)\), the volume constraints become \[Na^3=\dot{\tau}_1, \qquad \widehat N a^3=\dot{\tau}_2. \label{eq:BHT95flrw95conts95t}\tag{76}\] Unlike in BUG, the auxiliary densities \(\tau_i\) are dynamical objects. Their derivatives are not prescribed functions, and hence the relative volume element is not fixed a priori. The scalar therefore remains a genuine dynamical degree of freedom. The defining BHT condition is instead the constancy of the two multipliers, \[\lambda_1=\Lambda, \qquad \lambda_2=V_0, \label{eq:BHT95constants95cosmology}\tag{77}\] where \(\Lambda\) is the usual unimodular cosmological constant and \(V_0\) is a constant offset in the scalar sector.

4.4.2 BHT Background Expansion↩︎

For completeness, the Friedmann equation gives the algebraic constraint \[3\nu^2\widehat H^2 = \Lambda + \frac{1}{\nu} \left( \frac{1}{2}\phi'^2+V_0 \right), \label{eq:BHT95friedmann95constraint}\tag{78}\] though this has already been incorporated within the general matter-frame Raychaudhuri equation, which becomes \[\widehat H' = B\nu^2\widehat H\phi'\phi'' - \frac{\phi'^2}{2\nu} \left[ 1+BV_0+\frac{B}{2}\phi'^2 \right]. \label{eq:BHT95Hhat95ODE95phi95second}\tag{79}\] The biscalar equation follows by differentiating the first integral 51 . Since \(V_0\) is constant, this gives \[\phi'' = -3\widehat H\phi' \frac{ 1+B\left(V_0+\frac{1}{2}\phi'^2\right) }{ 1+B\left(V_0+\frac{3}{2}\phi'^2\right) }. \label{eq:BHT95phi95ODE}\tag{80}\] Together with \[a'=a\widehat H, \label{eq:BHT95a95ODE}\tag{81}\] equations 79 81 define the BHT background system.

4.4.3 Absence of Phantom Crossing in BHT↩︎

The possibility of phantom crossing in BHT may be determined directly from the matter-frame Raychaudhuri equation. In the FLRW background, we find that the numerator and denominator of the scalar speed become \[N = 1+B\left(V_0+\frac{1}{2}\phi'^2\right), \qquad D = 1+B\left(V_0+\frac{3}{2}\phi'^2\right).\] We then find that the scalar equation 80 takes the compact form \[\phi''=-3\widehat H\phi'\frac{N}{D}.\] Substituting this into 79 gives \[\widehat H' = -\phi'^2 A \left[ \frac{1}{2\nu} + \frac{3B\nu^2\widehat H^2}{D} \right]. \label{eq:BHT95Hprime95sign}\tag{82}\] Thus, on the branch free of ghosts and gradient instabilities, we have \(N,D>0\) (see Sec. 3.2.4). Since \(B>0\), \(\nu>0\), and \(\widehat H^2>0\), the bracket in 82 is positive. Hence \[\widehat H'\leq 0,\] with equality only on the frozen-scalar branch \(\phi'=0\) or at the boundary \(N=0\).

The matter-frame dark energy equation of state satisfies \[\widehat w_{\rm BHT}+1 = -\frac{2}{3}\frac{\widehat H'}{\widehat H^2}.\] Therefore, throughout the healthy rolling branch, \[\widehat w_{\rm BHT}+1\geq0.\] BHT can approach the de Sitter limit \(\widehat w_{\rm BHT}=-1\) as \(\phi'\to0\), but it cannot dynamically cross the phantom divide. Thus, unlike BUG, the BHT cosmology does not naturally realise a phantom crossing; it behaves as a more standard shift-symmetric k-essence theory.

5 Diffeomorphism Invariant Bimodular Gravity↩︎

So far we have made explicit the inequivalence between unimodularising bimetric gravity with the fixed-determinant constraints, BUG, and HT-type constraints, BHT. The former realises a genuinely novel bimodular theory, but one with broken diffeomorphism invariance. This prompts a natural question: can we construct a fully diffeomorphism-invariant completion that preserves the BUG balance law 21 and its associated kinematic constraint on \(X\) 12 ? In the remainder of this section we answer this question by developing a diffeomorphism-invariant formulation of bimodular gravity.

5.1 Motivation↩︎

For quantisation and for coupling to matter sectors in a standard covariant path integral, it is generally preferred to work with a formulation that is fully diffeomorphism invariant at the level of the action [10]. In the single-metric case, one may pass from UG to HT (or DUG) without concern, for they are classically equivalent and reproduce the same field equations as \(\mathrm{GR}\!+\!\Lambda\) [10], [28], [29]. In the present disformal bimetric setting, however, this equivalence is broken.

The general attempt to covariantise bimodular gravity by importing the HT/DUG formalism separately for each metric enforces both \(\lambda_1\) and \(\lambda_2\) to be strict integration constants. Conversely, in BUG, the individual multipliers remain dynamically correlated through the balance law. Treating the two HT sectors independently therefore overconstrains the multipliers—freezing both \(\lambda_{1,2}\)—and simultaneously underconstrains the relative volume element, leaving \(X\) as a fully dynamical variable. The biscalar then behaves as a standard, freely propagating degree of freedom with a genuine second-order evolution equation.

Hence, two formalisms that are equivalent in the single-metric case become on-shell inequivalent once two disformally related metrics are present. This motivates the search for a diffeomorphism-invariant reformulation that is tailored specifically to the bimodular structure which:

1. realises full diffeomorphism invariance off shell;

2. enforces the balance relation \(\nabla_\mu \lambda_1 = -\nu\nabla_\mu\lambda_2\) on shell;

3. reproduces the BUG kinematic constraint on \(\nu\) and subsequently on \(X=X(x;\phi)\);

4. reduces to BUG in an appropriate gauge, and hence is classically equivalent to BUG.

In what follows we present precisely such an action. One may regard it as a diffeomorphism-invariant completion of unimodularised bimetric scalar-tensor gravity, whereby full covariance is restored and the physical content of BUG is retained on shell.

5.2 Restoring Diffeomorphism Invariance↩︎

To this end, following the Stueckelberg covariantisation of generalised unimodular gravity [37] and the standard use of Stueckelberg fields to restore broken gauge symmetries [41], we introduce four scalar Stueckelberg fields \(Z^A(x)\), \(A=0,\ldots,3\), and define the Jacobian density \[\mathcal{J}[Z] \coloneq \det\!\left(\partial_\mu Z^A\right).\] Let \(\bar\epsilon_1(Z)\) and \(\bar\epsilon_2(Z)\) be two fixed volume densities on the internal \(Z^A\)-space, with \[\bar\epsilon_2(Z)=\bar\nu(Z)\bar\epsilon_1(Z).\] Their pullbacks to spacetime define the scalar densities \[\mathcal{E}_i[Z] \coloneq \bar\epsilon_i(Z)\mathcal{J}[Z], \qquad i=1,2.\] Equivalently, \(\mathcal{E}_i[Z]\) are the pullbacks of two fixed four-forms on the internal Stueckelberg space. The diffeomorphism-invariant BUG constraint sector is then \[S \supset -\int d^4x \left[ \lambda_1\left(\sqrt{-g}-\mathcal{E}_1[Z]\right) + \lambda_2\left(\sqrt{-\hat{g}}-\mathcal{E}_2[Z]\right) \right]. \label{eq:DBUG95action}\tag{83}\] This action is now fully diffeomorphism invariant.⁶

Variation with respect to the multipliers imposes \[\sqrt{-g}=\mathcal{E}_1[Z], \qquad \sqrt{-\hat{g}}=\mathcal{E}_2[Z].\] Taking the ratio gives \[\frac{\sqrt{-\hat{g}}}{\sqrt{-g}} = \frac{\mathcal{E}_2[Z]}{\mathcal{E}_1[Z]} = \nu(Z). \label{eq:stuck95rve}\tag{84}\] Thus the relative volume element is prescribed covariantly. In the unitary gauge \[Z^A=x^A,\] one has \(\mathcal{J}=1\), \(\mathcal{E}_i=\bar\epsilon_i(x)\), and hence \[\bar\nu(x) = \frac{\bar\epsilon_2(x)}{\bar\epsilon_1(x)} = \sqrt{1-2BX}.\] The usual BUG kinematic constraint is therefore recovered: \[X = \frac{1-\bar\nu^2(x)}{2B}.\]

It remains to check the multiplier equation. Let \(e^\mu{}_A\) denote the inverse Jacobian, satisfying \(e^\mu{}_A \partial_\mu Z^B=\delta_A{}^B\). Since \(\mathcal{E}_i[Z]\) are pullbacks of four-forms, their variations are total derivatives, \[\delta \mathcal{E}_i = \partial_\mu \left( \mathcal{E}_i e^\mu{}_A \delta Z^A \right).\] Variation of 83 with respect to \(Z^A\) therefore gives \[\mathcal{E}_1 e^\mu{}_A\partial_\mu\lambda_1 + \mathcal{E}_2 e^\mu{}_A\partial_\mu\lambda_2 =0,\] which, for \(\mathcal{J} \neq0\), is equivalent to \[\partial_\mu\lambda_1 + \nu(Z) \partial_\mu\lambda_2 = 0,\] where we have used 84 . In the unitary gauge, \(\nu(Z)=\bar \nu(x)\), and this becomes precisely the BUG balance law 21 . Consequently, the effective vacuum contribution \(\Pi(x)\) obeys 23 and is therefore not conserved for a varying relative volume element.

This construction is therefore to be understood as a Stueckelberg completion of Bimodular Gravity. The Stueckelberg fields restore full diffeomorphism invariance, whilst their equations of motion reproduce the BUG multiplier balance law and restricted dynamics \(X=X(x;\phi)\).

6 Conclusion and Discussion↩︎

We have unimodularised a simple bimetric scalar–tensor theory, in which gravity and matter couple to two disformally related metrics, by constraining the volume element of each. Although such unimodular formalisms dynamically coincide when only one metric is present, the procedure becomes inequivalent here. The fixed-determinant route prescribes the relative volume element \(\nu(x)\), which in turn locks the scalar kinetic term to a function of spacetime, \(X = X(x;\phi)\), and ties the two Lagrange multipliers \(\lambda_{1,2}\) into a single vacuum contribution \(\Pi(x)=\lambda_1(x) + \nu(x)\lambda_2(x)\) that runs as \(\nabla_\mu\Pi(x) \propto\nabla_\mu\nu(x)\) rather than remaining a rigid constant. We have termed this class of dynamics Bimodular Gravity. A unique cosmological result of bimodularity is a frame-dependent crossing of the phantom divide: whilst impossible in the Einstein-frame, an observer coupled to the matter metric would see \(w_{\rm DE}\) dynamically cross \(-1\), driven entirely by the evolving ratio between the two volumes. The alternative route, which retains full diffeomorphism invariance through Henneaux–Teitelboim (HT) constraints, does not realise this behaviour. Instead, it returns a more standard k-essence theory carrying a genuine cosmological constant and a constant scalar offset in the matter frame, with no phantom crossing. Finally, we have given a Stueckelberg completion that restores diffeomorphism invariance to Bimodular Gravity whilst reproducing its dynamics on shell.

Thus, Bimodular Gravity supplies a simple geometric route to the recently inferred phantom crossing, following from the evolving ratio between the two metrics that matter and gravity respectively see. Furthermore, the unimodular reinterpretation of \(\Lambda\) as an integration constant is carried one step further, with the constant being replaced by a contribution that evolves with the very volume element that drives the crossing, so that the two features share a common origin. Of course, \(\nu(x)\) is prescribed rather than derived, so the crossing is engineered by a phenomenological choice rather than predicted; and, as with ordinary unimodular gravity, nothing here speaks to the smallness of the vacuum energy itself.

It should be emphasised that the bimodular constraint 12 arises specifically from imposing two independent fixed-determinant constraints. If only one metric were constrained, i.e. if one unimodularised either \(g_{\mu\nu}\) or \(\hat{g}_{\mu\nu}\) alone, then the determinant of the remaining metric would be free to absorb variations in the relative volume element 5 , and the biscalar kinetic term would not be prescribed. Similarly, if mixed constraints were implemented, with one fixed-determinant constraint and one HT constraint, the biscalar dynamics would remain free because the relative volume element would not be completely fixed. In that case, the constancy of the HT multiplier would, through the balance law 21 , enforce the constancy of the remaining multiplier as well. The resulting theory would therefore fall on the side of a particular, but more standard, k-essence model, rather than Bimodular Gravity. Only when the two metric determinants are fixed, or equivalently fixed up to their covariant Stueckelberg completion, is the relative volume element prescribed and the bimodular kinetic constraint obtained.

It is also worth noting that the HT implementation admits a further, multiplier-free presentation. In [42], unimodularity is obtained through a mimetic-like metric redefinition involving a vector field of conformal weight four, yielding a Weyl-invariant higher-derivative vector–tensor theory which reduces, in gauge-invariant variables, to the HT formulation. In the bimetric setting, applying such a construction would therefore amount to a Weyl-invariant rewriting of the BHT-type theory rather than of BUG. The two integration constants would remain independently conserved, whilst the biscalar would remain fully dynamical.

The dynamics of Bimodular Gravity sits close to, but distinct from, Mimetic Gravity [32], [43] and Dust of Dark Energy [31]. In all three cases the scalar kinetic density is fixed algebraically, and hence the scalar sector is not described by an ordinary freely propagating k-essence mode. The similarities are nevertheless hierarchical. For constant \(B\) and constant \(\nu\), the bimodular constraint reduces to \(X=\mathrm{const}\), mimicking the mimetic constant-norm constraint, up to conventions. Allowing \(B=B(\phi)\), and more generally \(\nu(\phi)\), instead gives \(X=X(\phi)\), which is the kinematic structure underlying Dust of Dark Energy. Here, the pressure gradient is then aligned with the scalar flow, energy follows timelike geodesics, and the physical sound speed vanishes identically. Bimodular Gravity, however, corresponds to the further extension in which the prescribed relative volume element is a fixed spacetime density, \(\nu(x)\). Even with general \(B(\phi)\), the kinetic term is then locked to \(X=X(x;\phi)\), so space-time gradients of \(\nu\) source a non-zero acceleration of the scalar flow. Consequently, Bimodular Gravity is neither standard k-essence nor generically dusty dark energy. Indeed, its perturbations remain constrained rather than freely propagating, but the dust argument for an identically vanishing sound speed no longer applies.

It is likewise distinct from the mimetic k-essence construction [33]. There, the mimetic mixing of a scalar sector with a unimodular-type sector is shown to be on-shell equivalent to GR coupled to a k-essence scalar, whose overall scale appears as a conserved integration constant. Bimodular Gravity differs in both respects. Firstly, its scalar sector is not a freely propagating k-essence mode. Moreover, the vacuum contribution is inherited from two unimodular multipliers which, for generic \(\nu(x)\), are not individually conserved and cannot be absorbed into a single constant—\(\Pi(x)\) is instead a spacetime-dependent vacuum contribution.

Throughout, we have neglected the effects of matter, which must of course be included to take the cosmology seriously. The metric to which it couples is a genuine question, but the frame-dependent phantom crossing found here motivates a coupling to the matter metric \(\hat{g}_{\mu\nu}\), with \(S_{\rm m} = S_{\rm m}[\hat{g}_{\mu\nu},\psi]\). If matter then couples disformally to \(\hat{g}_{\mu\nu}\), the viability of the construction would be further subject to local and screening constraints which we do not address here. A fuller assessment must therefore await a treatment of perturbations and ultimately a confrontation with data, to which we hope to return. It would also be desirable to render the prescribed \(\nu(x)\) with a physical origin, so that the crossing is predicted rather than engineered. Finally, since part of the appeal of unimodular gravity in quantum cosmology lies in the emergence of a unimodular time variable [44], [45], the analogous notion of a bimodular time merits a dedicated Hamiltonian analysis, which would also clarify the precise propagating degrees of freedom. The Stueckelberg completion of Sec. 5 provides a natural starting point for such a treatment. But of course, we leave these questions for future work.

It brings us much pleasure in giving thanks to T. Baker, R. Terrazas-Santamaria, and A. Ridley for providing helpful comments and interesting conversations—in relation to this work and in general. We are also grateful to the anonymous referee, whose detailed and constructive comments substantially improved this work. We also thank T. Lawrence, M. Hatteea, and H. Wells for careful proofreading of the manuscript. J. Hallam is supported by an STFC studentship; J.Magueijo was partly supported by STFC Consolidated Grant ST/T000791/1. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising. Supporting research data are available on reasonable request from the corresponding authors.

References↩︎

---

1. Since, at fixed \(\phi\), \(\delta\hat{g}_{\mu\nu}=\delta g_{\mu\nu}\), one obtains \[\delta S_\phi = \frac{1}{2} \int d^4x \sqrt{- \hat{g}}\, \left[ \hat{T}^{\mu\nu} \right] \delta\hat{g}_{\mu\nu} = \frac{1}{2} \int d^4x \sqrt{- g} \, \left[ \nu \hat{T}^{\mu\nu} \right] \delta g_{\mu\nu}.\] Therefore \(T^{\mu\nu} = \nu \hat{T}^{\mu\nu}\) for the contravariant stress tensors, with the important caveat that \(\hat{T}^{\mu\nu}\) is raised using \(\hat{g}^{\mu\nu}\), not \(g^{\mu\nu}\). Equivalently, for the Einstein-frame covariant components, \(T_{\mu\nu} = \nu g_{\mu\kappa}g_{\nu\lambda}\hat{T}^{\kappa\lambda}\).↩︎

2. Since some components of \(\hat{\Box}\phi\) are also constrained by 16 , the kinematic constraint must be treated consistently in order for the dynamics to remain self-consistent [34].↩︎

3. The matter frame stress-energy tensor is found to be \[\hat{T}_{\mu\nu} = \big (\hat{X} - \lambda_2\big)\hat{g}_{\mu\nu} + \partial_\mu \phi \partial_\nu \phi.\]↩︎

4. The definition is deliberately dynamical, rather than tied to a particular off-shell presentation, because he same structure may be implemented by fixed density constraints, auxiliary vector densities, or equivalent covariant reformulations (see Sec. 5).↩︎

5. For the special case \(n=6\), the corresponding solution is instead \[E^2(a)=1+\Gamma\sqrt q \ln\frac{1}{a}. \label{eq:BUG95powerlaw95solution95n6}\tag{70}\] ↩︎

6. This essentially mirrors the Stueckelberg logic used in generalized unimodular gravity as k-essence [37], but not its dynamics. There the covariantised theory contains a genuine k-essence scalar; here the Stueckelberg fields restore diffeomorphism invariance whilst preserving the BUG constraint and multiplier balance law. In fact, BHT is closer to [37], since both contain a genuine k-essence-type scalar with standard perturbative sound speed. In BHT, however, the two HT multipliers remain separate constants, \(\Lambda=\lambda_1\) and \(V_0=\lambda_2\).↩︎