# The trimer paradox: the effect of stiff constraints on equilibrium
distributions in overdamped dynamics

Radost Waszkiewicz<EMAIL_ADDRESS>Institute of Theoretical
Physics, Faculty of Physics, University of Warsaw, L. Pasteura 5, 02-093
Warsaw, Poland  Maciej Lisicki<EMAIL_ADDRESS>Institute of Theoretical
Physics, Faculty of Physics, University of Warsaw, L. Pasteura 5, 02-093
Warsaw, Poland

###### Abstract

We reconsider the classical problem of a freely joined chain of Brownian
particles connected by elastic springs and study its conformational
probability distribution function in the overdamped regime in the limit of
infinite stiffness of constraints. We show that the well-known solution by
Fixman is missing a shape-related term, later alluded to but not computed by
Helfand. In our approach, the shape term, also termed zero-point energy,
arises explicitly from a careful treatment of the distributional limit. We
present a computationally feasible method of calculation of the shape term and
demonstrate its validity in a couple of examples.

##  I Introduction

Molecular Dynamics and Brownian Dynamics simulations are now standard tools
for detailed modelling of a plethora of molecular and mesoscopic systems,
where structural complexity poses a challenge to theoretical calculations.
Upon the introduction of a suitable coarse-graining scheme to represent the
composition of a given molecule, its conformational space is prescribed by
specifying intramolecular interactions, such as chemical bonds or
electrostatic forces, between the subunits. The nature of these interactions
endows molecular models with flexibility, which then in turn affects their
conformational variability, as well as their macroscopic and statistical
properties, such as rheology, diffusivity, and thermodynamics, and also
biological function.

A popular idealisation for a wide variety of systems is purely mechanical,
where the molecule is represented by beads connected with springs, an idea
which probably originated from Kramers [1], who built on concepts of Kuhn, and
termed it the pearl-necklace model. Variations of this idea are now the
cornerstone of polymer physics [2], but have also been used to describe other
classical many-body systems. Physically, the introduction of intramolecular
interactions leads to constrained dynamics with multiple time scales.
Typically, the characteristic time scales associated with internal vibrations
of molecular bonds are much shorter than those of translational and rotational
motion or bond angle dynamics. The disparity of time scale of relaxation of
hard (or stiff) degrees of freedom (such as bond lengths) in comparison to
soft degrees of freedom (such as bond angles) is a source of stiffness in the
problem and hinders fast numerical simulation of such systems [3].

In many applications, the dynamics of the hard degrees of freedom are of
secondary importance, and one method of circumventing this difficulty is to
treat them as fully constrained (i.e. rigid), thus eliminating the short
timescale, and enabling faster calculation. This is the basis of many
algorithms, such as SHAKE [4] and its further extensions. During the
development of such rigid models, a disconnect was discovered between the
equilibrium bond angle distributions in the rigid simulations and the
simulations with bonds modelled as very stiff springs.

The simplest example in which the complexity of the constrained dynamics can
be appreciated is the case of a flexible trimer: a hypothetical molecule
containing 3 subunits and 2 harmonic bonds, depicted in Fig. 2a. The angle
between the two bonds is denoted by $\psi$. In the classical textbook [3],
Frenkel and Smit present an apparently paradoxical result: the marginal
distribution of the bond angle in the case of rigid bonds
$\mathrm{d}p_{\text{rigid}}$ and the limiting distribution of spring-like
bonds where the spring stiffness $k$ is taken to infinity
$\mathrm{d}p_{k\to\infty}$ do not coincide. More precisely, they report

$\displaystyle\mathrm{d}p_{k\to\infty}$ $\displaystyle=$ $\displaystyle
N_{1}\sin\psi\,\,\mathrm{d}\psi$ (1) $\displaystyle\mathrm{d}p_{\text{rigid}}$
$\displaystyle=$ $\displaystyle
N_{2}\sin\psi\sqrt{{1-\frac{\cos^{2}\psi}{4}}}\,\,\mathrm{d}\psi$ (2)

for appropriate normalisation constants $N_{1}$ and $N_{2}$. Qualitatively,
distribution (1) corresponds to a spherically uniform distribution of the
second bond direction when the coordinates are aligned with the first bond,
while in Eq. (2) it seems that an additional force resists bond alignment. In
an early and widely cited work, Fixman [5] presented a derivation to argue
that the ratio of the two probability densities, Eqs. (1) and (2), can be
computed with the knowledge of the constraining surfaces alone and does not
require knowledge of the shape of the confining potential. This presupposes
that both of these distributions are well defined whenever the constraint is
given, which we argue is not the case.

In this contribution, we present a clear mathematical procedure leading to the
aforementioned limits in the overdamped regime when inertial effects can be
neglected. We show that the distribution given by Eq. (1) is not completely
determined by fixing the bond lengths. To correctly predict the limiting
distribution, knowledge about the nature of confinement is needed and its
structure survives in the final expression. Furthermore, we demonstrate that
Fixman’s expressions for the stiff limit miss a term, which is important even
in the case of harmonic springs.

##  II Stiff vs. rigid constraints

The first calculations of the equilibrium properties of flexible polymers can
be traced to the papers by Kramers [1], Gō and Sheraga [6] and Fixman [5].
Results of the latter paper were used essentially without changes in works
such as [7, 8, 9, 10]. Early papers focused on the treatment of constraints in
the Lagrangian and Hamiltonian pictures, both classical and quantum-
mechanical, by defining the soft and stiff degrees of freedom and integrating
out generalised coordinates corresponding to the hard variables. The limiting
distribution is then obtained for infinitely stiff constraints when stiff
variables are eliminated.

Fixman [5] presented a purely classical calculation of the limiting
distribution, correctly writing the action of the distribution on a test
function but erroneously assuming that the constraining potential depends only
on the hard coordinates. As a result, details of the confining potential
cannot appear in his final expressions, which contain only the projected
phase-space volume element terms. Gō and Sheraga [6] treated the problem
quantum-mechanically and correctly identified two contributions to the
equilibrium probability density of the soft coordinates: the projected volume
elements and the zero-point energy of the vibrational motion of the hard
coordinates. However, they then stated that the zero-energy contribution in
typical molecular systems varies only a little when changing the soft
coordinates and may be altogether neglected. We show in the following that
this is not necessarily the case in overdamped dynamics. By focusing on a
concrete classical two-dimensional example with one constraint, Helfand [11]
correctly identified the rigid-rod type distribution as a uniform distribution
on the constraining manifold and computed the stiff-spring limit in a similar
simple case. Unfortunately, these results are not expressed in a form
applicable to higher-dimensional examples with multiple constraints.

Due to this lack of generality, later works such as that of Hinch [8], cite
Helfand as a reference to an erroneous claim that the bond angle distribution
becomes non-uniform in the case of rigid-rod system. The pitfall lies in the
assumption that the distribution should be uniform under an action of a
rotation matrix on a single bond, presumably by analogy with a dimer. However,
there is a distinct difference between these two cases; In the case of a
dimer, the true physical symmetries of the equations of motion – rigid body
rotation of the whole system and translations – act transitively on the space
of all possible configurations, thus the equilibrium distribution can be
determined from symmetry arguments. More generally, for a rigid body the
rotation group $\mathrm{SO}(3)$ can be treated as a topological manifold,
resulting in a natural notion of uniformity, with the uniform distribution
given by the Haar measure. For the trimer, however, we have configurations
that are meaningfully distinct. There is no symmetry that changes the bond
angle $\psi$. While the equilibrium distribution for 6 out of the 7 degrees of
freedom can be determined from the six symmetries of translation and rotation,
the relative weight of configurations with different $\psi$ has to refer to
the physical problem, not just geometric considerations.

Hinch [8] approached the problem differently, by manipulating the Langevin
equations, which are notoriously difficult to handle, and erroneously assumed
that details of the confining potential are of little importance to the
equilibrium distribution of the constrained configurations. Van Kampen and
Lodder [12] noted that the approach of Helfand [11] is generally applicable to
constrained systems and comment that the discrepancy between the stiff-spring
and rigid-rod distributions in the case of the trimer molecule is due to the
fact that _the width of this gully is not the same everywhere_ [12]. However,
they did not quantify the effect of this width in a higher-dimensional case.
As a result, the knowledge of the influence of the shape of the constraining
potential appears to have been lost as multiple works (for example, [13, 14,
3]) mention only the volume element densities and omit the _zero point energy_
or, equivalently, the potential _shape_ terms from their descriptions of
constrained dynamics, widely using the result of Fixman [5] as reference.

In the following, we show a way to compute the limiting distribution by
integrating out the hard degrees of freedom in the overdamped limit. We
explicitly find the metric contribution arising from the transformation of
coordinates and the shape term that represents the details of the constraining
potential and which remains imprinted on the configurational distribution when
the hard coordinates are taken to be infinitely stiff.

##  III The stiff spring limit

Figure 1: A flexible polymer with $N$ degrees of freedom is modelled as beads
connected with springs. Conformations of the polymer are described by a vector
$\bm{q}$. The potential energy of the confining springs, $k^{2}W(\bm{q})$,
attains a minimum on a submanifold $\mathcal{M}$. We define an orthonormal
system of coordinates on a compact neighbourhood of $\mathcal{M}$, enclosed by
the red solid line, by introducing soft degrees of freedom $\bm{r}$ on
$\mathcal{M}$ and hard degrees of freedom $\bm{s}$, which are normal to
$\mathcal{M}$ and correspond to the confining springs. The hard degrees of
freedom will be integrated out to find the limiting distribution for stiff
springs.

We present a derivation of the limiting distribution in the following setting,
sketched in Fig. 1. Consider a molecule (or a polymer) with a handful of
subunits (atoms, beads, or monomers) and $N$ degrees of freedom in a heat bath
at constant temperature $T$, with $\beta=1/k_{B}T$. The conformation of the
molecule is described by an $N$-dimensional vector $\bm{q}$. We parameterise
conformations by Cartesian coordinates of monomers $q_{i}$ with $i\in 1\dots
N$. The potential energy of the molecule has two components: the conformation-
dependent energy $U(\bm{q})$ and the confining (springs) potential
$k^{2}W(\bm{q})$ where $k$ is a large parameter describing the spring
stiffness. Note that $k$ here would correspond to the square root of a
harmonic spring constant. The probability density $p_{k}$ of the equilibrium
distribution of the system configuration with respect to the Lebesgue measure
$\mathrm{d}\bm{q}$ is given by the Boltzmann distribution

$p_{k}=N(k)\exp\left\\{-\beta\left[U(\bm{q})+k^{2}W(\bm{q})\right]\right\\},$
(3)

with an appropriate normalising function $N(k)$. Suppose now that the spring
potential is non-negative, $W(\bm{q})\geq 0$, and attains a minimum for a
configuration lying on a smooth submanifold $\mathcal{M}$ of dimension $M<N$,
so that $W(\bm{q})=0\iff\bm{q}\in\mathcal{M}$. To determine the weak limit of
$p_{k}$, we consider integrals $I_{k}$ for a compactly supported smooth test
function $\phi(\bm{q})$ given by

$I_{k}=\int p_{k}(\bm{q})\phi(\bm{q})\mathrm{d}\bm{q}.$ (4)

Since $\mathcal{M}$ is smooth, we can use the tubular embedding theorem to
define new orthonormal coordinates
$\bm{q}\mapsto(\bm{r},\bm{s})=(r_{1},r_{2},\dots,r_{M},s_{1},s_{2},\dots)$ in
the vicinity of $\mathcal{M}$, where $\bm{q}\in\mathcal{M}\iff\bm{s}=0$. In
other words, coordinates $r_{i}$ represent the _soft_ degrees of freedom and
$s_{i}$ represent the _hard_ directions. Since the new coordinates are
orthonormal we know that the determinant of the transformation
$(q_{i})\mapsto(r_{i},s_{j})$ equals 1. Additionally, $I_{k}$ is finite by
compactness of the support of $\phi$ and continuity of $p_{k}$, so by Fubini’s
theorem we can replace the volume integral by an iterated one. Then, $I_{k}$
is given by

$I_{k}=\iint
p_{k}(\bm{r},\bm{s})\phi(\bm{r},\bm{s})\mathrm{d}\bm{s}\mathrm{d}\bm{r}$ (5)

We proceed by expanding $p_{k}$ in a Taylor series with respect to
$\bm{s}=(s_{1},s_{2},\dots)$ to second order as
$W(\bm{r},\bm{s})=\bm{s}^{T}\bm{\mathsf{H}}(\bm{r})\bm{s}$, where
$\bm{\mathsf{H}}(\bm{r})$ is the Hessian of $W$, evaluated at $\bm{s}=0$. Note
that the linear term vanishes, since at $\bm{s}=0$ the potential $W$ attains a
minimum. We can thus write

$I_{k}=\iint N(k)\exp\left\\{-\beta U(\bm{r},\bm{s})-\beta
k^{2}\left[\bm{s}^{T}\bm{\mathsf{H}}(\bm{r})\bm{s}+o(||\bm{s}||^{2})\right]\right\\}\phi(\bm{r},\bm{s})\mathrm{d}\bm{s}\mathrm{d}\bm{r}$
(6)

We now make a substitution $k\bm{s}=\bm{t}$ to arrive at

$I_{k}=\iint\frac{N(k)}{k}\exp\left\\{-\beta
U(\bm{r},k^{-1}\bm{t})-\beta\left[\bm{t}^{T}\bm{\mathsf{H}}(\bm{r})\bm{t}+||\bm{t}||^{2}\frac{o(k^{-2}||\bm{t}||^{2})}{k^{-2}||\bm{t}||^{2}}\right]\right\\}\phi(\bm{r},k^{-1}\bm{t})\mathrm{d}\bm{t}\mathrm{d}\bm{r}$
(7)

Whenever $W$ increases sufficiently fast with the distance from the
constraining manifold, we may take a limit inside the integral by the
dominated convergence theorem to arrive at the limiting value
$I=\lim_{k\to\infty}I_{k}$ given by

$I=\iint\left(\lim_{k\to\infty}\frac{N(k)}{k}\right)\exp\left[-\beta
U(\bm{r},\bm{0})-\beta\,\bm{t}^{T}\bm{\mathsf{H}}(\bm{r})\bm{t}\right]\phi(\bm{r},\bm{0})\mathrm{d}\bm{t}\mathrm{d}\bm{r}$
(8)

and, since $\bm{\mathsf{H}}$ is full rank, perform the integral over $\bm{t}$
to arrive at

