Periodic analogues of the Kerr solutions:
a numerical study

The configurations that we are looking for have three degrees of freedom that we will take to be the area $A$, the angular momentum $J$, and the period $L$, with $L$ linked to the physical separation between consecutive black holes. These parameters need to be incorporated into the boundary conditions of the main equations. The stationary and axisymmetric metrics are assumed to be in Weyl-Papapetrou form and are recalled in (2), (see for instance [12]). In this expression $\rho\in\mathbb{R}^{+}$ and $z\in\mathbb{R}$ are the Weyl-Papapetrou cylindrical coordinates and play a role in what follow. For Weyl-Papapetrou metrics, the Einstein equations reduce to find an axisymmetric harmonic map $(\omega(\rho,z),\eta(\rho,z)):\mathbb{R}^{3}\rightarrow\mathbb{H}^{2}$ from $\mathbb{R}^{3}$ into the hyperbolic space $\mathbb{H}^{2}$. Here $\eta$ is the norm squared of the axisymmetric Killing field and $\omega$ is the twist potential. All the metric components can be recovered from $\eta$ and $\omega$ after line-integrations. The full set of reduced equations are presented in section 2.2. The harmonic map equations are (3)-(4) and the metric components are (5) after the line integrations of (6) and (7).

As we look for metrics periodic in $z$, we restrict the analysis to the region $\{\rho\geq 0,-L/2\leq z\leq L/2\}$ keeping appropriate periodicity conditions on the top and bottom lines $\{\rho\geq 0,\ z=\pm L/2\}$. Here $L$ is the period mentioned earlier and is a free parameter. The harmonic map equations for $\omega$ and $\eta$ need to be supplied with boundary data on the border $\{\rho=0,-L/2\leq z\leq L/2\}$. This boundary contains the horizon $\mathcal{H}=\{\rho=0,\ -m\leq z\leq m\}$, and the two components of the axis ${\mathcal{A}}$ which is the complement of $\mathcal{H}$. The boundary conditions arise from two sides: from the natural conditions that $\eta$ and $\omega$ must verify on horizons and axes, and from fixing the values of the area $A$ and of the angular momentum $J$. The angular momentum $J$ is set by specifying Dirichlet data for $\omega$ on the axis, whereas the area $A$ is set by specifying the limit values of $\eta/\rho^{2}$ at the poles $z=\pm m$. The fact that one can incorporate $A$ and $J$ into the boundary conditions in this simple way is a great advantage of the formulation. The parameter $m$ that may seem to be a free parameter, is equal to the surface gravity $\kappa$, ($\kappa^{2}:=-\langle\nabla\partial_{t},\nabla\partial_{t}\rangle/2$), of the horizons times $A/4\pi$ and therefore can be fixed using the freedom in the definition of the stationary Killing field $\partial_{t}$, or the definition of time in (2). In this article we chose $\kappa=\kappa(A,J)$ equal to the surface gravity of the Kerr black holes for the given $A$ and $J$.

The full set of boundary conditions on $\{\rho=0,-L/2\leq z\leq L/2\}$ is discussed in section 2.4.

In practice, the harmonic map equations need to be solved numerically on a finite rectangle $[0,\rho_{\rm\small{MAX}}]\times[-L/2,L/2]$, and this adds the extra difficulty of finding natural boundary conditions also at $\rho=\rho_{\rm MAX}$, for $\omega$ and for $\eta$. On physical grounds one expects the solutions to become asymptotically independent of $z$ and indeed approaching the Lewis models of Stockum’s rotating cylinders that we recall in section 2.3. But the problem is that one does not know a priori which Lewis solution shows up for the given $A$, $J$ and $L$. If that information were known, then one could easily supply appropriate boundary conditions at $\rho_{\rm MAX}$. Now, all Lewis solutions have $\omega=wz$, so we make $\omega(\rho_{\rm MAX},z)=wz$. For $\eta$ however not such single choice is possible. We set a Neumann type of boundary condition for $\eta$ and for that we use that, on actual solutions, the Komar mass expression $M(\rho)$, in (34), is constant, so that $M(\rho_{\rm MAX})=M(\rho)$ easily relates $\partial_{\rho}\eta(\rho_{\rm MAX})$ to $\eta$ and $\omega$ at any $\rho<\rho_{\rm MAX}$. Then, to define the condition for $\eta$ we make use of this relation, equating $\partial_{\rho}\eta(\rho_{\rm MAX})$ to the average of the Komar mass expression in the bulk $0<\rho<\rho_{\rm MAX}$. We discuss this condition in section 2.4, and state it in (34). This peculiar Neumann condition for $\eta$ at $\rho_{\rm MAX}$ is what gave us the best numerical results in terms of the speed and the stability of our code.

