text:
# Simulating Cosmic Microwave Background maps in multi-connected spaces

## I Introduction

Future cosmic microwave background (CMB) experiments such as the MAP  and later the Planck satellites  will provide full sky maps of CMB anisotropies (up to the galactic cut). These datasets offer the opportunity to probe the topological properties of our universe. A series of tests to detect the topology, including the use of the angular power spectrum , the distribution of matched patterns such as circles , correlation of antipodal points  and non Gaussianity  have been proposed (see Refs.  for reviews of CMB methods to search for the topology). The study of the detectability of the topology by any of these methods first requires simulating maps with the topological signature for a large set of topologies. These maps will allow one to test the detection methods, estimate their run time, and, once all sources of noise are added, determine to what extent a given method detects the topological signal (in the same spirit as the investigation of the “crystallographic” methods based on galaxy catalogs ).

In a simply connected111Geometers and cosmologists often refer to simple-connectedness as the “trivial topology”. However, trivial topology has a different meaning in the context of point-set topology: in that formalism, the trivial topology is the smallest topology on a set $`X`$, namely the one in which the only open sets are the empty set and the entire set $`X`$., spatially homogeneous and isotropic universe, the angular correlation function depends only on the angle between the two directions and the coefficients $`a_\mathrm{}m`$ of the decomposition of the temperature fluctuation in spherical harmonics, which are uncorrelated for different sets of $`\mathrm{}`$ and $`m`$. Multi-connectedness breaks the global isotropy and sometimes the global homogeneity of the universe, except in projective space (see, e.g., Ref. ). Consequently, the CMB temperature angular correlation function will depend on the two directions of observation, not only on their relative angle, and possibly on the position of the observer. This induces correlations between the $`a_\mathrm{}m`$ of different $`\mathrm{}`$ and $`m`$. Such correlations are hidden when one considers only the angular correlation function and its coefficients, the so-called $`C_{\mathrm{}}`$, in a Legendre polynomial decomposition, because they pick up only the isotropic part of the information and are therefore a poor indicator of the topology. This work aims to detail the whole computation of the correlation matrix

$$C_\mathrm{}m^\mathrm{}^{}m^{}a_\mathrm{}ma_\mathrm{}^{}m^{}^{},$$
(1)
which encodes all the topological properties of the CMB, and from which one can compute the usual $`C_{\mathrm{}}`$, simulate maps, and so on.

The study of the detectability of a topological signal (if it exists) in forthcoming CMB datasets requires simulating high quality maps containing the topological signature for a wide class of topologies. Up to now, most CMB studies considered only compact Euclidean spaces  and some compact hyperbolic spaces , and focused mainly on the $`C_{\mathrm{}}`$. The approach developed in this article, and first introduced in Ref. , is well suited to simulate the required CMB maps in any topology once the eigenmodes of the Laplacian have been determined. It paves the way to the simulation of maps for a wide range of topologies, particularly spherical ones.

Recent measurements of the density parameter $`\mathrm{\Omega }`$ imply that the observable universe is “approximately flat”222The popular expression “flat universe” is misleading, because in General Relativity the “universe” is not three-, but four-dimensional and the Friedmann-Lemaître solutions are dynamical, so that the universe is not flat — only its spatial section may be flat or nearly so. In what follows, we will implicitly assume we are talking about the three-dimensional spacelike sections of the universe when talking about flat, hyperbolic, or spherical spaces/universes., perhaps with a slight curvature. The exact constraint on the total density parameter obtained from CMB experiments depends on the priors used during the data analysis. For example with a prior on the nature of the initial conditions, the Hubble parameter and the age of the universe, recent analysis of the DASI, BOOMERanG, MAXIMA and DMR data  lead to $`\mathrm{\Omega }=0.99\pm 0.12`$ at $`1\sigma `$ level, and to $`\mathrm{\Omega }=1.04\pm 0.05`$ at $`1\sigma `$ level if one takes into account only the DASI, BOOMERanG and CBI data. Including stronger priors can indeed sharpen the bound. For instance, including information respectively on large scale structure and on supernovae data leads to $`\mathrm{\Omega }=1.01_{0.06}^{+0.09}`$ and $`\mathrm{\Omega }=1.02_{0.08}^{+0.09}`$ at $`1\sigma `$ level while including both finally leads to $`\mathrm{\Omega }=1.00_{0.06}^{+0.10}`$. This has been recently improved by the Archeops balloon experiments  which get, with a prior on the Hubble constant, $`\mathrm{\Omega }=1.00_{0.02}^{+0.03}`$. In conclusion, it is fair to assert that current cosmological observations set a reliable bound $`0.9<\mathrm{\Omega }<1.1`$. These results are consistent with Friedmann-Lemaître universe models with spherical, flat or hyperbolic spatial sections. In the spherical and hyperbolic cases, $`\mathrm{\Omega }1`$ implies that the curvature radius must be larger than the horizon radius. In all three cases — spherical, flat, and hyperbolic — the universe may be simply connected or multiply connected.

The possibility of detecting the topology of a nearly flat universe was discussed in Ref. . It was noted that the chances of detecting a multiply connected topology are worst in a large hyperbolic universe. The reason is that the typical translation distance between a cosmic source and its nearest topological image seems to be on the order of the curvature radius, and that when $`\mathrm{\Omega }1`$ the distance to the last scattering surface is less than the half of that distance. See Refs.  for some studies on detectability of nearly flat hyperbolic universes. In a multiply connected flat universe the topology scale is completely independent of the horizon radius, because Euclidean geometry — unlike spherical and hyperbolic geometry — has no preferred scale and admits similarities. Note that in the Euclidean case, there are only ten compact topologies, which reduces and simplifies the analysis (in particular regarding the computation of the eigenmodes of the Laplacian). In a spherical universe the topology scale depends on the curvature radius, but, in contrast to the hyperbolic case, as the topology of a spherical manifold gets more complicated, the typical distance between two images of a single cosmic source decreases. No matter how close $`\mathrm{\Omega }`$ is to $`1`$, only a finite number of spherical topologies are excluded from detection. The particular case of the detectability of lens spaces was studied in Ref. , which also considers the detectability of hyperbolic topologies.

At present, the status of the constraints on the topology of the universe is still very preliminary. Regarding locally Euclidean spaces, it was shown on the basis of the COBE data that in the case of a vanishing cosmological constant the size of the fundamental domain of a 3-torus has to be larger than $`L4800h^1\text{Mpc}`$ , where the length $`L`$ is related to the smallest wavenumber $`2\pi /L`$ of the fundamental domain, which induces a suppression of fluctuations on scales beyond the size $`L`$ of the fundamental domain. This constraint does not exclude a toroidal universe since there can be up to eight copies of the fundamental cell within our horizon. This constraint relies mainly on the fact that the smallest wavenumber is $`2\pi /L`$, which induces a suppression of fluctuations on scales beyond the size of the fundamental domain. This result was generalized to all Euclidean manifolds in Ref. . A non-vanishing cosmological constant induces more power on large scales, via the integrated Sachs-Wolfe effect. For instance if $`\mathrm{\Omega }_\mathrm{\Lambda }=0.9`$ and $`\mathrm{\Omega }_\mathrm{m}=0.1`$, the constraint is relaxed to allow for 49 copies of the fundamental cell within our horizon. The constraint is also milder in the case of compact hyperbolic manifolds and it was shown  that the angular power spectrum was consistent with the COBE data on multipoles ranging from 2 to 20 for the Weeks and Thurston manifolds. Another approach, based on the method of images, was developed in . Only one spherical space using this method of images was considered in the literature, namely projective space . The tools developed in this article, as well as in our preceding works , will let us fill the gap regarding the simulation of CMB maps in spherical universes.

Technically, in standard relativistic cosmology, the universe is described by a Friedmann-Lemaître spacetime with locally isotropic and homogeneous spatial sections. These spatial sections can be defined as constant density or time hypersurfaces. With such a splitting, the equations of evolution of the cosmological perturbations, which give birth to the large scale structures of the universe, reduce to a set of coupled differential equations involving the Laplacian. This system is conveniently solved in Fourier space. In the case of a multiply connected universe, we visualize space as the quotient $`X/\mathrm{\Gamma }`$ of a simply connected space $`X`$ (which is just a $`3`$-sphere $`𝐒^3`$, a Euclidean space $`𝐑^3`$, or a hyperbolic space $`𝐇^3`$, depending on the curvature) by a group $`\mathrm{\Gamma }`$ of symmetries of $`X`$ that is discrete and fixed point free. The group $`\mathrm{\Gamma }`$ is called the holonomy group. To solve the evolution equations we must first determine the eigenmodes $`\mathrm{{\rm Y}}_k^{[\mathrm{\Gamma }]}`$ and eigenvalues $`\kappa _k^2`$ of the Laplacian on $`X/\mathrm{\Gamma }`$ through the generalized Helmholtz equation

$$\mathrm{\Delta }\mathrm{{\rm Y}}_k^{[\mathrm{\Gamma }]}=\kappa _k^2\mathrm{{\rm Y}}_k^{[\mathrm{\Gamma }]},$$
(2)
with

$$\kappa _k^2=k^2K,$$
(3)
where $`k`$ indexes the set of eigenmodes, the constant $`K`$ is positive, zero or negative according to whether the space is spherical, flat or hyperbolic, respectively,333We work here in comoving units, but when spatial section are not flat, the curvature of the comoving space is not normalized: it is not assumed to be $`+1`$ in the spherical case nor is it assumed to be $`1`$ in the hyperbolic case. Thus our eigenvalues may differ numerically from those found in the mathematical literature, although of course they agree up to a fixed constant multiple $`|K|^{\frac{1}{2}}`$. and the boundary conditions are compatible with the given topology. The Laplacian in Eq. (2) is defined as $`\mathrm{\Delta }D^iD_i`$, $`D_i`$ being the covariant derivative associated with the metric $`\gamma _{ij}`$ of the spatial sections ($`i`$, $`j=1`$, $`2`$, $`3`$). The eigenmodes of $`X/\mathrm{\Gamma }`$, on which any function on $`X/\mathrm{\Gamma }`$ can be developed, respect the boundary conditions imposed by the topology. That is, the eigenmodes of $`X/\mathrm{\Gamma }`$ correspond precisely to those eigenmodes of $`X`$ that are invariant under the action of the holonomy group $`\mathrm{\Gamma }`$. Thus any linear combination of such eigenmodes will satisfy, by construction, the required boundary conditions. In this way we visualize the space of eigenmodes of $`X/\mathrm{\Gamma }`$ as a subspace of the space of eigenmodes of $`X`$, namely the subspace that is invariant under the action of $`\mathrm{\Gamma }`$. The computational challenge is to find this $`\mathrm{\Gamma }`$-invariant subspace and construct an orthonormal basis for it. In the case of flat manifolds the eigenmodes of $`X/\mathrm{\Gamma }=𝐑^3/\mathrm{\Gamma }`$ can be found analytically. In the case of hyperbolic manifolds, many numerical investigations of the eigenmodes of $`X/\mathrm{\Gamma }=𝐇^3/\mathrm{\Gamma }`$ have been performed . In the case of spherical manifolds, the eigenmodes of $`X/\mathrm{\Gamma }=𝐒^3/\mathrm{\Gamma }`$ have been found analytically for lens and prism spaces  and otherwise can be found numerically .

In the following, we will develop the eigenmodes of $`X/\mathrm{\Gamma }`$ on the basis $`𝒴_{k\mathrm{}m}^{[X]}`$ of the eigenmodes of the universal covering space as

$$\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]}=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}𝒴_{k\mathrm{}m}^{[X]},$$
(4)
so that all the topological information is encoded in the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$, where $`s`$ labels the various eigenmodes sharing the same eigenvalue $`\kappa _k^2`$, both of which are discrete numbers444The spectrum is discrete when the space is compact, e.g. the torus or any spherical space. In a non-compact multi-connected space such as a cylinder it will have a continuous component.. The sum over $`\mathrm{}`$ runs from $`0`$ to infinity if the universal covering space is non-compact (i.e., hyperbolic or Euclidean). The $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$ coefficients can be determined analytically for Euclidean manifolds. In the case of a spherical manifold, the modes are discrete,

$$k=(\nu +1)\sqrt{K}$$
(5)
where $`\nu `$ is a nonnegative integer (the cases $`\nu =0`$ and $`\nu =1`$ correspond to pure gauge modes ), so that $`\kappa _k^2=\nu (\nu +2)K`$, and the sum over $`\mathrm{}`$ is finite and runs from 0 to $`\nu `$ since the multiplicity of a mode $`k`$ is at most its multiplicity in the universal cover, which is $`(\nu +1)^2`$. Among spherical spaces, the case of prism and lens spaces are the simplest since one can determine these coefficients analytically . The worst situation is that of compact hyperbolic manifolds, which is analogous to the Euclidean case since the universal covering $`𝐑^3`$ is not compact but for which no analytical forms are known for the eigenmodes. One then needs to rely on numerical computations (see, e.g., Refs. ).

Our preceding works provided a full classification of spherical manifolds  and developed methods to compute the eigenmodes of the Laplacian in them . Among all spherical manifolds, we were able to obtain analytically the spectrum of the Laplacian for lens and prism spaces  which form two countable families of spherical manifolds. The goal of the present article is to simulate CMB maps for these two families of spherical topologies, as well as for the Euclidean topologies.

