On matrix cosmology

 The BFSS/BMN matrix model (also called M-(atrix) theory) is a matrix quantum mechanics characetrized essentially by the Yang-Mills term in $d$ dimensions ($d=9$)  \cite{deWit:1988wriF,Banks:1996vhF,Berenstein:2002jqF}. The Euclidean action  of the bosonic truncation is given in terms of $N\times N$ Hermitian matrices $\Phi_i$ by
\begin{eqnarray}
S_{\rm Matrix}&=&\frac{1}{g^2}\int_0^{\beta}dt{\rm Tr}\bigg[\frac{1}{2}(D_t\Phi_i)^2-\frac{1}{4}[\Phi_i,\Phi_j]^2+\frac{\mu^2}{2}{\rm Tr}\Phi_a^2+\frac{2i\alpha}{3}\epsilon_{abc}\Phi_a\Phi_b\Phi_c\nonumber\\&+&\frac{\nu^2}{8}{\rm Tr}\Phi_A^2\bigg].\label{BFSS}
\end{eqnarray}
The index $i$ runs from $1$ to $d$, the index $a$ runs from $1$ to $3$, and the index $A$ runs from $4$ to $d$ where $d=9$ is the maximal dimension allowed for supersymmetry. The requirement of supersymmetry also imposes the relation $\mu=\nu=2\alpha/3$ between the parameters of the supersymmetric deformation operators. The parameter $\mu$ breakes rotational symmetry as $SO(9)\longrightarrow SO(3)\times SO(6)$. The value $\mu=0$ is the BFSS model whereas the one-parameter deformation $\mu\neq 0$ is the BMN model.

From the perspective of string theory the $N\times N$ Hermitian matrices $\Phi_i$ represent the $9$ non-commuting position coordinates of N D0-branes which are particles with Dirichlet boundary conditions, i.e. they are the locations of the ends of fundamental strings.

The vacua (zero-energy configurations) of the above matrix model are given by diagonal matrices (for $\mu=0$), the zero vacuum $X_i=0$ (for both $\mu=0$ and $\mu\neq 0$), and $SU(2)$ representations (for $\mu\neq 0$).

This model has a $U(N)$ gauge symmetry, i.e. $D_t=\partial_t-i[A_t,..]$ where $A_t$ is the $U(N)$ gauge field representing the $10$ dimensions. In the gauge $A_t=0$ we have to impose a constraint on the Hilbert space given by the Gauss's law which restricts quantum states to be singlets under $U(N)$. This constraint is given explicitly by
\begin{eqnarray}
[D_t\Phi_i,\Phi_i]\equiv 0.
\end{eqnarray}
Following the analysis of \cite{Maldacena:2018vsr} we will simply set $A_0=0$ without imposing the Gauss's constraint as it seems that gauge invariance does not impact significantly the gravity dual. Indeed, the gravity dual for the ungauged matrix model is found to be essentially the same as the one for the original gauged matrix model.

For  $d=9$ and $\mu=0$ the gravity dual of the above matrix model is a $10-$dimensional type IIA charged black hole  given by the near-horizon geometry \cite{Horowitz:1991cd}
\begin{eqnarray}
  ds^2=\alpha^{\prime}\bigg(-\frac{r^{7/2}f_0(r)}{\sqrt{\lambda d_0}}dt^2+\frac{\sqrt{\lambda d_0}}{r^{7/2}f_0(r)}dr^2+\frac{\sqrt{\lambda d_0}}{r^{3/2}}d\Omega_8^2\bigg).
