# Optimality of the Bussgang Linear MMSE Channel Estimator for MIMO Systems
with 1-Bit ADCs

Minhua Ding,1 Italo Atzeni,1 Antti Tölli,1 and A. Lee Swindlehurst2
This work was supported by the Research Council of Finland (336449 Profi6,
346208 6G Flagship, 348396 HIGH-6G, and 357504 EETCAMD), by the European
Commission (101095759 Hexa-X-II), and by the U. S. National Science Foundation
(CCF-2225575). 1Centre for Wireless Communications, University of Oulu,
Finland
2Department of Electrical Engineering and Computer Science, University of
California, Irvine, USA
E-mail: {minhua.ding, italo.atzeni<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we study the optimality of the Bussgang linear minimum mean
squared error (BLMMSE) channel estimator for multiple-input multiple-output
systems with 1-bit analog-to-digital converters. We compare the BLMMSE with
the optimal minimum mean squared error (MMSE) channel estimator, which is
generally non-linear, and we develop a novel framework based on the orthant
probability of a multivariate normal distribution to compute the MMSE channel
estimate. Then, we analyze the equivalence of the MMSE and BLMMSE channel
estimators under specific assumptions on the channel correlation or pilot
symbols. Interestingly, the BLMMSE channel estimator turns out to be optimal
in several specific cases. Our study culminates with the presentation of a
necessary and sufficient condition for the BLMMSE channel estimator to be
optimal.

###### Index Terms:

1-bit ADCs, channel estimation, massive MIMO.

##  I Introduction

Massive multiple-input multiple-output (MIMO) technology is a key enabler of
current fifth-generation cellular networks and is envisioned to assume an even
more prominent role in future wireless systems [1]. The large number of
antennas in massive MIMO motivates the study of the impact of low-complexity,
energy-efficient hardware on the system performance. Low-resolution and 1-bit
analog-to-digital converters (ADCs) are considered a promising solution to
decrease the complexity and power consumption of fully digital massive MIMO
systems [2, 3, 4, 5]. Since coherent MIMO communications require instantaneous
channel state information for the beamforming design, the subject of channel
estimation with 1-bit ADCs has received special attention in recent years [3,
4, 6]. As a chief example, the Bussgang linear minimum mean squared error
(BLMMSE) channel estimator proposed in [3] uses the second-order statistics
[7] of 1-bit quantized signals and its linear structure has inspired
subsequent signal processing methods applied to, e.g., data detection and
beamforming design.

In general, the optimal minimum mean squared error (MMSE) channel estimator is
a non-linear function of the 1-bit quantized signal received during the pilot
transmission. Nonetheless, under certain assumptions, the MMSE channel
estimator may assume a linear structure and thus be equivalent to its BLMMSE
counterpart; in this setting, the BLMMSE channel estimator is clearly optimal
in terms of mean squared error (MSE). Recently, it was revealed in [6] that
the MMSE and BLMMSE channel estimators are equivalent in single-input single-
output systems with a single pilot symbol or with specifically chosen pilot
symbols in the noiseless case. Given the widespread use of the BLMMSE channel
estimator, it is of great importance to analyze its optimality in general
scenarios, and this is precisely the goal of this work.

In this paper, we develop a new framework based on the orthant probability of
a multivariate normal (MVN) distribution to compute the MMSE channel estimate
for MIMO systems with 1-bit ADCs. Building on this, we identify specific cases
where the BLMMSE channel estimator is optimal. Then, we present a closed-form
expression for the non-linear MMSE channel estimator in a single-input
multiple-output (SIMO) system with three receive antennas. Although far from
the massive MIMO regime that motivates the use of 1-bit ADCs, studying these
specific cases paves the way for establishing a necessary and sufficient
condition for the optimality of the BLMMSE channel estimator.111A more
detailed analysis, complemented by additional specific cases and full
derivations, can be found in the longer version of this paper [8].

Notation: $\otimes$ denotes the Kronecker product; $|\cdot|$ represents the
determinant of a square matrix or the absolute value of a scalar;
$\mathrm{vec}(\cdot)$ denotes vectorization; ${\mathrm{tr}}(\cdot)$ represents
the trace; $\mathrm{sgn}(\cdot)$ denotes the sign function;
$\mathrm{Re}(\cdot)$ and $\mathrm{Im}(\cdot)$ represent the real and imaginary
parts, respectively; $j=\sqrt{-1}$ is the imaginary unit;
$\mathtt{diag}{{(}}\mathbf{z})$ is a diagonal matrix with the elements of
$\mathbf{z}$ on its diagonal; $\widetilde{\mathtt{diag}{{}}}(\mathbf{Z})$ is a
diagonal matrix with the diagonal elements of the square matrix $\mathbf{Z}$
on its diagonal; $\mathtt{DIAG}{{(}}\mathbf{Z}_{1},\ldots,\mathbf{Z}_{n})$ is
a block-diagonal matrix with blocks $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{n}$.

##  II System Model and Problem Statement

Consider a general MIMO system with $N_{T}$ transmit antennas and $N_{R}$
receive antennas, where $\mathbf{H}\in\mathbb{C}^{N_{R}\times N_{T}}$ denotes
the channel matrix between the transmitter and receiver. For the channel
estimation, $\tau$ pilot symbols per transmit antenna are transmitted, which
are collectively denoted by
$\mathbf{S}^{\mathrm{T}}=[\mathbf{s}_{1}\;\ldots\;\mathbf{s}_{\tau}]\in\mathbb{C}^{N_{T}\times\tau}$.
The received signal at the input of the ADCs is given by

$\displaystyle\mathbf{B}=\mathbf{H}\mathbf{S}^{\mathrm{T}}+\mathbf{N},$ (1)

where $\mathbf{N}\in\mathbb{C}^{N_{R}\times\tau}$ represents noise.
Furthermore, we introduce

$\displaystyle\mathbf{b}=\mathrm{vec}(\mathbf{B})=\mathbf{A}\mathbf{h}+\mathbf{n}\in\mathbb{C}^{\tau
N_{R}},$ (2)

with

$\displaystyle\mathbf{A}=\mathbf{S}\otimes\mathbf{I}_{N_{R}}\in\mathbb{C}^{(\tau
N_{R})\times(N_{T}N_{R})},$ (3)

$\mathbf{h}=\mathrm{vec}(\mathbf{H})\in\mathbb{C}^{N_{T}N_{R}}$, and
$\mathbf{n}=\mathrm{vec}(\mathbf{N})\in\mathbb{C}^{\tau N_{R}}$. We assume a
general correlated Rayleigh fading channel model with
$\mathbf{h}\sim\mathcal{CN}(\mathbf{0},\boldsymbol{\Sigma})$, where
$\boldsymbol{\Sigma}\in\mbox{$\mathbb{C}$}^{(N_{T}N_{R})\times(N_{T}N_{R})}$
is the channel covariance matrix. The probability density function for
$\mathbf{h}$ is expressed as

$\displaystyle
p(\mathbf{h})=\frac{1}{\pi^{N_{T}N_{R}}|\boldsymbol{\Sigma}|}e^{-\mathbf{h}^{\mathrm{H}}\boldsymbol{\Sigma}^{-1}\mathbf{h}}.$
(4)

Moreover, we assume
$\mathbf{n}\sim\mathcal{CN}(\mathbf{0},\sigma^{2}\mathbf{I}_{\tau N_{R}})$,
and the resulting signal-to-noise ratio (SNR) is
$\textrm{SNR}={\mathrm{tr}}(\mathbf{S\mathbf{S}^{\mathrm{H}}})/(\tau
N_{T}\sigma^{2})$, which reduces to $|s|^{2}/\sigma^{2}$ when $\tau=N_{T}=1$.

We denote the element-wise memoryless 1-bit quantization operation as

$\displaystyle\mathcal{Q}_{\textrm{1-bit}}(\cdot)=\mathrm{sgn}\big{(}\mathrm{Re}(\cdot)\big{)}+j\mathrm{sgn}\big{(}\mathrm{Im}(\cdot)\big{)}.$
(5)

The quantized received signal (at the output of the ADCs) is given by

$\displaystyle\mathbf{r}=\mathcal{Q}_{\textrm{1-bit}}(\mathbf{b})\in\mathbb{C}^{\tau
N_{R}}.$ (6)

Based on the quantized received signal $\mathbf{r}$, the MMSE estimate of
$\mathbf{h}$ is given by $\hat{\mathbf{h}}_{\rm
MMSE}=\mathbb{E}(\mathbf{h}|\mathbf{r})=\int_{\mathbb{C}^{N_{T}N_{R}}}\mathbf{h}p(\mathbf{h}|\mathbf{r}){\rm
d}\mathbf{h}$, with
$p(\mathbf{h}|\mathbf{r})=\mathbb{P}(\mathbf{r}|\mathbf{h})p(\mathbf{h})/\mathbb{P}(\mathbf{r})$.
Consequently, we have

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{\displaystyle\int_{\mathbb{C}^{N_{T}N_{R}}}\mathbf{h}\mathbb{P}(\mathbf{r}|\mathbf{h})p(\mathbf{h}){\rm
d}\mathbf{h}}{\mathbb{P}(\mathbf{r})},$ (7)

with

$\displaystyle\mathbb{P}(\mathbf{r})=\int_{\mathbb{C}^{N_{T}N_{R}}}\mathbb{P}(\mathbf{r}|\mathbf{h})p(\mathbf{h}){\rm
d}\mathbf{h}.$ (8)

An explicit expression for the solution to (7) is in general unavailable. In
the next section, we first develop a novel framework for computing (7). Then,
in Section IV, we study some specific cases and determine a general condition
under which the MMSE channel estimate is linear in $\mathbf{r}$, i.e., when
the BLMMSE channel estimator is optimal.

##  III A Framework for 1-Bit MMSE Channel Estimation

In this section, we analyze the structure of the MMSE channel estimator with
1-bit ADCs and reformulate the corresponding conditional mean expression into
a form that involves the orthant probability of a MVN distribution. This will
then be exploited to determine the optimality of the BLMMSE channel estimator.

Let us recall (2) and (5)–(6), and let
$\mathbf{r}=\mathbf{r}_{\mathsf{R}}+j\mathbf{r}_{\mathsf{I}}$, with

$\displaystyle\mathbf{r}_{\mathsf{R}}=\mathrm{Re}(\mathbf{r})=\mathrm{sgn}\left(\mathrm{Re}\left(\mathbf{b}\right)\right),\quad\mathbf{r}_{\mathsf{I}}=\mathrm{Im}(\mathbf{r})=\mathrm{sgn}\left(\mathrm{Im}\left(\mathbf{b}\right)\right).$
(9)

Then, we write the $k$th element of $\mathbf{r}_{\mathsf{R}}$ and
$\mathbf{r}_{\mathsf{I}}$ as $r_{\mathsf{R},k}=\mathrm{Re}(r_{k})$ and
$r_{\mathsf{I},k}=\mathrm{Im}(r_{k})$, respectively, with
$r_{\mathsf{R},k},r_{\mathsf{I},k}\in\\{\pm 1\\}$. Correspondingly, we write
the $k$th element of $\mathbf{r}$ as

$\displaystyle r_{k}$
$\displaystyle=r_{\mathsf{R},k}+jr_{\mathsf{I},k}=\mathrm{sgn}(\mathrm{Re}(b_{k}))+j\mathrm{sgn}(\mathrm{Im}(b_{k}))$
(10)
$\displaystyle=\mathrm{sgn}(\mathrm{Re}(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})+\mathrm{Re}(n_{k}))+j\mathrm{sgn}(\mathrm{Im}(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})+\mathrm{Im}(n_{k})),$
(11)