We first review briefly the physics of CMB anisotropies (Sec. II) mainly to explain how to take into account the topology (Sec. II.2) once the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$ are determined. We then detail the computation of these coefficients, focusing on the cases where it can be performed analytically, that is for flat spaces and lens and prism spaces. We then present results of numerical simulations (Sec. IV) as well as simulated maps. We discuss these maps qualitatively and confirm that the expected topological correlations (namely matching circles ) are indeed present. The effects of the integrated Sachs-Wolfe and Doppler terms, as well as the thickness of the last scattering surface, are discussed in order to give some insight into the possible detectability of these correlations. Figure 1 summarizes the different independent steps of the computation as well as their interplay.

## Notations

The local geometry of the universe is described by a Friedmann-Lemaître metric

$$\mathrm{d}s^2=c^2\mathrm{d}t^2+a^2(t)\left[\mathrm{d}\chi ^2+\text{s}_K^2(\chi )\mathrm{d}\mathrm{\Omega }^2\right],$$
(6)
where $`a(t)`$ is the scale factor, $`t`$ the cosmic time, $`\mathrm{d}\mathrm{\Omega }^2\mathrm{d}\theta ^2+\mathrm{sin}^2\theta \mathrm{d}\phi ^2`$ the infinitesimal solid angle, $`\chi `$ is the comoving radial distance, and where

$$\text{s}_K(\chi )=\{\begin{array}{cc}\mathrm{sinh}(\sqrt{|K|}\chi )/\sqrt{|K|}\hfill & (\mathrm{hyperbolic}\mathrm{case})\hfill \\ \chi \hfill & (\mathrm{flat}\mathrm{case})\hfill \\ \mathrm{sin}(\sqrt{K}\chi )/\sqrt{K}\hfill & (\mathrm{spherical}\mathrm{case})\hfill \end{array}$$
(7)
In the case of spherical and hyperbolic spatial sections, we also introduce the dimensionless coordinate

$$\overline{\chi }=\sqrt{|K|}\chi ,$$
(8)
which expresses radial distances in units of the curvature radius.

## II CMB anisotropies

The equations dictating the evolution of cosmological perturbations are local differential equations, and remain unchanged by the topology of the spatial sections. Indeed the goal of this section is not to derive the whole set of equations that has to be solved (see, e.g., Refs.  for reviews) but to explain the steps that must be changed to take into account the topology.

### II.1 Local physics of cosmological perturbations

The CMB is observed to a high accuracy as black-body radiation with temperature $`T_0=2.726\pm 0.002\text{K}`$ , almost independently of the direction. After accounting for the peculiar motion of the Sun and Earth, the CMB has tiny temperature fluctuations of order $`\delta T/T_010^5`$ that are usually decomposed in terms of spherical harmonics

$$\frac{\delta T}{T_0}(\theta ,\phi )=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_\mathrm{}m^{\mathrm{obs}}Y_{\mathrm{}}^m(\theta ,\phi ).$$
(9)
This relation can be inverted by using the orthonormality of the spherical harmonics to get

$$a_\mathrm{}m=\frac{\delta T}{T}Y_{\mathrm{}}^m{}_{}{}^{}\mathrm{sin}\theta \mathrm{d}\theta \mathrm{d}\phi .$$
(10)
The coefficients $`a_\mathrm{}m`$ obviously satisfy the conjugation relation

$$a_\mathrm{}m=(1)^ma_\mathrm{}m^{}.$$
(11)
The angular correlation function of these temperature anisotropies is observed on a $`2`$-sphere around us and can be decomposed on a basis of the Legendre polynomials $`P_{\mathrm{}}`$ as

$$\frac{\delta T}{T}(\widehat{\gamma }_1)\frac{\delta T}{T}(\widehat{\gamma }_2)_{\widehat{\gamma }_1.\widehat{\gamma }_2=\mathrm{cos}\theta _{12}}=C^{\mathrm{obs}}(\theta _{12})=\frac{1}{4\pi }\underset{\mathrm{}}{}(2\mathrm{}+1)C_{\mathrm{}}^{\mathrm{obs}}P_{\mathrm{}}(\mathrm{cos}\theta _{12}),$$
(12)
where the brackets stand for an average on the sky, i.e., on all pairs of directions $`(\widehat{\gamma }_1,\widehat{\gamma }_2)`$ subtending an angle $`\theta _{12}`$. The coefficients $`C_{\mathrm{}}^{\mathrm{obs}}`$ of the development of $`C^{\mathrm{obs}}(\theta _{12})`$ on the Legendre polynomials are thus given by

$$C_{\mathrm{}}^{\mathrm{obs}}=\frac{1}{2\mathrm{}+1}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_\mathrm{}m^{\mathrm{obs}}a_\mathrm{}m^{\mathrm{obs}}{}_{}{}^{}.$$
(13)
These $`C_{\mathrm{}}^{\mathrm{obs}}`$ can be seen as estimators of the variance of the $`a_\mathrm{}m`$555One can show that this is indeed the best estimator of their variance when the fluctuations are isotropic and Gaussian . and represent the rotationally invariant angular power spectrum. They have therefore to be compared to the values $`C_{\mathrm{}}`$ predicted by a given cosmological model, which is specified by (i) a model of structure formation which fixes the initial conditions for the perturbation (e.g., inflation, topological defects, etc), (ii) the geometry and matter content of the universe (via the cosmological parameters) and (iii) the topology of the universe.

For the case of a simply connected topology, the $`C_{\mathrm{}}`$ are usually computed as follows. First, one starts from a set of initial conditions given by an early universe scenario. Typically these are the initial conditions predicted in the framework of an inflationary model, although this is of no importance in the discussion that follows. Second, the modes of cosmological interest are evolved from an epoch before they enter into the Hubble radius till now. Third, one computes the variance for the $`a_\mathrm{}m`$. The details of this procedure are well known .

For simplicity let us sketch the case of a flat universe. The temperature fluctuation in a given direction of the sky can be related to (i) the eigenmodes $`\mathrm{exp}(i𝐤𝐱)`$ of the Laplacian (here $`𝐱`$ represents a vector in the usual Cartesian coordinate system) by a linear convolution operator $`O_k^{[𝐑^3]}\left(\mathrm{exp}(i𝐤𝐱)\right)`$ depending on the modulus $`k`$ only \[as well as on the universal cover (here, $`𝐑^3`$) and the cosmological parameters\], and (ii) a three-dimensional variable $`\widehat{e}_𝐤`$ related to the initial conditions, by the formula

$$\frac{\delta T}{T}(\theta ,\phi )=\frac{\mathrm{d}^3𝐤}{(2\pi )^{3/2}}O_k^{[𝐑^3]}\left(\mathrm{e}^{i𝐤𝐱}\right)\sqrt{𝒫_\varphi (k)}\widehat{e}_𝐤,$$
(14)
where $`𝒫_\varphi (k)`$ is the gravitational initial power spectrum, normalized so that $`𝒫_\varphi (k)k^3`$ for a Harrison-Zel’dovich spectrum, and where in most inflationary models the random variable $`\widehat{e}_𝐤`$ describes a Gaussian random field and satisfies

$$\widehat{e}_𝐤\widehat{e}_𝐤^{}^{}=\delta ^\mathrm{D}(𝐤𝐤^{})$$
(15)
with $`\widehat{e}_𝐤=\widehat{e}_𝐤^{}`$ when the space is simply connected. This relation stems from the fact that the temperature fluctuation is real, but it may be expressed differently for multi-connected spaces. Decomposing the exponential by means of Eq. (129) and using Eq. (127) allows to rewrite the temperature fluctuation as

$`{\displaystyle \frac{\delta T}{T}}(\theta ,\phi )`$ $`=`$ $`{\displaystyle \underset{\mathrm{},m}{}}i^{\mathrm{}}{\displaystyle }k^2\mathrm{d}k\sqrt{𝒫_\varphi (k)}O_k^{[𝐑^3]}\left(\sqrt{{\displaystyle \frac{2}{\pi }}}j_{\mathrm{}}(k\chi )\right)\left[{\displaystyle }\mathrm{d}\mathrm{\Omega }_𝐤Y_{\mathrm{}}^m{}_{}{}^{}(\theta _𝐤,\phi _𝐤)\widehat{e}_𝐤\right]Y_{\mathrm{}}^m(\theta ,\phi )`$ (16)
$`=`$ $`{\displaystyle \underset{\mathrm{},m}{}}i^{\mathrm{}}{\displaystyle k^2dk\sqrt{𝒫_\varphi (k)}O_k^{[𝐑^3]}\left(𝒴_{k\mathrm{}m}^{[𝐑^3]}\right)\widehat{e}_\mathrm{}m(k)},`$
where we have defined

$$\widehat{e}_\mathrm{}m(k)\mathrm{d}\mathrm{\Omega }_𝐤Y_{\mathrm{}}^m{}_{}{}^{}(\theta _𝐤,\phi _𝐤)\widehat{e}_𝐤,$$
(17)
which is the “average” of the random fields $`\widehat{e}_𝐤`$ over all the $`𝐤`$ of the same modulus. This quantity can therefore be identified as a two-dimensional Gaussian random variable satisfying $`\widehat{e}_\mathrm{}m(k)\widehat{e}_\mathrm{}^{}m^{}^{}(k^{})=\delta ^\mathrm{D}(kk^{})\delta _{\mathrm{}\mathrm{}^{}}\delta _{mm^{}}/k^2`$. It follows that the coefficients $`a_\mathrm{}m`$ take the general form

$$a_\mathrm{}m=i^{\mathrm{}}k^2dk\sqrt{𝒫_\varphi (k)}G_{\mathrm{}}(k)\widehat{e}_\mathrm{}m(k),$$
(18)
with

$$G_{\mathrm{}}(k)=O_k^{[𝐑^3]}\left(R_k\mathrm{}^{[𝐑^3]}\right),$$
(19)
and here the function $`G_{\mathrm{}}(k)`$ can be approximated by (see, e.g., )

$`G_{\mathrm{}}(k)`$ $`=`$ $`j_{\mathrm{}}(k(\eta _0\eta _{\mathrm{LSS}}))\left({\displaystyle \frac{\delta T}{T}}(k,\eta _{\mathrm{LSS}})+\mathrm{\Phi }(k,\eta _{\mathrm{LSS}})+\mathrm{\Psi }(k,\eta _{\mathrm{LSS}})\right)+j_{\mathrm{}}^{}(k(\eta _0\eta _{\mathrm{LSS}})){\displaystyle \frac{v_\mathrm{b}(k,\eta _{\mathrm{LSS}})}{k}}`$ (20)
$`+{\displaystyle _{\eta _{\mathrm{LSS}}}^{\eta _0}}j_{\mathrm{}}(k(\eta _0\eta ))\left(\dot{\mathrm{\Phi }}(k,\eta )+\dot{\mathrm{\Psi }}(k,\eta )\right)d\eta ,`$
which is indeed a linear convolution operator acting on $`R_k\mathrm{}^{[𝐑^3]}`$ as announced above. Here, $`\eta _{\mathrm{LSS}}`$ and $`\eta _0`$ are the conformal time at the last scattering epoch and today, respectively, $`j_{\mathrm{}}`$ is the spherical Bessel function of index $`\mathrm{}`$, $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ are the two Bardeen potentials, and $`v_\mathrm{b}`$ is the velocity divergence of the baryons. The only modification of note when one considers a non flat universe is that the $`j_{\mathrm{}}`$ are to be replaced by their analog for non flat geometries, the so-called ultraspherical Bessel functions , and for a closed universe the integral over $`k`$ is replaced by a discrete sum.

In conclusion, and without loss of generality, the temperature fluctuation can be decomposed as in Eq. (16) whatever the curvature of space. For a simply connected topology and Gaussian initial conditions the addition property (127) of the spherical harmonics imposes that

$$a_\mathrm{}ma_\mathrm{}^{}m^{}^{}=C_{\mathrm{}}\delta _{\mathrm{}\mathrm{}^{}}\delta _{mm^{}},$$
(21)
and therefore the $`C_{\mathrm{}}`$ coefficients encode all the information regarding the CMB anisotropies.

### II.2 Implementing the topology

The topology does not affect local physics, so the equations describing the evolution of the cosmological perturbations are left unchanged. As a consequence, quantities such as the Bardeen potentials $`\mathrm{\Phi }`$, $`\mathrm{\Psi }`$, etc., are computed in the same way as described above, and the operator $`O_k^{[X]}`$ is therefore the same. However, a change of topology translates into a change of the modes that can exist in the universe. In particular, the functions $`𝒴_{k\mathrm{}m}^{[X]}`$ are typically not well defined on the quotient space $`X/\mathrm{\Gamma }`$. Therefore, the only change that has to be performed is the substitution

$$𝒴_{k\mathrm{}m}^{[X]}\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]},$$
(22)
where the $`\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]}`$ form an orthonormal basis for the space of eigenmodes of the Laplacian on the given topology $`X/\mathrm{\Gamma }`$. One must then remember that the mode functions $`\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]}`$ can be decomposed uniquely by Eq. (4) and that the convolution operator $`O_k^{[X]}`$ is linear. When the multi-connected topology is compact, it follows that Eq. (16) will take the form

$$\frac{\delta T}{T}(\theta ,\phi )=\frac{(2\pi )^3}{V}\underset{k,s}{}O_k^{[X]}\left(\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]}\right)\sqrt{𝒫_\varphi (k)}\widehat{e}_𝐤,$$
(23)
where now $`\widehat{e}_𝐤`$ is a three-dimensional random variable which is related to the discrete mode $`𝐤`$. Equivalently one can write $`\widehat{e}_𝐤\widehat{e}_{ks}`$ where $`k`$ is the modulus of $`𝐤`$ and the index $`s`$ describes all the eigenmodes of the Laplacian for fixed modulus $`k`$ in the topology $`X/\mathrm{\Gamma }`$ of volume $`V`$. These random variables satisfy the normalization