$I=L\int|\bm{\mathsf{H}}(\bm{r})|^{-1/2}\exp\left[-\beta
U(\bm{r},\bm{0})\right]\phi(\bm{r},\bm{0})\mathrm{d}\bm{r}$ (9)

with the constant

$L=\lim_{k\to\infty}\frac{N(k)\sqrt{(2\pi)^{N-M}}}{k},$ (10)

being independent of the soft coordinates.

In a particular physical setting, finding the orthonormal coordinates
$(\bm{r},\bm{s})$ may be prohibitively difficult, but we can relax this strict
requirement by considering another parametrisation $\bm{w}$ such that
$\bm{q}=\bm{\zeta^{*}}(\bm{w})$, and that still separates _soft_ from _hard_
degrees of freedom by

$\bm{q}\in\mathcal{M}\iff\forall i>M:w_{i}=\xi_{i},$ (11)

for a set of constants $\xi_{i}$. This parametrisation is not necessarily
orthonormal. For convenience, we define a map
$\bm{\zeta}:\mathbb{R}^{M}\to\mathbb{R}^{N}$ such that
$\bm{\zeta}(w_{1},w_{2},\dots,w_{M})=\bm{\zeta^{*}}(w_{1},w_{2},\dots,w_{M},\xi_{M+1},\xi_{M+2},\dots,\xi_{N})$
which takes values of the _soft_ coordinates and returns points on
$\mathcal{M}$. We now change variables under the integral from $\bm{r}$ to
$\bm{w}$ to arrive at

$I=L\int\frac{\left|\bm{\mathsf{J}}^{T}\bm{\mathsf{J}}\right|^{1/2}}{\left|\bm{\mathsf{H}}(\bm{\zeta})\right|^{1/2}}\exp\left[-\beta
U(\bm{\zeta})\right]\phi(\bm{\zeta})\mathrm{d}w_{1}\dots\mathrm{d}w_{M},$ (12)

where $\bm{\mathsf{J}}$ is the Jacobian of $\bm{\zeta}(\bm{w})$. Note that
$\bm{\zeta}:\mathbb{R}^{M}\to\mathbb{R}^{N}$, so the Jacobian is not a square
matrix. We can express the above using the Dirac-$\delta$ distributions.
Taking

$\delta(\bm{w}^{\prime}-\bm{\xi}^{\prime})=\prod_{i=M+1}^{N}\delta(w_{i}-\xi_{i}),$
(13)

the limiting integral is recast as

$I=L\int\frac{\left|\bm{\mathsf{J}}^{T}\bm{\mathsf{J}}\right|^{1/2}}{|\bm{\mathsf{H}}|^{1/2}}\exp\left(-\beta
U\right)\phi(\bm{\zeta^{*}}(\bm{w}))\delta(\bm{w}^{\prime}-\bm{\xi}^{\prime})\,\mathrm{d}\bm{w}$
(14)

for an appropriate constant $L$, and since $\phi$ is arbitrary we arrive at
$p_{\infty}=\lim_{k\to\infty}p_{k}$ such that

$\mathrm{d}p_{\infty}=L\frac{|\bm{\mathsf{J}}^{T}\bm{\mathsf{J}}|^{1/2}}{\left|\bm{\mathsf{H}}\right|^{1/2}}\exp\left[-\beta
U(\bm{\zeta^{*}}(\bm{w}))\right]\delta(\bm{w}^{\prime}-\bm{\xi}^{\prime})\
\mathrm{d}\bm{w},$ (15)

in the weak sense. We have thus arrived at the central result of this
manuscript.

The limiting distribution comprises two important factors. First, the metric
term $|\bm{\mathsf{J}}^{T}\bm{\mathsf{J}}|$ describing the projection of the
surface element of $\mathcal{M}$ onto $\mathbb{R}^{M}$ in the parametrisation
$\bm{\zeta}(\bm{w})$; This term was correctly computed by Fixman [5]. Second,
the _shape_ term $|\bm{\mathsf{H}}|$, also called the zero-point energy in the
quantum-mechanical setting, which was often missing in the derivations.

###  III.1 Computation of the shape term

A careful choice of parametrisation allows for computation of the metric term
with relative ease. However, the Hessian term $|\bm{\mathsf{H}}|$ may be more
difficult to compute. In this Section, we present a feasible approach to its
computation.

Let $\bm{\mathsf{B}}$ be the full Hessian of $W$ at a point $\bm{q}$,

$\bm{\mathsf{B}}_{ij}=\frac{\partial^{2}}{\partial q_{i}\partial
q_{j}}W(\bm{q}).$ (16)

Since the eigenvectors of $\bm{\mathsf{B}}$ are orthogonal and zero
eigenvectors lie inside the tangent space $T\mathcal{M}$, we can compute
$|\bm{\mathsf{H}}|$ from the product of the non-zero eigenvalues of
$\bm{\mathsf{B}}$. We can use the knowledge of the tangent space of
$\mathcal{M}$ to help us compute the eigenvalues of $\bm{\mathsf{B}}$. If our
confining function is of the form $W(\bm{q})=\sum_{i}P_{i}^{2}$, we can find
vectors that lie in the normal space $\mathcal{M}^{\perp}$ [15, proposition
8.15] by computing

$\bm{\mathsf{A}}_{ji}=\frac{\partial}{\partial q_{j}}P_{i}$ (17)

The vectors $\\{\bm{\mathsf{A}}_{j1},\ldots,\bm{\mathsf{A}}_{j,N-M}\\}$ are
not pairwise orthogonal but are orthogonal to $\mathcal{M}$, and thus can form
a basis of the normal space at $\bm{q}$. We can write the eigenvalue problem
using an arbitrary vector $b_{i}$ as

$\lambda\bm{\mathsf{A}}_{ki}b_{i}=\bm{\mathsf{B}}_{kj}\bm{\mathsf{A}}_{ji}b_{i}.$
(18)

This is a system of $N$ equations with $N-M$ unknowns. We can eliminate
redundant equations by contracting each side with $\bm{\mathsf{A}}_{ij}$ to
finally get

$\lambda(\bm{\mathsf{A}}^{T}\bm{\mathsf{A}})\bm{b}=(\bm{\mathsf{A}}^{T}\bm{\mathsf{B}}\bm{\mathsf{A}})\bm{b}$
(19)

and thus the product of eigenvalues in this problem is simply

$|\bm{\mathsf{H}}|=\frac{|\bm{\mathsf{A}}^{T}\bm{\mathsf{B}}\bm{\mathsf{A}}|}{|\bm{\mathsf{A}}^{T}\bm{\mathsf{A}}|}.$
(20)

We see that there are two terms in this expression. The denominator
$|\bm{\mathsf{A}}^{T}\bm{\mathsf{A}}|$ measures the angles between the
gradients of the constraining functions, while the numerator
$|\bm{\mathsf{A}}^{T}\bm{\mathsf{B}}\bm{\mathsf{A}}|$ measures the shape and
strength of the confining field in these directions.

##  IV Calculation for the trimer problem

Figure 2: (a) Parametrisation used for the trimer uses the position of the
central particle ($b$), bond angle $\psi$ and bond lengths $l_{a}$ and
$l_{b}$, with identical equilibrium spring lengths $l_{0}$. Rigid-body
rotation around the particle $b$ is omitted for clarity. (b) Motion along the
symmetric (top) and antisymmetric eigenvectors of the hessian
$\bm{\mathsf{H}}$ of the harmonic confinement function. Note that the central
bead must move for the overall motion to be orthogonal to the constraining
manifold $\mathcal{M}$. (c) Equilibrium probability density function (PDF) for
the bond angle $\psi$ from the Langevin numerical simulation. Solid line is
the theoretical prediction of Eq. (32).

In the aforementioned case of a trimer, depicted in Fig. 2(a), we have 3 beads
located at $\bm{r}_{a},\bm{r}_{b}$ and $\bm{r}_{c}$ with bond extensions
$P_{a},P_{c}$ given by

$\displaystyle P_{a}$ $\displaystyle=$
$\displaystyle|\bm{r}_{a}-\bm{r}_{b}|-l_{0}$ (21) $\displaystyle P_{c}$
$\displaystyle=$ $\displaystyle|\bm{r}_{c}-\bm{r}_{b}|-l_{0},$ (22)

$l_{0}$ being the equilibrium bond length, and the confining potential given
by

$W(\bm{q})=P_{a}^{2}+P_{c}^{2}.$ (23)

We now pick the parametrisation $\bm{q}=\bm{\zeta^{*}}(\bm{w})$ with
$\bm{w}=[x,y,z,\alpha,\beta,\gamma,\psi,l_{a},l_{c}]^{T}$. Cartesian
coordinates $(x,y,z)$ describe the position of the central bead ($b$),
$(\alpha,\beta,\gamma)$ are the Euler angles that encode the direction of bead
$a$ seen from $b$, and $\psi$ is the bond angle between the edges $ab$ and
$bc$. The positions of beads are thus given by

$\displaystyle\bm{r}_{a}$ $\displaystyle=$
$\displaystyle[x,y,z]^{T}+l_{a}\bm{\mathsf{E}}_{\alpha\beta\gamma}[1,0,0]^{T},$
(24) $\displaystyle\bm{r}_{b}$ $\displaystyle=$ $\displaystyle[x,y,z]^{T},$
(25) $\displaystyle\bm{r}_{c}$ $\displaystyle=$
$\displaystyle[x,y,z]^{T}+l_{c}\bm{\mathsf{E}}_{\alpha\beta\gamma}[\cos\psi,\sin\psi,0]^{T}.$
(26)

We defined the rotation matrix $\bm{\mathsf{E}}_{\alpha\beta\gamma}$ in terms
of the Euler angles in the Appendix. In addition, we note that
$P_{a}=l_{a}-l_{0}$ and $P_{c}=l_{c}-l_{0}$. We compute the metric term as

$|\bm{\mathsf{J}}^{T}\bm{\mathsf{J}}|=\frac{1}{2}l_{0}^{8}\sin^{2}\beta\sin^{2}\psi\left(7-\cos
2\psi\right),$ (27)

from the matrix $\bm{\mathsf{J}}$ listed in the Appendix. Next, we compute the
shape term using

$\displaystyle\bm{\mathsf{A}}^{T}\bm{\mathsf{A}}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}2&\cos\psi\\\ \cos\psi&2\end{bmatrix},$ (28)
$\displaystyle\bm{\mathsf{A}}^{T}\bm{\mathsf{B}}\bm{\mathsf{A}}$
$\displaystyle=$ $\displaystyle\begin{bmatrix}9+\cos 2\psi&8\cos\psi\\\
8\cos\psi&9+\cos 2\psi\end{bmatrix},$ (29)

thus finding

$|\bm{\mathsf{H}}|=\frac{|\bm{\mathsf{A}}^{T}\bm{\mathsf{B}}\bm{\mathsf{A}}|}{|\bm{\mathsf{A}}^{T}\bm{\mathsf{A}}|}=2(7-\cos
2\psi).$ (30)

We therefore have

$\frac{|\bm{\mathsf{J}}^{T}\bm{\mathsf{J}}|^{1/2}}{|\bm{\mathsf{H}}|^{1/2}}=\frac{l_{0}^{4}}{2}\sin\beta\sin\psi,$
(31)

and for the final distribution, we find

$\mathrm{d}p_{\infty}=N\sin\beta\sin\psi\,\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}\alpha\,\mathrm{d}\beta\,\mathrm{d}\gamma\,\mathrm{d}\psi,$
(32)

where $N$ is a normalisation constant. Therefore, the marginal density of the
bond angle is proportional to $\sin\psi$ as expected. However, this is simply
not the case of a uniform distribution––that distribution corresponds to
$|\bm{\mathsf{H}}|=\text{const}$. This can be examined in greater detail by
looking at the eigenvectors of $\bm{\mathsf{H}}$, visualised in Fig. 2(b).
There are two eigenvectors which correspond to symmetric and antisymmetric
motion of the terminal beads. Counterintuitively, the vectors in
$\mathcal{M}^{\perp}$ change the orientation of at least one bond when the
bond length changes.

The simplification lies in the fortuitous cancellation of the $(7-\cos 2\psi)$
term, rather than its absence. This is best highlighted by an example where no
such cancellation occurs, which we present in the next Section.

We additionally corroborate this result with a Brownian Dynamics simulation in
Python using the package pychastic [16]. Choosing $k_{B}T$ as the energy
scale, $l_{0}$ as the distance scale and $l_{0}^{2}/D$, as the timescale, with
$D$ being the diffusion coefficient, the corresponding stochastic BD equation
takes the form

$\mathrm{d}\bm{r}_{i}=-\bm{\nabla}(U+k^{2}W)\mathrm{d}t+\sqrt{2}\,\mathrm{d}\bm{\mathcal{B}}_{i},$
(33)

with $\bm{\mathcal{B}}_{i}$ denoting the standard Wiener process. We have
performed $n_{\text{traj}}=4000$ simulations with $k=35$ (in units of
$k_{B}T/l_{0}$) for a time up to $t_{max}=10$ with time step $\delta
t=10^{-5}$. The simulation results in Fig. 2(c) complement the prediction of
Eq. (32). We use the same code for a less studied example of a tetramer in the
following.

##  V Cyclic tetramer

Figure 3: (a) The shape of a cyclic tetramer is specified with the bead $b$
taken as reference point. The directions of sections $ab$ and $bc$ span a
rhombus, which is used to define the position of bead $d$. (b) Equilibrium
probability density function (PDF) for the bond angle $\psi$ from the
numerical simulation. Solid line is the theoretical prediction of Eq. (47).

A more complex, and yet fully tractable analytically, is the case of a cyclic
tetramer, a molecule of 4 beads joined into a quadrilateral. For the sake of
brevity, we discuss the dynamics in two dimensions. The calculation in 3D is
completely analogous but involves much longer intermediate expressions due to
the additional rotational degrees of freedom and a dihedral angle. We pick the
parametrisation where the 4 beads are located at
$\bm{r}_{a},\bm{r}_{b},\bm{r}_{c}$ and $\bm{r}_{d}$, with their positions
given by

$\displaystyle\bm{r}_{a}$ $\displaystyle=$
$\displaystyle[x,y]^{T}+l_{a}\bm{\mathsf{E}}_{\alpha}[1,0]^{T},$ (34)
$\displaystyle\bm{r}_{b}$ $\displaystyle=$ $\displaystyle[x,y]^{T},$ (35)
$\displaystyle\bm{r}_{c}$ $\displaystyle=$
$\displaystyle[x,y]^{T}+l_{c}\bm{\mathsf{E}}_{\alpha}[\cos\psi,\sin\psi]^{T},$
(36) $\displaystyle\bm{r}_{d}$ $\displaystyle=$
$\displaystyle[x,y]^{T}+\bm{\mathsf{E}}_{\alpha}\left(\left(1+\delta\right)\bm{v}_{\delta}+\epsilon\,\bm{v}_{\epsilon}\right),$
(37)