where $\mathbf{a}_{k}^{\mathrm{T}}$ denotes the $k$th row of $\mathbf{A}$,
$\mathrm{Re}(n_{k})$ and $\mathrm{Im}(n_{k})$ are independent and identically
distributed (i.i.d.) $\mathcal{N}(0,\frac{\sigma^{2}}{2})$ random variables,
and

$\displaystyle(\mathrm{Re}(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h}))^{2}+(\mathrm{Im}(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h}))^{2}=\mathbf{h}^{\mathrm{H}}\mathbf{a}_{k}\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h}.$
(12)

It can be readily shown that

$\displaystyle\mathbb{P}(r_{\mathsf{R},k}|\mathbf{h})$
$\displaystyle=\int\limits_{0}^{\infty}\frac{\exp\left\\{-\frac{1}{\sigma^{2}}\left(x_{k}-r_{\mathsf{R},k}\mathrm{Re}(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})\right)^{2}\right\\}}{\sqrt{\pi\sigma^{2}}}{\rm
d}x_{k},$ (13) $\displaystyle\mathbb{P}(r_{\mathsf{I},k}|\mathbf{h})$
$\displaystyle=\int\limits_{0}^{\infty}\frac{\exp\left\\{-\frac{1}{\sigma^{2}}\left(y_{k}-r_{\mathsf{I},k}\mathrm{Im}(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})\right)^{2}\right\\}}{\sqrt{\pi\sigma^{2}}}{\rm
d}y_{k}.$ (14)

Based on the distribution of the noise, we have
$\mathbb{P}(r_{k}|\mathbf{h})=\mathbb{P}(r_{\mathsf{R},k}|\mathbf{h})\mathbb{P}(r_{\mathsf{I},k}|\mathbf{h})$
and, hence,

$\displaystyle\mathbb{P}(\mathbf{r}|\mathbf{h})=\prod_{k=1}^{\tau
N_{R}}\mathbb{P}(r_{k}|\mathbf{h}).$ (15)

Define $\mathbf{x}=[x_{1}\;\ldots\;x_{\tau N_{R}}]^{\mathrm{T}}$,
$\mathbf{y}=[y_{1}\;\ldots\;y_{\tau N_{R}}]^{\mathrm{T}}$, and

$\displaystyle\boldsymbol{\Lambda}_{\mathsf{R}}=\mathtt{diag}{{(}}\mathbf{r}_{\mathsf{R}}),\quad\boldsymbol{\Lambda}_{\mathsf{I}}=\mathtt{diag}{{(}}\mathbf{r}_{\mathsf{I}}).$
(16)

Clearly, we have
$\boldsymbol{\Lambda}_{\mathsf{R}}\boldsymbol{\Lambda}_{\mathsf{R}}=\boldsymbol{\Lambda}_{\mathsf{I}}\boldsymbol{\Lambda}_{\mathsf{I}}=\mathbf{I}_{\tau
N_{R}}$. Furthermore, it can be shown that (15) is equivalent to

$\displaystyle\mathbb{P}(\mathbf{r}|\mathbf{h})=\int_{\mathbb{R}_{+}^{\tau
N_{R}}}\int_{\mathbb{R}_{+}^{\tau
N_{R}}}\\!\\!\frac{\exp\left\\{-\frac{\|\boldsymbol{\Lambda}_{\mathsf{R}}\mathbf{x}+j\boldsymbol{\Lambda}_{\mathsf{I}}\mathbf{y}-\mathbf{A}\mathbf{h}\|^{2}}{\sigma^{2}}\right\\}}{(\pi\sigma^{2})^{\tau
N_{R}}}{\rm d}\mathbf{x}{\rm d}\mathbf{y}.$ (17)

Now, introduce the covariance matrix of $\mathbf{b}$ in (2), i.e.,

$\displaystyle\boldsymbol{\Omega}_{\rm
b}=\mathbf{A}\boldsymbol{\Sigma}\mathbf{A}^{\mathrm{H}}+\sigma^{2}\mathbf{I}_{\tau
N_{R}}\in\mathbb{C}^{(\tau N_{R})\times(\tau N_{R})},$ (18)

and define

$\displaystyle\mathbf{D}_{\mathsf{R}}$
$\displaystyle=\mathrm{Re}\left(\boldsymbol{\Omega}_{\rm
b}^{-1}\right)=\mathbf{D}_{\mathsf{R}}^{\mathrm{T}}\in\mathbb{R}^{(\tau
N_{R})\times(\tau N_{R})},$ (19) $\displaystyle\mathbf{D}_{\mathsf{I}}$
$\displaystyle=\mathrm{Im}\left(\boldsymbol{\Omega}_{\rm
b}^{-1}\right)=-\mathbf{D}_{\mathsf{I}}^{\mathrm{T}}\in\mathbb{R}^{(\tau
N_{R})\times(\tau N_{R})}.$ (20)

Based on (4) and (17), we rewrite (8) as