$$\widehat{e}_𝐤\widehat{e}_𝐤^{}^{}=\frac{V}{(2\pi )^3}\delta _{kk^{}}\delta _{ss^{}}.$$
(24)
For a given value of the wavenumber $`k`$, there are fewer eigenmodes in the multi-connected case, so that $`s`$ has to be seen as a “subset” of the set $`\{\mathrm{},m\}`$.

By inserting the expansion of $`\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]}`$ in terms of the covering space eigenmodes, as given by Eq. (4), we obtain

$$\frac{\delta T}{T}(\theta ,\phi )=\frac{(2\pi )^3}{V}\underset{k,s}{}\underset{\mathrm{},m}{}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}O_k^{[X]}\left(𝒴_{k\mathrm{}m}^{[X]}\right)\sqrt{𝒫_\varphi (k)}\widehat{e}_𝐤.$$
(25)
It follows that the $`a_\mathrm{}m`$, seen as random variables, are given by

$$a_\mathrm{}m=\frac{(2\pi )^3}{V}\underset{k}{}\sqrt{𝒫_\varphi (k)}O_k^{[X]}\left(R_k\mathrm{}^{[X]}\right)\underset{s}{}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}\widehat{e}_𝐤.$$
(26)
Note that the sum over $`s`$ is analogous to the sum over angles defining the two-dimensional random variable $`\widehat{e}_\mathrm{}m`$ in Eq. (16). Since the $`a_\mathrm{}m`$ are linear functions of the initial three-dimensional random variables, they are still Gaussian distributed but they are not independent anymore (as explained before, this is the consequence of the breakdown of global isotropy and/or homogeneity). The correlation between the coefficients $`a_\mathrm{}m`$ is given by

$`a_\mathrm{}ma_\mathrm{}^{}m^{}^{}={\displaystyle \frac{(2\pi )^3}{V}}{\displaystyle \underset{k}{}}𝒫_\varphi (k)O_k^{[X]}\left(R_k\mathrm{}^{[X]}\right)O_k^{[X]}\left(R_k\mathrm{}^{}^{[X]}\right){\displaystyle \underset{s}{}}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}\xi _{k\mathrm{}^{}m^{}}^{[\mathrm{\Gamma }]s}.`$ (27)
Clearly these correlations can have non-zero off-diagonal terms, reflecting the global anisotropy induced by the multi-connected topology, so that Eq. (21) no longer holds and the observational consequences of a given topology on the CMB anisotropies are given by the correlation matrix (1). This means in particular that for fixed $`\mathrm{}`$, the $`a_\mathrm{}m`$ might not have the same variance, although they all follow Gaussian statistics as long as the initial conditions do. This translates into an apparent non Gaussianity in the sense that the $`C_{\mathrm{}}`$ will not follow the usual $`\chi ^2`$ distribution. Strictly speaking, this is not a signature of non Gaussianity but of anisotropy.

Note also that the correlation matrix (27) is not rotation invariant. It will explicitly depend on the orientation of the manifold with respect of the coordinate system. However, knowing how the spherical harmonics transform under a rotation allows us to compute the correlation matrix under any other orientation of the coordinate system. To finish let us note that one can define the usual $`C_{\mathrm{}}`$ coefficients in any topology by the formula

$$C_{\mathrm{}}\frac{1}{2\mathrm{}+1}\underset{m}{}C_\mathrm{}m^\mathrm{}m,$$
(28)
which is easily shown to be rotationally invariant.

## III Eigenmodes of multi-connected spaces

### III.1 Flat spaces

The purpose of this subsection is to compute in detail the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$ in the case of a cubic 3-torus of comoving size $`L`$, referred to as $`T_1`$. The method generalized easily to any compact flat manifold.

In the case at hand, the topology implies the “quantization” of the allowed wave vectors

$$𝐤=\frac{2\pi }{L}\left(n_1\widehat{𝐱}+n_2\widehat{𝐲}+n_3\widehat{𝐳}\right)=\frac{2\pi }{L}𝐧=k\widehat{𝐧},$$
(29)
with

$`𝐧`$ $``$ $`(n_1,n_2,n_3),`$ (30)
$`n`$ $``$ $`\sqrt{𝐧𝐧},`$ (31)
$`\widehat{𝐧}`$ $``$ $`(n_1,n_2,n_3)/n,`$ (32)
$`k`$ $``$ $`\sqrt{𝐤𝐤}={\displaystyle \frac{2\pi }{L}}n,`$ (33)
so that the label $`s`$ can be chosen to be the “unit” integer vector $`\widehat{𝐧}`$. The label multiplicity $`\mathrm{mult}(k)\mathrm{card}\{s\}`$ of a mode $`k`$ is given in this case by the number of representations of $`n^2`$ by 3 squares, i.e., by $`\mathrm{card}\{\widehat{𝐧}\}`$ (see Fig. 2 below). The corresponding normalized eigenmodes in Cartesian coordinates are thus simply given by

$$\mathrm{{\rm Y}}_{k\widehat{𝐧}}^{[T_1]}(𝐫)=\frac{\mathrm{e}^{i𝐤𝐫}}{(2\pi )^{3/2}},$$
(34)
with $`𝐤`$ given by Eq. (29). Using the decomposition (129) of the exponential and plugging in the closure relation (127), one gets that

$$\xi _{k\mathrm{}m}^{[T_1]s}=\xi _{2\pi n/L\mathrm{}m}^{[T_1]\widehat{𝐧}}=i^{\mathrm{}}Y_{\mathrm{}}^m{}_{}{}^{}(\widehat{𝐧}),$$
(35)
where $`\widehat{𝐧}`$ can also be defined by the two spherical angles $`(\theta _𝐧`$, $`\phi _𝐧)`$ which are explicitly given by

$`\mathrm{tan}\theta _𝐧`$ $`=`$ $`{\displaystyle \frac{\sqrt{n_1^2+n_2^2}}{n_3}},`$ (36)
$`\mathrm{tan}\phi _𝐧`$ $`=`$ $`{\displaystyle \frac{n_2}{n_1}}.`$ (37)
This expression could also have been obtained by simply considering the Fourier transform of $`𝒴_{k\mathrm{}m}^{[𝐑^3]}`$ as given by Eq. (128). One can check that the normalization of the basis $`\mathrm{{\rm Y}}_{k\widehat{𝐧}}^{[T_1]}`$ \[i.e., $`\mathrm{{\rm Y}}_{k\widehat{𝐧}}^{[T_1]}\mathrm{{\rm Y}}_{k^{}\widehat{𝐧}^{}}^{[T_1]}{}_{}{}^{}\mathrm{d}_{}^{3}𝐱=\delta ^\mathrm{D}(𝐤𝐤^{})`$\] implies that the coefficients $`\xi _{k\mathrm{}m}^{[T_1]\widehat{𝐧}}`$ satisfy the closure relation

$$\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\xi _{k\mathrm{}m}^{[T_1]\widehat{𝐧}}\xi _{k^{}\mathrm{}m}^{[T_1]\widehat{𝐧}^{}}=\delta ^\mathrm{D}(\mathrm{cos}\theta _𝐧\mathrm{cos}\theta _𝐧^{})\delta ^\mathrm{D}(\phi _𝐧\phi _𝐧^{})\delta _{kk^{}}.$$
(38)
The Dirac distribution in the above expression can be shown, by using either Eqs. (126,127) or Eq. (130) alone, to be

$$\delta ^\mathrm{D}(\mathrm{cos}\theta _𝐧\mathrm{cos}\theta _𝐧^{})\delta ^\mathrm{D}(\phi _𝐧\phi _𝐧^{})=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\frac{2\mathrm{}+1}{4\pi }P_{\mathrm{}}(\widehat{𝐧}.\widehat{𝐧}^{}).$$
(39)
¿From these results, we deduce that the correlation matrix of the $`a_\mathrm{}m`$ is given by

$$C_\mathrm{}m^\mathrm{}^{}m^{}=\left(\frac{2\pi }{L}\right)^3\underset{n}{}i^{\mathrm{}\mathrm{}^{}}𝒫_\varphi \left(\frac{2\pi n}{L}\right)G_{\mathrm{}}\left(\frac{2\pi n}{L}\right)G_{\mathrm{}^{}}\left(\frac{2\pi n}{L}\right)\underset{\widehat{𝐧}}{}Y_{\mathrm{}}^m{}_{}{}^{}(\widehat{𝐧})Y_{\mathrm{}^{}}^m^{}(\widehat{𝐧}).$$
(40)
Using Eq. (127) and the fact that $`_{n,\widehat{𝐧}}=_𝐧`$, the $`C_{\mathrm{}}`$ coefficients are simply

$$C_{\mathrm{}}=\frac{1}{2\pi }\left(\frac{2\pi }{L}\right)^3\underset{𝐧}{}𝒫_\varphi \left(\frac{2\pi n}{L}\right)G_{\mathrm{}}^2\left(\frac{2\pi n}{L}\right),$$
(41)
which was used in many earlier works  on the influence of topology on the CMB.

Note that spherical harmonics satisfy the following symmetry relation

$$Y_{\mathrm{}}^m(\theta ,\phi )=\mathrm{e}^{2im\phi }Y_{\mathrm{}}^m(\theta ,\phi ).$$
(42)
Due to the symmetry of the torus with respect to the $`y=0`$ plane, in the sum over $`\widehat{𝐧}`$ in Eq. (40) a term $`(\theta _𝐧,\phi _𝐧)`$ will always be associated to a term $`(\theta _𝐧,\phi _𝐧)`$ leading to a term of the form $`4\mathrm{cos}\left((m^{}m)\phi _𝐧\right)Y_{\mathrm{}}^m{}_{}{}^{}(\theta _𝐧,0)Y_{\mathrm{}^{}}^m^{}(\theta _𝐧,0)`$, the only exception being the term arising when $`n_2=0`$, which is real. ¿From this result, one easily shows that the correlation matrix satisfies

$`C_\mathrm{}m^\mathrm{}^{}m^{}𝐑.`$ (43)
Along similar lines, the relations

$`Y_{\mathrm{}}^m(\theta ,\pi )`$ $`=`$ $`\text{e}^{im\phi }Y_{\mathrm{}}^m(\theta ,0),`$ (44)
$`Y_{\mathrm{}}^m(\pi \theta ,\phi )`$ $`=`$ $`(1)^{\mathrm{}+m}Y_{\mathrm{}}^m(\theta ,\phi ),`$ (45)
imply

$$C_\mathrm{}m^\mathrm{}^{}m^{}=\frac{1}{4}\left[1+(1)^{mm^{}}\right]\left[1+(1)^{\mathrm{}\mathrm{}^{}}\right]C_\mathrm{}m^\mathrm{}^{}m^{},$$
(46)
so that

$$C_\mathrm{}m^\mathrm{}^{}m^{}0mm^{}0\mathrm{mod}(2)\mathrm{and}\mathrm{}\mathrm{}^{}0\mathrm{mod}(2).$$
(47)
Furthermore Eqs.(43,124) imply

$$C_\mathrm{}m^\mathrm{}^{}m^{}=C_\mathrm{}m^\mathrm{}^{}m^{}.$$
(48)
Let us emphasize that these properties of the correlation matrix still hold even if the torus is not cubic. However, a cubic torus is invariant under a $`\pi /2`$-rotation about the $`z`$ axis, so if $`(n_1,n_2,n_3)`$ corresponds to a wavenumber then so does $`(n_2,n_1,n_3)`$, and one has

$$C_\mathrm{}m^\mathrm{}^{}m^{}0mm^{}0\mathrm{mod}(4).$$
(49)
### III.2 Spherical spaces

The goal of this section is to recall the basic analytical results concerning the lens and prism spaces (Sec. III.2.2). In these spaces, the eigenmodes and eigenvalues of the Laplacian operator can be determined analytically using toroidal coordinates . CMB computations use spherical coordinates, so we must perform a change of coordinates and a change of basis (detailed in Appendix B). Fortunately, this can also be achieved analytically to compute the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$.

#### III.2.1 Generalities

In our preceding article , we presented in a pedestrian way the complete classification of three-dimensional spherical topologies and we described how to compute their holonomy transformations.

The isometry group of the $`3`$-sphere is $`\mathrm{SO}(4)`$. Every isometry in $`\mathrm{SO}(4)`$ can be decomposed as the product of a right-handed and a left-handed Clifford translation, and the factorization is unique up to simultaneous multiplication of both factors by -1. Furthermore, the space $`𝐒^3`$ itself enjoys a group structure as the set $`𝒮^3`$ of unit length quaternions. Each right-handed (resp. left-handed) Clifford translation corresponds to left (resp. right) quaternion multiplication of $`𝒮^3`$, so the group of right-handed (resp. left-handed) Clifford translations is isomorphic to $`𝒮^3`$. It follows that $`\mathrm{SO}(4)`$ is isomorphic to $`𝒮^3\times 𝒮^3/\{\pm (\mathrm{𝟏},\mathrm{𝟏})\}`$ and thus the classification of the subgroups of $`\mathrm{SO}(4)`$ can be deduced from the classification of subgroups of $`𝒮^3`$. There is a two-to-one homomorphism from $`𝒮^3`$ to $`\mathrm{SO}(3)`$; the finite subgroups of $`\mathrm{SO}(3)`$ are the cyclic, dihedral, tetrahedral, octahedral and icosahedral groups, so the finite subgroups of $`𝒮^3`$ are their lifts, namely