where $\bm{\mathsf{E}}_{\alpha}$ is a rotation matrix defined in the Appendix,
and the vectors $\bm{v}_{\delta,\epsilon}$, written as

$\displaystyle\bm{v}_{\delta}$ $\displaystyle=$
$\displaystyle\frac{[l_{a}+l_{c}\cos\psi,l_{c}\sin\psi]^{T}}{2},$ (39)
$\displaystyle\bm{v}_{\epsilon}$ $\displaystyle=$
$\displaystyle\frac{[l_{a}-l_{c}\cos\psi,-\sin\psi]^{T}}{2},$ (40)

point along the diagonals of the rhombus spanned by the edges $ab$ and $bc$.
Finally, $\epsilon$ and $\delta$ encode the position of the bead located at
$\bm{r}_{d}$. The four constraining surfaces are given by

$\displaystyle P_{a}$ $\displaystyle=$
$\displaystyle|\bm{r}_{a}-\bm{r}_{b}|-l_{0}$ (41) $\displaystyle P_{b}$
$\displaystyle=$ $\displaystyle|\bm{r}_{b}-\bm{r}_{c}|-l_{0}$ (42)
$\displaystyle P_{c}$ $\displaystyle=$
$\displaystyle|\bm{r}_{c}-\bm{r}_{d}|-l_{0}$ (43) $\displaystyle P_{d}$
$\displaystyle=$ $\displaystyle|\bm{r}_{d}-\bm{r}_{a}|-l_{0}$ (44)

and $\mathcal{M}$ is given by $l_{a}=l_{c}=l_{0},\delta=1,\epsilon=0$.

The metric term is simply constant with
$|\bm{\mathsf{J}}^{T}\bm{\mathsf{J}}|=16l_{0}^{4}$, while the shape term is
non-trivial. Denoting $\cos\psi=c_{\psi}$ for brevity, the shape factors are
given by

$\displaystyle\bm{\mathsf{A}}^{T}\bm{\mathsf{A}}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}2&c_{\psi}&0&-c_{\psi}\\\
c_{\psi}&2&-c_{\psi}&0\\\ 0&-c_{\psi}&2&c_{\psi}\\\ -c_{\psi}&0&c_{\psi}&2\\\
\end{bmatrix},$ (45)
$\displaystyle\bm{\mathsf{A}}^{T}\bm{\mathsf{B}}\bm{\mathsf{A}}$
$\displaystyle=$ $\displaystyle
4\begin{bmatrix}c_{\psi}^{2}+2&2c_{\psi}&-c_{\psi}^{2}&-2c_{\psi}\\\
2c_{\psi}&c_{\psi}^{2}+2&-2c_{\psi}&-c_{\psi}^{2}\\\
-c_{\psi}^{2}&-2c_{\psi}&c_{\psi}^{2}+2&2c_{\psi}\\\
-2c_{\psi}&-c_{\psi}^{2}&2c_{\psi}&c_{\psi}^{2}+2\\\ \end{bmatrix},$ (46)

yielding $|\bm{\mathsf{H}}|=256\sin^{2}\psi$. As a result, we get that

$\mathrm{d}p_{\infty}\propto\frac{1}{\sin\psi}\mathrm{d}\psi,$ (47)

which is not non-normalizable near 0 and $\pi$. Physically, as we approach the
stiff spring limit, the tetramer spends more and more time in the folded
state, and unless some repulsive potential is added, the molecule in the stiff
spring limit looks predominantly like a trimer, alternating between beads
$(a,c)$ and $(b,d)$ coinciding in space.

This surprising result is also seen in the numerical simulation. In this case,
a strong repulsive potential was added to prevent bead overlaps for distances
less than $0.05l_{0}$. The final bond angle distribution shown in Figure 3(b)
coincides well with $0.15/\sin\psi$.

##  VI Conclusions

In this work, we have provided a procedure by which one can find the
equilibrium probability probability distribution function of the
configurations of a bead-spring model of a polymer in the overdamped regime in
the limit of infinite stiffness of the bonds. In particular, we have
demonstrated that the shape of the confining spring potential persists in the
limiting distribution. By rephrasing the problem as a mathematically rigorous
limit, we have shown that the classical expression given by Fixman [5], later
reproduced in numerous works and books, for example in Ref. [3], is missing a
critical term describing the shape of the confining potential,
$|\bm{\mathsf{H}}|$.

This shows that the reasoning presented by Frenkel & Smit [3] cannot be
universally applied. First, the dynamics and the distribution of a polymer
connected with (stiff or soft) springs is independent of the mass of the
monomers in the overdamped regime, and thus the limiting distribution cannot
depend on the mass of the monomers. Second, putting aside the question of
existence of the rigid rod distribution (or the lack thereof), the supposition
that the stiff spring limit should always be the same, regardless of the
details of the confining potential is clearly not true, and harmonic
potentials will lead to a different outcome than springs that realise uniform
confinement around the constraining manifold, as previously remarked, e.g., by
van Kampen & Lodder [12], although without a general mathematical formulation.
The difference can be arbitrarily large, as shown in the tetramer example.
Frenkel & Smit remark in their book that for the bond length constraints of
the type most often used in Molecular Dynamics simulations, the effect of hard
constraints on the distribution functions seems to be relatively small [3,
Chapter 15.1.1]. We believe this statement to be misleading. In the discussed
case of the trimer molecule, the effect might indeed be small, but we have
shown that for the tetramer system, which is still a relatively simple setup,
the effect is pronounced. To provide a practical insight, we propose a
straightforward technique for the analysis of the limiting distributions once
the confining potential is specified in detail.

###### Acknowledgements.