$\displaystyle\mathbb{P}(\mathbf{r})$
$\displaystyle=\displaystyle\int_{\mathbb{C}^{N_{T}N_{R}}}\mathbb{P}(\mathbf{r}|\mathbf{h})\frac{e^{-\mathbf{h}^{\mathrm{H}}\boldsymbol{\Sigma}^{-1}\mathbf{h}}}{\pi^{N_{T}N_{R}}|\boldsymbol{\Sigma}|}{\rm
d}\mathbf{h}$ (21) $\displaystyle=\frac{1}{\pi^{\tau
N_{R}}|\boldsymbol{\Omega}_{\rm b}|}\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm d}\mathbf{z},$
(22)

with

$\displaystyle\mathbf{C}$
$\displaystyle=\begin{bmatrix}\boldsymbol{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{R}}\boldsymbol{\Lambda}_{\mathsf{R}}&\boldsymbol{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{I}}^{\mathrm{T}}\boldsymbol{\Lambda}_{\mathsf{I}}\\\
\boldsymbol{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{I}}\boldsymbol{\Lambda}_{\mathsf{R}}&\boldsymbol{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{R}}\boldsymbol{\Lambda}_{\mathsf{I}}\end{bmatrix}\in\mathbb{R}^{(2\tau
N_{R})\times(2\tau N_{R})},$ (23) $\displaystyle\mathbf{z}$
$\displaystyle=\left[\mathbf{x}^{\mathrm{T}}\;\;\mathbf{y}^{\mathrm{T}}\right]^{\mathrm{T}}\in\mathbb{R}^{2\tau
N_{R}}.$ (24)

Following similar steps, we rewrite the numerator of (7) as

$\displaystyle\displaystyle\int_{\mathbb{C}^{N_{T}N_{R}}}\mathbf{h}\mathbb{P}(\mathbf{r}|\mathbf{h})p(\mathbf{h}){\rm
d}\mathbf{h}$
$\displaystyle=\frac{\boldsymbol{\Sigma}\mathbf{A}^{\mathrm{H}}\boldsymbol{\Omega}_{\rm
b}^{-1}\left[\boldsymbol{\Lambda}_{\mathsf{R}}\;\;j\boldsymbol{\Lambda}_{\mathsf{I}}\right]}{\pi^{\tau
N_{R}}|\boldsymbol{\Omega}_{\rm b}|}$ $\displaystyle\phantom{=}\
\cdot\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\mathbf{z}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm
d}\mathbf{z}.$ (25)

Finally, based on (22) and (25), we express (7) as

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\boldsymbol{\Sigma}\mathbf{A}^{\mathrm{H}}\boldsymbol{\Omega}_{\rm
b}^{-1}\left[\boldsymbol{\Lambda}_{\mathsf{R}}\;\;j\boldsymbol{\Lambda}_{\mathsf{I}}\right]\left(\\!\frac{\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\mathbf{z}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm
d}\mathbf{z}}{\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}e^{-{\mathbf{z}}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm
d}{\mathbf{z}}}\\!\right).$ (26)

The fraction in (26) denotes the mean of a truncated MVN distribution with
$\mathbb{R}_{+}^{2\tau N_{R}}$ as its support. Moreover, since $\frac{{\rm
d}}{{\rm
d}\mathbf{z}}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}=-2\mathbf{Cz}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}$,
we have

$\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\mathbf{z}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm
d}\mathbf{z}=\frac{1}{2}\mathbf{C}^{-1}\boldsymbol{\alpha},$ (27)

where the $k$th element of $\boldsymbol{\alpha}\in\mathbb{R}^{2\tau N_{R}}$ is
given by

$\displaystyle\alpha_{k}=\int_{\mathbb{R}_{+}^{2\tau
N_{R}-1}}e^{-\mathbf{z}_{-k}^{\mathrm{T}}\mathbf{C}_{-k}\mathbf{z}_{-k}}{\rm
d}\mathbf{z}_{-k},$ (28)

with $\mathbf{z}_{-k}=[z_{1}\;\ldots\;z_{k-1}\;z_{k+1}\;\ldots\;z_{2\tau
N_{R}}]^{\mathrm{T}}\in\mathbb{R}^{2\tau N_{R}-1}$ and where
$\mathbf{C}_{-k}\in\mathbb{R}^{(2\tau N_{R}-1)\times(2\tau N_{R}-1)}$ is
obtained from $\mathbf{C}$ by deleting its $k$th row and column. Hence, (26)
is equivalent to

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{\boldsymbol{\Sigma}\mathbf{A}^{\mathrm{H}}\boldsymbol{\Omega}_{\rm
b}^{-1}\left[\boldsymbol{\Lambda}_{\mathsf{R}}\;\;j\boldsymbol{\Lambda}_{\mathsf{I}}\right]\mathbf{C}^{-1}\boldsymbol{\alpha}}{\displaystyle
2\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm d}\mathbf{z}}.$
(29)

######  Remark 1.

Each off-diagonal element of $\mathbf{C}$ is affected by exactly two different
elements from
$[\mathbf{r}_{\mathsf{R}}^{\mathrm{T}}\;\;\mathbf{r}_{\mathsf{I}}^{\mathrm{T}}]^{\mathrm{T}}$.
The same comment also applies to $\mathbf{C}_{-k}$. This fact is instrumental
for deriving computationally efficient expressions of $\hat{\mathbf{h}}_{\rm
MMSE}$ under specific assumptions on the channel correlation or pilot symbols
[8].

The expression in (29) is directly related to the orthant probability of a
real-valued MVN distribution, defined as [9]

$\displaystyle\mathcal{P}(\boldsymbol{\Psi})=\frac{1}{(2\pi)^{\frac{L}{2}}\left|\boldsymbol{\Psi}\right|^{\frac{1}{2}}}\displaystyle\int_{\mathbb{R}_{+}^{L}}e^{-\frac{1}{2}\mathbf{u}^{\mathrm{T}}\boldsymbol{\Psi}^{-1}\mathbf{u}}{\rm
d}\mathbf{u},$ (30)

where $\boldsymbol{\Psi}\in\mbox{$\mathbb{R}$}^{L\times L}$ denotes the
covariance matrix. Based on (30), we rewrite (28) as

$\displaystyle\alpha_{k}=\left(2\pi\right)^{\frac{2\tau
N_{R}-1}{2}}\left|\frac{1}{2}\mathbf{C}_{-k}^{-1}\right|^{\frac{1}{2}}\mathcal{P}\left(\frac{1}{2}\mathbf{C}_{-k}^{-1}\right)$
(31)

and the integral in the denominator of (29) as

$\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm
d}\mathbf{z}=\left(2\pi\right)^{\tau
N_{R}}\left|\frac{1}{2}\mathbf{C}^{-1}\right|^{\frac{1}{2}}\mathcal{P}\left(\frac{1}{2}\mathbf{C}^{-1}\right).$
(32)

It can also be shown that

$\displaystyle|\mathbf{C}_{-k}^{-1}|/|\mathbf{C}^{-1}|=|\mathbf{C}|/|\mathbf{C}_{-k}|=1/[\mathbf{C}^{-1}]_{kk}.$
(33)

Based on (31)–(33), we rewrite (29) as

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{\boldsymbol{\Sigma}\mathbf{A}^{\mathrm{H}}\boldsymbol{\Omega}_{\rm
b}^{-1}\left[\boldsymbol{\Lambda}_{\mathsf{R}}\;\;j\boldsymbol{\Lambda}_{\mathsf{I}}\right]\mathbf{C}^{-1}\left(\widetilde{\mathtt{diag}{{}}}(\mathbf{C}^{-1})\right)^{-\frac{1}{2}}\mathbf{g}}{2\sqrt{\pi}\mathcal{P}\left(\frac{1}{2}\mathbf{C}^{-1}\right)},$
(34)

where the $k$th element of $\mathbf{g}\in\mbox{$\mathbb{R}$}^{2\tau N_{R}}$ is
given by $g_{k}=\mathcal{P}\left(\frac{1}{2}\mathbf{C}_{-k}^{-1}\right)$.

######  Remark 2.