To find solutions for the harmonic map system, we run a harmonic map heat flow with certain initial data, and look for the stationary solutions in the long time. We can show that the numerical solutions for two different values of $\rho$, say $1\ll\rho_{\rm MAX}^{1}<\rho_{\rm MAX}^{2}$, but the same values of $A$, $J$ and $L$, overlap on the smaller rectangle $[0,\rho^{1}_{\rm MAX}]\times[-L/2,L/2]$. This shows that the solutions constructed indeed depend on $A$, $J$ and $L$ only and the boundary condition at $\rho_{MAX}$ does not introduce any spurious new degree of freedom, but just drives the harmonic map heat flow to settle with the right asymptotic for the given $A$, $J$ and $L$.

The black hole configurations that we are looking for will be in Weyl-Papapetrou form,

$$g=-Vdt^{2}+2Wdtd\phi+\eta d\phi^{2}+\frac{e^{2\gamma}}{\eta}(d\rho^{2}+dz^{2}),$$

(2)

where $(t,\rho,z,\phi)\in\mathbb{R}\times\mathbb{R}^{+}\times\mathbb{R}\times S^{1}$ and where the components $V,W,\eta$ and $\gamma$ depend only on $(\rho,z)$. We require $V,W,\eta$ and $\gamma$ to be $z$-periodic with period $L$, and, to prevent struts on the axis, we demand in addition $V,W,\eta$ and $\gamma$ to be symmetric with respect to the reflection $z\rightarrow-z+L$.

If we denote by $\omega$ the twist potential of $\partial_{\varphi}$, then $\eta(\rho,z)$ and $\omega(\rho,z)$ satisfy the harmonic map equations (Ernst equations),

	$\displaystyle\Delta_{\rho}\eta$	$\displaystyle=\frac{|\partial\eta|^{2}-|\partial\omega|^{2}}{\eta},$		(3)
	$\displaystyle\Delta_{\rho}\omega$	$\displaystyle=2\langle\partial\omega,\frac{\partial\eta}{\eta}\rangle,$		(4)

where the inner product $\langle\ ,\ \rangle$ is with respect to the flat metric $d\rho^{2}+dz^{2}$, and $\Delta_{\rho}=\partial^{2}_{\rho}+\frac{1}{\rho}\partial_{\rho}+\partial^{2}_{z}$ is the $3$-dim Laplacian in cylindrical coordinates $(\rho,z,\phi)$ under axial symmetry. The metric components $W$ and $V$ follow from the identities,

$$W=\eta\Omega,\quad V=\frac{\rho^{2}-W^{2}}{\eta},$$

(5)

after finding the angular velocity function $\Omega$ by line integration of,

$\displaystyle\Omega_{z}=\rho\frac{\omega_{\rho}}{\eta^{2}},\quad\Omega_{\rho}=-\rho\frac{\omega_{z}}{\eta^{2}}.$

(6)

The exponent $\gamma$ is found after line integration of,

$\displaystyle\gamma_{z}=\frac{\rho}{2\eta^{2}}(\eta_{\rho}\eta_{z}+\omega_{\rho}\omega_{z}),\quad\gamma_{\rho}=\frac{\rho}{4\eta^{2}}(\eta_{\rho}^{2}-\eta_{z}^{2}+\omega_{\rho}^{2}-\omega_{z}^{2}).$