* the cyclic groups $`Z_n`$ of order $`n`$,
* the binary dihedral groups $`D_m^{}`$ of order $`4m`$, $`m2`$,
* the binary tetrahedral group $`T^{}`$ of order 24,
* the binary octahedral group $`O^{}`$ of order 48,
* the binary icosahedral group $`I^{}`$ of order 120,

where a binary group is the two-fold cover of the corresponding plain group.

¿From this classification, it can be shown that there are three categories of spherical 3-manifolds.

* The single action manifolds are those for which a subgroup $`R`$ (resp. $`L`$) of $`𝒮^3`$ acts as pure right-handed (resp. pure left-handed) Clifford translations. They are thus the simplest spherical manifolds and can all be written as $`𝐒^3/\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }=Z_n,D_m^{},T^{},O^{},I^{}`$.
* The double-action manifolds are those for which subgroups $`R`$ and $`L`$ of $`𝒮^3`$ act simultaneously as right- and left-handed Clifford translations, and every element of $`R`$ occurs with every element of $`L`$. There are obtained for the groups $`\mathrm{\Gamma }=\mathrm{\Gamma }_1\times \mathrm{\Gamma }_2`$ with $`(\mathrm{\Gamma }_1,\mathrm{\Gamma }_2)=(Z_m,Z_n)`$, $`(D_m^{},Z_n)`$, $`(T^{},Z_n)`$, $`(O^{},Z_n)`$, $`(I^{},Z_n)`$, with $`\mathrm{gcd}(m,n)=1`$, $`\mathrm{gcd}(4m,n)=1`$, $`\mathrm{gcd}(24,n)=1`$, $`\mathrm{gcd}(48,n)=1`$, and $`\mathrm{gcd}(120,n)=1`$, respectively.
* The linked-action manifolds are similar to the double action manifolds, except that each element of $`R`$ occurs with only some of the elements of $`L`$.

The classification of these manifolds is summarized in Fig. 8 of Ref. .

#### III.2.2 Lens and prism spaces

In this article, we focus on prism spaces $`𝐒^3/D_m^{}`$ and lens spaces $`L(p,q)`$ . The latter are defined by identifying the lower surface of a lens-shaped solid to the upper surface with a rotation of angle $`2\pi q/p`$ for $`p`$ and $`q`$ relatively prime integers with $`0<q<p`$. Furthermore, we may restrict our attention to $`0<qp/2`$ because for values of $`q`$ in the range $`p/2<q<p`$ the twist $`2\pi q/p`$ is the same as $`2\pi (pq)/p`$, thus $`L(p,q)`$ is the mirror image of $`L(p,pq)`$. Lens spaces can be single action, double action, or linked action; Fig. 9 of Ref.  summarizes their classification.

The eigenmodes and eigenvalues of prism and lens spaces can be obtained analytically by working in toroidal coordinates . Starting from Cartesian coordinates, in which the equation for the $`3`$-sphere is $`x^2+y^2+z^2+w^2=1`$, the toroidal coordinates $`(\chi _\mathrm{T},\theta _\mathrm{T},\phi _\mathrm{T})`$ are defined via the equations

$`x`$ $`=`$ $`\mathrm{cos}\chi _\mathrm{T}\mathrm{cos}\theta _\mathrm{T},`$ (50)
$`y`$ $`=`$ $`\mathrm{cos}\chi _\mathrm{T}\mathrm{sin}\theta _\mathrm{T},`$ (51)
$`z`$ $`=`$ $`\mathrm{sin}\chi _\mathrm{T}\mathrm{cos}\phi _\mathrm{T},`$ (52)
$`w`$ $`=`$ $`\mathrm{sin}\chi _\mathrm{T}\mathrm{sin}\phi _\mathrm{T},`$ (53)
with

$`0`$ $``$ $`\chi _\mathrm{T}\pi /2,`$ (54)
$`0`$ $``$ $`\theta _\mathrm{T}2\pi ,`$ (55)
$`0`$ $``$ $`\phi _\mathrm{T}2\pi .`$ (56)
Ref.  gives the eigenmodes of $`𝐒^3`$ explicitly as

$$𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}=B_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}\mathrm{cos}^{|\mathrm{}_\mathrm{T}|}\chi _\mathrm{T}\mathrm{sin}^{|m_\mathrm{T}|}\chi _\mathrm{T}P_d^{(|m_\mathrm{T}|,|\mathrm{}_\mathrm{T}|)}(\mathrm{cos}2\chi _\mathrm{T})f(|\mathrm{}_\mathrm{T}|\theta _\mathrm{T})f(|m_\mathrm{T}|\phi _\mathrm{T}),$$
(57)
where $`\nu `$ is the integer parameterizing $`k=(\nu +1)\sqrt{K}`$ as in Eqn. (5), $`P_d^{(|m_\mathrm{T}|,|\mathrm{}_\mathrm{T}|)}`$ is the Jacobi polynomial, and $`f`$ stands for the cosine (resp. sine) function when $`\mathrm{}_\mathrm{T}`$ or $`m_\mathrm{T}0`$ (resp. $`\mathrm{}_\mathrm{T}`$ or $`m_\mathrm{T}<0`$). For each value of $`\nu `$, the indices $`\mathrm{}_\mathrm{T}`$ and $`m_\mathrm{T}`$ range over all integers satisfying

$`|\mathrm{}_\mathrm{T}|+|m_\mathrm{T}|`$ $``$ $`\nu \mathrm{and}`$ (58)
$`|\mathrm{}_\mathrm{T}|+|m_\mathrm{T}|`$ $``$ $`\nu \mathrm{mod}(2),`$ (59)
and for convenience we define

$$d=\frac{1}{2}(\nu |\mathrm{}_\mathrm{T}||m_\mathrm{T}|).$$
(60)
The normalization coefficients $`B_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}`$ are given by666Note the factor $`1/(\nu +1)`$ which differs from Ref.  due to a different choice of normalization.

$$B_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}=\frac{\sigma _\mathrm{}_\mathrm{T}\sigma _{m_\mathrm{T}}}{\pi (\nu +1)}\sqrt{\frac{2(\nu +1)d!(|\mathrm{}_\mathrm{T}|+|m_\mathrm{T}|+d)!}{(|\mathrm{}_\mathrm{T}|+d)!(|m_\mathrm{T}|+d)!}}$$
(61)
with $`\sigma _i=1/\sqrt{2}`$ if $`i=0`$ and $`\sigma _i=1`$ otherwise.

Using these definitions, Ref.  shows that for lens spaces the explicit set of coefficients $`\eta _{s\mathrm{}_\mathrm{T}m_\mathrm{T}}^{[\mathrm{\Gamma }]\nu }`$ such that

$$\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]}=\underset{\mathrm{}_\mathrm{T},m_\mathrm{T}}{}\eta _{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}^{[\mathrm{\Gamma }]s}𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}(\chi _\mathrm{T},\theta _\mathrm{T},\phi _\mathrm{T}),$$
(62)
can be obtained as follows:

## Theorem 1: lens spaces

The eigenspace of the Laplacian on the lens space $`L(p,q)`$ has an orthonormal basis that, when lifted to $`Z_p`$-invariant eigenmodes of the $`3`$-sphere, comprises those eigenmodes in the left column for which the corresponding condition in the right column is satisfied, subject to the restriction that an eigenmode $`𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}`$ exists if and only if the integers $`\nu `$, $`\mathrm{}_\mathrm{T}`$, and $`m_\mathrm{T}`$ satisfy $`|\mathrm{}_\mathrm{T}|+|m_\mathrm{T}|\nu `$ and $`|\mathrm{}_\mathrm{T}|+|m_\mathrm{T}|\nu \mathrm{mod}(2)`$.

| basis vectors |  | condition |
| --- | --- | --- |
| $`𝒬_{\nu \mathrm{\hspace{0.17em}0\hspace{0.17em}0}}`$ |  | always |
| $`𝒬_{\nu \mathrm{}_\mathrm{T}\mathrm{\hspace{0.17em}0}}`$ |  | $`\mathrm{}_\mathrm{T}0\mathrm{mod}(p)`$ |
| $`𝒬_{\nu \mathrm{\hspace{0.17em}0}m_\mathrm{T}}`$ |  | $`qm_\mathrm{T}0\mathrm{mod}(p)`$ |
| $`(𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}+𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}})/\sqrt{2},(𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}})/\sqrt{2}`$ |  | $`\mathrm{}_\mathrm{T}qm_\mathrm{T}\mathrm{mod}(p)`$ |
| $`(𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}})/\sqrt{2},(𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}+𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}})/\sqrt{2}`$ |  | $`\mathrm{}_\mathrm{T}qm_\mathrm{T}\mathrm{mod}(p)`$ |

An analogous theorem was demonstrated for prism spaces and can be found in Ref. .

Unfortunately, for practical purposes the eigenmodes of the lens and prism spaces are needed in spherical coordinates, while they are most easily obtained in toroidal coordinates. As explained in the Introduction, one needs the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$ of the decomposition (4). Since $`𝒴_{k\mathrm{}m}^{[𝐒^3]}`$ and $`𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}`$ are two orthogonal bases of dimension $`(\nu +1)^2`$, all of whose elements have the same norm, there is an orthogonal transformation taking one to the other

$$𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}=\underset{\mathrm{},m}{}\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}𝒴_{k\mathrm{}m}^{[𝐒^3]}.$$
(63)
The “transpose” of this transformation $`\alpha `$ takes a given eigenmode’s $`𝒬`$-based coefficients $`\eta `$ to its $`𝒴`$-based coefficients $`\xi `$:

$$\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}=\underset{\mathrm{}_\mathrm{T}=\nu }{\overset{\nu }{}}\underset{m_\mathrm{T}=\nu }{\overset{\nu }{}}\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}\eta _{s\mathrm{}_\mathrm{T}m_\mathrm{T}}^{[\mathrm{\Gamma }]\nu }.$$
(64)
The orthonormality of the basis $`\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]}`$ implies that the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$ satisfy

$$\underset{\mathrm{}=0}{\overset{\nu }{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s^{}}=\delta _{ss^{}}.$$
(65)
This relation is simpler than the closure relation (38) obtained in the flat case because, for a given $`k`$, the space of modes is finite-dimensional. The computation of the coefficients $`\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}`$ appears in Appendix B. Because both the $`𝒬`$ basis and the $`𝒴`$ basis are orthonormal, the transformation $`\alpha `$ is orthogonal:

$$\underset{\mathrm{}=0}{\overset{\nu }{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}^{}m_\mathrm{T}^{}}^{}=\delta _{\mathrm{}_\mathrm{T}\mathrm{}_\mathrm{T}^{}}\delta _{m_\mathrm{T}m_\mathrm{T}^{}}\epsilon _\nu (\mathrm{}_\mathrm{T},m_\mathrm{T})\epsilon _\nu (\mathrm{}_\mathrm{T}^{},m_\mathrm{T}^{})$$
(66)
where $`\epsilon _\nu (\mathrm{}_\mathrm{T},m_\mathrm{T})=1`$ if the conditions (58,59) are satisfied and 0 otherwise.

With these coefficients, the CMB computation goes as in the flat case, except for the fact that some integrals have to be replaced by discrete sums. One easily gets that

$$\frac{\delta T}{T}(\theta ,\phi )=\frac{(2\pi )^3}{V}\underset{\nu =2}{\overset{\mathrm{}}{}}\underset{s}{}\underset{\mathrm{},m}{}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}O_k^{[𝐒^3]}\left(𝒴_{k\mathrm{}m}^{[𝐒^3]}\right)\sqrt{𝒫_\varphi (k)}\widehat{e}_𝐤,$$
(67)
so that

$$a_\mathrm{}m=\frac{(2\pi )^3}{V}\underset{\nu =2}{\overset{\mathrm{}}{}}\sqrt{𝒫_\varphi (k)}O_k^{[𝐒^3]}\left(R_k\mathrm{}^{[𝐒^3]}\right)\underset{s}{}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}\widehat{e}_𝐤,$$
(68)
and

$$C_\mathrm{}m^\mathrm{}^{}m^{}=\frac{(2\pi )^3}{V}\underset{\nu =2}{\overset{\mathrm{}}{}}𝒫_\varphi (k)G_{\mathrm{}}(k)G_{\mathrm{}^{}}(k)\underset{s}{}\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}\xi _{k\mathrm{}^{}m^{}}^{[\mathrm{\Gamma }]s},$$
(69)
as first obtained in Ref. . Note also that, following Refs. , a scale invariant spectrum will in that case be defined as

$$𝒫_\varphi (k)\frac{1}{k(k^2K)}.$$
(70)
To finish, let us discuss the properties of the random variable $`\widehat{e}_𝐤`$. Since the eigenmodes in toroidal coordinates, $`𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}`$, and the coefficients $`\eta _{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}^{[\mathrm{\Gamma }]s}`$ are real valued, it follows from Eq. (C2) that

$$\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}^{}=(1)^m\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}.$$
(71)
It follows that

$$\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}{}_{}{}^{}=(1)^m\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s},$$
(72)
whatever $`s`$ and thus that the eigenmodes $`\mathrm{{\rm Y}}_{ks}^{[\mathrm{\Gamma }]}`$ are real-valued. This implies that $`\widehat{e}_𝐤`$ is a real random variable, contrary to the preceding example of the torus.