The channel covariance matrix $\boldsymbol{\Sigma}$ and the pilot matrix
$\mathbf{S}$, which are embedded in $\boldsymbol{\Omega}_{\rm b}$ in (18) and
in $\mathbf{C}$ in (23) through (19)–(20), determine the complexity of
computing $\hat{\mathbf{h}}_{\rm MMSE}$ in (34). The number of transmit
antennas $N_{T}$ is not involved in the computation of the orthant
probability.

In [10], a method based on the characteristic function (CF), i.e., the Fourier
transform with sign reversal, was used to express the orthant probability in
(30). Based on [9, 10], we have

$\displaystyle\mathcal{P}(\boldsymbol{\Psi})=\begin{cases}\frac{1}{4}+\frac{\arcsin(\psi_{12})}{2\pi}&\textrm{for
}L=2,\\\
\frac{1}{8}+\frac{\sum_{i=1}^{2}\sum_{k=i+1}^{3}\arcsin(\psi_{ik})}{4\pi}&\textrm{for
}L=3,\end{cases}$ (35)

where $\boldsymbol{\Psi}$ is standardized. In the Appendix, we extend the CF-
based method to compute the numerator of the fraction in (26), thereby
providing an alternative to (27)–(28).

##  IV Optimality of the BLMMSE Channel Estimator

In this section, building on the above framework and utilizing the
observations in Remark 2, we present three specific cases where the MMSE
channel estimate is linear in $\mathbf{r}$, i.e., the BLMMSE channel estimator
is optimal (Sections IV-A, IV-B, and IV-C), along with one specific case where
it is not (Section IV-D). Then, we determine a necessary and sufficient
condition for the optimality of the BLMMSE channel estimator (Section IV-E).
To facilitate the subsequent exposition, we first express the BLMMSE channel
estimate as [3]

$\displaystyle\hat{\mathbf{h}}_{\rm BLMMSE}$
$\displaystyle=\frac{\sqrt{\pi}}{2}\boldsymbol{\Sigma}\mathbf{A}^{\mathrm{H}}\mathbf{D}_{\Omega}^{-\frac{1}{2}}\Big{(}\arcsin\big{(}\mathbf{D}_{\Omega}^{-\frac{1}{2}}\mathrm{Re}\left(\boldsymbol{\Omega}_{\rm
b}\right)\mathbf{D}_{\Omega}^{-\frac{1}{2}}\big{)}$ $\displaystyle\phantom{=}\
+j\arcsin\big{(}\mathbf{D}_{\Omega}^{-\frac{1}{2}}\mathrm{Im}\left(\boldsymbol{\Omega}_{\rm
b}\right)\mathbf{D}_{\Omega}^{-\frac{1}{2}}\big{)}\Big{)}^{-1}\mathbf{r},$
(36)

with
$\mathbf{D}_{\Omega}=\widetilde{\mathtt{diag}{{}}}(\boldsymbol{\Omega}_{\rm
b})$ and $\boldsymbol{\Omega}_{\rm b}$ defined in (18).

###  IV-A Linear Case: Uncorrelated Channel and Unitary Pilot Matrix

Let $\tau=N_{T}$, and assume $\boldsymbol{\Sigma}=\mathbf{I}_{N_{T}N_{R}}$ and
$\mathbf{S}\mathbf{S}^{\mathrm{H}}=\eta\mathbf{I}_{\tau}$, with $\eta>0$.
Then,

$\displaystyle\mathbf{D}_{\mathsf{R}}=\frac{1}{\eta+\sigma^{2}}\mathbf{I}_{\tau
N_{R}},\quad\mathbf{D}_{\mathsf{I}}=\mathbf{0},\quad\mathbf{C}=\frac{1}{\eta+\sigma^{2}}\mathbf{I}_{2\tau
N_{R}}.$ (37)

Further algebraic manipulations based on (34) and (36) lead to

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{\left(\mathbf{S}^{\mathrm{H}}\otimes\mathbf{I}_{N_{R}}\right)\mathbf{r}}{\sqrt{\pi(\eta+\sigma^{2})}}=\hat{\mathbf{h}}_{\rm
BLMMSE}.$ (38)

This specific case is illustrated in Fig. 1 for different numbers of transmit
antennas. The result in (38) encompasses the uncorrelated SIMO system
considered in [6] as a special case and can be further generalized using the
proposed framework, as shown next.

###  IV-B Linear Case: Transmit-Only Channel Correlation

$-10$$0$$10$$20$$30$$0.35$$0.4$$0.45$$0.5$$0.55$$0.6$SNR [dB]MSE per
antenna$N_{T}=8$ (uncorrelated)$N_{T}=16$ (uncorrelated)$N_{T}=32$
(uncorrelated)$N_{T}=8$ (correlated)$N_{T}=16$ (correlated)$N_{T}=32$
(correlated) Figure 1: MSE per antenna versus SNR for the specific cases of
uncorrelated channel and unitary pilot matrix described in Section IV-A
(uncorrelated) and transmit-only channel correlation described in Section IV-B
(correlated). The lines and markers denote the analytical results and
simulations, respectively.

Let $\tau=N_{T}$ and assume $\mathbf{H}=\mathbf{H}_{\rm
w}\boldsymbol{\Sigma}_{\rm TX}^{\frac{1}{2}}$, where $\boldsymbol{\Sigma}_{\rm
TX}\in\mbox{$\mathbb{C}$}^{N_{T}\times N_{T}}$ represents the transmit channel
correlation matrix and where the elements of $\mathbf{H}_{\rm
w}\in\mbox{$\mathbb{C}$}^{N_{R}\times N_{T}}$ are i.i.d. $\mathcal{CN}(0,1)$
random variables. Then, we have $\boldsymbol{\Sigma}=\boldsymbol{\Sigma}_{\rm
TX}\otimes\mathbf{I}_{N_{R}}$. Now, let $\boldsymbol{\Sigma}_{\rm
TX}=\mathbf{U}\boldsymbol{\Xi}\mathbf{U}^{\mathrm{H}}$ denote the eigenvalue
decomposition of $\boldsymbol{\Sigma}_{\rm TX}$, with
$\boldsymbol{\Xi}=\mathtt{diag}{{}}(\left[\xi_{i}\right]_{i=1}^{N_{T}})$ and
$\xi_{i}>0$. Utilizing the insights from Remark 2, we set
$\mathbf{S}=\sqrt{\eta}\mathbf{U}^{\mathrm{H}}$, which yields

$\displaystyle\mathbf{D}_{\mathsf{R}}=\mathtt{diag}{{\left(\left[\frac{1}{\eta\xi_{i}+\sigma^{2}}\right]_{i=1}^{N_{T}}\right)}}\otimes\mathbf{I}_{N_{R}},\quad\mathbf{D}_{\mathsf{I}}=\mathbf{0}.$
(39)

Furthermore,

$\displaystyle\mathbf{C}=\mathtt{DIAG}{{(}}\mathbf{D}_{\mathsf{R}},\mathbf{D}_{\mathsf{R}})$
(40)

is also diagonal and the integrals in (26) are over scalar variables, which
allows us to write

$\displaystyle\hat{\mathbf{h}}_{\rm MMSE}$
$\displaystyle=\left(\left(\mathbf{U}\;\mathtt{diag}{{}}\left(\left[\frac{\xi_{i}\sqrt{\eta}}{\sqrt{\eta\xi_{i}+\sigma^{2}}}\right]_{i=1}^{N_{T}}\right)\right)\otimes\mathbf{I}_{N_{R}}\right)\frac{\mathbf{r}}{\sqrt{\pi}}$
(41) $\displaystyle=\hat{\mathbf{h}}_{\rm BLMMSE}.$ (42)

This specific case is included in Fig. 1 for different numbers of transmit
antennas, where the transmit channel correlation coefficients are given by
$[\boldsymbol{\Sigma}_{\rm TX}]_{ik}=\textsf{J}_{0}(2\pi(k-i)\Delta\gamma_{\rm
max}\cos(\theta))e^{-j2\pi(k-i)\Delta\sin(\theta)}$, where
$\textsf{J}_{0}(\cdot)$ denotes the zeroth-order Bessel function of the first
kind and where we have set $\gamma_{\rm max}=0.1$ (maximum angle spread),
$\theta=\frac{\pi}{6}$ (angle of arrival), and $\Delta=0.5$ wavelengths
(antenna spacing) [11]. Clearly, (38) represents a special case of (42) with
$\boldsymbol{\Sigma}_{\rm TX}=\mathbf{I}_{N_{T}}$. If
$\boldsymbol{\Sigma}_{\rm TX}$ is circulant, then a discrete Fourier transform
matrix can be used as the pilot matrix.