(7)

Thus, the Einstein equations are solved in a ladder-like fashion: first solve the harmonic map equations (3)-(4), then solve the quadratures (6) and (7), and finally get from them the other coefficients $V$ and $W$ of the metric (2).

The Lewis solutions [10] (see also [11]) are cylindrically symmetric (i.e. independent on $\varphi$ and $z$) rotating stationary vacuum solutions. Some of the solutions extend to infinity and some do not. The possible forms for $\eta$ are the following,

	$\displaystyle{\rm(I\pm)}:\quad$	$\displaystyle\eta=\rho\frac{|w|}{a}\sin(\pm a\ln(\rho)+b),\quad a>0,\ b\in\mathbb{R},$		(8)
	$\displaystyle{\rm(II\pm)}:\quad$	$\displaystyle\eta=\rho|w|(\pm\ln(\rho)+b),\quad b\in\mathbb{R},$		(9)
	$\displaystyle{\rm(III\pm)}:\quad$	$\displaystyle\eta=\rho\frac{|w|}{a}\sinh(\pm a\ln(\rho)+b),\quad a>0,\ b\in\mathbb{R},$		(10)

and the twist potential is always,

$$\omega=wz,\quad w\neq 0.$$

(11)

The solutions extending to infinity are (9) and (10) with the $+$ sign (note that they are positive after some $\rho$). For the case (10,+), the metric components $V$, $W$, $\eta$ and $e^{2\gamma}/\eta$ get the form,

	$$\displaystyle V=\frac{2a}{|w|}e^{-b}\rho^{1-a},\quad W=\text{s}(w)e^{-b}\rho^{1-a},$$		(12)
	$$\displaystyle\eta=\frac{|w|}{2a}(e^{b}\rho^{1+a}-e^{-b}\rho^{1-a}),\quad\frac{e^{2\gamma}}{\eta}=c\rho^{(a^{2}-1)/2},$$		(13)

where $\text{s}(w)=w/|w|$ is the sign of $w$ and $b\in\mathbb{R},c>0$ and $a>0$ are free parameters. The angular velocity function $\Omega$ is,

$$\Omega=\frac{2a}{w}\frac{e^{-b}\rho^{-a}}{(e^{b}\rho^{a}-e^{-b}\rho^{-a})},$$

(14)

and note that $\Omega$ is set to be zero at infinity. Not all of the solutions can model the asymptotic of coaxial arrays of black holes like the ones we are considering. Concretely, if $a\geq 1$ the lapse function for the hypersurface $t=0$ tends to zero as infinity something that can be ruled out by the maximum principle on the lapse equation taking into account that the lapse is zero on the horizons. As $a$ tends to zero, $\eta$ in (III+) degenerates into (II+). The asymptotic models are therefore (III+) with $0<a<1$ and (II+).

The metrics of the models (III+) are Kasner to leading order (except for the cross term $-2\text{s}(w)e^{-b}\rho^{1-a}dtd\phi$). Recalling that the Kasner metrics have the form,

$$g_{K}=-e^{-b}\rho^{\alpha}dt^{2}+e^{b}\rho^{2-\alpha}d\phi^{2}+c\rho^{\alpha^{2}/2-\alpha}\left(d\rho^{2}+dz^{2}\right),$$

(15)

we see that in case (III+) we have $\alpha=1-a$. Thus, the Kasner exponent $\alpha$ of the models (III+) is restricted to vary in $(0,1)$.

From now it will be more convenient to work with $\sigma=\ln\eta-2\ln\rho$ (this choice has some advantages, see [13, 14, 12]). In terms of $(\sigma,\omega)$ the harmonic map equations (3) and (4) become,

	$\displaystyle\Delta_{\rho}\sigma$	$\displaystyle=-\frac{e^{-2\sigma}|\partial\omega|^{2}}{\rho^{4}},$		(19)
	$\displaystyle\Delta_{\rho}\omega$	$\displaystyle=4\frac{\langle\partial\omega,\partial\rho\rangle}{\rho}+2\langle\partial\omega,\partial\sigma\rangle.$		(20)