## IV Numerical computations

### IV.1 Implementation

The correlation matrix for $`\mathrm{}\mathrm{}_{\mathrm{max}}`$ has $`\mathrm{}_{\mathrm{max}}^4`$ coefficients. However, the parity and symmetry relations

$`C_\mathrm{}m^\mathrm{}^{}m^{}`$ $`=`$ $`(1)^{mm^{}}C_\mathrm{}m^\mathrm{}^{}m^{}^{}`$ (73)
$`C_\mathrm{}m^\mathrm{}^{}m^{}`$ $`=`$ $`C_\mathrm{}^{}m^{}^\mathrm{}m^{}`$ (74)
reduce the problem to computing only a quarter of them. Then, for a given topology, symmetries can further reduce the number of coefficients to compute. For example, with a cubic torus, Eqs. (47,49) insure that only one coefficient out of eight is nonzero, and the symmetries (48) also give the coefficients when one changes the sign of both $`m`$ and $`m^{}`$. This leaves only $`\mathrm{}_{\mathrm{max}}^4/64`$ coefficients to compute. For example, a COBE scale map ($`\mathrm{}_{\mathrm{max}}30`$) requires $`12500`$ coefficients, while a Planck scale map ($`\mathrm{}_{\mathrm{max}}1500`$) requires $`8\times 10^{10}`$ coefficients.

Each coefficients is computed using Eq. (40), which involves a sum over all the wavemodes $`𝐤`$. For a given resolution $`\mathrm{}_{\mathrm{max}}`$ the modulus of the largest wavemode is given by $`k_{\mathrm{max}}3\mathrm{}_{\mathrm{max}}/\eta _0`$. Moreover, the density of wavemodes is proportional to the size of the torus, so that we have $`O(\mathrm{}_{\mathrm{max}}^3L^3)`$ modes. Therefore, the computational time and the memory requirement scale as $`\mathrm{}_{\mathrm{max}}^7L^3`$ and $`\mathrm{}_{\mathrm{max}}^3L^3`$, respectively. This is obviously a serious limitation of our algorithm. For example, computing the correlation matrix for a COBE scale map on a relatively small torus ($`L=2R_H`$, where $`R_H`$ is the Hubble radius) takes around $`10`$ hours on a $`900\text{MHz}`$ CPU and allocates $`60`$ Megabytes of memory.

When the topology is simply connected, it is a well-known fact that the $`C_{\mathrm{}}`$ are in general a smooth function of the multipole $`\mathrm{}`$. This reduces the computational time because for $`\mathrm{}_{\mathrm{max}}=1500`$ one needs to compute only $`50`$ coefficients. For the topologies we have studied, we did not find any evidence for a smooth structure of the correlation matrix, at least at the relatively large scales we considered. At which scale one can reliably approximate the correlation matrix by its isotropic diagonal part (the $`C_{\mathrm{}}`$) remains an open question.

Also, if one wants to simulate CMB maps from the correlation matrix, one needs to diagonalize it. This procedure can also take a lot of time because it is an $`\mathrm{}_{\mathrm{max}}^6`$ process. For the case of the torus, however, this problem is not serious as the symmetries of the torus insure that the matrix is block diagonal, with eight blocks if the torus is cubic or four blocks otherwise.

Strictly speaking, one does not need the correlation matrix to compute maps. One can do it directly by using Eq. (18). This amounts to performing a realization of the three-dimensional random field describing the cosmological perturbations, and projecting it onto the sphere. In this case, one has only $`\mathrm{}_{\mathrm{max}}`$ coefficients to compute (the $`a_\mathrm{}m`$) instead of the correlation matrix, so that the memory requirements are roughly the same (one only needs to store the value of the random field for each mode), but the computational time scales as $`\mathrm{}_{\mathrm{max}}^5L^3`$. In this case, computing maps for a cubic torus of size $`L=2R_H`$ till $`\mathrm{}_{\mathrm{max}}=120`$ takes $`3`$ hours on a $`1.7\text{GHz}`$ CPU and allocates 300 Megabytes of memory.

In the case of spherical spaces, the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$ must also be computed numerically. This involves determining both the coefficients $`\eta _{s\mathrm{}_\mathrm{T}m_\mathrm{T}}^{[\mathrm{\Gamma }]\nu }`$ and $`\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}`$. This computation can be reduced by taking into account their symmetries, as described in Appendix B, which imply that $`\nu |m|+|\mathrm{}_\mathrm{T}|`$ is odd, $`|m|=|m_\mathrm{T}|`$, $`|\mathrm{}_\mathrm{T}|+|m_\mathrm{T}|\nu `$, and $`|\mathrm{}_\mathrm{T}|+|m_\mathrm{T}|\nu \mathrm{mod}(2)`$. The computation can be performed analytically with the use of symbolic computation software such as Mathematica. In the case of the lens space $`L(17,5)`$, the computation up to $`\nu =43`$ and $`\nu =55`$ takes $`3`$ and $`12`$ hours respectively on a $`1\text{Ghz}`$ CPU with a negligible amount of memory.

### IV.2 Expected results

The three main effects that are expected on a CMB map computed in a multi-connected topology are (i) the appearance of $`\mathrm{}`$$`\mathrm{}^{}`$ and $`m`$$`m^{}`$ correlations reflecting the breakdown of global isotropy, (ii) the existence of a cutoff in the CMB angular power spectrum on large angular scales (low $`\mathrm{}`$), and (iii) the existence of pattern correlations such as pairs of circles where the temperature fluctuations are strongly correlated as they represent the intersection of the last scattering surface with itself and therefore show the temperature of the same emission region from different directions. Note that the effects (ii) and (iii) will show up only if the topological scale is smaller than the radius of the last scattering surface while the effect (i) may be present even if the topological scale is a bit larger than the diameter of the last scattering surface.

So far, the main constraints that have been given on multi-connected topologies come from the absence of a cutoff at large angular scales in the COBE spectrum. This gives strong constraints on the minimal size of the topology as the cutoff is given by the angular size of the torus projected on the last scattering surface. However, as previously discussed, this cutoff in the “true” temperature fluctuations can be compensated, at least partially, by an integrated Sachs-Wolfe effect which arises, for example, when the cosmological constant is large.

The third, and up-to-now never computed, effect of a multi-connected topology is the appearance of pairs of circles which are correlated in temperature. This correlation, however, is not perfect. It would be perfect if the temperature fluctuation were a pure scalar function on the last scattering surface $`2`$-sphere around us which would be the case only if (i) the temperature anisotropies were given only by the Sachs-Wolfe effect \[first term of Eq. (20)\] and (ii) the last scattering surface were infinitely thin.

It is well-known that the temperature fluctuations observed in a given direction are in fact a combination of several effects: first, one has the intrinsic temperature fluctuations of the emitting region, which is eventually affected by a gravitational redshift. These two contributions form the so-called Sachs-Wolfe effect \[first term in the righthand side of Eq. (20)\]. Second, if the emission region is not at rest with respect to the observer, one will observe some apparent temperature fluctuations which in fact result from a Doppler shift \[second term in the righthand side of Eq. (20)\]. Third, several events can alter the photons energy and trajectory while travelling toward us. In particular, they can be slightly disturbed from their trajectory (lensing) and, more importantly, they can exchange energy when they cross time-varying potential wells. This last effect is usually referred to as the integrated Sachs-Wolfe effect \[third term in the righthand side of Eq. (20), see Fig. 9 below\]. Obviously, the Sachs-Wolfe effect is a scalar quantity that depends only on the emission region. Therefore, it should be the same whatever the direction of observation. By contrast, the Doppler effect will explicitly depend on the direction of observation. If one observes two directions which correspond to the same point of the last scattering surface and which form a small angle, then one expects that the Doppler contribution will be almost the same. If the matching points are $`90`$ degrees from each other, then one expects on average no correlation at all, whereas the Doppler effect between two antipodal points will become anti-correlated. Finally, since photons originating from the same emission region but observed from different directions will travel through different regions of space, they will undergo different integrated Sachs-Wolfe effects, so that no significant correlations are expected from this effect, which is therefore considered a noise term for our purposes.

Actually, one aim of this work is precisely to compute the typical amount of correlation one can expect on pairs of circles. Note that this correlation is likely to depend on scale: on large scales, one should be annoyed by the late integrated Sachs-Wolfe effect; between the Sachs-Wolfe plateau and the first Doppler peak (and, at a lesser extent, at every dip between two Doppler peaks), the Doppler effect dominates; at the first peak, there is usually (especially when the matter content is low) a significant contribution of the early integrated Sachs-Wolfe effect (see Fig. 9); at very small scales, one feels the finite width of the last scattering surface (see below), etc. Also, as explained above, the relative position of the circles will play a role because of the Doppler contribution. It is therefore interesting to look at the best way to find matching circles on a realistic CMB map. We leave this important point to future work .

The second (and probably less important) effect that reduces the correlation between the circles is the finite width of the last scattering surface. As far as we know, this effect has not yet been carefully analyzed. It plays a role when one looks at fluctuations on scales smaller than the projected width of the last scattering surface. In this case when looking in a given direction, one picks up fluctuations which are situated “on one side” of the last scattering surface, but for pairs of circles, one sees opposite sides of the last scattering surface. On larger scales, the effect is negligible as one averages temperature fluctuations on regions much larger than the thickness of the last scattering surface.

## V Results

We now outline some of the results we have already obtained from our simulations. The main aim of this Section is to provide a series of tests to check our simulations. A more detailed analysis of the structure of the correlation matrix $`C_\mathrm{}m^\mathrm{}^{}m^{}`$ as well as a search for accurate tests to detect the topology are left for future work .

### V.1 Flat case: cubic torus

In all the simulations we performed, we have considered a flat $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, a Hubble parameter of $`H_0100h\text{km}\text{s}^1\text{Mpc}^1`$ with $`h=0.62`$, a baryon density $`\omega _\mathrm{b}\mathrm{\Omega }_\mathrm{b}h^2=0.019`$ and a spectral index $`n_\mathrm{S}=1`$. With this choice of cosmological parameters, the Hubble radius is $`R_H4.8\text{Gpc}`$, the “horizon” radius (under the hypothesis of a radiation dominated universe at early times) is $`R_h15.6\text{Gpc}`$, and the radius of the last scattering surface is $`R_{\mathrm{LSS}}15.3\text{Gpc}`$. The volume of the observable universe is therefore $`V_{\mathrm{obs}}15\times 10^3\text{Gpc}^3`$.

Let us first compare the $`C_{\mathrm{}}`$ in the simply connected topology to the $`C_{\mathrm{}}`$ in a torus (Fig. 3). As expected, we see a cutoff at some angular size which corresponds to the angular size of the torus on the last scattering surface. This corresponds to the multipole

$$\mathrm{}_c\frac{2\pi R_{\mathrm{LSS}}}{L},$$
(75)
where $`L`$ is the length of the cubic torus’ fundamental domain . Note that even when the torus is larger than the size of the observable universe, the spectrum exhibits a loss of power on large scales. This is because the Harrison-Zel’dovich spectrum exhibits a significant amount of power at large scales (by definition, it is scale invariant), and in practice, the modes that contribute to the quadrupole of the CMB anisotropies can be as large as ten times the size of the observable universe (the exact number depends mostly on the spectral index and on the amplitude of the integrated Sachs-Wolfe term). Therefore, this leaves hope to detect the topology “beyond the horizon” where the circles method would fail.

It is not easy to predict the amplitude of the power at scales larger than the cutoff because it depends mostly on the amplitude of the integrated Sachs-Wolfe effect, which is difficult to estimate even when the topology is multiconnected. Another consequence of a multiconnected topology are oscillations in the spectrum. These come both from the fact that there is a sharp cutoff in the spectrum (which causes oscillations in Fourier/Legendre space) and that the spectrum is “spiky” on large scales. Should we consider a simply connected universe with a cutoff at some scales, then the corresponding $`C_{\mathrm{}}`$ would be less irregular. Finally, note that on small angular scales, the spectrum tends to behave as in the simply connected case, but computing this scale remains an open problem at the moment.

So far, we have considered only the $`C_{\mathrm{}}`$, which represent only some average of the diagonal part of the correlation matrix. The true diagonal part of the correlation matrix is given by the $`C_\mathrm{}m^\mathrm{}m`$ which represent the variance of the $`a_\mathrm{}m`$. An example of their behavior is shown in Fig. 4 for the same topologies as in Fig. 3. Several features appear on this figure. First and most importantly, the dispersion in the variance of the $`a_\mathrm{}m`$ at fixed $`\mathrm{}`$ is very large. It appears that it is maximal at the cutoff scale $`\mathrm{}_c`$, where the dispersion in the variance of the $`a_\mathrm{}m`$ can be as large as two orders of magnitude. This dispersion slowly decays at larger multipoles, where one “tends” (in the sense of observable quantities) towards the simply connected case, and surprisingly also decays at scales larger than the cutoff. With the hypothesis of a multi-connected universe, this dispersion would be (incorrectly) interpreted as non Gaussianity. Since at present no non Gaussianity or anisotropy was observed in the data, this allows new constraints of the size of the fundamental domain.