###  IV-C Linear Case: SIMO System with Real-Valued Channel Correlation
($N_{R}=2$)

Assume a correlated SIMO system with a single complex-valued pilot symbol $s$,
i.e., $\tau=N_{T}=1$, so that the received signal in (2) becomes

$\displaystyle\mathbf{b}=\mathbf{h}s+\mathbf{n}\in\mbox{$\mathbb{C}$}^{N_{R}}.$
(43)

Let $\boldsymbol{\Sigma}$ be real-valued, as for example in the exponentially
correlated channel model [12]. Then, we have

$\displaystyle\mathbf{D}_{\mathsf{R}}=\boldsymbol{\Omega}_{\rm
b}^{-1}=(|s|^{2}\boldsymbol{\Sigma}+\sigma^{2}\mathbf{I}_{N_{R}})^{-1},\quad\mathbf{D}_{\mathsf{I}}=\mathbf{0},$
(44)

and

$\displaystyle\mathbf{C}=\mathtt{DIAG}{{(}}\boldsymbol{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{R}}\boldsymbol{\Lambda}_{\mathsf{R}},\boldsymbol{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{R}}\boldsymbol{\Lambda}_{\mathsf{I}}).$
(45)

When $N_{R}=2$, the diagonal blocks in (45) are $2\times 2$ matrices. Finally,
assuming a standardized $\boldsymbol{\Sigma}$ with channel correlation
coefficient $[\boldsymbol{\Sigma}]_{12}=\rho_{12}$ and defining
$\beta_{12}=\frac{\rho_{12}|s|^{2}}{|s|^{2}+\sigma^{2}}$, we obtain

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{s^{*}\boldsymbol{\Sigma}\mathbf{T}^{-1}\mathbf{r}}{\sqrt{\pi(|s|^{2}+\sigma^{2})}}=\hat{\mathbf{h}}_{\rm
BLMMSE}.$ (46)

with $[\mathbf{T}]_{12}=[\mathbf{T}]_{21}=\frac{2}{\pi}\arcsin(\beta_{12})$
and $[\mathbf{T}]_{11}=[\mathbf{T}]_{22}=1$.

This specific case is included in Fig. 2, where the channel correlation
coefficient between antennas $i$ and $k$ is given by
$[\boldsymbol{\Sigma}]_{ik}=\rho^{|i-k|}$ [12].

###  IV-D Non-Linear Case: SIMO System with Real-Valued Channel Correlation
($N_{R}=3$)

$-10$$0$$10$$20$$30$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$\rho=0.35$$\rho=0.65$$\rho=0.95$SNR
[dB]MSE per antenna$N_{R}=2$ (MMSE = BLMMSE)$N_{R}=3$ (MMSE)$N_{R}=3$ (BLMMSE)
Figure 2: MSE per antenna versus SNR for the specific cases of SIMO with real-
valued channel correlation and $N_{R}=2$ described in Section IV-C and
$N_{R}=3$ described in Section IV-D.

Here, we consider the same SIMO system described in Section IV-C but with
$N_{R}=3$. In this case, (43)–(45) remain valid. Without loss of generality,
assume a standardized $\boldsymbol{\Sigma}$ with channel correlation
coefficient between antennas $i$ and $k$ denoted by
$[\boldsymbol{\Sigma}]_{ik}=\rho_{ik}$. Further define
$\beta_{ik}=\frac{\rho_{ik}|s|^{2}}{|s|^{2}+\sigma^{2}}$. Based on (26) or
(34) and using (35), we obtain the following result. The proof is omitted due
to space limitations and can be obtained from [8].

######  Theorem 1.

Consider the SIMO system with $\tau=N_{T}=1$ described by (43). Assume real-
valued channel correlation and $N_{R}=3$. Then, we have

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{s^{*}\boldsymbol{\Sigma}}{2\sqrt{\pi\left(|s|^{2}+\sigma^{2}\right)}}\left(\frac{\mathbf{v}_{\rm
3R}}{P_{\rm 3R}}+j\frac{\mathbf{v}_{\rm 3I}}{P_{\rm 3I}}\right),$ (47)

with

$\displaystyle\mathbf{v}_{\rm 3R}$
$\displaystyle=\begin{bmatrix}\frac{r_{\mathsf{R},1}}{4}+\frac{r_{\mathsf{R},1}r_{\mathsf{R},2}r_{\mathsf{R},3}}{2\pi}\arcsin\left(\frac{\beta_{23}-\beta_{12}\beta_{13}}{\sqrt{1-\beta_{12}^{2}}\sqrt{1-\beta_{13}^{2}}}\right)\\\
\frac{r_{\mathsf{R},2}}{4}+\frac{r_{\mathsf{R},1}r_{\mathsf{R},2}r_{\mathsf{R},3}}{2\pi}\arcsin\left(\frac{\beta_{13}-\beta_{12}\beta_{23}}{\sqrt{1-\beta_{12}^{2}}\sqrt{1-\beta_{23}^{2}}}\right)\\\
\frac{r_{\mathsf{R},3}}{4}+\frac{r_{\mathsf{R},1}r_{\mathsf{R},2}r_{\mathsf{R},3}}{2\pi}\arcsin\left(\frac{\beta_{12}-\beta_{13}\beta_{23}}{\sqrt{1-\beta_{13}^{2}}\sqrt{1-\beta_{23}^{2}}}\right)\end{bmatrix},$
(48) $\displaystyle P_{\rm 3R}$
$\displaystyle=\frac{1}{8}+\frac{\sum_{i=1}^{2}\sum_{k=i+1}^{3}r_{\mathsf{R},i}r_{\mathsf{R},k}\arcsin(\beta_{ik})}{4\pi},$
(49)

and where $\mathbf{v}_{\rm 3I}$ (resp. $P_{\rm 3I}$) is the same as
$\mathbf{v}_{\rm 3R}$ (resp. $P_{\rm 3R}$) except that $r_{\mathsf{R},i}$ in
(48) (resp. (49)) is replaced by $r_{\mathsf{I},i}$.

Under the assumptions of Theorem 1, from (36), we have

$\displaystyle\hat{\mathbf{h}}_{\rm
BLMMSE}=\frac{s^{*}\boldsymbol{\Sigma}\widetilde{\mathbf{T}}^{-1}\mathbf{r}}{\sqrt{\pi\left(|s|^{2}+\sigma^{2}\right)}},$
(50)

with $[\widetilde{\mathbf{T}}]_{ik}=\frac{2}{\pi}\arcsin(\beta_{ik})$ for
$i\neq k$ and $[\widetilde{\mathbf{T}}]_{ii}=1$.

######  Remark 3.

Note that (47) is non-linear in $\mathbf{r}$ and thus not equivalent to (50),
unless one of the following holds: $\rho_{23}=\rho_{13}=0$,
$\rho_{13}=\rho_{12}=0$, or $\rho_{12}=\rho_{23}=0$.

######  Remark 4.

Theorem 1 provides the simplest closed-form expression for the MMSE channel
estimate when it is non-linear in $\mathbf{r}$.

This specific case is included in Fig. 2, which shows that the BLMMSE channel
estimator is nearly optimal in cases with weak channel correlation.

###  IV-E Necessary and Sufficient Optimality Condition

Finally, we provide a necessary and sufficient condition for the optimality of
the BLMMSE channel estimator. The proof builds on the intuition obtained from
Sections III, IV-C, and IV-D; the details are omitted due to space limitations
and can be found in [8].

######  Theorem 2.

The BLMMSE channel estimator is optimal if and only if the off-diagonal
elements of $\mathbf{C}$ in (23) satisfy $[\mathbf{C}]_{il}\neq 0$ for at most
one value of $l$ other than $l=i$, i.e., an element $z_{i}$ of the random
vector $\mathbf{z}$ defined in (24) is correlated with at most one other
element $z_{l}$, for $l\neq i$.

It is straightforward to observe that $\mathbf{C}$ in (37), (40), and (45)
with $N_{R}=2$ satisfies the condition of Theorem 2 and, therefore, the BLMMSE
channel estimator is optimal. On the other hand, (45) with $N_{R}=3$ does not
satisfy the condition of Theorem 2 except for the special cases mentioned in
Remark 3.