The boundary conditions for this system will be the same as the boundary condition for the harmonic map heat flow, and are as follows.

The boundary conditions for $\omega$ are,

	$\displaystyle\omega\big{|}_{\mathcal{A}_{+}}=4J,\quad\omega\big{|}_{\mathcal{A}_{-}}=-4J,$		(21)
	$\displaystyle\partial_{\rho}\omega\big{|}_{\mathcal{H}}=0.$		(22)

The first, Dirichlet condition, comes from fixing the angular momentum of the horizons to $J$. The second, Neumann condition, can be obtained by demanding the spacetime smoothness of $\Omega=W/\eta$ that is linked to $\omega$ by (6).

The boundary conditions for $\sigma$ are,

	$\displaystyle\partial_{\rho}\sigma\big{|}_{\mathcal{A}\setminus\partial\mathcal{H}}=0,$		(23)
	$\displaystyle\partial_{\rho}(\sigma+2\ln\rho)\big{|}_{\mathcal{H}}=0,$		(24)
	$\displaystyle\lim_{(\rho,z)\rightarrow(0,\pm m)}\frac{\sigma}{\sigma_{0}}=1,$		(25)

where $\sigma_{0}$ is the reference Kerr solution given $A,J$ and $m(A,J)$, where

$$m=\sqrt{\frac{A}{16\pi}}\frac{1-\left(8\pi J/A\right)^{2}}{\sqrt{1+\left(8\pi J/A\right)^{2}}}.$$

(26)

This enforces the area of the solution to be actually $A$, noting that the area and temperature are given by

	$\displaystyle A$	$\displaystyle=2\pi\int_{-m}^{m}e^{\gamma}dz,$
	$\displaystyle\kappa$	$\displaystyle=e^{-\gamma}$

The boundary condition for $\sigma$ is delicate and we motivate it as follows. Given a solution to (3)-(4) one can get a Smarr type of expression for the Komar mass (per black hole) for the Killing field $\partial_{t}$ by integrating the Komar form on any 2-torus $t=0$ and $\rho={\rm const}$ of the spacetime, of course after identifying $z=-L/2$ to $z=L/2$. The expression is,

$$M(\rho)=\frac{1}{4}\int_{-L/2}^{L/2}(-\rho\partial_{\rho}\sigma+\Omega\partial_{z}\omega)dz,$$

(28)

which is in fact $\rho-$independent.

For the Lewis solutions (III+) we have $\Omega\rightarrow 0$ as $\rho\rightarrow\infty$, $\rho\partial_{\rho}\sigma\rightarrow a-1$ as $\rho\rightarrow\infty$, $\partial_{z}\omega=w$, so we obtain,

$$M=\frac{L}{4}(1-a).$$

(29)

Thus, the Kasner exponent $\alpha=1-a$ is equal to $4M/L$.

The Lewis solution (III+) on the finite rectangle $[0,\rho_{\rm MAX}]\times[-L/2,L/2]$ satisfies the decay

$$\partial_{\rho}\sigma\big{|}_{\rho_{\rm MAX}}=\frac{4M}{L\rho_{\rm MAX}}+o(1/\rho_{\rm MAX}),\quad\rho_{\rm MAX}\gg 1$$

(30)

where $M$ is the value of the Smarr mass, still unknown.

Now, since we are evolving a harmonic map heat flow, at any time $\tau\geq 0$ the integral in (28) associated to $(\sigma,\omega)$ at $\tau$ is not necessarily constant as a function of $\rho$, and even the function $\Omega$ is not well defined since the integrability conditions (6) do not necessarily hold. Therefore, we have to give a definition for $M(\rho,\tau)$ so that we can use the approximation (30) as boundary condition. By integrating by parts the last term in (28), and using (6), we have a candidate