We did not find any convenient way to represent the off-diagonal terms of the correlation matrix. We therefore switch to showing and analyzing some realizations corresponding to the numerically computed correlation matrix. In what follows, we have fixed the size of the torus to $`L=2R_H`$. We therefore have $`N16`$ copies of the torus in the observable universe, and from Eq. (75) this corresponds to $`\mathrm{}_c9`$, as can be checked in Fig. 3. Although such a model is now excluded by the data , we analyze it in detail mostly for pedagogical purposes (and also because the computing time scales as $`L^3`$). In the case where the torus, or more generally the fundamental domain, is smaller than the last scattering surface, one expects to see pairs of circles where the temperature is correlated . Seeing these circles at their expected position is therefore the most crucial test of the procedure outlined in Secs. II,III,IV. As already announced, the aim here is not to derive a detailed procedure to detect these circles, but to check our algorithm and to explore some of the properties of these matching circles. Here, we have $`2R_{\mathrm{LSS}}/L=3.17\sqrt{10}`$. One therefore expects to have $`61`$ pairs of circles and $`12`$ pairs of points having correlated temperature777The circles correspond to the intersection of the last scattering surface with translates of the form $`L(n_1,n_2,n_3)`$, where $`(n_1,n_2,n_3)=(0,0,1)`$, $`(0,0,2)`$, $`(0,0,3)`$, $`(0,1,1)`$, $`(0,1,2)`$, $`(0,2,2)`$, $`(1,1,1)`$, $`(1,1,2)`$, and $`(1,2,2)`$ plus all permutations and sign changes among each triplet $`(n_1,n_2,n_3)`$. The pairs of points correspond to the case where the intersection between the last scattering surface and its translate reduces to almost a single point, as is the case for $`(n_1,n_2,n_3)=(0,1,3)`$ and its permutations and sign changes.. In order to see the circles, it is convenient to show the last scattering surface as a sphere seen from the outside and to look at its intersection with itself after a translation of the form $`L(n_1,n_2,n_3)`$, as shown in Figs. 58.

In the last Section, we pointed out that the correlation would depend on the amount of Sachs-Wolfe, Doppler and integrated Sachs-Wolfe effects. The decomposition of the temperature anisotropies both in the simply connected topology and in the toroidal cases are shown in Fig. 9. Note that we show only the relative amplitude of these effects and not their cross correlation.

We now turn to the correlation between pairs of circles, as introduced in .

If the topology is not known in advance, the relative position between matching circles can be arbitrary, so that in general the search for circles is a six parameter problem: two parameters for the center of the first circle, two more for the center of the second circle, one for their common radius, and one for their relative phase (i.e. the twist with which they are identified). In the case of a torus, the circles sit directly opposite each other on the sky (eliminating two parameters) and there is no twist (eliminating another parameter), so the problem reduces to a three parameters search888Another way to see that in a torus it’s a three parameter search is to visualize the situation in the universal covering space. Place one copy of the last scattering surface with its center at the origin, and imagine a translated copy with its center at some point $`(x,y,z)`$. Each choice of $`(x,y,z)`$ uniquely determines a circle of intersection (assuming $`0<x^2+y^2+z^2<R_{\mathrm{LSS}}^2`$), and conversely each pair of circles arises from exactly two points $`(x,y,z)`$ and $`(x,y,z)`$ and no others. Thus the point $`(x,y,z)`$ serves to parameterize the circle search.. We are not going to perform such a study, but rather focus on some features of matching circles in a toroidal universe.

A simple estimator for the correlation between pairs of circles that are horizontal with respect to the coordinate system is obviously

$$C(\theta _1,\theta _2)\frac{1}{2\pi }\frac{2\mathrm{\Theta }(\theta _1,\phi )\mathrm{\Theta }(\theta _2,\phi )}{\mathrm{\Theta }^2(\theta _1,\phi )+\mathrm{\Theta }^2(\theta _2,\phi )}d\phi ,$$
(76)
where we have set $`\mathrm{\Theta }\delta T/T`$. If the temperature fluctuations on a pair of circles are completely uncorrelated, then on average $`C=0`$. If they are completely correlated, then $`C=1`$, and if they are anti correlated, then $`C=1`$. Each $`a_\mathrm{}m`$ can be seen either as the theoretical expectation for a given model, or an observed quantity that can be measured from the CMB sky. In the first case, it represents a feature that one can expect from a given model, and in the second case it represents an estimator of some features predicted by the topology. Here, we shall concentrate on the observed $`C(\theta _1,\theta _2)`$ that we compute from simulated maps, first to check the validity of our procedure to compute the correlation matrix $`C_\mathrm{}m^\mathrm{}^{}m^{}`$, and second to convince us that it is possible to see the presence of matching circles using simple techniques (although we do not pretend that this method is optimal).

In principle, two matching circles have the same angular diameter, so that only the case $`\theta _2=\theta _1`$ is relevant, but we have chosen to leave $`\theta _2`$ as a free parameter to see to what extent uncorrelated circles might happen to seem correlated by chance.

We first show in Fig. 10 a few examples of the observed function $`C(\theta _1,\theta _2)`$ for a simply connected universe. As expected, the correlation is quite large when $`\theta _1\theta _2`$ because the circles are near each other and the real space correlation function $`C(\theta )`$ (the Legendre transform of the $`C_{\mathrm{}}`$) is not 0 when $`\theta 0`$. With our normalization of $`C(\theta _1,\theta _2)`$, one has

$$\underset{\theta _2\theta _1}{lim}C(\theta _1,\theta _2)=1,$$
(77)
which will of course remain valid when the topology is multi-connected. When the separation between $`\theta _1`$ and $`\theta _2`$ is large, one can neglect the correlation between the two circles, and the main contribution to $`C`$ comes from statistical fluctuations: it is always possible that two circles exhibit similar temperature patterns by chance. The variance of these statistical fluctuations is probably given by the number of independent pixels on the map and therefore by a combination of the scale at which the power spectrum is large and of the resolution of the map (here, $`\mathrm{}_{\mathrm{max}}=30`$). In any case, the amplitude of the largest statistical fluctuations of $`C(\theta _1,\theta _2)`$ gives an idea of the amplitude of the signal needed to detect a multiconnected topology999The signal threshold could therefore be reduced by performing the same analysis on a higher resolution map, but because the search for circles is in general a six-parameter problem, it might be necessary to search low resolution maps first to find likely candidates, and then search higher resolution maps to confirm them.. For the maps we have generated, the correlation reaches $`30\%`$ for a few pairs of falsely matched circles.

Figs. 11, 12, 13 and 14 show contour plots for several realizations of the correlation matrix in a torus universe. The torus aligns naturally with the coordinate system, so one expects correlated pairs of circles at

$$\theta _1=\theta _2=\pm \mathrm{arcsin}\left(\frac{nL}{2R_{\mathrm{LSS}}}\right),$$
(78)
for each positive integer $`n`$ such that the arcsine exists. For our choice of cosmological parameters we have $`R_{\mathrm{LSS}}=3.17`$ and $`L=2`$ (in units of the Hubble radius), giving

$$\theta _1=\theta _2=\pm 18^{},\pm 39^{},\pm 71^{}.$$
(79)
For these values of $`\theta _1`$ and $`\theta _2`$, one expects a perfect correlation for the Sachs-Wolfe contribution:

$$C_{\mathrm{SW}}=1.$$
(80)
This formula holds both when one considers $`C_{\mathrm{SW}}`$ as an ensemble average and when one considers a given realization of the density field since in both cases it follows from the fact that one sees the same region from different directions. These correlations appear clearly in Fig. 11 which considers only the Sachs-Wolfe contribution. In this case one would have even expected perfect correlations for the values of $`\theta _1`$ and $`\theta _2`$ given in (79). This is not what we have, but the reason for this is easy to understand: imposing $`\mathrm{\Theta }(\theta _1,\phi )=\mathrm{\Theta }(\theta _1,\phi )`$ in real space induces in Legendre space correlations at arbitrary large multipoles $`\mathrm{}`$. Here, for computational reasons, we were forced to truncate the correlation matrix at a rather low value of $`\mathrm{}`$, so the matching is significant but not perfect. It would presumably increase in higher resolution maps.

If one considers the Doppler contribution to the CMB anisotropies, the situation is somewhat different. As announced above, the correlation between two circles depends on their relative angle. More precisely, it is given by

$$C_{\mathrm{Dop}}=\frac{2(\widehat{𝐧}_1\widehat{𝐧})(\widehat{𝐧}_2\widehat{𝐧})}{(\widehat{𝐧}_1\widehat{𝐧})^2+(\widehat{𝐧}_2\widehat{𝐧})^2}_{\widehat{𝐧}},$$
(81)
where $`\widehat{𝐧}_1`$ and $`\widehat{𝐧}_2`$ are two constant unit vectors spanning an angle $`|\theta _1\theta _2|`$, and where the brackets denote an average over all the directions of the unit vector $`\widehat{𝐧}`$. After some manipulations, one obtains

$$C_{\mathrm{Dop}}=\mathrm{tan}\left(\frac{\pi }{4}\frac{|\theta _1\theta _2|}{2}\right).$$
(82)
Again, for the same reason as for the Sachs-Wolfe contribution, this formula holds both if one considers $`C_{\mathrm{Dop}}`$ as an ensemble average or if we consider a given realization of the density field. One recovers as expected that the correlation is $`1`$, $`0`$, $`1`$ for $`|\theta _1\theta _2|=0^{}`$, $`90^{}`$, $`180^{}`$, respectively. For the values of the angle given in Eq. (79), one obtains $`C=0.51`$, $`0.07`$, $`0.49`$, respectively. These are the results that we obtain qualitatively in Fig. 12, where no correlation at all is seen for the circles at $`\pm 39^{}`$, and positive (resp. negative) correlation is seen for the circles at $`\pm 71^{}`$ (resp. $`\pm 18^{}`$).

Finally, the correlations due to the integrated Sachs-Wolfe effect are shown in Fig. 13. As expected, no particular correlation is seen for the values of $`\theta _1`$ and $`\theta _2`$ of Eq. (79). The contour plots are however quite different from those of the Sachs-Wolfe and Doppler contributions. The reason is twofold. First, most of the power lies at the smallest multipoles. This translates into the fact that the contours are broad in the sense that they do not vary a lot on small intervals of $`\theta _1`$ and $`\theta _2`$. Second, the fact that most of the power is at large scales implies that a very small number of modes contributes to it (since we see a finite region of the universe), so there is a large cosmic variance that makes the statistical uncertainty very large (thus serendipitously similar temperature patterns on two unrelated circles are easily achieved here).

Combining all the contributions to the CMB anisotropies allows one to simulate realizations of the exact $`C(\theta _1,\theta _2)`$ as shown in Fig. 14. Since the Sachs-Wolfe contribution is dominant, the spikes are still clearly visible at their expected positions, but appear less prominent than in Fig. 11. As expected, it seems that the circles at $`\theta _{1,2}=\pm 39^{}`$ are slightly less correlated than the other two pairs because their Doppler contribution is not correlated, but this deserves a more careful analysis.

### V.2 Spherical case: lens spaces

Among the spherical spaces, the procedure presented above can be applied most easily to lens and prism spaces, because their eigenmodes are known explicitly. The eigenmodes are known analytically in toroidal coordinates (see Section III.2), and Appendix B shows how to convert them to spherical coordinates. In this section, we present some sample maps exhibiting the matching circles to demonstrate that the whole computational chain (computation of the modes and implementation in a CMB code) is working. A complete and detailed study, along the same lines as the study done for the cubic torus in the previous section, will be presented in a follow-up article.

As explained in Ref. , because our universe is almost flat, observational methods such as the circles method will typically detect only a cyclic subgroup of the holonomy group, so the universe “looks like a lens space” no matter what its true topology is. It follows that lens spaces are particularly interesting to capture the observational properties of multi-connected spherical spaces. In particular, we showed  that a cyclic factor $`Z_p`$ creates matching circles in the CMB only when $`\mathrm{\Omega }1>1/p^2`$ and that the second factor, if it exists, is in general undetectable.

Let us emphasize some differences with the torus case. First, concerning the eigenmodes, let us take the example of a lens space $`L(p,1)`$ of order $`p`$. For $`p=1`$ it reduces to $`𝐒^3`$ and for $`p=2`$ it reduces to projective space; more generally the index plays a role analogous to the size, $`L`$, of the torus in Euclidean space. The first non-zero eigenvalue is always $`\nu =2`$ and has a multiplicity $`3`$ for $`p>2`$ and $`9`$ otherwise. This constancy of the first eigenvalue contrasts sharply with the case of a cubic torus, for which the smallest eigenvalue scales as $`L^1`$. It can be understood by realizing that when $`p`$ increases the space is becoming smaller only in one direction and remains large in perpendicular directions.

The lens spaces $`L(p,1)`$ are globally homogeneous (like the torus) so that the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$ do not depend on the observer’s position (i.e., they are the same no matter where in the space you choose the basepoint). Thus neither the correlation matrix $`C_\mathrm{}m^\mathrm{}^{}m^{}`$ nor the positions of the matching circles depend on the basepoint. Unfortunately, this is not the case for a general lens space $`L(p,q)`$. For a general lens space, the coefficients $`\xi _{k\mathrm{}m}^{[\mathrm{\Gamma }]s}`$, the correlation matrix $`C_\mathrm{}m^\mathrm{}^{}m^{}`$ and the positions of the matching circles all depend on the observer’s position. For instance, the “canonical” choice of coordinates used in Sec. III.2.2 for the toroidal coordinate system puts the preferred symmetry axes in the $`(xy)`$ and $`(zw)`$ directions (i.e., the axis are the intersection of $`𝐒^3`$ with the $`(x,y)`$-plane and the $`(z,w)`$-plane, respectively, in four-dimensional Euclidean space). From a cosmological point of view, this is a poor choice, because the observer’s translated images are “atypically close”. For example in $`L(12,5)`$, which has cyclic factors $`Z_3`$ and $`Z_4`$, a generic observer will see three lines of four images each, but a nongeneric observer sitting on a symmetry axis will see a single line of $`12`$ images.