Motivated by [10], we provide an alternative to (27)–(28) for computing the
numerator of the fraction in (26). To this end, let $p(\mathbf{u})$ correspond
to $\mathcal{N}(\mathbf{0},\boldsymbol{\Psi})$ and write

$\displaystyle\displaystyle\int_{\mathbb{R}_{+}^{L}}\mathbf{u}p(\mathbf{u}){\rm
d}\mathbf{u}=\displaystyle\int_{\mathbb{R}_{+}^{L}}\frac{\mathbf{u}e^{-\frac{1}{2}\mathbf{u}^{\mathrm{T}}\boldsymbol{\Psi}^{-1}\mathbf{u}}}{(2\pi)^{\frac{L}{2}}\left|\boldsymbol{\Psi}\right|^{\frac{1}{2}}}{\rm
d}\mathbf{u},$ (51)

whose $i$th element is given by

$\displaystyle\mathcal{I}_{i}=\displaystyle\int_{\mathbb{R}_{+}^{L}}u_{i}p(\mathbf{u}){\rm
d}\mathbf{u}=\displaystyle\int_{\mathbb{R}^{L}}\underbrace{u_{i}\mathcal{U}(u_{i})\prod_{\begin{subarray}{c}k=1\\\
k\neq
i\end{subarray}}^{L}\mathcal{U}(u_{k})}_{\varphi_{i}(u_{1},\ldots,u_{L})=\varphi_{i}(\mathbf{u})}p(\mathbf{u}){\rm
d}\mathbf{u},$ (52)

where $\mathcal{U}(\cdot)$ represents the Heaviside step function. Let
$\boldsymbol{\omega}=[\omega_{1}\ldots\omega_{L}]^{\mathrm{T}}$ and express
the Fourier transform of $\varphi_{i}(\mathbf{u})$ as

$\displaystyle\mathcal{F}_{\varphi_{i}}(\boldsymbol{\omega})=\left(j\pi\delta^{\prime}(\omega_{i})-\frac{1}{\omega_{i}^{2}}\right)\prod_{\begin{subarray}{c}k=1\\\
k\neq
i\end{subarray}}^{L}\left(\pi\delta(\omega_{k})+\frac{1}{j\omega_{k}}\right),$
(53)

where $\delta(\cdot)$ is the Dirac delta function and $\delta^{\prime}(\cdot)$
its derivative. Meanwhile, the Fourier transform of $p(\mathbf{u})$ is given
by
$e^{-\frac{1}{2}\boldsymbol{\omega}^{\mathrm{T}}\boldsymbol{\Psi}\boldsymbol{\omega}}$.
As in [10], we use the Parseval–Plancherel identity to rewrite (52) as

$\displaystyle\mathcal{I}_{i}=\frac{1}{(2\pi)^{L}}\displaystyle\int_{\mathbb{R}^{L}}\mathcal{F}_{\varphi_{i}}(\boldsymbol{\omega})e^{-\frac{1}{2}\boldsymbol{\omega}^{\mathrm{T}}\boldsymbol{\Psi}\boldsymbol{\omega}}{\rm
d}\boldsymbol{\omega}.$ (54)

Combining (51), (53)–(54), and [10, Eq. (6)] provides a method based solely on
the CF to compute the fraction in (26).

As an example, consider the case of $L=2$ and assume a standardized
$\boldsymbol{\Psi}$ with off-diagonal element $\psi_{12}$. Using (54), it can
be shown that $\mathcal{I}_{1}=\mathcal{I}_{2}=(1+\psi_{12})/(2\sqrt{2\pi})$.
The case of $L=3$ can be readily developed following similar steps.

## References

  * [1] N. Rajatheva, I. Atzeni, E. Björnson _et al._ , “White paper on broadband connectivity in 6G,” 2020. [Online]. Available: https://oulurepo.oulu.fi/handle/10024/36799
  * [2] S. Jacobsson, G. Durisi, M. Coldrey, U. Gustavsson, and C. Studer, “Throughput analysis of massive MIMO uplink with low-resolution ADCs,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 6, pp. 4038–4051, 2017.
  * [3] Y. Li, C. Tao, G. Seco-Granados, A. Mezghani, A. L. Swindlehurst, and L. Liu, “Channel estimation and performance analysis of one-bit massive MIMO systems,” _IEEE Trans. Signal Process._ , vol. 65, no. 15, pp. 4075–4089, 2017.
  * [4] I. Atzeni and A. Tölli, “Channel estimation and data detection analysis of massive MIMO with 1-bit ADCs,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 6, pp. 3850–3867, 2022.
  * [5] I. Atzeni, A. Tölli, and G. Durisi, “Low-resolution massive MIMO under hardware power consumption constraints,” in _Proc. Asilomar Conf. Signals, Syst., and Comput. (ASILOMAR)_ , 2021.
  * [6] B. Fesl, M. Koller, and W. Utschick, “On the mean square error optimal estimator in one-bit quantized systems,” _IEEE Trans. Signal Process._ , vol. 71, pp. 1968–1980, 2023.
  * [7] J. J. Bussgang, “Crosscorrelation functions of amplitude-distorted Gaussian signals,” _Res. Lab. Electron., MIT, Cambridge, MA, U.S.A._ , 1952.
  * [8] M. Ding, I. Atzeni, A. Tölli, and A. L. Swindlehurst, “On the optimal MMSE channel estimation for one-bit quantized MIMO systems,” 2024. [Online]. Available: https://arxiv.org/pdf/2404.05536.pdf
  * [9] I. G. Abrahamson, “Orthant probabilities of the quadrivariate normal distribution,” _Ann. Math. Statist._ , vol. 35, no. 4, pp. 1685–1703, 1964\.
  * [10] D. R. Childs, “Reduction of the multivariate normal integral to characteristic form,” _Biometrika_ , vol. 54, no. 1–2, pp. 293–300, 1967.
  * [11] A. Abdi and M. Kaveh, “A space-time correlation model for multielement antenna systems in mobile fading channels,” _IEEE J. Sel. Areas Commun._ , vol. 20, no. 3, pp. 550–560, 2002.
  * [12] M. Chiani, M. Z. Win, and A. Zanella, “On the capacity of spatially correlated MIMO Rayleigh-fading channels,” _IEEE Trans. Inf. Theory_ , vol. 49, no. 10, pp. 2363–2371, 2003.


11institutetext: Departamento de Física, Universidad Técnica Federico Santa
María, Casilla 110-V, Valparaíso, Chile 22institutetext: Centro Científico
Tecnológico de Valparaíso, Casilla 110-V, Valparaíso, Chile

# Polynomial affine gravity in $3+1$ dimensions

José Perdiguero 11<EMAIL_ADDRESS>Oscar Castillo-Felisola 1122
<EMAIL_ADDRESS>
###### Abstract

The polynomial affine gravity is an alternative model of gravity whose
fundamental field is the affine connection, and it is invariant under the
complete group of diffeomorphisms. In $3+1$ dimensions the field equations
generalise those of Einstein–Hilbert, the coupling constants are
dimensionless, the action has a finite numbers of term, and although the
action does not involve a (fundamental) metric, some metric tensor fields
might _emerge_ from the connection. Provided a cosmological ansatz, the
properties of diverse cosmological models are discussed.

##  1 Introduction

The relativistic equations for gravitational interactions were introduced in
1915 by Einstein and Hilbert independently hilbert15_grund_physik ;
einstein15_grund_allgem_relat_und_anwen , cementing the idea of geometrisation
of gravity. The foundation of General Relativity is Riemannian geometry, which
is determined by the metric tensor field. Such choice endows the metric with a
double role: (i) it is the object that serves to _measure distances_ , and
(ii) it is the mediator of gravitational interactions.