$$M(\rho,\tau)=2\Omega(\rho,\tau)J\bigg{|}_{z=-L/2}-\frac{1}{4}\int_{\rho}(\rho\partial_{\rho}\sigma+\frac{\rho}{\eta^{2}}\partial_{\rho}\omega\omega)dz,$$

(31)

where $\Omega(\rho,\tau)$ is prescribed to be

$$\Omega(\rho,\tau)=\int_{\rho_{\rm MAX}}^{\rho}\frac{\rho^{\prime}}{\eta^{2}}\partial_{z}\omega d\rho^{\prime},$$

(32)

with the functions evaluated at $z=-L/2$. If $(\sigma,\omega)$ converges to a solution, then $M(\rho,\tau)$ converges to a constant function. Taking this into account, a natural dynamical condition arise if we take the $\rho$-average of $M(\rho,\tau)$ over the numerical domain, at each time step, and impose

$$M(\rho_{\rm MAX})=\overline{M}(\tau),$$

(33)

which in particular links the values of $M$ near the axis with those values far away.¹¹1In general we have a weighted average, which range from being completely uniform on the interval $(0,\rho_{\rm MAX})$ to associating more weight to the region near the axis and cutting off the region past certain $\rho_{0}<\rho_{\rm MAX}$. Different weights were explore in the numerical simulations. Then, our boundary condition for $\sigma$ at $\rho_{\rm MAX}$ is given by

$$\partial_{\rho}\sigma\big{|}_{\rho_{\rm MAX}}=-\frac{4\overline{M}(\tau)}{L\rho_{\rm MAX}}.$$

(34)

Observe that in this condition there is no explicit intervention of the total angular momentum, and there is also no reference to any specific asymptotic model. This is an important characteristic of this boundary condition: it will be the same for any periodic configuration, as long as the solutions become $z-$independent asymptotically.

It remains to fix the boundary conditions for the integration of the quadratures (6) and (16). For $\Omega$, we impose

$$\Omega\big{|}_{\rho_{\rm MAX}}=0,$$

(35)

which is the natural condition for the asymptotic models, and also consistent with the definition of $M$. The boundary condition on $q$ is just that it vanishes at just one point in $\mathcal{A}_{0}$. As it was presented above, the symmetries on $\sigma$ and $\omega$ and the integrability equations (16) imply this same condition holds at any other point of any axis component $\mathcal{A}_{i}$, preventing struts in between black holes.

We need some initial condition at $\tau=0$ of the heat flow, which we will call seed. In particular, it is desirable that the seed contains the prescribed singular behaviour of the solutions at the horizons, so that we can define a well-posed numerical problem without singularities. We decompose $\sigma$ and $\omega$ as follows: we split them as a sum of known solutions to the non-periodic problem plus a perturbation $\bar{\sigma},\bar{\omega}$. In the case a single horizon per period, this sum is constructed in the same fashion as the function $\sigma_{MKN}$ was defined: following [13], let $\sigma_{0}(\rho,z)$ and $\omega_{0}(\rho,z)$ be the solutions to the asymptotically flat Kerr black hole with momentum $J$ and horizon $\mathcal{H}_{0}$, and define

$$\begin{split}\sigma=&\sigma_{0}+\sigma_{r}+\bar{\sigma}\\ \omega=&\omega_{0}+\omega_{r}+\bar{\omega},\end{split}$$

(36)

where,

	$\displaystyle\sigma_{r}(\rho,z)$	$\displaystyle=C+\sum_{n=1}^{\infty}\left(\sigma_{0}(\rho,z-nL,J)+\sigma_{0}(\rho,z+nL,J)-\frac{4M}{nL}\right)$
	$\displaystyle\omega_{r}(\rho,z)$	$\displaystyle=\sum_{n=1}^{\infty}\left(\omega_{0}(\rho,z-nL,J)+\omega_{0}(\rho,z+nL-D/2,J)\right).,$