For a globally homogeneous space $`L(p,1)`$, the closest topological image is located at a distance $`\overline{\chi }=2\pi /p`$, so the topology is detectable just so $`p`$ isn’t too small. More precisely, the topology is potentially detectable if and only if $`\pi /p<\overline{\chi }_{\mathrm{LSS}}`$. For example, this implies that the topology is detectable for all $`p>10`$ if $`\mathrm{\Omega }_{\mathrm{tot}}110^2`$.

Circles match differently in a homogeneous lens space than in a tours. In a torus the circles match straight across because the holonomies are all pure translations. In a homogeneous lens space, by contrast, the holonomies are Clifford translations, and so the matching circles are still diametrically opposite but match with a twist that’s a multiple of $`2\pi /p`$, because Clifford translations twist and translate of the same amount.

Fig. 15 shows a CMB map with resolution $`\mathrm{}_{\mathrm{max}}=30`$ for the lens space $`L(21,1)`$ considering the Sachs-Wolfe term only. A far more detailed discussion about CMB anisotropies in lens spaces will appear elsewhere .

## VI Discussion and conclusions

This articles describes the implementation of topology in CMB codes and gives explicitly the required tools to perform such an implementation in flat and spherical spaces. As emphasized in the Introduction, these two cases are observationally the most relevant for an almost flat universe.

Examples of simulated maps were given in the two cases. Here we presented only low resolution maps due to the computational time limitation but higher resolution maps will be presented elsewhere. It was checked that the expected topological correlations (the matched circles) were present, confirming the quality of our simulations.

Our method relies on the computation of the correlation matrix of the coefficients of the decomposition of the temperature fluctuation in spherical harmonics. This matrix encodes all the topological information. We emphasize that, due to the breakdown of global isotropy, this matrix is not purely diagonal. This also offers a working example to construct tests for the detection of deviation from global isotropy.

We have illustrated the influence of different effects that will tend to blur these patterns and affect the perfect circle matching, namely the Doppler effect and the integrated Sachs-Wolfe effect. We also considered the effect of the thickness of the last scattering surface, but found it to be negligible on the scales considered here. A more detailed quantitative analysis of these effects on the detectability of the topological signal is left for future studies .

A complete investigation of the detectability of the topology in coming CMB data requires the construction of reliable simulation tools. Besides the quantification of the amplitude of the effects cited above, one would also need to include all other observational effects such as instrumental noise, foreground contamination, etc. The present work paves the way to all these essential studies.

## Acknowledgments

We want to thank Nabila Aghanim, Francis Bernardeau, Francois Bouchet, Gilles Esposito-Farèse, Simon Prunet for providing computational support, Dorian Goldfeld for his help in evaluating some sums of Appendix B, and Jean-Pierre Luminet for very careful reading of this manuscript and numerous discussions. J.W. thanks the MacArthur Foundation for its support. J.-P.U. thanks the University of Barcelona for hospitality while a part of this work was performed. Part of this work was achieved while A.R. was a CMBnet postdoc in the Département de Physique Théorique of Geneva University.

## Appendix A Eigenmodes of constant curvature three-dimensional spaces

This appendix follows the work by Abbott and Schaeffer  and Harrison  and borrows heavily from Appendix A of Ref. . It summarizes, without proof, the explicit forms of the scalar harmonic functions solutions of the Helmholtz equation (2).

It is convenient to factor the eigenfunctions into radial and angular functions as

$$𝒴_{k\mathrm{}m}^{[X]}(\chi ,\theta ,\phi )=R_k\mathrm{}^{[X]}(\chi )Y_{\mathrm{}}^m(\theta ,\phi ),$$
(83)
with $`Y_{\mathrm{}}^m(\theta ,\phi )`$ being the spherical harmonics. The associated eigenvalues are $`\kappa _k^2=k^2K`$, with

$`K<0`$ $``$ $`k[0,\mathrm{}[\text{or}ik[0,\sqrt{|K|}],`$ (84)
$`K=0`$ $``$ $`k[0,\mathrm{}[,`$ (85)
$`K>0`$ $``$ $`k=(\nu +1)\sqrt{K},\nu 𝐍.`$ (86)
With the normalization

$$𝒴_{k\mathrm{}m}^{[X]}(\chi ,\theta ,\phi )𝒴_{k^{}\mathrm{}^{}m^{}}^{[X]}{}_{}{}^{}(\chi ,\theta ,\phi )\text{s}_K^2(\chi )d\chi d\mathrm{\Omega }=\frac{1}{k^2}\delta ^\mathrm{D}(kk^{})\delta _{\mathrm{}\mathrm{}^{}}\delta _{mm^{}},$$
(87)
where $`\text{s}_K(\chi )`$ is defined in Eq. (7), the normalized radial functions take the form

$`R_k\mathrm{}^{[𝐇^3]}(\chi )`$ $`=`$ $`\left({\displaystyle \frac{N_k\mathrm{}}{k\text{s}_K(\chi )}}\right)^{1/2}P_{1/2+i\omega }^{1/2\mathrm{}}\left(\mathrm{cosh}(\sqrt{|K|}\chi )\right),`$ (88)
$`R_k\mathrm{}^{[𝐑^3]}(\chi )`$ $`=`$ $`\left({\displaystyle \frac{2}{\pi }}\right)^{1/2}j_{\mathrm{}}(k\chi ),`$ (89)
$`R_k\mathrm{}^{[𝐒^3]}(\chi )`$ $`=`$ $`\left({\displaystyle \frac{M_k\mathrm{}}{k\text{s}_K(\chi )}}\right)^{1/2}P_{1/2+\nu }^{1/2\mathrm{}}\left(\mathrm{cos}(\sqrt{K}\chi )\right)`$ (90)
with

$$\omega =k/\sqrt{|K|},\nu =k/\sqrt{K}1.$$
(91)
In the case of spatially hyperbolic spaces, this normalization is valid only for sub-curvature modes, and for the super-curvature modes ($`ik[0,\sqrt{|K|}]`$) the radial function is obtained by analytic continuation (see Ref.  for details). The two numerical coefficients are given by

$`N_k\mathrm{}`$ $``$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left(\omega ^2+n^2\right),`$ (92)
$`M_k\mathrm{}`$ $``$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left((\nu +1)^2n^2\right),M_k\mathrm{}=0\text{if}\mathrm{}>\nu .`$ (93)
For any function, we can perform the mode decomposition

$$f(𝐱)=\underset{\mathrm{},m}{}\left\{\begin{array}{c}k^2dk\hfill \\ K^{3/2}_{\nu =2}^{\mathrm{}}(\nu +1)^2\hfill \end{array}\right\}f_{k\mathrm{}m}𝒴_{k\mathrm{}m}^{[X]}f_{k\mathrm{}m}=f(𝐱)𝒴_{k\mathrm{}m}^{[X]}{}_{}{}^{}\sqrt{\gamma }\mathrm{d}^3𝐱,$$
(94)
choosing the sum or the integral according to whether the universal covering space is compact or not. The symbol $`\sqrt{\gamma }`$ stands for the square root of the determinant of the spatial metric. In the case of spatially hyperbolic spaces, the super-curvature modes add a term to this mode expansion, namely $`_0^1|k^2|\mathrm{d}(ik)_{\mathrm{},m}f_{k\mathrm{}m}𝒴_{k\mathrm{}m}^{[𝐇^3]}`$; see Ref.  for details.

In the spherical case, one can however find a solution of the Helmholtz equation (2) which does not involve Legendre functions. The radial part of the Helmholtz equation reduces, after setting $`𝒴_{k\mathrm{}m}^{[𝐒^3]}=R_k\mathrm{}^{[𝐒^3]}(\chi )Y_{\mathrm{}}^m(\theta ,\phi )`$, to

$$\frac{1}{\text{s}_K^2(\chi )}\frac{\mathrm{d}}{\mathrm{d}\chi }\left(\text{s}_K^2(\chi )\frac{\mathrm{d}}{\mathrm{d}\chi }R_k\mathrm{}^{[𝐒^3]}\right)+\left[(k^2K)\frac{\mathrm{}(\mathrm{}+1)}{\text{s}_K^2(\chi )}\right]R_k\mathrm{}^{[𝐒^3]}=0.$$
(95)
It is obviously much more convenient to work in a coordinate system where the curvature $`K`$ reduces to $`1`$. In terms of the dimensionless radial variable $`\overline{\chi }`$ defined in Eq. (8), the Helmholtz equation then reduces to

$$\frac{1}{\mathrm{sin}^2\overline{\chi }}\frac{\mathrm{d}}{\mathrm{d}\overline{\chi }}\left(\mathrm{sin}^2\overline{\chi }\frac{\mathrm{d}}{\mathrm{d}\overline{\chi }}R_k\mathrm{}^{[𝐒^3]}\right)+\left[\nu (\nu +2)\frac{\mathrm{}(\mathrm{}+1)}{\mathrm{sin}^2\overline{\chi }}\right]R_k\mathrm{}^{[𝐒^3]}=0.$$
(96)
Note that this is a second order equation and that only one of the two independent solutions is well behaved at the origin, so the radial functions are completely determined once the normalization has been chosen. After setting $`R_k\mathrm{}^{[𝐒^3]}=(\mathrm{sin}\overline{\chi })^{\mathrm{}}f_\nu \mathrm{}`$, it can be checked that it reduces to Eq. (133), the solution of which is simply given in terms of ultraspherical Gegenbauer polynomials as $`f_\nu \mathrm{}=A_\nu \mathrm{}C_\nu \mathrm{}^{\mathrm{}+1}`$. The normalization condition (87) implies, using the integral relation (134), that

$$A_\nu \mathrm{}=\frac{2^{\mathrm{}+1/2}}{\nu +1}\mathrm{}!\sqrt{\frac{\nu +1}{\pi }}\sqrt{\frac{(\nu \mathrm{})!}{(\nu +\mathrm{}+1)!}}.$$
(97)
Expressing the spherical harmonics in terms of Gegenbauer polynomials by means of Eq. (131), one ends up with an expression of the eigenmodes in terms of Gegenbauer polynomials only as

$$𝒴_{k\mathrm{}m}^{[𝐒^3]}=A_\nu \mathrm{}\stackrel{~}{A}_\mathrm{}m\mathrm{sin}^{|m|}\theta \mathrm{sin}^{\mathrm{}}\overline{\chi }C_\nu \mathrm{}^{\mathrm{}+1}(\mathrm{cos}\overline{\chi })C_{\mathrm{}|m|}^{|m|+1/2}(\mathrm{cos}\theta )\mathrm{e}^{im\phi },$$
(98)
with

$$\stackrel{~}{A}_\mathrm{}m=\zeta _m(2|m|1)!!\sqrt{\frac{(2\mathrm{}+1)}{4\pi }\frac{(\mathrm{}|m|)!}{(\mathrm{}+|m|)!}}$$
(99)
with $`\zeta _m`$ given by Eq. (132) and we used the notation $`(2|m|1)!!=(2|m|1)(2|m|3)\mathrm{}1`$.

## Appendix B Change of basis between toroidal and spherical coordinates

Sec. III.2.2 found the eigenmodes of lens and prism spaces in toroidal coordinates and converted them to spherical coordinates. In this appendix, we give the expression for the matrix $`\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}`$ necessary to perform the change of basis.

The spherical coordinate system, as used in Eq. (6), is related to the embedding of the $`3`$-sphere $`𝐒^3`$ in four-dimensional Euclidean space by

$`x`$ $`=`$ $`\mathrm{cos}\overline{\chi },`$ (100)
$`y`$ $`=`$ $`\mathrm{sin}\overline{\chi }\mathrm{cos}\theta ,`$ (101)
$`z`$ $`=`$ $`\mathrm{sin}\overline{\chi }\mathrm{sin}\theta \mathrm{cos}\phi ,`$ (102)
$`w`$ $`=`$ $`\mathrm{sin}\overline{\chi }\mathrm{sin}\theta \mathrm{sin}\phi ,`$ (103)
with

$`0\overline{\chi }\pi ,`$ (104)
$`0\theta \pi ,`$ (105)
$`0\phi 2\pi .`$ (106)
The (complex) coefficients $`\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}`$ characterizing the change of basis are defined by

$$\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}=(\nu +1)^2𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}(\chi _\mathrm{T},\theta _\mathrm{T},\phi _\mathrm{T})𝒴_{k\mathrm{}m}^{[𝐒^3]}{}_{}{}^{}(\overline{\chi },\theta ,\phi )(\mathrm{sin}\overline{\chi })^2\mathrm{d}\overline{\chi }\mathrm{sin}\theta \mathrm{d}\theta \mathrm{d}\phi .$$
(107)
In this expression the integer $`\mathrm{}`$ ranges from $`0`$ to $`\nu `$ and $`m`$ ranges from $`\mathrm{}`$ to $`\mathrm{}`$, while $`\mathrm{}_\mathrm{T}`$ and $`m_\mathrm{T}`$ range from $`\nu `$ to $`\nu `$. To compute this integral, one needs:

1. to express the eigenmodes $`𝒬_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}`$ as polynomials in $`x`$, $`y`$, $`z`$ and $`w`$ (see Ref. ), and replace these rectangular coordinates by their expressions (100),
2. use the relations