However, the development of the _Erlangen Program_ (initiated by Klein
klein08_compar_review_recen_resear_geomet  and concluded by Cartan
ivey03_cartan_begin ) yields the idea of non-Riemannian spaces
eisenhart27_non_rieman , sometimes called _metric-affine_ spaces
hehl95_metric_affin_gauge_theor_gravit , inducing several proposal
generalising the model of Einstein–Hilbert. For example, Weyl considered an
affine connection which was not compatible with the metric weyl18_reine_infin
; weyl18_gravit_und_elekt , Cartan considered a non-symmetric connection
cartan23_sur_les_connex_affin_et ; cartan24_sur_les_connex_affin_et ;
cartan25_sur_les_connex_affin_et , Veblen proposed a projective version of
General Relativity veblen33_projec_gener_relat , and Einstein and Eddington
proposed an affine formulation of General Relativity
einstein23_zur_affin_feldt ; eddington23 , where the metric does not play the
role of mediating the interactions.

From those scenarios, the affine formulations of gravity are the more general,
since in principle the space do not require a (pre-existing) metric structure.
The polynomial affine model of gravity castillo-
felisola15_polyn_model_purel_affin_gravit ; castillo-
felisola18_einst_gravit_from_polyn_affin_model  was built up as an attempt to
obtain the less constrained model for gravitational interactions.

In this letter we briefly discuss the properties of the cosmological models
found on the polynomial affine model of gravity complementing the results from
Ref. castillo-felisola21_aspec_polyn_affin_model_gravit_three ; castillo-
felisola22_polyn_affin_model_gravit .

##  2 The model

As we mentioned, the polynomial affine model of gravity is a field theory for
the affine connection, $\hat{\Gamma}_{{\mu}}{}^{{\lambda}}{}_{{\nu}}$, and the
most general action can be built up using their irreducible components,
$\Gamma_{{\mu}}{}^{{\lambda}}{}_{{\nu}}$, ${\mathcal{A}}_{\mu}$ and
$\mathcal{B}_{{\mu}}{}^{{\lambda}}{}_{{\nu}}$, defined by the relation
$\hat{\Gamma}_{{\mu}}{}^{{\lambda}}{}_{{\nu}}=\Gamma_{{\mu}}{}^{{\lambda}}{}_{{\nu}}+\mathcal{B}_{{\mu}}{}^{{\lambda}}{}_{{\nu}}+\delta^{\lambda}_{[\mu}{\mathcal{A}}_{\nu]}$,
and the volume element
$\mathop{\mathrm{d}V}\nolimits^{\alpha\beta\gamma\delta}=\mathrm{d}x^{\alpha}\mathop{\scriptstyle\wedge}\mathrm{d}x^{\beta}\mathop{\scriptstyle\wedge}\mathrm{d}x^{\gamma}\mathop{\scriptstyle\wedge}\mathrm{d}x^{\delta}$
(See Ref. castillo-felisola20_emerg_metric_geodes_analy_cosmol  for details).

Up to boundary and topological terms, the action is written as castillo-
felisola20_emerg_metric_geodes_analy_cosmol ,

$S=\int\mathop{\mathrm{d}V}\nolimits^{\alpha\beta\gamma\delta}\biggl{[}B_{1}\mathcal{R}_{{\mu\nu}}{}^{{\mu}}{}_{{\rho}}\mathcal{B}_{{\alpha}}{}^{{\nu}}{}_{{\beta}}\mathcal{B}_{{\gamma}}{}^{{\rho}}{}_{{\delta}}+B_{2}\mathcal{R}_{{\alpha\beta}}{}^{{\mu}}{}_{{\rho}}\mathcal{B}_{{\gamma}}{}^{{\nu}}{}_{{\delta}}\mathcal{B}_{{\mu}}{}^{{\rho}}{}_{{\nu}}+B_{3}\mathcal{R}_{{\mu\nu}}{}^{{\mu}}{}_{{\alpha}}\mathcal{B}_{{\beta}}{}^{{\nu}}{}_{{\gamma}}{\mathcal{A}}_{\delta}+B_{4}\mathcal{R}_{{\alpha\beta}}{}^{{\sigma}}{}_{{\rho}}\mathcal{B}_{{\gamma}}{}^{{\rho}}{}_{{\delta}}{\mathcal{A}}_{\sigma}+B_{5}\mathcal{R}_{{\alpha\beta}}{}^{{\rho}}{}_{{\rho}}\mathcal{B}_{{\gamma}}{}^{{\sigma}}{}_{{\delta}}{\mathcal{A}}_{\sigma}+C_{1}\mathcal{R}_{{\mu\alpha}}{}^{{\mu}}{}_{{\nu}}\nabla_{\beta}\mathcal{B}_{{\gamma}}{}^{{\nu}}{}_{{\delta}}+C_{2}\mathcal{R}_{{\alpha\beta}}{}^{{\rho}}{}_{{\rho}}\nabla_{\sigma}\mathcal{B}_{{\gamma}}{}^{{\sigma}}{}_{{\delta}}+D_{1}\mathcal{B}_{{\nu}}{}^{{\mu}}{}_{{\lambda}}\mathcal{B}_{{\mu}}{}^{{\nu}}{}_{{\alpha}}\nabla_{\beta}\mathcal{B}_{{\gamma}}{}^{{\lambda}}{}_{{\delta}}+D_{2}\mathcal{B}_{{\alpha}}{}^{{\mu}}{}_{{\beta}}\mathcal{B}_{{\mu}}{}^{{\lambda}}{}_{{\nu}}\nabla_{\lambda}\mathcal{B}_{{\gamma}}{}^{{\nu}}{}_{{\delta}}+D_{3}\mathcal{B}_{{\alpha}}{}^{{\mu}}{}_{{\nu}}\mathcal{B}_{{\beta}}{}^{{\lambda}}{}_{{\gamma}}\nabla_{\lambda}\mathcal{B}_{{\mu}}{}^{{\nu}}{}_{{\delta}}+D_{4}\mathcal{B}_{{\alpha}}{}^{{\lambda}}{}_{{\beta}}\mathcal{B}_{{\gamma}}{}^{{\sigma}}{}_{{\delta}}\nabla_{\lambda}{\mathcal{A}}_{\sigma}+D_{5}\mathcal{B}_{{\alpha}}{}^{{\lambda}}{}_{{\beta}}{\mathcal{A}}_{\sigma}\nabla_{\lambda}\mathcal{B}_{{\gamma}}{}^{{\sigma}}{}_{{\delta}}+D_{6}\mathcal{B}_{{\alpha}}{}^{{\lambda}}{}_{{\beta}}{\mathcal{A}}_{\gamma}\nabla_{\lambda}{\mathcal{A}}_{\delta}+D_{7}\mathcal{B}_{{\alpha}}{}^{{\lambda}}{}_{{\beta}}{\mathcal{A}}_{\lambda}\nabla_{\gamma}{\mathcal{A}}_{\delta}+E_{1}\nabla_{\rho}\mathcal{B}_{{\alpha}}{}^{{\rho}}{}_{{\beta}}\nabla_{\sigma}\mathcal{B}_{{\gamma}}{}^{{\sigma}}{}_{{\delta}}+E_{2}\nabla_{\rho}\mathcal{B}_{{\alpha}}{}^{{\rho}}{}_{{\beta}}\nabla_{\gamma}{\mathcal{A}}_{\delta}+F_{1}\mathcal{B}_{{\alpha}}{}^{{\mu}}{}_{{\beta}}\mathcal{B}_{{\gamma}}{}^{{\sigma}}{}_{{\delta}}\mathcal{B}_{{\mu}}{}^{{\lambda}}{}_{{\rho}}\mathcal{B}_{{\sigma}}{}^{{\rho}}{}_{{\lambda}}+F_{2}\mathcal{B}_{{\alpha}}{}^{{\mu}}{}_{{\beta}}\mathcal{B}_{{\gamma}}{}^{{\nu}}{}_{{\lambda}}\mathcal{B}_{{\delta}}{}^{{\lambda}}{}_{{\rho}}\mathcal{B}_{{\mu}}{}^{{\rho}}{}_{{\nu}}+F_{3}\mathcal{B}_{{\nu}}{}^{{\mu}}{}_{{\lambda}}\mathcal{B}_{{\mu}}{}^{{\nu}}{}_{{\alpha}}\mathcal{B}_{{\beta}}{}^{{\lambda}}{}_{{\gamma}}{\mathcal{A}}_{\delta}+F_{4}\mathcal{B}_{{\alpha}}{}^{{\mu}}{}_{{\beta}}\mathcal{B}_{{\gamma}}{}^{{\nu}}{}_{{\delta}}{\mathcal{A}}_{\mu}{\mathcal{A}}_{\nu}\biggr{]}.$