The Authors would like to thank Bartłomiej Lewandowski for his insightful
comments regarding convergence of distributions. The work of ML and RW was
supported by the National Science Centre of Poland (FundRef DOI:
http://dx.doi.org/10.13039/501100004281) grant Sonata to ML no.
2018/31/D/ST3/02408.

##  Appendix A Auxiliary matrices used in calculations

Using the notation $s_{\alpha}=\sin\alpha$, $c_{\alpha}=\cos\alpha$ we can
compactly express the matrices used. The 2D rotation matrix is given by

$\bm{\mathsf{E}}_{\alpha}=\begin{bmatrix}c_{\alpha}&-s_{\alpha}\\\
s_{\alpha}&c_{\alpha}\end{bmatrix}.$ (48)

The 3D rotation matrix is given by

$\bm{\mathsf{E}}_{\alpha\beta\gamma}=\begin{bmatrix}c_{\alpha}c_{\beta}c_{\gamma}-s_{\alpha}s_{\gamma}&\
-c_{\alpha}c_{\beta}s_{\gamma}-c_{\gamma}s_{\alpha}&\ c_{\alpha}s_{\beta}\\\
c_{\beta}c_{\gamma}s_{\alpha}+c_{\alpha}s_{\gamma}&\
c_{\alpha}c_{\gamma}-c_{\beta}s_{\alpha}s_{\gamma}&\ s_{\alpha}s_{\beta}\\\
-c_{\gamma}s_{\beta}&s_{\beta}s_{\gamma}&c_{\beta}\\\ \end{bmatrix}.$ (49)

For our paremetrisation of the trimer, the Jacobian $\bm{\mathsf{J}}$ can be
expressed as a block matrix with $3\times 4$ blocks $\bm{\mathsf{Q}}_{i}$

$\bm{\mathsf{J}}=\begin{bmatrix}\bm{\mathsf{1}}&\bm{\mathsf{Q}}_{1}\\\
\bm{\mathsf{1}}&\bm{\mathsf{Q}}_{2}\\\ \bm{\mathsf{1}}&\bm{\mathsf{Q}}_{3}\\\
\end{bmatrix}.$ (50)

Writing $\gamma+\psi=\theta$ for brevity, the blocks $\bm{\mathsf{Q}}_{i}$ are
given by

$\bm{\mathsf{Q}}_{1}^{T}=l_{0}\begin{bmatrix}-s_{\alpha}c_{\beta}c_{\gamma}-c_{\alpha}s_{\gamma}&c_{\alpha}c_{\beta}c_{\gamma}-s_{\alpha}s_{\gamma}&0\\\
-c_{\alpha}s_{\beta}c_{\gamma}&-s_{\alpha}s_{\beta}c_{\gamma}&-c_{\beta}c_{\gamma}\\\
-s_{\alpha}c_{\gamma}-c_{\alpha}c_{\beta}s_{\gamma}&c_{\alpha}c_{\gamma}-s_{\alpha}c_{\beta}s_{\gamma}&s_{\beta}s_{\gamma}\\\
0&0&0\end{bmatrix}$ (51)

$\bm{\mathsf{Q}}_{2}^{T}=\begin{bmatrix}0&0&0\\\ 0&0&0\\\ 0&0&0\\\ 0&0&0\\\
\end{bmatrix}$ (52)

.

$\bm{\mathsf{Q}}_{3}^{T}=l_{0}\begin{bmatrix}-s_{\alpha}c_{\beta}c_{\theta}-c_{\alpha}s_{\theta}&c_{\alpha}c_{\beta}c_{\theta}-s_{\alpha}s_{\theta}&0\\\
-c_{\alpha}s_{\beta}c_{\theta}&-s_{\alpha}s_{\beta}c_{\theta}&-c_{\beta}c_{\theta}\\\
-s_{\alpha}c_{\theta}-c_{\alpha}c_{\beta}s_{\theta}&c_{\alpha}c_{\theta}-s_{\alpha}c_{\beta}s_{\theta}&s_{\beta}s_{\theta}\\\
-s_{\alpha}c_{\theta}-c_{\alpha}c_{\beta}s_{\theta}&c_{\alpha}c_{\theta}-s_{\alpha}c_{\beta}s_{\theta}&s_{\beta}s_{\theta}\\\
\end{bmatrix}$ (53)

## References

  * Kramers [1946] H. A. Kramers, The behavior of macromolecules in inhomogeneous flow, The Journal of Chemical Physics 14, 415 (1946).
  * Flory [1953] P. Flory,  _Principles of Polymer Chemistry_ , Baker lectures 1948 (Cornell University Press, 1953).
  * Frenkel and Smit [2023] D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ (Elsevier, 2023).
  * Ryckaert  _et al._ [1977] J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen, Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes, Journal of Computational Physics 23, 327 (1977).
  * Fixman [1974] M. Fixman, Classical statistical mechanics of constraints: a theorem and application to polymers., Proceedings of the National Academy of Sciences of the United States of America 71, 3050 (1974).
  * Gō and Scheraga [1969] N. Gō and H. A. Scheraga, Analysis of the Contribution of Internal Vibrations to the Statistical Weights of Equilibrium Conformations of Macromolecules, The Journal of Chemical Physics 51, 4751 (1969).
  * Pear and Weiner [1979] M. R. Pear and J. H. Weiner, Brownian dynamics study of a polymer chain of linked rigid bodies, The Journal of Chemical Physics 71, 212 (1979).
  * Hinch [1994] E. Hinch, Brownian motion with stiff bonds and rigid constraints, Journal of Fluid Mechanics 271, 219 (1994).
  * den Otter and Briels [1998] W. K. den Otter and W. J. Briels, The calculation of free-energy differences by constrained molecular-dynamics simulations, The Journal of Chemical Physics 109, 4139 (1998).
  * Liu [1989] T. W. Liu, Flexible polymer chain dynamics and rheological properties in steady flows, The Journal of Chemical Physics 90, 5826 (1989).
  * Helfand [1979] E. Helfand, Flexible vs rigid constraints in statistical mechanics, The Journal of Chemical Physics 71, 5000 (1979).
  * van Kampen and Lodder [1984] N. van Kampen and J. Lodder, Constraints, American Journal of Physics 52, 419 (1984).
  * Vitalis and Pappu [2009] A. Vitalis and R. V. Pappu, Methods for monte carlo simulations of biomacromolecules, Annual reports in computational chemistry 5, 49 (2009).
  * Allen [2004] M. P. Allen, Introduction to molecular dynamics simulation, in  _Computational soft matter: from synthetic polymers to proteins_ , Vol. 23, edited by N. Attig, K. Binder, H. Grubmüller, and K. Kremer (NIC Series Jülich, Germany, 2004) pp. 1–28.
  * Lee [2003] J. M. Lee,  _Introduction to Smooth Manifolds_ (Springer New York, 2003).
  * Waszkiewicz  _et al._ [2023] R. Waszkiewicz, M. Bartczak, K. Kolasa, and M. Lisicki, Pychastic: Precise brownian dynamics using Taylor-Itō integrators in Python, SciPost Physics Codebases , 011 (2023).


New Mexico Tech (September 26, 2023)

MOND via Matrix General Relativity

Ivan Avramidi and Roberto Niardi

Department of Mathematics

New Mexico Institute of Mining and Technology

Socorro, NM 87801, USA

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>

MOND theory has arisen as a promising alternative to dark matter in explaining
the collection of discrepancies that constitute the so-called missing mass
problem. The MOND paradigm is briefly reviewed. It is shown that MOND theory
can be incorporated in the framework of Matrix General Relativity. In
particular, we demonstrate that Matrix General Relativity contains MOND as a
particular case.

##  1 Introduction

The theory of gravity, which Sir Isaac Newton famously introduced in the 17th
century and Albert Einstein’s General Relativity later developed in the early
20th century, has served as a cornerstone of our knowledge of the fundamental
forces that shape the world for many years. Numerous phenomena, ranging from
the motion of celestial bodies to the behavior of matter on Earth, have been
successfully described by these theories. However, a rising body of evidence
indicates that adjustments to our existing theories of gravity may be required
as our scientific knowledge has grown and our observational abilities have
increased.

Here we list some of the problems that constitute the so-called missing mass
problem:

  * •

In 1932, Oort [21] found that the total number of visible stars in the solar
neighborhood was insufficient to fully explain the vertical motions of stars
in the Milky Way disk. The observed star vertical oscillations could not be
explained entirely by the luminous matter. This experimental fact became known
as Oort discrepancy.

  * •

Concomitantly, Zwicky [28] noted that in order for galaxies within galaxy
clusters to form a bound system over a significant portion of cosmic time,
their velocity dispersion should have been way lower than measured.

  * •

Galactic disk stability was one of the earliest recent signs of the necessity
for dark matter. Cold, self-gravitating disks were shown to be extremely
unstable in early simulations [22]. A potential well was necessary for the
stability of the system. This could be provided by a halo of dark matter, as
Ostriker & Peebles [22] proposed in 1973. This would prevent the growth of
these instabilities, which would have destroyed the whole system on short a
time scale, compared to a Hubble time.

  * •

Bosma [8] and Rubin [23] established that the rotation curves of spiral
galaxies are approximately flat: according to Newtonian gravity, when the
majority of its mass is enclosed, a system should have a rotation curve that
diminishes in a Keplerian fashion, similar to the solar system: $v\sim
r^{-1/2}$. Nonetheless, rotation curves for spiral galaxies tended to remain
roughly flat with increasing radius instead: $v\sim{\rm constant}$. A second
time, the proposal by Ostriker and Peebles of a dark matter halo would have
been the solution to this mystery.

  * •

Further anomalies come from the observation of dwarf spheroidal Milky Way
satellite galaxies [25, 26]. By analyzing the motion of individual stars that
constitute these galaxies, the masses of these systems can be inferred, and it
is way beyond the mass visible in luminous stars.

  * •

Another piece of evidence comes from larger scales. It is the so called
“timing argument” in the Local Group [13]. At this time, the Milky Way and
Andromeda are approaching one another. But it is widely believed that the
distance between the material constituting these galaxies was initially
increasing, as a consequence of the Big Bang expansion. In order for the two
to have overcome the initial expansion, a gravitating mass greater than their
average mass has to be located between them. The mass-to-light ratio which is
necessary to explain the approaching velocity currently measured at the
present time is way beyond the mass-to-light ratio of the stars themselves,
which have the same mass-to-light ratio of our Sun, as order of magnitude [6].

  * •

Cluster of galaxies exhibit mass discrepancies, too. In particular, this is
observed in three different ways, which give similar estimates for the amount
of missing mass. The first way is through the measurement of the redshifts of
individual cluster members; from these, velocity dispersions can be inferred
and consequently mass-to-light ratio can be estimated [5]. The second way is
through mapping the temperature and emission of the diffuse intra-cluster gas,
which emits in the X-rays; an estimate of the mass is obtained through the
equation of hydrostatic equilibrium [12]. The third way is through the lensing
produced by these cluster of galaxies [15]. In all cases, the amount of mass
needed to explain these phenomena is way beyond the amount actually observed.

  * •

Further evidence comes from cosmology: by comparing the rate of growth of
large scale structures and the estimate of the baryonic mass density
$\Omega_{b}$ coming from big bang nucleosynthesis, we deduce that
$\Omega_{gravitating}>\Omega_{barionic}$ [10, 27, 9]. Therefore, whatever dark
matter is, it has to be non-baryonic.

  * •

Another cosmological constraint comes from the CMB: from the temperature
fluctuations at the time of decoupling between photons and baryons, we infer
that the Universe was very homogeneous at the beginning, approximately to one
in $10^{5}$. Nonetheless, the Universe today presents remarkable density
contrasts between various astrophysical objects, and the emptiness of
intergalactic space. Gravity is the only interaction we know that could be
responsible for the formation of such structures. But the baryon density from
BBN can only produce a growth factor of $\sim 10^{2}$ in a Hubble time [24],
while the observed one is roughly $10^{5}$. In this context, dark matter would
be the ideal candidate to enhance the growth rate. Moreover, this dark matter
must not interact with electromagnetic radiation, because otherwise its
contribution on the CMB could not be distinguished from the one produced by
baryons.

All the experimental fact listed above (and many others) have prompted
scientists to think that those anomalies were caused by not taking into
account a new kind of matter with peculiar features: dark matter. Nonetheless,
so far dark matter has always escaped detection.

Another path for the explanation of these anomalies is admitting that our
understanding of gravity is still incomplete, and a modification of General
Relativity is necessary. The Modified Newtonian Dynamics (MOND) theory has
emerged as a promising alternative framework that modifies the laws of gravity
in order to explain the observed dynamics without invoking the presence of
dark matter.

This paper explores the incorporation of MOND theory into the framework of
Matrix General Relativity [1, 2, 3]. By introducing noncommutativity into the
geometry, it becomes possible to encode additional physical information.

We begin by reviewing the main concepts of MOND theory. We then explore the
possible ways in which the noncommutative framework can be used to reformulate
MOND, leading to a modified gravitational theory that incorporates
noncommutativity as a fundamental aspect of the metric tensor. We examine the
consequences of this modification for the behavior of gravity.

##  2 MOND Theory

To fix notation, we work in the three-dimensional Euclidean space
$\mathbb{R}^{3}$ with the standard Euclidean metric $\delta_{ij}$ and
Cartesian coordinates $x^{i}$, where the Latin indices are running over
$i=1,2,3$. For the relativistic theory we use the units with $c=1$ and the
coordinates $(x^{\mu})=(t,x^{i})$, where the Greek indices run over
$\mu=0,1,2,3$ so that $x^{0}=t$; the Minkowski metric is
$\eta_{\mu\nu}=\mathrm{diag\,}(-1,1,1,1)$. The volume element is denoted by
$dx=dx^{1}\,dx^{2}\,dx^{3}$. Also, summation over repeated indices is assumed
if not specified otherwise. The partial derivatives with respect to $x^{\mu}$
are denoted by $\partial_{\mu}$ and the Laplacian by
$\Delta=\nabla_{i}\nabla^{i}$; the norm of a vector $A_{i}$ is denoted by
$|A|=\sqrt{A_{i}A^{i}}$.

The basic equation of the Newtonian gravity is the Poisson equation

$\Delta\varphi=4\pi G\rho,$ (2.1)

where $\varphi$ is the gravitational potential, $\rho$ is the mass density and
$G$ is the gravitational constant, which is obtained from the action

$S=S_{\rm Newton}+S_{\rm matter},$ (2.2)

with

$\displaystyle S_{\rm Newton}$ $\displaystyle=$ $\displaystyle-\frac{1}{8\pi
G}\int dt\int\limits_{\Omega}dx\;|\nabla\varphi|^{2},$ (2.3) $\displaystyle
S_{\rm matter}$ $\displaystyle=$ $\displaystyle-\int
dt\int\limits_{\Omega}dx\;\rho\varphi,$ (2.4)

where $\Omega$ is a region in space.

The idea of the Modified Newtonian Dynamics (MOND) theory of gravity suggested
by Milgrom [16] is to modify the Newtonian dynamics in the weak-acceleration
regime, $a<a_{0}$, with a universal acceleration constant

$a_{0}\approx 10^{-8}\;\mathrm{cm}\;\mathrm{sec}^{-2}.$ (2.5)

This amounts to modifying the Poisson equation by [7]

$\nabla_{i}\left[\mu\left(\frac{|\nabla\varphi|}{a_{0}}\right)\nabla_{i}\varphi\right]=4\pi
G\rho,$ (2.6)

where $\mu(x)$ is a function that is supposed to have the properties

$\mu(x)\sim\left\\{\begin{array}[]{ll}1&\mbox{ for }x>>1,\\\\[5.0pt] x&\mbox{
for }x<<1.\end{array}\right.$ (2.7)

This equation is obtained from the following action

$S_{\rm MOND}=-\frac{a_{0}^{2}}{8\pi G}\int
dt\int\limits_{\Omega}dx\;F\left(\frac{|\nabla\varphi|^{2}}{a^{2}_{0}}\right),$
(2.8)

where the function $F(z)$ satisfies

$F(z)\sim\left\\{\begin{array}[]{ll}z&\mbox{ for }z>>1,\\\\[12.0pt]
\frac{2}{3}z^{3/2}&\mbox{ for }z<<1.\end{array}\right.$ (2.9)

The function $\mu(x)$ is then determined by $\mu(x)=F^{\prime}(x^{2})$.

An alternative modification called QUMOND [19] requires the introduction of
two potential functions

$S_{\rm QUMOND}=\frac{1}{8\pi G}\int
dt\int\limits_{\Omega}dx\;\left\\{-2\nabla_{i}\varphi\nabla_{i}\psi+a_{0}^{2}Q\left(\frac{|\nabla\psi|^{2}}{a^{2}_{0}}\right)\right\\},$
(2.10)

where the function $Q(z)$ satisfies

$Q(z)\sim\left\\{\begin{array}[]{ll}z&\mbox{ for }z>>1,\\\\[12.0pt]
\frac{4}{3}z^{3/4}&\mbox{ for }z<<1.\end{array}\right.$ (2.11)

This leads to two equations

$\displaystyle\Delta\psi$ $\displaystyle=$ $\displaystyle 4\pi G\rho,$ (2.12)
$\displaystyle\Delta\varphi$ $\displaystyle=$
$\displaystyle\nabla_{i}\left[\nu\left(\frac{|\nabla\psi|}{a_{0}}\right)\nabla_{i}\psi\right],$
(2.13)

where $\nu(x)=Q^{\prime}(x^{2})$. By choosing new variable

$\phi=\varphi-\psi$ (2.14)

this action takes the form

$S_{\rm QUMOND}=\frac{1}{8\pi G}\int
dt\int\limits_{\Omega}dx\;\left\\{-|\nabla\varphi|^{2}+|\nabla\phi|^{2}+a_{0}^{2}H\left(\frac{|\nabla(\varphi-\phi)|^{2}}{a^{2}_{0}}\right)\right\\},$
(2.15)

where

$H(z)=Q(z)-z.$ (2.16)

Milgrom also proposed a bimetric relativistic MOND theory (called BIMOND) [17,
18] with two metrics, $g_{\mu\nu}$ and $\hat{g}_{\mu\nu}$, and two types of
matter, $\rho$ and $\hat{\rho}$ (called the twin matter). Let $M$ be a region
in spacetime with a spacelike boundary $\partial M$. Let $x^{\mu}$ be the
spacetime coordinates and $y^{i}$ be the coordinates on the boundary. Let
$\gamma_{ij}$ and $\hat{\gamma}_{ij}$ be the induced metrics on the boundary.
Let $g=-\det g_{\mu\nu}$, $\hat{g}=-\det\hat{g}_{\mu\nu}$, and
$\gamma=\det\gamma_{ij}$, $\hat{\gamma}=\det\hat{\gamma}_{ij}$. Let $R$ and
$\hat{R}$ be the scalar curvatures of the metrics $g_{\mu\nu}$ and
$\hat{g}_{\mu\nu}$, $K$ and $\hat{K}$ are the corresponding extrinsic
curvatures of the boundary. Let $\Gamma^{\alpha}{}_{\mu\beta}$ and
$\hat{\Gamma}^{\alpha}{}_{\mu\beta}$ be the Christoffel symbols of the metrics
$g_{\mu\nu}$ and $\hat{g}_{\mu\nu}$,

$C^{\alpha}{}_{\mu\beta}=\Gamma^{\alpha}{}_{\mu\beta}-\hat{\Gamma}^{\alpha}{}_{\mu\beta},$
(2.17)

and

$X=g^{\mu\nu}\left(C^{\alpha}{}_{\mu\beta}C^{\beta}{}_{\nu\alpha}-C^{\alpha}{}_{\mu\nu}C^{\beta}{}_{\beta\alpha}\right).$
(2.18)

The action of BIMOND theory has the form

$\displaystyle S_{\rm BIMOND}$ $\displaystyle=$ $\displaystyle\frac{1}{16\pi
G}\int\limits_{M}dx\left\\{\beta
g^{1/2}R+\alpha\hat{g}^{1/2}\hat{R}-2a_{0}^{2}g^{1/4}\hat{g}^{1/4}f\left(\frac{X}{a_{0}^{2}}\right)\right\\}{}$
(2.19) $\displaystyle+\frac{1}{8\pi G}\int\limits_{\partial
M}dy\;\left(\beta\gamma^{1/2}K+\alpha\hat{\gamma}^{1/2}\hat{K}\right),$

where $dx=dx^{0}\,dx^{1}\,dx^{2}\,dx^{3}$, $dy=dy^{1}\,dy^{2}\,dy^{3}$,
$\alpha$ and $\beta$ are some parameters, and $f(x)$ is a function that has
the properties

$f(x)\sim\left\\{\begin{array}[]{ll}\mathrm{const\,}&\mbox{ for
}x>>1,\\\\[12.0pt] \frac{4}{3}x^{3/4}&\mbox{ for }x<<1.\end{array}\right.$
(2.20)

In the static weak field limit with $g_{0i}=0$ and

$\displaystyle g_{00}=-1-2\varphi,\qquad\hat{g}_{00}=-1-2\hat{\varphi},$
(2.21)

we have

$X=|\nabla\psi|^{2},$ (2.22)

where

$\psi=\varphi-\hat{\varphi}.$ (2.23)

Therefore, the field equations become

$\displaystyle\beta\Delta\varphi$ $\displaystyle=$ $\displaystyle 4\pi
G\rho+\nabla_{i}\left[f^{\prime}\left(\frac{|\nabla\psi|^{2}}{a_{0}^{2}}\right)\nabla_{i}\psi\right],$
(2.24) $\displaystyle\alpha\Delta\hat{\varphi}$ $\displaystyle=$
$\displaystyle 4\pi
G\hat{\rho}-\nabla_{i}\left[f^{\prime}\left(\frac{|\nabla\psi|^{2}}{a_{0}^{2}}\right)\nabla_{i}\psi\right].$
(2.25)

By combining these equations we get a linear equation

$\Delta(\beta\varphi+\alpha\hat{\varphi})=4\pi G(\rho+\hat{\rho});$ (2.26)

In a particular case, $\beta=1$, $\alpha=-1$ this equation determines the
function $\psi$

$\Delta\psi=4\pi G(\rho+\hat{\rho}),$ (2.27)

and the functions $\varphi$ and $\hat{\varphi}$ are determined by the above
equations. This is equivalent to non-relativistic QUMOND theory.

##  3 Non-commutative Metric

An extended theory of gravity called Matrix General Relativity (MGR) was
developed in our papers [1, 2, 3]. In this approach all matter fields are
supposed to be vector-valued in a complex $N$-dimensional vector space and
gravity is described by a $N\times N$ matrix-valued self-adjoint symmetric
2-tensor field $a^{\mu\nu}$. We briefly review this approach (with some
modifications) and apply it to construct a phenomenological MOND model which
is, in fact, a generalization of Milgrom’s BIMOND theory.

The main idea of this approach is to replace the ingredients of Riemannian
geometry, like the metric, connection and curvature, by the real elements of a
noncommutative complex associative algebra $L(\mathcal{H})$ of operators in a
Hilbert space $\mathcal{H}$. We define the symmetrized multilinear product of
self-adjoint operators

$\langle\langle A_{1}\cdots A_{k}\rangle\rangle=\frac{1}{k!}\sum_{\varphi\in
S_{k}}A_{\varphi(1)}\cdots A_{\varphi(k)},$ (3.1)

where the summation goes over the permutation group $S_{k}$.

It is worth pointing out that this notation comes with a caveat. Consider two
matrices $A$, $B$, and their product $C$. We have

$\langle\langle C\rangle\rangle=C\neq\langle\langle
AB\rangle\rangle=\frac{1}{2}(AB+BA).$ (3.2)

This happens because the expression of the symmetrized product depends on the
number of factors, therefore the operations of taking the symmetrized product
and substituting a matrix with a product of matrices do not commute.

Obviously, the symmetrized product is completely symmetric, that is,

$\langle\langle A_{1}\cdots A_{i}\cdots A_{j}\cdots
A_{k}\rangle\rangle=\langle\langle A_{1}\cdots A_{j}\cdots A_{i}\cdots
A_{k}\rangle\rangle.$ (3.3)

In particular, this defines a bilinear product

$A\circ B=\langle\langle AB\rangle\rangle=\frac{1}{2}(AB+BA);$ (3.4)

it is easy to see that it is commutative but not associative. Further, we
define a linear functional (called the trace) on the algebra,

$\mathrm{tr\,}:L(\mathcal{H})\to\mathbb{R},$ (3.5)

which has the cyclic property

$\mathrm{tr\,}(A_{1}A_{2}\cdots A_{k})=\mathrm{tr\,}(A_{2}\cdots A_{k}A_{1}).$
(3.6)

Most of the time we will think of the Hilbert space $\mathcal{H}$ as a finite-
dimensional vector space $V$ and the algebra $L(\mathcal{H})$ simply as a
matrix algebra $GL(V,\mathbb{C})$. So, the trace will be just the usual trace
of a finite-dimensional matrix. Notice that the multilinear product satisfies
the Leibnitz rule

$D\langle\langle A_{1}\cdots A_{k}\rangle\rangle=\sum_{j=1}^{n}\langle\langle
A_{1}\cdots(DA_{j})\cdots A_{k}\rangle\rangle,$ (3.7)

where $D$ is a derivative.

Let $A^{\mu}{}_{\nu}$ be a matrix-valued tensor. It is easy to see that the
tensor

$T_{\nu_{1}\dots\nu_{n}}=\varepsilon_{\mu_{1}\dots\mu_{n}}\langle\langle
A^{\mu_{1}}{}_{\nu_{1}}\cdots A^{\mu_{n}}{}_{\nu_{n}}\rangle\rangle,$ (3.8)

where $\varepsilon_{\mu_{1}\dots\mu_{n}}$ is the Levi-Civita symbol, is
completely anti-symmetric and, therefore, it defines a matrix that we call the
symmetrized determinant

$\varepsilon_{\mu_{1}\dots\mu_{n}}\langle\langle A^{\mu_{1}}{}_{\nu_{1}}\cdots
A^{\mu_{n}}{}_{\nu_{n}}\rangle\rangle=\varepsilon_{\nu_{1}\dots\nu_{n}}\;\mathrm{sdet\,}A^{\mu}{}_{\nu},$
(3.9)

that is,

$\displaystyle\mathrm{sdet\,}A^{\mu}{}_{\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{n!}\varepsilon_{\mu_{1}\dots\mu_{n}}\varepsilon^{\nu_{1}\dots\nu_{n}}\langle\langle
A^{\mu_{1}}{}_{\nu_{1}}\cdots A^{\mu_{n}}{}_{\nu_{n}}\rangle\rangle{}$ (3.10)
$\displaystyle=$ $\displaystyle\langle\langle A^{\mu_{1}}{}_{[\mu_{1}}\cdots
A^{\mu_{n}}{}_{\mu_{n}]}\rangle\rangle\,.$

This is proved by using the property

$\varepsilon_{\mu_{1}\dots\mu_{n}}\varepsilon^{\nu_{1}\dots\nu_{n}}=n!\delta^{\nu_{1}}{}_{[\mu_{1}}\cdots\delta^{\nu_{n}}{}_{\mu_{n}]}\,.$
(3.11)

Similarly we define

$\displaystyle\mathrm{sdet\,}(B^{\mu\nu})=\frac{1}{n!}\varepsilon_{\mu_{1}\dots\mu_{n}}\varepsilon_{\nu_{1}\dots\nu_{n}}\langle\langle
B^{\mu_{1}\nu_{1}}\cdots B^{\mu_{n}\nu_{n}}\rangle\rangle\,,$ (3.12)
$\displaystyle\mathrm{sdet\,}(C_{\mu\nu})=\frac{1}{n!}\varepsilon^{\mu_{1}\dots\mu_{n}}\varepsilon^{\nu_{1}\dots\nu_{n}}\langle\langle
C_{\mu_{1}\nu_{1}}\cdots C_{\mu_{n}\nu_{n}}\rangle\rangle\,.$ (3.13)

Notice that the symmetrized determinant of a $(1,1)$ tensor is invariant under
diffeomorphisms, while the determinant of tensors of type $(2,0)$ and $(0,2)$
define scalar densities of the weight $-2$ and $2$ correspondingly. Also,
under the transformation

$A^{\prime\mu}{}_{\nu}=UA^{\mu}{}_{\nu}U^{-1}$ (3.14)

it transforms accordingly

$\mathrm{sdet\,}A^{\prime\mu}{}_{\nu}=U(\mathrm{sdet\,}A^{\mu}{}_{\nu})U^{-1}\,.$
(3.15)

Further, it is easy to see that the symmetrized determinant is multiplicative
for any two commuting tensors, that is, if

$C_{\mu\nu}B^{\alpha\beta}=B^{\alpha\beta}C_{\mu\nu},$ (3.16)

then for a tensor $A^{\mu}{}_{\nu}=B^{\mu\alpha}C_{\alpha\nu}$ there holds

$\mathrm{sdet\,}A=\mathrm{sdet\,}B\;\mathrm{sdet\,}C.$ (3.17)

###### Lemma 1

There holds

$\displaystyle\mathrm{sdet\,}(I+A)$ $\displaystyle=$ $\displaystyle
I+\sum_{k=1}^{n}\langle\langle A^{\mu_{1}}{}_{[\mu_{1}}\cdots
A^{\mu_{k}}{}_{\mu_{k}]}\rangle\rangle$ $\displaystyle=$ $\displaystyle
I+A+\frac{1}{2}A^{2}-\frac{1}{2}\langle\langle
A^{\mu}{}_{\nu}A^{\nu}{}_{\mu}\rangle\rangle+\cdots+\mathrm{sdet\,}A,{}$

where

$A=A^{\mu}{}_{\mu}\,.$ (3.19)

Let $a^{\mu\nu}$ be a self-adjoint symmetric matrix valued tensor, that is,

$(a^{\mu\nu})^{*}=a^{\mu\nu},\qquad a^{\mu\nu}=a^{\nu\mu}.$ (3.20)

Let $b_{\mu\nu}$ be another matrix-values tensor defined by

$a^{\mu\nu}b_{\nu\alpha}=b_{\alpha\nu}a^{\nu\mu}=\delta^{\mu}{}_{\alpha}I;$
(3.21)

this tensor is neither self-adjoint nor symmetric but rather satisfies

$(b_{\mu\nu})^{*}=b_{\nu\mu}.$ (3.22)

Let $g_{\mu\nu}$ be a pseudo-Riemannian metric and $h^{\mu\nu}$ be some self-
adjoint symmetric matrix valued tensor defined by We define

$\displaystyle a^{\mu\nu}$ $\displaystyle=$ $\displaystyle
g^{\mu\nu}I+h^{\mu\nu},$ (3.23)

where $I$ is the unit matrix. We do not assume that the matrix valued tensor
$h^{\mu\nu}$ is traceless.

We raise the indices with the metric $g_{\mu\nu}$. Then the tensor
$b_{\mu\nu}$ can be found in perturbation theory in the tensor $h$

$b_{\mu\nu}=g_{\mu\nu}I+\sum_{k=1}^{\infty}(-1)^{k}h_{\mu}{}^{\beta_{1}}h_{\beta_{1}}{}^{\beta_{2}}\cdots
h_{\beta_{k-2}}{}^{\beta_{k-1}}h_{\beta_{k-1}\nu}.$ (3.24)

We compute the symmetrized determinant of the tensor $a^{\mu\nu}$

$\mathrm{sdet\,}a^{\mu\nu}=-g^{-1}(I+Q),$ (3.25)

where $g=-\det g_{\mu\nu}$ and

$Q=\sum_{k=1}^{n}\langle\langle h^{\mu_{1}}{}_{[\mu_{1}}\cdots
h^{\mu_{k}}{}_{\mu_{k}]}\rangle\rangle;$ (3.26)

this enables us to define a self-adjoint matrix-valued density $\rho$ of
weight $1$ by

$\rho=(-\mathrm{sdet\,}a^{\mu\nu})^{-1/2}=g^{1/2}\Phi,$ (3.27)

where

$\displaystyle\Phi$ $\displaystyle=$ $\displaystyle(I+Q)^{-1/2}{}$ (3.28)
$\displaystyle=$ $\displaystyle
I+\sum_{k=1}^{\infty}(-1)^{k}\frac{(2k-1)!}{k!(k-1)!2^{2k-1}}Q^{k}.$

##  4 Non-commutative Connection and Curvature

Let $\nabla_{\mu}$ be the covariant derivative defined with respect to the
Levi-Civita connection of the metric $g_{\mu\nu}$ so that

$\displaystyle\nabla_{\alpha}a^{\mu\nu}$ $\displaystyle=$
$\displaystyle\nabla_{\alpha}h^{\mu\nu}.$ (4.1)

Let $\varphi^{\mu}$ and $\psi_{\mu}$ be vector-valued tensors. We define the
non-commutative covariant derivative by

$\displaystyle\mathcal{D}_{\mu}\varphi^{\alpha}=\partial_{\mu}\varphi^{\alpha}+\mathcal{A}^{\alpha}{}_{\beta\mu}\varphi^{\beta},$
(4.2)
$\displaystyle\mathcal{D}_{\mu}\psi_{\alpha}=\partial_{\mu}\psi_{\alpha}-\mathcal{A}^{\beta}{}_{\alpha\mu}\psi_{\beta},$
(4.3)

where $\mathcal{A}^{\alpha}{}_{\beta\mu}$ is a matrix valued connection. We
impose the following conditions on the connection.

  1. 1.

For the operator $\mathcal{D}_{\mu}$ to be well defined it should transform as
connection under the diffeomorphisms $x^{\prime\mu}=x^{\prime\mu}(x)$,

$\displaystyle\mathcal{A}^{\prime\varepsilon}{}_{\gamma\delta}$
$\displaystyle=$ $\displaystyle
A^{\varepsilon}{}_{\theta}B^{\xi}{}_{\gamma}B^{\phi}{}_{\delta}\left(\mathcal{A}^{\theta}{}_{\xi\phi}-E^{\theta}{}_{\xi\phi}\right),$
(4.4)

where

$\displaystyle A^{\mu}{}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{\partial x^{\prime\mu}}{\partial x^{\alpha}},$ (4.5)
$\displaystyle B^{\mu}{}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{\partial x^{\mu}}{\partial x^{\prime\alpha}},$ (4.6)
$\displaystyle E^{\alpha}{}_{\kappa\sigma}$ $\displaystyle=$ $\displaystyle
B^{\alpha}{}_{\lambda}C^{\lambda}{}_{\kappa\sigma},$ (4.7) $\displaystyle
C^{\mu}{}_{\alpha\sigma}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}x^{\prime\mu}}{\partial x^{\sigma}\partial
x^{\alpha}}.$ (4.8)

Therefore, the connection should have the form

$\mathcal{A}^{\alpha}{}_{\mu\nu}=\Gamma^{\alpha}{}_{\mu\nu}I+\theta^{\alpha}{}_{\mu\nu},$
(4.9)

where $\theta^{\alpha}{}_{\mu\nu}$ is a matrix-valued tensor. Then the non-
commutative covariant derivative takes the form

$\displaystyle\mathcal{D}_{\mu}\varphi^{\alpha}=\nabla_{\mu}\varphi^{\alpha}+\theta^{\alpha}{}_{\beta\mu}\varphi^{\beta},$
(4.10)
$\displaystyle\mathcal{D}_{\mu}\psi_{\alpha}=\nabla_{\mu}\psi_{\alpha}-\theta^{\beta}{}_{\alpha\mu}\psi_{\beta}.$
(4.11)

  2. 2.

Let $\langle\;,\;\rangle$ be the inner product in the vector space
$\mathcal{H}$. Then we require

$\langle\mathcal{D}_{\mu}\psi_{\alpha},\varphi^{\alpha}\rangle+\langle\psi_{\alpha},\mathcal{D}_{\mu}\varphi^{\alpha}\rangle=\partial_{\mu}\langle\psi_{\alpha},\varphi^{\alpha}\rangle.$
(4.12)

Then the connection must be self-adjoint,

$\mathcal{A}^{\alpha}{}_{\mu\nu}^{*}=\mathcal{A}^{\alpha}{}_{\mu\nu},$ (4.13)

and, therefore, the tensor $\theta$ is also self-adjoint

$\theta^{\alpha}{}_{\mu\nu}^{*}=\theta^{\alpha}{}_{\mu\nu}.$ (4.14)

  3. 3.

We also require the connection to be torsion-free

$\mathcal{A}^{\alpha}{}_{\mu\nu}=\mathcal{A}^{\alpha}{}_{\nu\mu},$ (4.15)

therefore,

$\theta^{\alpha}{}_{\mu\nu}=\theta^{\alpha}{}_{\nu\mu}.$ (4.16)

Now we need to fix the connection, i.e. to relate it somehow to the tensors
$a^{\mu\nu}$ and $b_{\mu\nu}$ and their first derivatives. The derivatives of
the tensor $a^{\mu\nu}$ transform under diffeomorphism
$x^{\prime}=x^{\prime}(x)$ as

$\frac{\partial a^{\prime\mu\nu}}{\partial
x^{\prime\rho}}=A^{\mu}{}_{\alpha}A^{\nu}{}_{\beta}B^{\sigma}{}_{\rho}\left\\{\frac{\partial
a^{\alpha\beta}}{\partial
x^{\sigma}}+E^{\alpha}{}_{\kappa\sigma}a^{\kappa\beta}+E^{\beta}{}_{\kappa\sigma}a^{\alpha\kappa}\right\\}.$
(4.17)

To fix the connection we should impose an invariant compatibility condition.
We will impose the compatibility condition in the form

$W^{\alpha\beta}{}_{\mu}=0,$ (4.18)

where

$W^{\alpha\beta}{}_{\mu}=\partial_{\mu}a^{\alpha\beta}+c_{1}{\cal
A}^{\alpha}{}_{\lambda\mu}a^{\lambda\beta}+c_{2}{\cal
A}^{\beta}{}_{\lambda\mu}a^{\alpha\lambda}+c_{3}a^{\lambda\beta}{\cal
A}^{\alpha}{}_{\lambda\mu}+c_{4}a^{\alpha\lambda}{\cal
A}^{\beta}{}_{\lambda\mu}.$ (4.19)

By using the transformation rule of the connection one can show that this
quantity transforms as a tensor if

$c_{1}+c_{3}=1.\,\qquad c_{2}+c_{4}=1.$ (4.20)

Then it can be written in the form

$\displaystyle W^{\alpha\beta}{}_{\mu}$ $\displaystyle=$
$\displaystyle\nabla_{\mu}a^{\alpha\beta}+c_{1}\theta^{\alpha}{}_{\lambda\mu}a^{\lambda\beta}+c_{2}\theta^{\beta}{}_{\lambda\mu}a^{\alpha\lambda}+(1-c_{1})a^{\lambda\beta}\theta^{\alpha}{}_{\lambda\mu}+(1-c_{2})a^{\alpha\lambda}\theta^{\beta}{}_{\lambda\mu}{}$
(4.21) $\displaystyle=$
$\displaystyle\nabla_{\mu}a^{\alpha\beta}+c_{1}R^{\alpha\beta}{}_{\mu}+c_{2}R^{\beta\alpha}{}_{\mu}+(1-c_{1})L^{\alpha\beta}{}_{\mu}+(1-c_{2})L^{\beta\alpha}{}_{\mu},$

where

$\displaystyle R^{\alpha\beta}{}_{\mu}$ $\displaystyle=$
$\displaystyle\theta^{\alpha}{}_{\lambda\mu}a^{\beta\lambda},$ (4.22)
$\displaystyle L^{\alpha\beta}{}_{\mu}$ $\displaystyle=$ $\displaystyle
a^{\beta\lambda}\theta^{\alpha}{}_{\lambda\mu}.$ (4.23)

Then the compatibility condition has the form

$\nabla_{\mu}a^{\alpha\beta}+c_{1}R^{\alpha\beta}{}_{\mu}+c_{2}R^{\beta\alpha}{}_{\mu}+(1-c_{1})L^{\alpha\beta}{}_{\mu}+(1-c_{2})L^{\beta\alpha}{}_{\mu}=0.$
(4.24)

By antisymmetrizing over the indices $\alpha$ and $\beta$ we get

$(c_{1}-c_{2})\left(R^{\alpha\beta}{}_{\mu}-R^{\beta\alpha}{}_{\mu}-L^{\alpha\beta}{}_{\mu}+L^{\beta\alpha}{}_{\mu}\right)=0.$
(4.25)

By taking the anti-adjoint part of eq. (4.24) we get

$\displaystyle(2c_{1}-1)(R^{\alpha\beta}{}_{\mu}-L^{\alpha\beta}{}_{\mu})+(2c_{2}-1)(R^{\beta\alpha}{}_{\mu}-L^{\beta\alpha}{}_{\mu})=0.$
(4.26)

We rewrite these equations in the form

$\displaystyle(c_{1}-c_{2})\left(S^{\alpha\beta}{}_{\mu}-S^{\beta\alpha}{}_{\mu}\right)$
$\displaystyle=$ $\displaystyle 0,$ (4.27)
$\displaystyle(2c_{1}-1)S^{\alpha\beta}{}_{\mu}+(2c_{2}-1)S^{\beta\alpha}{}_{\mu}$
$\displaystyle=$ $\displaystyle 0,$ (4.28)

where

$S^{\alpha\beta}{}_{\mu}=\frac{1}{2}\left(L^{\alpha\beta}{}_{\mu}-R^{\alpha\beta}{}_{\mu}\right)=\frac{1}{2}[a^{\beta\lambda},\theta^{\alpha}{}_{\lambda\mu}].$
(4.29)

By decomposing the tensor $S$ in symmetric and anti-symmetric parts in the
first two indices we get

$\displaystyle(c_{1}-c_{2})S^{[\alpha\beta]}{}_{\mu}$ $\displaystyle=$
$\displaystyle 0,$ (4.30)
$\displaystyle(c_{1}+c_{2}-1)S^{(\alpha\beta)}{}_{\mu}$ $\displaystyle=$
$\displaystyle 0.$ (4.31)

Therefore, if $c_{1}\neq c_{2}$ then the tensor $S$ is symmetric in the first
two indices

$S^{[\alpha\beta]}{}_{\mu}=0,$ (4.32)

and if $c_{1}+c_{2}\neq 1$ then the tensor $S$ is anti-symmetric in the first
two indices

$S^{(\alpha\beta)}{}_{\mu}=0.$ (4.33)

Thus, if both of these condistions are satisfied then the tensor $S$ vanishes,

$S^{\alpha\beta}{}_{\mu}=0.$ (4.34)

These conditions puts severe algebraic constraints on the algebra. Therefore,
we will choose $c_{1}$ and $c_{2}$ in such a way to avoid additional algebraic
constraints coming from the compatibility conditions.

Let us try to use eq. (4.24) with $c_{1}=c_{2}=1$ or $c_{1}=c_{2}=0$.

Case I: In the case

$c_{1}=c_{2}=1$ (4.35)

the compatibility condition looks like

$\nabla_{\mu}a^{\alpha\beta}+R^{\alpha\beta}{}_{\mu}+R^{\beta\alpha}{}_{\mu}=0$
(4.36)

or

$K^{\mu\alpha\beta}+X^{\alpha\mu\beta}+X^{\beta\mu\alpha}=0,$ (4.37)

where

$\displaystyle K^{\mu\alpha\beta}$ $\displaystyle=$ $\displaystyle
a^{\mu\gamma}\nabla_{\gamma}a^{\alpha\beta},$ (4.38) $\displaystyle
X^{\alpha\beta\gamma}$ $\displaystyle=$ $\displaystyle
a^{\beta\lambda}\theta^{\alpha}{}_{\lambda\mu}a^{\mu\gamma}$ (4.39)

This equation can be solved in a closed form. The solution has the form

$X^{\alpha\mu\beta}=-\frac{1}{2}\left(K^{\alpha\mu\beta}-K^{\mu\alpha\beta}-K^{\beta\alpha\mu}\right)+Y^{\alpha\mu\beta},$
(4.40)

where $Y$ is an arbitrary tensor which is antisymmetric in the first and the
third index,

$Y^{\alpha\beta\mu}=-Y^{\mu\beta\alpha}.$ (4.41)

Case II: In the case

$c_{1}=c_{2}=0$ (4.42)

the compatibility condition looks like

$\nabla_{\mu}a^{\alpha\beta}+L^{\alpha\beta}{}_{\mu}+L^{\beta\alpha}{}_{\mu}=0$
(4.43)

or

$N^{\mu\alpha\beta}+X^{\alpha\mu\beta}+X^{\beta\mu\alpha}=0,$ (4.44)

where

$N^{\mu\alpha\beta}=(\nabla_{\gamma}a^{\alpha\beta})a^{\mu\gamma}.$ (4.45)

This equation can be solved in a closed form. The solution has the form

$X^{\alpha\mu\beta}=-\frac{1}{2}\left(N^{\alpha\mu\beta}-N^{\mu\alpha\beta}-N^{\beta\alpha\mu}\right)+Z^{\alpha\mu\beta},$
(4.46)

where $Z$ is an arbitrary tensor which is antisymmetric in the first and the
third index,

$Z^{\alpha\beta\mu}=-Z^{\mu\beta\alpha}$ (4.47)

However, the tensors $K$ and $N$ are not self-adjoint and, therefore, the
connection is also not self-adjoint.

Case III: We might try the compatibility condition with the noncommutative
metric, that is,

$\displaystyle\mathcal{D}_{\mu}(a^{\alpha\beta}\psi_{\beta})=a^{\alpha\beta}\mathcal{D}_{\mu}\psi_{\beta}.$
(4.48)

This condition is satisfied if

$\nabla_{\mu}a^{\alpha\beta}+\theta^{\alpha}{}_{\lambda\mu}a^{\lambda\beta}+a^{\alpha\lambda}\theta^{\beta}{}_{\lambda\mu}=0,$
(4.49)

which is equivalent to

$\nabla_{\mu}b_{\alpha\beta}-\theta^{\lambda}{}_{\alpha\mu}b_{\lambda\beta}-b_{\alpha\lambda}\theta^{\lambda}{}_{\beta\mu}=0,$
(4.50)

and guarantees the compatibility condition

$\displaystyle\mathcal{D}_{\mu}(b_{\alpha\beta}\varphi^{\beta})=b_{\alpha\beta}\mathcal{D}_{\mu}\varphi^{\beta}.$
(4.51)

This compatibility condition (4.49) corresponds to

$c_{1}=1,\qquad c_{2}=0.$ (4.52)

In this case the compatibility conditions reads

$\nabla_{\mu}a^{\alpha\beta}+R^{\alpha\beta}{}_{\mu}+L^{\beta\alpha}{}_{\mu}=0,$
(4.53)

and it cannot be solved in a closed form. Let

$\displaystyle B^{\alpha\beta}{}_{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(L^{\alpha\beta}{}_{\mu}+R^{\alpha\beta}{}_{\mu}\right)=a^{\beta\lambda}\circ\theta^{\alpha}{}_{\lambda\mu},$
(4.54)

Then we have

$\displaystyle B^{(\alpha\beta)}{}_{\mu}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\nabla_{\mu}a^{\alpha\beta},$ (4.56) $\displaystyle
S^{[\alpha\beta]}{}_{\mu}$ $\displaystyle=$ $\displaystyle 0,$ (4.57)

which means

$[a^{\beta\lambda},\theta^{\alpha}{}_{\lambda\mu}]=[a^{\alpha\lambda},\theta^{\beta}{}_{\lambda\mu}].$
(4.58)

These constraints cannot be satisfied even in perturbation theory for a
generic matrix algebra.

Case IV: That is why, we choose

$c_{1}=c_{2}=\frac{1}{2}.$ (4.59)

In this case

$W^{\alpha\beta}{}_{\mu}=\partial_{\mu}a^{\alpha\beta}+{\cal
A}^{\alpha}{}_{\lambda\mu}\circ a^{\lambda\beta}+{\cal
A}^{\beta}{}_{\lambda\mu}\circ a^{\alpha\lambda}$ (4.60)

and the compatibility condition takes the form

$\displaystyle\nabla_{\mu}a^{\alpha\beta}+\theta^{\alpha}{}_{\lambda\mu}\circ
a^{\lambda\beta}+\theta^{\beta}{}_{\lambda\mu}\circ a^{\alpha\lambda}=0.$
(4.61)

Therefore, the compatibility condition fixes the symmetric part of the tensor
$B$,

$B^{(\alpha\beta)}{}_{\mu}=-\frac{1}{2}\nabla_{\mu}a^{\alpha\beta},$ (4.62)

but leaves the antisymmetric part of the tensor $B$ and the tensor $S$
arbitrary. Therefore, there are no extra algebraic constraints.

We solve the compatibility condition (4.61) in perturbation theory (indices
are raised with the metric $g$). We have

$\theta^{\alpha\beta}{}_{\mu}+\theta^{\beta\alpha}{}_{\mu}=-\nabla_{\mu}h^{\alpha\beta}-\theta^{\alpha}{}_{\lambda\mu}\circ
h^{\beta\lambda}-\theta^{\beta}{}_{\lambda\mu}\circ h^{\alpha\lambda},$ (4.63)

or

$\theta_{\alpha\beta\mu}+\theta_{\beta\alpha\mu}=-P_{\alpha\beta\mu},$ (4.64)

where

$P_{\alpha\beta\mu}=\nabla_{\mu}h_{\alpha\beta}+\theta_{\alpha\lambda\mu}\circ
h^{\lambda}{}_{\beta}+\theta_{\beta\lambda\mu}\circ h^{\lambda}{}_{\alpha}.$
(4.65)

Recall that $\theta_{\mu\nu\alpha}$ is symmetric in the last two indices.

###### Lemma 2

Let $T_{\alpha\beta\mu}$ be a tensor satisfying the conditions

$\displaystyle T_{\alpha\beta\mu}$ $\displaystyle=$ $\displaystyle
T_{\alpha\mu\beta},$ (4.66) $\displaystyle
T_{\alpha\beta\mu}+T_{\beta\alpha\mu}$ $\displaystyle=$ $\displaystyle
V_{\alpha\beta\mu}.$ (4.67)

Then

$T_{\alpha\beta\mu}=\frac{1}{2}\left(V_{\alpha\beta\mu}+V_{\alpha\mu\beta}-V_{\beta\mu\alpha}\right)$
(4.68)

Proof: We have

$\displaystyle T_{\alpha\beta\mu}+T_{\beta\alpha\mu}$ $\displaystyle=$
$\displaystyle V_{\alpha\beta\mu},$ (4.69) $\displaystyle
T_{\mu\alpha\beta}+T_{\alpha\mu\beta}$ $\displaystyle=$ $\displaystyle
V_{\alpha\mu\beta},$ (4.70) $\displaystyle-
T_{\mu\beta\alpha}-T_{\beta\mu\alpha}$ $\displaystyle=$ $\displaystyle-
V_{\beta\mu\alpha}.$ (4.71)

By adding these equations we obtain the result.

The equation (4.64) can be solved in terms of the tensor $P$,

$\theta_{\alpha\beta\mu}=\frac{1}{2}\left(P_{\beta\mu\alpha}-P_{\alpha\mu\beta}-P_{\alpha\beta\mu}\right),$
(4.72)

that is,

$\displaystyle\theta_{\alpha\beta\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Biggl{\\{}\nabla_{\alpha}h_{\mu\beta}-\nabla_{\mu}h_{\alpha\beta}-\nabla_{\beta}h_{\alpha\mu}$
$\displaystyle+h^{\lambda}{}_{\beta}\circ(\theta_{\mu\lambda\alpha}-\theta_{\alpha\lambda\mu})+h^{\lambda}{}_{\mu}\circ(\theta_{\beta\lambda\alpha}-\theta_{\alpha\lambda\beta})-h^{\lambda}{}_{\alpha}\circ(\theta_{\mu\lambda\beta}+\theta_{\beta\lambda\mu})\Biggr{\\}}.{}$

This equation can be solved in perturbation theory. By rescaling
$h\mapsto\varkappa h$ and expanding the tensor $\theta$ in the power series

$\theta=\sum_{k=1}^{\infty}\varkappa^{k}\theta_{(k)},$ (4.74)

we obtain in the first order

$\theta_{(1),\alpha\beta\mu}=\frac{1}{2}\left(\nabla_{\alpha}h_{\mu\beta}-\nabla_{\mu}h_{\alpha\beta}-\nabla_{\beta}h_{\alpha\mu}\right),$
(4.75)

and for the higher orders, $k\geq 2$, we get the recurrence relations

$\displaystyle\theta_{(k)\alpha\beta\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Biggl{\\{}h^{\lambda}{}_{\beta}\circ(\theta_{(k-1)\mu\lambda\alpha}-\theta_{\alpha\lambda\mu})+h^{\lambda}{}_{\mu}\circ(\theta_{(k-1)\beta\lambda\alpha}-\theta_{\alpha\lambda\beta}){}$
(4.76)
$\displaystyle-h^{\lambda}{}_{\alpha}\circ(\theta_{(k-1)\mu\lambda\beta}+\theta_{(k-1)\beta\lambda\mu})\Biggr{\\}}.$

In particular, in the second order we get

$\theta_{(2),\alpha\beta\mu}=\frac{1}{2}\Biggl{\\{}h^{\lambda}{}_{\alpha}\circ\nabla_{\lambda}h_{\mu\beta}+h^{\lambda}{}_{\beta}\circ(\nabla_{\mu}h_{\alpha\lambda}-\nabla_{\alpha}h_{\mu\lambda})+h^{\lambda}{}_{\mu}\circ(\nabla_{\beta}h_{\alpha\lambda}-\nabla_{\alpha}h_{\beta\lambda})\Biggr{\\}}.$
(4.77)

It is important to emphasize that this defines the tensor $\theta$ uniquely.

We compute the contraction of the connection

$\mathcal{A}^{\alpha}{}_{\alpha\mu}=-\omega_{,\mu}+\theta_{\mu},$ (4.78)

where $\omega=-\log g^{1/2}$, and

$\theta_{\mu}=\theta^{\alpha}{}_{\alpha\mu}.$ (4.79)

In the first order we have simply

$\theta_{(1),\mu}=-\frac{1}{2}\nabla_{\mu}h,$ (4.80)

where

$h=h^{\mu}{}_{\mu}.$ (4.81)

For the second order we obtain

$\theta_{(2),\mu}=\frac{1}{2}h^{\alpha\beta}\circ\nabla_{\mu}h_{\alpha\beta}=\frac{1}{4}\nabla_{\mu}(h^{\alpha\beta}\circ
h_{\alpha\beta})$ (4.82)

Similarly, we expand the symmetrized determinant,
$\rho=(-\mathrm{sdet\,}a^{\mu\nu})^{-1/2}$, of the noncommutative metric more
precisely, the matrix $\Phi=g^{-1/2}\rho$ defined by eqs. (3.27), (3.28),

$\Phi=I+\sum_{k=1}^{\infty}\varkappa^{k}\Phi_{(k)}.$ (4.83)

By using (3.28) we get

$\displaystyle\Phi_{1}$ $\displaystyle=$ $\displaystyle-\frac{1}{2}h,$ (4.84)
$\displaystyle\Phi_{2}$ $\displaystyle=$ $\displaystyle\frac{1}{8}h\circ
h+\frac{1}{4}h_{\alpha\beta}\circ h^{\alpha\beta}.$ (4.85)

Therefore,

$\nabla_{\mu}\Phi_{1}=\theta_{(1),\mu}=-\frac{1}{2}\nabla_{\mu}h$ (4.86)

and

$\displaystyle\nabla_{\mu}\Phi_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4}h\circ\nabla_{\mu}h+\frac{1}{4}\nabla_{\mu}(h_{\alpha\beta}\circ
h^{\alpha\beta}){}$ (4.87) $\displaystyle=$
$\displaystyle\theta_{(1),\mu}\circ\Phi_{1}+\theta_{(2),\mu}.$

Therefore, the matrix $\Phi$ satisfies the equation

$\nabla_{\mu}\Phi-\theta_{\mu}\circ\Phi=O(\varkappa^{3})$ (4.88)

up to the third order.

Now we can define the curvature (noncommutative Riemann tensor) as follows

$\displaystyle\mathcal{R}^{\alpha}{}_{\beta\mu\nu}$ $\displaystyle=$
$\displaystyle\partial_{\mu}\mathcal{A}^{\alpha}{}_{\beta\nu}-\partial_{\nu}\mathcal{A}^{\alpha}{}_{\beta\mu}+\mathcal{A}^{\alpha}{}_{\gamma\mu}\circ\mathcal{A}^{\gamma}{}_{\beta\nu}-\mathcal{A}^{\alpha}{}_{\gamma\nu}\circ\mathcal{A}^{\gamma}{}_{\beta\mu}{}$
(4.89) $\displaystyle=$ $\displaystyle
R^{\alpha}{}_{\beta\mu\nu}I+\nabla_{\mu}\theta^{\alpha}{}_{\beta\nu}-\nabla_{\nu}\theta^{\alpha}{}_{\beta\mu}+\theta^{\alpha}{}_{\gamma\mu}\circ\theta^{\gamma}{}_{\beta\nu}-\theta^{\alpha}{}_{\gamma\nu}\circ\theta^{\gamma}{}_{\beta\mu}.$

This tensor is antisymmetric in the last two indices and is self-adjoint. It
naturally leads to the unique noncommutative Ricci tensor

$\displaystyle\mathcal{R}_{\beta\nu}=\mathcal{R}^{\mu}{}_{\beta\mu\nu}$
$\displaystyle=$
$\displaystyle\partial_{\mu}\mathcal{A}^{\mu}{}_{\beta\nu}-\partial_{\nu}\mathcal{A}_{\beta}+\mathcal{A}_{\gamma}\circ\mathcal{A}^{\gamma}{}_{\beta\nu}-\mathcal{A}^{\mu}{}_{\gamma\nu}\circ\mathcal{A}^{\gamma}{}_{\beta\mu}{}$
(4.90) $\displaystyle=$ $\displaystyle
R_{\beta\nu}I+\nabla_{\mu}\theta^{\mu}{}_{\beta\nu}-\nabla_{\nu}\theta_{\beta}+\theta_{\gamma}\circ\theta^{\gamma}{}_{\beta\nu}-\theta^{\mu}{}_{\gamma\nu}\circ\theta^{\gamma}{}_{\beta\mu}.$

This tensor is self-adjoint but not symmetric. Next, we define the scalar
curvature by

$\displaystyle\mathcal{R}$ $\displaystyle=$ $\displaystyle
a^{\beta\nu}\circ\mathcal{R}_{\beta\nu}{}$ (4.91) $\displaystyle=$
$\displaystyle
a^{\beta\nu}\circ(R_{\beta\nu}I+\nabla_{\mu}\theta^{\mu}{}_{\beta\nu}-\nabla_{\nu}\theta_{\beta})+\Theta,$

where

$\Theta=a^{\beta\nu}\circ\left(\theta_{\gamma}\circ\theta^{\gamma}{}_{\nu\beta}-\theta^{\mu}{}_{\gamma\nu}\circ\theta^{\gamma}{}_{\mu\beta}\right);$
(4.92)

which is also self-adjoint.

##  5 Noncommutative Lagrangian

To construct the action of the noncommutative gravity we need to construct
invariant scalar densities from the noncommutative metric and its derivatives.
We can use the scalar curvature to define a real scalar density

$\mathcal{L}=\frac{1}{N}\mathrm{tr\,}\left(\rho\;\mathcal{R}\right);$ (5.1)

more explicitly, this invariant has the form

$\mathcal{L}=\frac{1}{N}g^{1/2}\mathrm{tr\,}(\Phi\circ
a^{\beta\nu})\Biggl{\\{}R_{\beta\nu}I+\nabla_{\mu}\theta^{\mu}{}_{\beta\nu}-\nabla_{\nu}\theta_{\beta}+\theta_{\gamma}\circ\theta^{\gamma}{}_{\beta\nu}-\theta^{\mu}{}_{\gamma\nu}\circ\theta^{\gamma}{}_{\beta\mu}\Biggr{\\}}.$
(5.2)

We note that the commutators of the connection

$W^{\alpha\beta}{}_{\mu\nu\rho\sigma}=[\mathcal{A}^{\alpha}{}_{\mu\nu},\mathcal{A}^{\beta}{}_{\rho\sigma}]=[\theta^{\alpha}{}_{\mu\nu},\theta^{\beta}{}_{\rho\sigma}]$
(5.3)

as well as the commutators of the connection with the noncommutative metric
$a^{\mu\nu}$,

$V^{\mu\nu\alpha}{}_{\rho\sigma}=[a^{\mu\nu},\mathcal{A}^{\alpha}{}_{\rho\sigma}]=[h^{\mu\nu},\theta^{\alpha}{}_{\rho\sigma}],$
(5.4)

transform as tensors. These tensors are anti-self-adjoint. One can use these
tensors to construct real scalars; one needs to add terms linear in $W$ and
quadratic in $V$

$\mathcal{L}\sim\frac{1}{N}g^{1/2}\mathrm{tr\,}\Phi\Biggl{\\{}ia^{\mu\nu}W^{\alpha\beta}{}_{\mu\nu\alpha\beta}+\dots+b_{\alpha\beta}V^{\mu\nu\alpha}{}_{\rho\sigma}V^{\rho\sigma\beta}{}_{\mu\nu}+\dots\Biggr{\\}}$
(5.5)

In particular, we define a natural self-adjoint symmetric tensor is

$\mathcal{F}_{\mu\nu}=iW^{\alpha\beta}{}_{\mu\nu\alpha\beta}=i[\mathcal{A}^{\alpha}{}_{\mu\nu},\mathcal{A}_{\alpha}]=i[\theta^{\alpha}{}_{\mu\nu},\theta_{\alpha}];$
(5.6)

it defines the anti-self-adjoint scalar

$\mathcal{F}=ia^{\mu\nu}\circ[\theta^{\alpha}{}_{\mu\nu},\theta_{\alpha}],$
(5.7)

which can be used to form a real non-commutative action term

$\mathcal{L}=\frac{1}{N}\mathrm{tr\,}\left(\rho\mathcal{F}\right);$ (5.8)

more explicitly

$\mathcal{L}=\frac{1}{N}g^{1/2}\mathrm{tr\,}(\Phi\circ
a^{\mu\nu})i[\theta^{\alpha}{}_{\mu\nu},\theta_{\alpha}].$ (5.9)

We could also use double commutators to construct tensors

$[\mathcal{A}^{\lambda}{}_{\phi\chi},[\mathcal{A}^{\alpha}{}_{\mu\nu},\mathcal{A}^{\beta}{}_{\rho\sigma}]]=[\theta^{\lambda}{}_{\phi\chi},[\theta^{\alpha}{}_{\mu\nu},\theta^{\beta}{}_{\rho\sigma}]]$
(5.10)

However, these tensors are of odd rank and it is impossible to construct
scalars out of them. Also, they are cubic in the derivatives, $(\nabla
h)^{3}$.

##  6 MOND-like Noncommutative Gravity

We assume that all matter tensor fields are vector valued and the covariant
derivatives are defined with the noncommutative connection
$\mathcal{A}^{\alpha}{}_{\mu\nu}=\Gamma^{\alpha}{}_{\mu\nu}I+\theta^{\alpha}{}_{\mu\nu}$
(which is defined in terms of the noncommutative metric
$a^{\mu\nu}=g^{\mu\nu}I+h^{\mu\nu}$). So, the action of the matter fields
$\psi$ is generalized as follows

$S_{\rm mat}(\psi,g)\mapsto S_{\rm mat}(\psi,a)=S_{\rm mat}(\psi,g,h).$ (6.1)

The noncommutative MOND gravitational action is defined by

$S_{\rm ncgrav}(a)=\frac{1}{16\pi
G}\int\limits_{M}dx\;\frac{1}{N}\mathrm{tr\,}\rho\Biggl{\\{}\mathcal{R}+\lambda\mathcal{F}-2a^{2}_{0}f\left(\frac{\Theta}{a_{0}^{2}}\right)\Biggr{\\}},$
(6.2)

where $G$ is the Newton constant, $\lambda$ is a coupling constant and $f(x)$
is the MOND function (2.20). This function is supposed to have the following
properties

$f(x)\sim\left\\{\begin{array}[]{ll}\displaystyle\frac{\Lambda}{a^{2}_{0}}&\mbox{
for }x>>1,\\\\[15.0pt] \displaystyle\frac{4}{3}x^{3/4}&\mbox{ for
}x<<1,\end{array}\right.$ (6.3)

where $\Lambda$ is the cosmological constant and $a_{0}$ is the MOND parameter
(2.5).

This is a generalization of the bimetric MOND theory (BIMOND) (2.19). In
particular, we consider the Abelian case, when the matrix $a^{\mu\nu}$ is two-
dimensional, $N=2$, and diagonal, that is,

$\displaystyle a^{\mu\nu}$ $\displaystyle=$ $\displaystyle
g^{\mu\nu}I+h^{\mu\nu}T,$ (6.4)

where $g^{\mu\nu}$ is the (pseudo)-Riemannian metric, $h^{\mu\nu}$ is a
symmetric tensor and $T$ is the diagonal matrix

$T=\left(\begin{array}[]{cc}1&0\\\ 0&-1\\\ \end{array}\right).$ (6.5)

This can be rewritten in the form

$\displaystyle a^{\mu\nu}$ $\displaystyle=$ $\displaystyle
g^{\mu\nu}_{+}P_{+}+g^{\mu\nu}_{-}P_{-},$ (6.6)

where $g^{\mu\nu}_{\pm}$ are the symmetric tensors defined by

$g^{\mu\nu}_{+}=g^{\mu\nu}+h^{\mu\nu},\qquad
g^{\mu\nu}_{-}=g^{\mu\nu}-h^{\mu\nu},$ (6.7)

and $P_{+},P_{-}$ are the projections

$P_{+}=\left(\begin{array}[]{cc}1&0\\\ 0&0\\\ \end{array}\right),\qquad
P_{-}=\left(\begin{array}[]{cc}0&0\\\ 0&1\\\ \end{array}\right).$ (6.8)

In this case all geometric quantities diagonalize

$\displaystyle b_{\mu\nu}$ $\displaystyle=$ $\displaystyle
g^{+}_{\mu\nu}P_{+}+g^{-}_{\mu\nu}P_{-},$ (6.9) $\displaystyle\rho$
$\displaystyle=$ $\displaystyle g^{1/2}_{+}P_{+}+g^{1/2}_{-}P_{-},$ (6.10)
$\displaystyle\mathcal{A}^{\alpha}{}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\Gamma_{+}^{\alpha}{}_{\mu\nu}P_{+}+\Gamma_{-}^{\alpha}{}_{\mu\nu}P_{-},$
(6.11) $\displaystyle\mathcal{R}^{\alpha}{}_{\beta\mu\nu}$ $\displaystyle=$
$\displaystyle
R_{+}^{\alpha}{}_{\beta\mu\nu}P_{+}+R_{-}^{\alpha}{}_{\beta\mu\nu}P_{-}.$
(6.12)

The metric $g^{\mu\nu}$ is related to the tensors $g^{\mu\nu}_{\pm}$ by

$g^{\mu\nu}=\frac{1}{2}(g^{\mu\nu}_{+}+g^{\mu\nu}_{-}),$ (6.13)

and the tensor $\theta^{\alpha}{}_{\mu\nu}$ is then defined by

$\theta^{\alpha}{}_{\mu\nu}=\theta_{+}^{\alpha}{}_{\mu\nu}P_{+}+\theta_{-}^{\alpha}{}_{\mu\nu}P_{-},$
(6.14)

where

$\theta_{\pm}^{\alpha}{}_{\mu\nu}=\Gamma_{\pm}^{\alpha}{}_{\mu\nu}-\Gamma^{\alpha}{}_{\mu\nu},$
(6.15)

$\Gamma^{\alpha}{}_{\mu\nu}$ are Christoffel symbols for the metric
$g_{\mu\nu}$ and $\Gamma^{\alpha}_{\pm}{}_{\mu\nu}$ are Christoffel symbols
for the metric $g^{\pm}_{\mu\nu}$. Therefore, we obtain the quantity (4.92)

$\Theta=\Theta_{+}P_{+}+\Theta_{-}P_{-},$ (6.16)

where

$\Theta_{\pm}=g_{\pm}^{\beta\nu}\left(\theta^{\pm}_{\gamma}\theta_{\pm}^{\gamma}{}_{\nu\beta}-\theta_{\pm}^{\mu}{}_{\gamma\nu}\theta_{\pm}^{\gamma}{}_{\mu\beta}\right).$
(6.17)

Therefore, the action (6.2) becomes

$\displaystyle S_{\rm ncgrav}(a)$ $\displaystyle=$
$\displaystyle\frac{1}{16\pi
G}\int\limits_{M}dx\;\frac{1}{2}\Biggl{\\{}g_{+}^{1/2}R_{+}+g^{1/2}_{-}R_{-}{}$
(6.18)
$\displaystyle-2a^{2}_{0}g_{+}^{1/2}f\left(\frac{\Theta_{+}}{a_{0}^{2}}\right)-2a^{2}_{0}g_{-}^{1/2}f\left(\frac{\Theta_{-}}{a_{0}^{2}}\right)\Biggr{\\}},$

where $g_{\pm}=-\det g^{\pm}_{\mu\nu}$ and $R_{\pm}$ are the scalar curvatures
of the metrics $g^{\pm}_{\mu\nu}$. This is a generalization of the
relativistic bimetric BIMOND action (2.19) that should be carefully studied.

##  7 Static Weak Field Limit

In the static weak field limit we have

$\displaystyle g^{\mu\nu}_{\pm}$ $\displaystyle=$
$\displaystyle\eta^{\mu\nu}+2\varphi_{\pm}\delta^{\mu\nu},$ (7.1)
$\displaystyle g^{\mu\nu}$ $\displaystyle=$
$\displaystyle\eta^{\mu\nu}+2\varphi\delta^{\mu\nu},$ (7.2)

where $\eta^{\mu\nu}=\mathrm{diag\,}(-1,1,1,1)$ is the Minkowski metric, and

$\varphi=\frac{1}{2}(\varphi_{+}+\varphi_{-}),$ (7.3)

so that, up to the second order in $\varphi$

$g_{\pm}^{1/2}=1-2\varphi_{\pm}.$ (7.4)

For each metric we compute the Christoffel symbols in the first order

$\displaystyle\Gamma^{0}{}_{00}$ $\displaystyle=$
$\displaystyle\Gamma^{0}{}_{ij}=\Gamma^{i}{}_{j0}=0,$ (7.5)
$\displaystyle\Gamma^{0}{}_{0i}$ $\displaystyle=$
$\displaystyle\Gamma^{i}{}_{00}=\varphi_{,i},$ (7.6)
$\displaystyle\Gamma^{i}{}_{jk}$ $\displaystyle=$
$\displaystyle\varphi_{,i}\delta_{jk}-\varphi_{,j}\delta_{ik}-\varphi_{,k}\delta_{ij}.$
(7.7)

Here, as usual, comma denotes the partial derivative,
$\varphi_{,i}=\partial_{i}\varphi$.

Next, we compute the tensor $\theta^{\mu}{}_{\alpha\beta}$, in the first
order,

$\displaystyle\theta_{\pm}^{0}{}_{00}$ $\displaystyle=$
$\displaystyle\theta_{\pm}^{0}{}_{ij}=\theta_{\pm}^{i}{}_{j0}=0,$ (7.8)
$\displaystyle\theta_{\pm}^{0}{}_{0i}$ $\displaystyle=$
$\displaystyle\theta_{+}^{i}{}_{00}=\pm\frac{1}{2}\psi_{,i},$ (7.9)
$\displaystyle\theta_{\pm}^{i}{}_{jk}$ $\displaystyle=$
$\displaystyle\pm\frac{1}{2}\left(\psi_{,i}\delta_{jk}-\psi_{,j}\delta_{ik}-\psi_{,k}\delta_{ij}\right),$
(7.10)

where

$\psi=\varphi_{+}-\varphi_{-}.$ (7.11)

Then we compute the tensor $\theta_{\mu}$; we get

$\displaystyle\theta^{\pm}_{0}$ $\displaystyle=$ $\displaystyle 0,$ (7.12)
$\displaystyle\theta^{\pm}_{i}$ $\displaystyle=$
$\displaystyle\mp\frac{3}{2}\psi_{,i}.$ (7.13)

Finally, we compute the scalars $\Theta_{\pm}$ in the second order,

$\Theta_{\pm}=\frac{1}{2}\psi_{,i}\psi_{,i}=\frac{1}{2}|\nabla\psi|^{2}.$
(7.14)

We also need the Einstein-Hilbert action in the second order (after
integration by parts)

$\int dx\;g^{1/2}R=-\int dx\;2|\nabla\varphi|^{2}$ (7.15)

By using the above formulas we obtain the static weak field action

$S_{\rm ncgrav}(a)=-\frac{1}{8\pi
G}\int\limits_{M}dx\;\Biggl{\\{}|\nabla\varphi_{+}|^{2}+|\nabla\varphi_{-}|^{2}+a^{2}_{0}f\left(\frac{|\nabla\psi|^{2}}{2a_{0}^{2}}\right)\Biggr{\\}},$
(7.16)

which is very similar to the QUMOND model (2.15). Assuming the matter action

$S_{\rm
matter}=\int\limits_{M}dx\;\left\\{-\rho_{+}\varphi_{+}-\rho_{-}\varphi_{-}\right\\},$
(7.17)

we obtain two coupled equations of motion

$\displaystyle\Delta\varphi_{+}$ $\displaystyle=$ $\displaystyle 4\pi
G\rho_{+}-\frac{1}{2}\nabla_{i}\left[f^{\prime}\left(\frac{|\nabla\psi|^{2}}{2a_{0}^{2}}\right)\nabla_{i}\psi\right],$
(7.18) $\displaystyle\Delta\varphi_{-}$ $\displaystyle=$ $\displaystyle 4\pi
G\rho_{-}+\frac{1}{2}\nabla_{i}\left[f^{\prime}\left(\frac{|\nabla\psi|^{2}}{2a_{0}^{2}}\right)\nabla_{i}\psi\right].$
(7.19)

By subtracting and summing these equations we get

$\displaystyle\Delta\varphi$ $\displaystyle=$ $\displaystyle 4\pi G\rho,$
(7.20) $\displaystyle\Delta\psi$ $\displaystyle=$ $\displaystyle 4\pi
G(\rho_{+}-\rho_{-})-\nabla_{i}\left[f^{\prime}\left(\frac{|\nabla\psi|^{2}}{2a_{0}^{2}}\right)\nabla_{i}\psi\right],$
(7.21)

where

$\rho=\frac{1}{2}(\rho_{+}+\rho_{-}).$ (7.22)

For a symmetric matter, $\rho_{+}=\rho_{-}=\rho$, these equations coincide
with the QUMOND model (2.15).

##  8 Conclusions

We have shown that the non-commutative framework of Matrix General Relativity
can be used to reformulate MOND. In particular, MOND appears to be a special
case of Matrix General Relativity in which the non-commutative metric tensor
is two-dimensional and diagonal. Further investigation is needed to tackle the
more general case in which there are off-diagonal terms and/or the rank of the
non-commutative metric tensor is greater than two.

## References

  * [1] I. G. Avramidi, A noncommutative deformation of general relativity, Phys. Lett. B, 576 (2003) 195-198;
  * [2] I. G. Avramidi, Matrix general relativity: a new look at old problems, Class. Quant. Grav., 21 (2004) 103-120;
  * [3] I. G. Avramidi, Gauged gravity via spectral asymptotics of non-Laplace type operators, J. High Energy Phys., 07 (2004) 030;
  * [4] I. G. Avramidi, Dirac operator in matrix geometry, Int. J. Geom. Methods in Mod. Phys., 2 (2005) 227-264;
  * [5] N. A. Bahcall, L. M. Lubin and V. Dorman, Where is the dark matter?, Astrophys. J., 447 (1995), L81–L84.
  * [6] E. F. Bell, D. H. McIntosh, N. Katz and M. D. Weinberg, A first estimate of the baryonic mass function of galaxies, Astrophys. J., 585 (2003), L117–L120.
  * [7] J. Bekenstein and M. Milgrom, Does the missing mass problem signal the breakdown of Newtonian gravity?, Astrophys. J., 286, (1984) 7–14
  * [8] A. Bosma, 21-cm line studies of spiral galaxies. The distribution and kinematics of neutral hydrogen in spiral galaxies of various morphological types, Astron. J 86 (1981), 1825–1846.
  * [9] C. J. Copi, D. N. Schramm and M. S. Turner, Big-Bang nucleosynthesis and the baryon density of the Universe, Science, 267 (1995), 192-199.
  * [10] M. Davis, The cosmological matter density, Phys. Rep., 333 (2000), 147-165.
  * [11] B. Famaey and S. S. McGaugh, Modified Newtonian dynamics (MOND): observational phenomenology and relativistic extensions, Living Rev. Relativity, 15, (2012), 10
  * [12] S. Giodini, D. Pierini, A. Finoguenov, G. W. Pratt and the COSMOS Collaboration, Stellar and total baryon mass fractions in groups and clusters since redshift 1, Astrophys. J., 703 (2009), 982-993.
  * [13] F. D. Kahn and L. Woltjer, Intergalactic matter and the galaxy, Astrophys. J. 130 (1959), 705.
  * [14] J. S. Kalirai1, R. L. Beaton, M. C. Geha, K. M. Gilbert, P. Guhathakurta, E. N. Kirby, S. R. Majewski, J. C. Ostheimer, R. J. Patterson, and J. Wolf, The SPLASH Survey: Internal kinematics, chemical abundances, and masses of the Andromeda I, II, III, VII, X, and XIV dwarf spheroidal galaxies, Astrophys. J. 711 (2010), 671–692.
  * [15] J. P. Kneib, R. S. Ellis, I. Smail, W. J. Couch and R. M. Sharples, Hubble space telescope observations of the lensing cluster Abell 2218, Astrophys. J., 471 (1996), 643.
  * [16] M. Milgrom, A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis, Astrophys. J., 270, (1983) 365–370.
  * [17] M. Milgrom, Bimetric MOND gravity, Phys. Rev. D, 80, (2009) 123536.
  * [18] M. Milgrom, MD or DM? Modified dynamics at low accelerations vs dark matter, in Quarks, Strings and the Cosmos - H´ector Rubinstein Memorial Symposium, Stockholm, Sweden, August 9 – 11, 2010, Proceedings of Science PoS(HRMS)033, (SISSA, Trieste, 2010)
  * [19] M. Milgrom, Quasi-linear formulation of MOND, Mon. Not. R. Astron. Soc., 403, (2010) 886–895.
  * [20] M. Milgrom, Broader view of bimetric MOND, Phys. Rev. D 106 (2022), 084010
  * [21] J. H. Oort, The force exerted by the stellar system in the direction perpendicular to the galactic plane and some related problems, Bull. Astron. Inst. Netherlands 6 (1932), 249–287.
  * [22] J.P. Ostriker and P.J.E. Peebles, A Numerical study of the stability of flattened galaxies: or, can cold galaxies survive?, Astrophys. J. 186 (1973), 467–480.
  * [23] V. C. Rubin, J. W. K. Ford, N. Thonnard and D. Burstein, Rotational properties of 23 SB galaxies, Astrophys. J. 261 (1982), 439–456.
  * [24] J. Silk, Fluctuations in the primordial fireball, Nature, 215 (1967), 1155–1156.
  * [25] J.D. Simon and M. Geha, The kinematics of the ultra-faint Milky Way satellites: solving the missing satellite problem, Astrophys. J. 670 (2007), 313-331.
  * [26] M. G. Walker, M. Mateo and E. W. Olszewski, Stellar velocities in the Carina, Fornax, Sculptor, and Sextans dSph Galaxies: data from the Magellan/MMFS Survey, Astron. J. 137 (2009), 3100–3108.
  * [27] T. P. Walker, G. Steigman, H. S. Kang, D. M. Schramm and K. A. Olive, Primordial nucleosynthesis redux, Astrophys. J., 376 (1991), 51-69.
  * [28] F. Zwicky, Die Rotverschiebung von extragalaktischen Nebeln, Helv. Phys. Acta 6 (1932), 110–127.