I'm solving the Einstein equation assuming a cylindrical symmetry and found something interesting which I never saw elsewhere. I now feel that I may have found the Melvin magnetic universe solution in another form but need a confirmation (sorry for the lengthy exposition below). Assume the following metric (a cylindrical universe): \begin{equation}\tag{1} ds^2 = dt^2 - a^2(t) g^2(r)(dx^2 + dy^2) - b^2(t) \, dz^2, \end{equation} where $a(t)$, $b(t)$ and $g(r)$ are unknown functions. Of course: $r =\sqrt{x^2 + y^2}$ and we could use $dx^2 + dy^2 = dr^2 + r^2 d\varphi^2$. I want an homogeneous and isotropic spatial section in the plane $x \, y$, so the function $g(r)$ isn't completely arbitrary. Calculating the Riemann tensor give these components (here, $i, j, k = 1, 2$ only, for the coordinates $x$ and $y$): \begin{align} R^0_{i 0 j} &= a \ddot{a} \, g^2 \, \delta_{ij}, \\[12pt] R^0_{303} &= b \ddot{b}, \\[12pt] R^i_{33j} &= -\, \frac{\dot{a} \dot{b} b}{a} \, \delta_{ij}, \\[12pt] R^i_{jkl} &= (\dot{a}^2 g^2 + \mathcal{Q}(r))(\delta_{ki} \, \delta_{lj} - \delta_{li} \, \delta_{kj}) + \text{anisotropic terms}. \end{align} The anisotropic terms cancel if \begin{equation}\tag{2} \frac{g^{\prime \prime}}{g} = \frac{g^{\prime}}{g \, r} + 2 \frac{g^{\prime 2}}{g^2}, \end{equation} which implies the solution \begin{equation}\tag{3} g(r) = \frac{A}{1 + B \, r^2}, \end{equation} where $A$ and $B$ are constants (they can be absorbed into the coordinates, so I'll take $A = 1$). Then \begin{equation}\tag{4} \mathcal{Q}(r) = 4 B \, g^2. \end{equation} Note that if $B$ is positive, then the whole $x \, y$ plane has a finite area: $\mathcal{A}_{x y} = \frac{\pi}{B}$ (that plane is compact)! The components of the Einstein tensor could be calculated without troubles (again, $i, j = 1, 2$ only): \begin{align} G_{00} &= -\, \Big( 2 \frac{\dot{a} \dot{b}}{a b} + \frac{\dot{a}^2}{a^2} + \frac{4 B}{a^2} \Big), \tag{5} \\[12pt] G_{ij} &= -\, \Big( \frac{\ddot{a}}{a} + \frac{\ddot{b}}{b} + \frac{\dot{a} \dot{b}}{a b} \Big) \, g_{ij}, \tag{6} \\[12pt] G_{33} &= -\, \Big( 2 \frac{\ddot{a}}{a} + \frac{\dot{a}^2}{a^2} + \frac{4 B}{a^2} \Big) \, g_{33}. \tag{7} \end{align} I then need an anisotropic (i.e. cylindrical) perfect fluid stress-tensor, of components \begin{equation}\tag{8} T_{\mu \nu} = (\rho + \sigma) \, u_{\mu} u_{\nu} + (p - \sigma) \, n_{\mu} n_{\nu} - \sigma \, g_{\mu \nu}, \end{equation} where $p \equiv p_{\parallel}$ and $\sigma \equiv p_{\perp}$ are the pressure along the $z$ axis and in the plane $x \, y$, respectively. The fluid four-velocity is $u^{\mu} = (1, 0, 0, 0)$ (such that $u_{\mu} \, u^{\mu} = 1$) and the symmetry axis is represented by the spacelike four-vector $n^{\mu} = (0, 0, 0, 1/b(t))$ (such that $n_{\mu}\, n^{\mu} = -1$). The local conservation of fluid $\nabla_{\mu} T^{\mu \nu} = 0$ implies the following: \begin{equation}\tag{9} \dot{\rho} + \Big( 2 \frac{\dot{a}}{a} + \frac{\dot{b}}{b} \Big) \rho + \frac{\dot{b}}{b} p + 2 \frac{\dot{a}}{a} \sigma = 0. \end{equation} Take note that a purely magnetic field could be cast into (8): \begin{equation} T_{\mu \nu} = F_{\mu \lambda} F^{\lambda}_{\; \; \nu}+ \frac{1}{4} g_{\mu \nu} F_{\lambda \kappa} F^{\lambda \kappa}, \end{equation} with $F_{12}$ the only component. We then have $p = -\, \rho$ (negative pressure = tension along the magnetic field lines) and $\sigma = + \rho$, such that the trace of (8) cancels: $g^{\mu \nu} T_{\mu \nu} = \rho - p - 2 \sigma = 0$. Equ. (9) then implies the standard energy density of radiation in cosmology: \begin{equation}\tag{10} \rho(t) = \frac{\text{cste}}{a^4(t)}. \end{equation} Putting all these into Einstein equation $G_{\mu \nu} = -\, 8 \pi G \, T_{\mu \nu}$ give the following equations to be solved: \begin{align} 2 \frac{\dot{a} \dot{b}}{a b} + \frac{\dot{a}^2}{a^2} + \frac{4 B}{a^2} = 8 \pi G \rho, \tag{11} \\[12pt] \frac{\ddot{a}}{a} + \frac{\ddot{b}}{b} + \frac{\dot{a} \dot{b}}{a b} = -\, 8 \pi G \sigma = -\, 8 \pi G \rho, \tag{12} \\[12pt] 2 \frac{\ddot{a}}{a} + \frac{\dot{a}^2}{a^2} + \frac{4 B}{a^2} = -\, 8 \pi G p = +\, 8 \pi G \rho. \tag{13} \end{align} If the constant $B$ is positive, solving these give $a(t) = \text{cste}$ (so $\rho = \text{cste}$ !) and $b(t) = \sin{\omega t}$ (for $0 < \omega t < \pi$, where $\omega = \sqrt{8 \pi G \rho}$).
I'm perplexed by this solution, which appears to describe a magnetic universe, homogeneous, isotropic and static in the plane $x \, y$, and expanding/contracting along the $z$ axis.
Is this the Melvin's magnetic universe in another form? Usually, Melvin's universe is described by the following static metric: \begin{equation}\tag{14} ds^2 = \Phi^2(r) (dt^2 - dr^2 - dz^2) - \frac{1}{\Phi^2(r)} \, r^2 d\varphi^2, \end{equation} where $\Phi(r) = 1 + C \, r^2$ (this function is similar to (3) above).
So how can I verify that metric (1) is equivalent to metric (14)? \begin{equation}\tag{15} ds^2 = d\tilde{t}^2 - \frac{1}{(1 + B \, \tilde{r}^2)^2}(d\tilde{r}^2 + \tilde{r}^2 \, d\varphi^2) - \sin^2{\omega \tilde{t}} \, d\tilde{z}^2. \end{equation}