$$\mathrm{cos}^{|\mathrm{}_\mathrm{T}|}\chi _\mathrm{T}\{\begin{array}{c}\mathrm{cos}(|\mathrm{}_\mathrm{T}|\theta )\hfill \\ \mathrm{sin}(|\mathrm{}_\mathrm{T}|\theta )\hfill \end{array}=\{\begin{array}{cc}\mathrm{}\left[(\mathrm{cos}\overline{\chi }+i\mathrm{sin}\overline{\chi }\mathrm{cos}\theta )^{|\mathrm{}_\mathrm{T}|}\right]\hfill & \mathrm{}_\mathrm{T}0,\hfill \\ \mathrm{}\left[(\mathrm{cos}\overline{\chi }+i\mathrm{sin}\overline{\chi }\mathrm{cos}\theta )^{|\mathrm{}_\mathrm{T}|}\right]\hfill & \mathrm{}_\mathrm{T}<0,\hfill \end{array}$$
(108)
and

$$\mathrm{sin}^{|m_\mathrm{T}|}\chi _\mathrm{T}\{\begin{array}{c}\mathrm{sin}(|m_\mathrm{T}|\phi )\hfill \\ \mathrm{cos}(|m_\mathrm{T}|\phi )\hfill \end{array}=(\mathrm{sin}\overline{\chi }\mathrm{sin}\theta )^{|m_\mathrm{T}|}\{\begin{array}{cc}\mathrm{cos}(|m_\mathrm{T}|\phi )\hfill & m_\mathrm{T}0,\hfill \\ \mathrm{sin}(|m_\mathrm{T}|\phi )\hfill & m_\mathrm{T}<0,\hfill \end{array}$$
(109)
3. develop the Jacobi polynomial appearing in Eq. (57) by using Eq. (133) and with

$$\mathrm{cos}2\chi _\mathrm{T}1=2\mathrm{sin}^2\overline{\chi }\mathrm{sin}^2\theta .$$
(110)

This leads, after an easy integration on $`\phi `$, to the somewhat heavy expressions involving two sums \[arising from the development of the Jacobi polynomials and the power in Eq. (108)\] and an integral over $`\theta `$ and $`\overline{\chi }`$,

$$\alpha _{\nu \mathrm{}m\mathrm{}_\mathrm{T}m_\mathrm{T}}=(\nu +1)^2\pi A_\nu \mathrm{}\stackrel{~}{A}_\mathrm{}mB_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}C_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}\{\begin{array}{c}\mathrm{}\hfill \\ \mathrm{}\hfill \end{array}\left[𝒥(\nu ;\mathrm{},m;\mathrm{}_\mathrm{T},m_\mathrm{T})\right]\times \{\begin{array}{c}\left(\delta _{m|m_\mathrm{T}|}+\delta _{m|m_\mathrm{T}|}\right)\hfill \\ i\left(\delta _{m|m_\mathrm{T}|}\delta _{m|m_\mathrm{T}|}\right)\hfill \end{array}.$$
(111)
Here, the first and second line of the first brace are for $`\mathrm{}_\mathrm{T}0`$ and $`\mathrm{}_\mathrm{T}<0`$, respectively \[see Eq. (108)\], the first and second line of the second brace are for $`m_\mathrm{T}0`$ and $`m_\mathrm{T}<0`$, respectively \[see Eq. (109)\], and the numerical coefficient $`C_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}`$ is given by

$$C_{\nu \mathrm{}_\mathrm{T}m_\mathrm{T}}=\frac{\mathrm{\Gamma }(|m_\mathrm{T}|+d+1)}{d!\mathrm{\Gamma }(|m_\mathrm{T}|+|\mathrm{}_\mathrm{T}|+d+1)}.$$
(112)
The function $`𝒥`$ is explicitly given by

$`𝒥(\nu ;\mathrm{},m;\mathrm{}_\mathrm{T},m_\mathrm{T})`$ $`=`$ $`{\displaystyle \underset{q=0}{\overset{d}{}}}\left(\begin{array}{c}d\\ q\end{array}\right)(1)^q{\displaystyle \frac{\mathrm{\Gamma }(|m_\mathrm{T}|+|\mathrm{}_\mathrm{T}|+d+q+1)}{\mathrm{\Gamma }(|m_\mathrm{T}|+q+1)}}`$ (119)
$`\times {\displaystyle \underset{j=0}{\overset{|\mathrm{}_\mathrm{T}|}{}}}i^j\left(\begin{array}{c}|\mathrm{}_\mathrm{T}|\\ j\end{array}\right)I(j,2q+|m|+|m_\mathrm{T}|;\mathrm{}|m|,|m|+{\displaystyle \frac{1}{2}})`$
$`\times I(|\mathrm{}_\mathrm{T}|j,2q+j+\mathrm{}+|m_\mathrm{T}|+1;\nu \mathrm{},\mathrm{}+1).`$
Note that the index $`j`$ is even when $`\mathrm{}_\mathrm{T}0`$ and odd when $`\mathrm{}_\mathrm{T}<0`$. This quantity involves only two sums, once the quantity $`I(p,q;n,\alpha )`$, defined by

$$I(p,q;n,\alpha )_1^1x^p(1x^2)^{q/2}C_n^\alpha (x)dx,$$
(120)
is known. Using the expression (136) for the Gegenbauer polynomials in term of hypergeometric functions and the integral (138), it can be shown that

$`I(p,q;2m,\alpha )`$ $`=`$ $`{\displaystyle \frac{(1)^m}{2(\alpha +m)}}{\displaystyle \frac{B(q/2+1,(p+1)/2)}{B(\alpha ,m+1)}}\left[1+(1)^p\right]{}_{3}{}^{}F_{2}^{}(m,m+\alpha ,{\displaystyle \frac{p+1}{2}};{\displaystyle \frac{1}{2}},{\displaystyle \frac{q+p+3}{2}};1),`$ (121)
$`I(p,q;2m+1,\alpha )`$ $`=`$ $`(1)^m{\displaystyle \frac{B(q/2+1,(p+2)/2)}{B(\alpha ,m+1)}}\left[1(1)^p\right]{}_{3}{}^{}F_{2}^{}(m,m+\alpha +1,{\displaystyle \frac{p+2}{2}};{\displaystyle \frac{3}{2}},{\displaystyle \frac{q+p+4}{2}};1),`$ (122)
where $`B`$ is the Euler Beta function. It follows directly from these expressions that $`𝒥(\nu ;\mathrm{},m;\mathrm{}_\mathrm{T},m_\mathrm{T})=0`$ when $`\nu |m|+|\mathrm{}_\mathrm{T}|`$ is odd.

## Appendix C Some properties of some special functions

This appendix gathers some useful relations used in the article, to make the article more self-contained.

The spherical harmonics $`Y_{\mathrm{}}^m`$ are related to the associated Legendre polynomials $`P_{\mathrm{}}^m`$ by (see Eq. (5.2.1) of Ref. )

$$Y_{\mathrm{}}^m(\theta ,\phi )=\sqrt{\frac{2\mathrm{}+1}{4\pi }\frac{(\mathrm{}m)!}{(\mathrm{}+m)!}}P_{\mathrm{}}^m(\mathrm{cos}\theta )\mathrm{e}^{im\phi }.$$
(123)
They satisfy the conjugation relation (Eq. (5.4.1) of Ref. )

$$Y_{\mathrm{}}^m{}_{}{}^{}(\theta ,\phi )=(1)^mY_{\mathrm{}}^m(\theta ,\phi )=Y_{\mathrm{}}^m(\theta ,\phi ),$$
(124)
the normalization (Eq. (5.6.1) of Ref. )

$$_0^{2\pi }\mathrm{d}\phi _0^\pi \mathrm{sin}\theta \mathrm{d}\theta Y_{\mathrm{}}^m{}_{}{}^{}(\theta ,\phi )Y_{\mathrm{}^{}}^m^{}(\theta ,\phi )=\delta _{\mathrm{}\mathrm{}^{}}\delta _{mm^{}},$$
(125)
the closure relation (Eq. (5.2.2) of Ref. )

$$\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}Y_{\mathrm{}}^m(\theta ,\phi )Y_{\mathrm{}}^m{}_{}{}^{}(\theta ^{},\phi ^{})=\delta ^\mathrm{D}(\mathrm{cos}\theta \mathrm{cos}\theta ^{})\delta ^\mathrm{D}(\phi \phi ^{}),$$
(126)
and the addition theorem (Eq. (5.17.2.9) of Ref. )

$$\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}Y_{\mathrm{}}^m(\theta ,\phi )Y_{\mathrm{}}^m{}_{}{}^{}(\theta ^{},\phi ^{})=\frac{2\mathrm{}+1}{4\pi }P_{\mathrm{}}(\mathrm{cos}\alpha ),$$
(127)
where $`\alpha `$ is the angle between the two directions $`(\theta ,\phi )`$ and $`(\theta ^{},\phi ^{})`$, and $`P_{\mathrm{}}`$ is the Legendre polynomial. The Fourier transform of the spherical harmonics is given by (Eq. (5.9.2.6) of Ref. )

$$_0^{\mathrm{}}r^2dr_0^{2\pi }d\phi _0^\pi \mathrm{sin}\theta \mathrm{d}\theta \frac{\mathrm{e}^{i𝐤.𝐫}}{(2\pi )^{3/2}}j_{\mathrm{}}(k^{}r)Y_{\mathrm{}}^m(\theta ,\phi )=\sqrt{\frac{2}{\pi }}i^{\mathrm{}}\frac{\delta ^\mathrm{D}(k^{}k)}{k^2}Y_{\mathrm{}}^m(\theta _𝐤,\phi _𝐤),$$
(128)
where $`j_{\mathrm{}}`$ is a spherical Bessel function, from which it follows that (Eq. (5.17.3.14) of Ref. )

$$\mathrm{e}^{i𝐤.𝐫}=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}(2\mathrm{}+1)i^{\mathrm{}}j_{\mathrm{}}(kr)P_{\mathrm{}}(\mathrm{cos}\theta _{𝐤,𝐫}),$$
(129)
and (Eq. (5.17.4.18) of Ref. )

$$\delta ^\mathrm{D}(𝐫_1𝐫_2)=\frac{\delta ^\mathrm{D}(r_1r_2)}{r_1^2}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\frac{2\mathrm{}+1}{4\pi }P_{\mathrm{}}(\mathrm{cos}\theta _{12}).$$
(130)
The spherical harmonics can also be expressed in terms of Gegenbauer polynomials as (Eq. (5.2.6.39c) of Ref. )

$$Y_{\mathrm{}}^m(\theta ,\phi )=\zeta _m\mathrm{e}^{im\phi }\sqrt{\frac{2\mathrm{}+1}{4\pi }}\sqrt{\frac{(\mathrm{}|m|)!}{(\mathrm{}+|m|)!}}(2|m|1)!!(\mathrm{sin}\theta )^{|m|}C_{\mathrm{}|m|}^{|m|+1/2}(\mathrm{cos}\theta )$$
(131)
with $`\zeta _m`$ defined by

$$\zeta _m=\{\begin{array}{cc}(1)^m\hfill & m>0\hfill \\ 1\hfill & m0\hfill \end{array}.$$
(132)
The ultraspherical (or Gegenbauer) polynomials $`C_n^\alpha `$ are solutions of the differential equation (Eq. (22.6.5) of Ref. )

$$(1x^2)y^{\prime \prime }(2\alpha +1)y^{}+n(n+2\alpha )y=0,$$
(133)
and they satisfy the normalization condition (Eq. (7.313) of Ref. )

$$_1^1(1x^2)^{\alpha 1/2}\left[C_n^\alpha (x)\right]^2dx=\frac{\pi 2^{12\alpha }\mathrm{\Gamma }(2\alpha +n)}{n!(n+\alpha )\left[\mathrm{\Gamma }(\alpha )\right]^2},$$
(134)
if $`\mathrm{}(\alpha )>1/2`$.

The Jacobi polynomials are given by (Eq. (22.3.2) of Ref. )

$$P_n^{(\alpha ,\beta )}(x)=\frac{\mathrm{\Gamma }(\alpha +n+1)}{n!\mathrm{\Gamma }(\alpha +\beta +n+1)}\underset{j=0}{\overset{n}{}}\left(\begin{array}{c}n\\ j\end{array}\right)\frac{\mathrm{\Gamma }(\alpha +\beta +n+j+1)}{2^j\mathrm{\Gamma }(\alpha +j+1)}(x1)^j,$$
(135)
under the conditions $`\alpha >1`$ and $`\beta >1`$. Interestingly, the Gegenbauer polynomials can be expressed in terms of hypergeometric functions as (Eqs. (8.932.2,8.932.3) of Ref. )

$`C_{2n}^\lambda (x)`$ $`=`$ $`{\displaystyle \frac{(1)^n}{(\lambda +n)B(\lambda ,n+1)}}F(n,\lambda +n;1/2;x^2),`$ (136)
$`C_{2n+1}^\lambda (x)`$ $`=`$ $`{\displaystyle \frac{(1)^n}{B(\lambda ,n+1)}}(2x)F(n,\lambda +n+1;3/2;x^2),`$ (137)
which satisfies the integral property (Eqs. (7.513) of Ref. )

$$_0^1x^{s1}(1x^2)^\nu F(n,a;b;x^2)dx=\frac{1}{2}B(\nu +1,s/2){}_{3}{}^{}F_{2}^{}(n,a,s/2;b,\nu +1+s/2;1)$$
(138)
if $`\mathrm{}(s)>0`$ and $`\mathrm{}(\nu )>1`$.