One key feature of this purely affine formulation is that the number of terms
that can go to the action is limited, due to the restriction of the _index
structural analysis_ (for more details see Ref. castillo-
felisola18_einst_gravit_from_polyn_affin_model ; castillo-
felisola20_emerg_metric_geodes_analy_cosmol ). This feature contrast strongly
with other gravitational theories. For example the Starobinsky model adds a
correction $\alpha\mathcal{R}^{2}$ to the Einstein–Hilbert action, but there
is no constraint on the number of possible correction that can go the action.
Whereas in Polynomial Affine Gravity one has constraints coming from index
structural analysis. This is what we call the _rigidity_ of the model.

The field equations for the model can be obtained by the standard process of
optimisation, but with the subtlety of preserving the symmetries of the
fields, e.g. the $\mathcal{B}$-field is traceless. Since the action depends on
the first derivatives of the fields, in principle there is no need of
considering affine analogous to the Gibbons–Hawking–York term.

In order to solve the field equations, it is customary to start from an
ansatz. For the sake of simplicity, we consider an affine connection which is
compatible with the cosmological principle, which is parameterised as follows,

$\displaystyle{\Gamma_{t}{}^{t}{}_{t}}$ $\displaystyle=f(t),$
$\displaystyle\Gamma_{i}{}^{t}{}_{j}$ $\displaystyle=g(t)S_{ij},$
$\displaystyle\Gamma_{t}{}^{i}{}_{j}$ $\displaystyle=h(t)\delta^{i}_{j},$
$\displaystyle\Gamma_{i}{}^{j}{}_{k}$ $\displaystyle=\gamma_{i}{}^{j}{}_{k},$
$\displaystyle\mathcal{B}_{\theta}{}^{r}{}_{\varphi}$
$\displaystyle=\psi(t)r^{2}\sin\theta\sqrt{1-\kappa r^{2}},$
$\displaystyle\mathcal{B}_{r}{}^{\theta}{}_{\varphi}$
$\displaystyle=\frac{\psi(t)\sin\theta}{\sqrt{1-\kappa r^{2}}},$
$\displaystyle\mathcal{B}_{r}{}^{\varphi}{}_{\theta}$
$\displaystyle=\frac{\psi(t)}{\sqrt{1-\kappa r^{2}}\sin\theta},$
$\displaystyle\mathcal{A}_{t}$ $\displaystyle=\eta(t),$

where $S_{ij}$ and $\gamma_{i}{}^{j}{}_{k}$ are the three dimensional rank two
symmetric tensor and connection compatible with the cosmological symmetries.
The function $f$ can be set to zero by a suitable re-parameterisation of $t$.
Given such choice, the function $f$ will not be present in the dynamic
equations.

Using the above ansatz, the cosmological field equations are given by

$0=\left(B_{3}(2\kappa+gh+\dot{g})+2B_{4}(gh-\dot{g})+2D_{6}\eta
g-2F_{3}\psi^{2}\right)\psi$ $0=\left(\left(B_{3}-2B_{4}\right)\eta\psi-
C_{1}(2\eta\psi-\dot{\psi})\right)g$ $0=\left(B_{3}+2B_{4}\right)\eta
g\psi+2C_{1}\left(\kappa\psi+4gh\psi-g\dot{\psi}-\psi\dot{g}\right)+2\psi^{3}\left(-D_{1}+2D_{2}-D_{3}\right)$
$0=B_{3}\left((h\psi-\dot{\psi})\eta-\psi\dot{\eta}\right)-2B_{4}\left((-h\psi-\dot{\psi})\eta-\psi\dot{\eta}\right)+C_{1}\left(4h^{2}\psi+2\psi\dot{h}-\ddot{\psi}\right)+D_{6}\eta^{2}\psi$
$0=B_{3}(2\kappa+gh+\dot{g})\eta+2B_{4}(gh-\dot{g})\eta+C_{1}(2\kappa
h+4gh^{2}+2g\dot{h}-\ddot{g})+6h\psi^{2}(-D_{1}+2D_{2}-D_{3})+D_{6}\eta^{2}g-6F_{3}\eta\psi^{2}$

##  3 Final remarks

Despite the complexity of the field equations, solving the first couple of
equations provide four scenarios (given their factorisation). However, the
remaining equations from those scenarios are still very loose, making very
difficult to extract useful information from them. This push us to consider
other alternatives.

Another approach is to ask for compatibility between the connection and the
emergent metric tensor structure coming from either the Ricci tensor
$\mathcal{R}_{\mu\nu}$ or the torsion-descent Popłaswki metric
$\mathcal{P}_{\mu\nu}=\mathcal{A}_{\mu}\mathcal{A}_{\nu}+\mathcal{B}_{\mu}{}^{\lambda}{}_{\sigma}\mathcal{B}_{\nu}{}^{\sigma}{}_{\lambda}$
poplawski14_affin_theor_gravit . This additional condition together with the
field equations allows us to find analytically solutions.

In the case of the Popłaswki, we have an _affine scale factor_ given by the
$\psi(t)$ function, which serves as analogue of the scale factor $a(t)$ from
the Friedman–Robertson–Walker metric. In that sense, we are able to obtain
solutions as (Anti) de Sitter, where the parameter playing the role of
cosmological constant comes from an integration constant. This type of
solution is considered to be an eternal inflationary scenario, for
$\ddot{a}>0$. It might be argued that in this purely affine formulation of
gravity we can have a _geometric inflation_ , in the sense that our model does
not include matter.

Two of the most consolidated inflationary scenarios are the power-law scale
factor and the Starobinsky model. Those models require the introduction of a
self-interacting scalar field. A relevant question is whether or not
phenomenological viable scenarios of _geometric inflation_ can be found in the
context of polynomial affine gravity, even in the absence of an energy-
momentum tensor. That question is part of the problems that drive our
research.

## References

  * (1) D. Hilbert, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 1915, 395 (1915)
  * (2) A. Einstein, Sitzungsber. Preuss. Akad. Wiss. 315, 778 (1915)
  * (3) F.C. Klein (2008), 0807.3161, http://arxiv.org/abs/0807.3161v1
  * (4) T.A. Ivey, J.M. Landsberg, _Cartan for Beginners_ , Vol. 61 of _Graduate Studies in Mathematics_ (AMS, 2003)
  * (5) L. Eisenhart, _Non-Riemannian geometry_ (AMS, 1927)
  * (6) F.W. Hehl, J.D. McCrea, E.W. Mielke, Y. Ne’eman, Phys. Rep. 258, 1 (1995), gr-qc/9402012
  * (7) H. Weyl, Mathematische Zeitschrift 2, 384 (1918)
  * (8) H. Weyl, Sitzungsber. Königlish Preuss. Akad. Wiss. p. 465 (1918)
  * (9) E. Cartan, Ann. Ec. Norm. Super. 40, 325 (1923)
  * (10) E. Cartan, Ann. Ec. Norm. Super. 41, 1 (1924)
  * (11) E. Cartan, Ann. Ec. Norm. Super. 42, 17 (1925)
  * (12) O. Veblen, _Projective General Relativity_ (Springer, 1933)
  * (13) A. Einstein, Sitzungsber. Preuss. Akad. Wiss. pp. 137–140 (1923)
  * (14) A.S. Eddington, _The mathematical theory of relativity_ (CUP, 1923)
  * (15) O. Castillo-Felisola, A. Skirzewski, Rev. Mex. Fis. 61, 421 (2015), 1410.6183
  * (16) O. Castillo-Felisola, A. Skirzewski, Class. Quant. Grav. 35, 055012 (2018), 1505.04634
  * (17) O. Castillo-Felisola, O. Orellana, J. Perdiguero, F. Ramírez, A. Skirzewski, A.R. Zerwekh, Eur. Phys. J. C 82, 8 (2022), 2107.07209
  * (18) O. Castillo-Felisola, B. Grez, O. Orellana, J. Perdiguero, F. Ramirez, A. Skirzewski, A.R. Zerwekh, Universe 8, 68 (2022)
  * (19) O. Castillo-Felisola, J. Perdiguero, O. Orellana, A.R. Zerwekh, Class. Quant. Grav. 37, 075013 (2020), 1908.06654
  * (20) N.J. Popławski, Gen. Rel. Grav. 46, 1625 (2014), 1203.0294