where $C$ is a constant such that $\sigma_{r}\mid_{\partial\mathcal{H}}=0$, that is, its value at the poles is zero. Here the constants $\frac{4M}{nL}$ are needed for the series to be convergent (since asymptotically, each term goes as $-\frac{2M}{\sqrt{(x-nL)^{2}+\rho^{2}}}$ and therefore we need to cancel out this divergent term, as in [6]). In our actual numerical calculations we use a cut-off value $N_{d}$ for $n$, which we take $N_{d}\gg 1$, which can be thought of as the number of “domains” we stack on both the top and below the central domain.

We expect $(\bar{\sigma},\bar{\omega})$ to be regular throughout the evolution. By inserting the decomposition (36) into (19) and (20), and using the fact that $(\sigma_{0},\omega_{0})$ is a solution, we obtain the evolution equations for $\bar{\sigma}$ and $\bar{\omega}$,

$$\begin{split}\dot{\bar{\sigma}}&=\Delta\bar{\sigma}+\Delta\sigma_{r}+\frac{e^{-2\sigma_{0}}|\partial\omega_{0}|}{\rho^{4}}\left(e^{-2(\bar{\sigma}+\sigma_{r})}-1\right)\\ &\quad+\frac{e^{-2(\sigma_{0}+\sigma_{r}+\bar{\sigma})}}{\rho^{4}}\left(|\partial\omega_{r}|^{2}+|\partial\bar{\omega}|^{2}+2\left(\partial_{i}\omega_{r}\partial^{i}\omega_{0}+\partial_{i}\omega_{r}\partial^{i}\bar{\omega}+\partial_{i}\bar{\omega}\partial^{i}\omega_{0}\right)\right),\end{split}$$

(37)

$$\begin{split}\dot{\bar{\omega}}&=\Delta\bar{\omega}+\Delta\omega_{r}-\frac{4}{\rho}\left(\partial_{\rho}\omega_{r}+\partial_{\rho}\bar{\omega}\right)-2\left(\partial_{i}\omega_{0}\partial^{i}\sigma_{r}+\partial_{i}\omega_{0}\partial^{i}\bar{\sigma}\right.\\ &\quad\left.+\partial_{i}\omega_{r}\partial^{i}\sigma_{0}+\partial_{i}\omega_{r}\partial^{i}\sigma_{r}+\partial_{i}\omega_{r}\partial^{i}\bar{\sigma}+\partial_{i}\bar{\omega}\partial^{i}\sigma_{0}+\partial_{i}\bar{\omega}\partial^{i}\sigma_{r}+\partial_{i}\bar{\omega}\partial^{i}\bar{\sigma}\right).\end{split}$$

(38)

These are the equations that we solve numerically, with the following boundary conditions, which can be read off from the conditions for the total functions $(\omega,\sigma)$,

$$\bar{\omega}(\rho_{\rm MAX})=0,\quad\partial_{\rho}\bar{\omega}\mid_{\mathcal{H}}=0,\quad\bar{\omega}\mid_{\mathcal{A}}=0,$$

(39)

since $\omega_{0}+\omega_{r}$ already satisfies the asymptotic linear behaviour for $\omega$, and

$$\partial_{\rho}\bar{\omega}\mid_{(\mathcal{A}\cup\mathcal{H})\setminus\partial\mathcal{H}}=0,\quad\bar{\sigma}\mid_{\partial\mathcal{H}}=0,$$

(40)

with the asymptotic condition

$$\partial_{\rho}\bar{\sigma}\big{|}_{\rho_{\rm MAX}}=-\frac{4\overline{M(\rho)}^{\rho}}{L\rho_{\rm MAX}}-\overline{\partial_{\rho}(\sigma_{0}+\sigma_{r})}^{z}\big{|}_{\rho_{\rm MAX}}$$

(41)

The last term is not $z-$independent, since the series are truncated, and therefore we take its average on $z$, denoting by $\overline{X}^{x}$ the average along $x$ coordinate of the variable $X$. We will call $\beta$ the dinamical quantity given by the right hand side of equations (41).