\end{eqnarray}
In this equation $\alpha^{'}=l_s^2$ is the inverse of the string tension,  $\lambda$ is the t'Hooft coupling given by $\lambda=g^2N$ and $f_0$ is the Schwarzschild radial coefficient in $d=9$ dimension given explicitly by $f_0(r)=1-r_0^7/r^7$.

The Hamiltonian in real time corresponding to the Euclidean action (\ref{BFSS}) is given by
\begin{eqnarray}
  H_{\rm Matrix}=Tr\bigg[\frac{1}{2}P_iP_i-\frac{1}{4g^2}[\Phi_i,\Phi_j]^2+\frac{\mu^2}{2g^2}{\rm Tr}\Phi_a^2+\frac{2i\alpha}{3g^2}\epsilon_{abc}\Phi_a\Phi_b\Phi_c+\frac{\nu^2}{8g^2}{\rm Tr}\Phi_A^2\bigg].\label{Ham}
    \end{eqnarray}
The Hamilton equations of motion derived from the Hamiltonian (\ref{Ham}) are given straightforwardly by
\begin{eqnarray}
  \frac{d({P}_i)_{\alpha\beta}}{dt}=-\frac{\partial H_{\rm Matrix}}{\partial (\Phi_i)_{\beta\alpha}}\equiv(F_i)_{\alpha\beta}\nonumber\\
   \frac{d({\Phi}_i)_{\alpha\beta}}{dt}=\frac{\partial H_{\rm Matrix}}{\partial (P_i)_{\beta\alpha}}\equiv ({P}_i)_{\alpha\beta}.
\end{eqnarray}
The force $F_i$ is given explicitly by
\begin{eqnarray}
  -F_i=\frac{1}{g^2}[\Phi_j,[\Phi_j,\Phi_i]]+\frac{\mu^2}{g^2}\Phi_a\delta_{i,a}+\frac{i\alpha}{g^2}\epsilon_{abc}[\Phi_b,\Phi_c]\delta_{i,a}+\frac{\nu^2}{g^2}\Phi_A\delta_{i,A}.
\end{eqnarray}
The second and third terms are non-zero only if $i$ equals $1-3$ whereas the fourth term is non-zero if $i$ equals $4-9$. Equivalently, the equations of motion read
\begin{eqnarray}
 \frac{d^2\Phi_i}{dt^2}+\frac{1}{g^2}[\Phi_j,[\Phi_j,\Phi_i]]=-\frac{\mu^2}{g^2}\Phi_a\delta_{i,a}-\frac{i\alpha}{g^2}\epsilon_{abc}[\Phi_b,\Phi_c]\delta_{i,a}-\frac{\nu^2}{g^2}\Phi_A\delta_{i,A}.\label{Heq}
\end{eqnarray}
Matrix cosmology starts like any other cosmology from the cosmological principle, i.e. with the homothetic ansatz \cite{Alvarez:1997fy,Hoppe:1997gr,Freedman:2004xg}
\begin{eqnarray}
 \Phi_a(t)=a(t)J_a~,~\Phi_A=b(t)J_A~,~[J_a,J_b]=i\epsilon_{abc}J_c~,~[J_A,J_B]=0~,~[J_a,J_A]=0.
\end{eqnarray}
Again the lowercase index $a$ runs over the values $1,2,3$ while the uppercase index $A$ runs over the values $4-9$. By substituting we get for  $i=a$ the two equations
\begin{eqnarray}
[J_b,[J_b,J_a]]+\theta J_a=0~,~\theta=-2.
\end{eqnarray}
\begin{eqnarray}
\ddot{a}-\Theta a^3=\Lambda_0 a+\lambda^{\prime} a^2~,~\Theta=\frac{\theta}{g^2}~,~\Lambda_0=-\frac{\mu^2}{g^2}~,~\lambda^{\prime}=\frac{2\alpha}{g^2}.
\end{eqnarray}
The first equation is the analogue of the equation governing central configurations in Newtonian cosmology, i.e. it represents essentially rotational invariance. See \cite{Freedman:2004xg} for a brief discussion and an adequate list of references. The second equation above is
precisely Raychaudhuri's equation which takes also the form
\begin{eqnarray}
  \frac{\ddot{a}}{a^3}-\frac{\Lambda_0}{a^2}=\frac{\lambda^{\prime}}{a}+\Theta.
\end{eqnarray}
This gives the scale factor of an expanding universe with a cosmological constant $\Lambda_0$ (a quadratic potential) together with two extra exotic sources of the energy density given by  $\lambda^{\prime}$ (a cubic potential, Chern-Simons term) and $\Theta$ (a quartic potential, Yang-Mills term) and. The Hamilton-Jacobi first integral of the above equation is immediately given by the Friedmann equation
\begin{eqnarray}
  \big(\frac{\dot{a}}{a}\big)^2+\frac{k}{a^2}=\Lambda_0+\frac{2}{3}\lambda^{'}a+\frac{1}{2}\Theta a^2.
\end{eqnarray}
The constant of integration $k$ plays the role of the spatial scalar curvature. The first term on the right-hand side is a constant and thus $\Lambda_0$ represents a cosmological constant (dark energy). However, due to the sign $\Lambda_0<0$ we have here the behavior of an anti-de Sitter space. The last term on the right-hand side is also negative ($\Theta<0$) and thus seems to be unstable which is certainly not correct as the original model is perfectly stable.  In fact, this term is also equivalent to a cosmological constant. Indeed, the Yang-Mills term at large number of dimensions $d\longrightarrow\infty$ can be replaced by a harmonic oscillator \cite{Mandal:2009vzF,Mandal:2011hbF}
\begin{eqnarray}
-\frac{1}{4}{\rm Tr}[\Phi_i,\Phi_j]^2\longrightarrow \frac{1}{2}m^2{\rm Tr}\Phi_i^2~,~m=d^{1/3}.
\end{eqnarray}
The above Friedmann equation becomes then
\begin{eqnarray}
  \big(\frac{\dot{a}}{a}\big)^2+\frac{k}{a^2}=\Lambda+\frac{2}{3}\lambda^{'}a~,~\Lambda=-\frac{\mu^2+m^2}{g^2}<0.
\end{eqnarray}
The effective cosmological constant is again negative. For zero spatial curvature, i.e. for $k=0$ the solution of this equation takes the form (with $A=a^{5/2}$ and $\epsilon={3\Lambda}/{2\lambda^{'}}$)
\begin{eqnarray}
\int \frac{dA}{(1+\frac{\epsilon}{A^{2/5}})^{1/2}}=\frac{5}{2}\sqrt{\frac{2\lambda^{'}}{3}}t.
\end{eqnarray}
We seek an expanding solution for which $a(t)\longrightarrow \infty$  as a power law when $t\longrightarrow \infty$. This allows us to expand the root square in the above equation and we end up with the following leading behavior (with $\lambda\equiv g^2N= 1$)
\begin{eqnarray}
a(t)=(\frac{25\lambda^{'}}{6})^{1/5}t^{2/5}=(\frac{25\alpha}{3g^2})^{1/5}t^{2/5}=(\frac{25N\alpha}{3})^{1/5}t^{2/5}.
\end{eqnarray}
This approximation is also valid in the regime
\begin{eqnarray}
\alpha>>\frac{3}{4}(\mu^2+m^2).
\end{eqnarray}
Finally, by substituting $i=A$ in (\ref{Heq}) we get instead the two equations
\begin{eqnarray}
[J_B,[J_B,J_A]]+\hat{\Theta} J_A=0~,~.\ddot{b}-\hat{\Theta} b^3=\hat{\lambda} b~,~\hat{\lambda}=-\frac{\nu^2}{g^2}.
\end{eqnarray}
We have clearly $\hat{\Theta}=0$ and $\hat{\lambda}<0$ and hence $b=b(t)$ is an oscillating function in time. In other words, the matrix coordinates $X_4$,...,$X_9$ are oscillating directions.

In conclusion, we have three (or in fact two) expanding directions corresponding to those of the emergent geometry in the BFSS/IKKT mode and six oscillatory directions corresponding to the remaining flat directions.