Of course this construction is closely related to the situation of Section 8.
In that case the partition function of ${\cal N}=4$ dyons in a string
compactification on $\mathrm{K3}\times{\mathbb{T}}^{2}$ is proportional to the
generating function of the elliptic genera of the symmetric products of K3s,
via (8.4). In this Section we will focus on elliptic genera from a different
perspective.

The elliptic genus (9.1) has an interesting equivariant (or twisted)
generalization. Given an element $g\in\mathbb{M}_{24}$, we have

$\mathscr{E}_{g}(\tau,z)={\rm tr}_{{\cal
H}_{RR}}\left(g\,(-1)^{J_{0}+\overline{J}_{0}}\,y^{J_{0}}\,q^{L_{0}-\frac{c}{24}}\,\overline{q}^{\overline{L}_{0}-\frac{\overline{c}}{24}}\right)\
.$ (9.4)

In this case a version of (9.2) holds, where the dimension of the
representations are replaced by ${\rm tr}_{R_{n}}\,g$ (e similarly for the
virtual representations)555We have been a bit imprecise with the definitions
since K3 sigma models have a moduli space. Implicitly we are always choosing a
point in the CFT moduli space where the subgroup of $\mathbb{M}_{24}$
generated by $g$ is a symmetry of the Hilbert space.. In particular we can
rewrite (9.4) as

$\mathscr{E}_{g}(\tau,z)=\frac{\theta_{1}^{2}(\tau,z)}{\eta^{3}(\tau)}\left(\chi_{g}\,\mu(\tau,z)+H_{g}(\tau)\right)\,,$
(9.5)

where $\chi_{g}={\rm tr}_{R}\,g$ is the character of the defining
representation $R_{n}$ of $\mathbb{M}_{24}$ and

$\mu(\tau,z)=\frac{{\,{\rm
i}\,}y^{1/2}}{\theta_{1}(\tau,z)}\sum_{n\in{\mathbb{Z}}}\,(-1)^{n}\frac{y^{n}\,q^{n(n+1)/2}}{1-y\,q^{n}}\,,$
(9.6)

is an Appell-Lerch sum. The Mathieu McKay-Thompson series $H_{g}(\tau)$ are
weight $\frac{1}{2}$ mock modular forms, with shadow
$\chi_{g}\,\eta(\tau)^{3}$, and can be expanded as

$H_{g}(\tau)=q^{-\frac{1}{8}}\left(-2+\sum_{n=1}^{\infty}q^{n}{\rm
tr}_{K_{n}}g\right)\,,$ (9.7)

where $K=\bigoplus_{n=1}^{\infty}\,K_{n}$ is an infinite dimensional
$\mathbb{M}_{24}$ module [59, 60, 61, 62]. When $g={1\\!\\!1}$, (9.7) is the
generating function of the degeneracies of massive irreducible representations
of the worldsheet ${\cal N}=4$ conformal algebra, as they appear in (9.1) and
(9.2).

One can see that $\mathscr{E}_{g}(\tau,z)$ only depends on the conjugacy class
$[g]$ of $g\in\mathbb{M}_{24}$. The conjugacy classes lead to 21 distinct
elliptic genera; they are labelled by $g=\ $1A, 2A, 2B, 3A, 3B, 4A, 4B, 4C,
5A, 6A, 6B, 7AB, 8A, 10A, 11A, 12A, 12B, 14AB, 15AB, 21AB, 23AB. The notation
is standard and used in [57]. For each conjugacy class $[g]$, Table 2 of [57]
gives $\chi_{g}$ and the non universal ingredients needed to construct the
function $H_{g}(\tau)$.

The equivariant elliptic genus has the Fourier expansion

$\mathscr{E}_{g}(\tau,z)=\sum_{n\geq 0,l\in{\mathbb{Z}}}\ c_{g}(n,l)\
q^{n}\,y^{l}\ .$ (9.8)

Due to modularity one has $c_{g}(n,l)=c_{g}(4n-l)$. These coefficients satisfy
$c_{g}(k)=\mathrm{str}_{\hat{K}_{k}}\ g$, where $\hat{K}_{k}$ is the degree
$k$ component of a certain infinite dimensional ${\mathbb{Z}}$-graded
$\mathbb{M}_{24}$ module, obtained from $K$, and whose precise form is known
but not relevant for the present discussion. The coefficients $c_{g}(n)$ thus
have a precise relation to the irreducible representations of
$\mathbb{M}_{24}$. From the determination of the exact coefficients, their
distribution and asymptotic formulae, one gathers essential informations about
the representation theory of $\mathbb{M}_{24}$ and its realization as a vertex
operator algebra.

Now we want to address the problem from a different angle and use our
formalism to understand topological features of the distributions of the
coefficients of the equivariant elliptic genera $\mathscr{E}_{g}(\tau,z)$.

We construct a point cloud ${\mathsf{X}}^{g}$ whose elements have the form
$\mathsf{x}=(\log|c_{g}(n,l)|,n,l)$ for each conjugacy class associated to an
element $g\in\mathbb{M}_{24}$. Using the results of [57] we construct
explicitly $\mathscr{E}_{g}(\tau,z)$ and compute the first non-vanishing
coefficients up to $q^{20}$. Again the coefficients are growing very fast, and
the logarithm of the degeneracies makes it easier to capture the essential
features using persistent topology. In all cases the point cloud consists of
259 points. For each $[g]$ listed as above, we construct the
$\mathbb{N}$-persistence modules
$H_{i}\left(\mathsf{VR}_{\epsilon}({\mathsf{X}}^{g});{\mathbb{Z}}_{2}\right)$
and compute their barcodes.

We collect all the results in Appendix A. All the conjugacy classes present
similar topological features. It is quite likely that the explanation for this
fact is modularity, which as we have seen in the previous Section produces
strong constraints on the allowed persistent homology classes. We stress that
we are not suggesting any relation whatsoever between the coefficients of
different elliptic genera; apparently modularity has a subtle way of showing
up in the topological characterization of point clouds. It would be very
interesting to pursue this issue further, if one general grounds given the
information of modularity one can predict persistent features in the
distribution of barcodes. In this paper we will have nothing more to say on
the issue.

Taking at face value that modularity constrains the persistent homology
groups, there are still interesting more subtle feature in the barcode
spectra. Indeed, even at the qualitative level, the barcode distributions
associated with the elliptic genera $\mathscr{E}_{g}(\tau,z)$, fall into two
families.

The first family contains the conjugacy classes: 2B, 3B, 4A, 4C, 6B, 10A, 12A,
12B, 21AB; the second family the remaining classes. What singles out the first
class is that, despite having a similar form, the barcodes for the zeroth and
first homology are consistently more short-lived than those in the second
family. In the Betti number zero distribution, the persistent classes all die
for values of the proximity parameter $\epsilon\sim 2$. In degree one the
features are less clear-cut, the barcode distribution tends to be centered
around smaller values of $\epsilon$, and most persistent classes disappear
before $\epsilon\sim 3$. Another striking feature is the number of simplices
involved in the computation of persistent homology, which is not apparent just
by looking at the Figures in Appendix A. All the classes in the second family,
namely 1A, 2A, 3A, 4B, 5A, 6A, 7AB, 8A, 11A, 14AB, 15AB and 23AB generate a
number of simplices roughly comparable, around 65000-68000. On the other hand
the number of simplices generated in the Vietoris-Rips complex for the first
family fluctuates wildly, roughly in the range 38000-135000. This indications
are very qualitative, but from the behavior of the persistent homology it is
natural to wonder if there is any particular property that distinguish these
families.

Indeed there is. The above classes in the first family are precisely the only
conjugacy classes in $\mathbb{M}_{24}$ which do not have a representative in
the subgroup $\mathbb{M}_{23}$ of $\mathbb{M}_{24}$. Some classes of
$\mathbb{M}_{23}$ are singled out by Mukai’s theorem, which states that given
a finite subgroup of the group of symplectic automorphisms of a K3 surface, it
can be embedded into $\mathbb{M}_{23}$. The other classes in $\mathbb{M}_{23}$
do not have a geometrical interpretation.

At the technical level there is a simple explanation of the appearance of
these two families. The above classes in the first family all appear to have
$\chi_{g}=0$, which therefore acts as a relative shift of the coefficients
with respect to the elliptic genera of the second family. What is however
interesting here is that, with a bit of hindsight, such behavior could have
been predicted just by looking at the equivariant elliptic genera, for example
computed via conformal field theory, without any knowledge of the
representation theory of Mathieu’s groups. At a very qualitative level such
features are captured by persistent homology.

##  10 Discussion

In this paper we have taken the perspective of persistent homology to analyze
certain enumerative BPS invariants which arise in some physical and number
theoretical problems. The focus of this paper has been the comparison between
the topological features which characterize different distributions. We have
done so in different ways, by studying the same theory in different limits or
chambers, or by studying different theories in a similar situation. The
underlying theme has been that the distributions of supersymmetric states are
rich in topological information, which often can be clearly traced back to
physical properties of the system under consideration. This information is
qualitative, and concerns the overall structure of the space of supersymmetric
states. In a sense it is a measure of its topological complexity.

On a more practical level, this note has taken to the task of exemplifying the
uses of new methods based on persistent homology to study physical problems in
string and field theory; as well as calling the attention of the computational
topology community on the wealth of enumerative and number theoretical
datasets which arise in string/field theories. From a physical perspective, to
have a meaningful enumerative problem it is necessary to resort to protected
quantities, in this case supersymmetric states.

In the first part of the paper we have applied the tools of topological data
analysis to the wall-crossing phenomenon in ${\cal N}=2$ theories. The indices
of BPS protected quantities jump at walls of marginal stability and the effect
of the jump can be in principle computed using a wall-crossing formula. In
Sections 4 and 6 we have studied how the topological properties of the BPS
spectra change upon crossing certain walls. We have done so for a quantum
field theory, namely SU(3) ${\cal N}=2$ pure Yang-Mills, which is the simple
example of a theory containing wild chambers, and for string theory on the
conifold. Both theories exhibit physically interesting phenomena and we have
shown their impact on the persistent homology of the BPS spectra. In SU(3)
Yang-Mills we have analyzed the difference at the topological level between
the spectra in the weak coupling and wild chamber, where there is an
exponential growth of states and higher spin particles appear. At the level of
the barcodes the transition is striking, going from an ordered simple pattern
to an irregular (“wild”) distribution. This provides a qualitative measure of
the complexity of the spectrum in wild chambers. In the case of the conifold
the transitions are less apparent. In this case the physical features, for
example the fact that physical states are realized as very different objects
in the large radius geometrical chamber or in the noncommutative crepant
resolution chamber, or the presence of conjugation walls where the core of the
bound states changes, appear as rather small fluctuations in the BPS barcodes.
On the other hand the conifold is known to be a very simple geometry, where
all the BPS generating functions are known in closed form. It would be
interesting to carry on this analysis for more complicated non-compact
threefolds.

In Sections 7, 8 and 9 we have considered BPS spectra in compact Calabi-Yaus.
Interestingly the series of Donaldson-Thomas invariants corresponding to black
hole microstates in the large radius limit of one parameter compact
threefolds, all exhibit very similar persistent features. We interpret this
fact as another sense in which the microstate counting of black holes is
universal. Note that in general when talking about black hole universality,
one refers to the large charge limit, where the supergravity approximation is
reliable. In this limit the entropy obeys the area law and the subleading
corrections are expressed in terms of certain higher derivative terms [43].
Here we are looking at the first few hundreds of microstates; these data are
enough to extrapolate the large order behavior, which agrees with the
supergravity expectations [42]. Therefore some universal behavior is somewhat
expected; it is interesting to see it arise at the level of the barcodes.

Sections 8 and 9 are concerned with another aspect of the enumerative BPS
problems, more precisely in a few cases where the degeneracies of microstates
arise as coefficients of known functions. In this cases the relevant functions
have modular properties, and we have focused on the interplay between
modularity and topology. Indeed we have seen experimentally in a series of
examples that modularity strongly constrains the shape of the barcode
distributions. For example we have seen in Section 9 that the barcodes
associated with different elliptic genera present clear similarities. Note
that these functions are very different, in the sense that there is no
relation between the coefficients of the Fourier expansions of different
elliptic genera. On the other hand their topological features are quite
similar and we interpret this fact as a consequence of modularity, although we
don’t have any a priori argument based on topology. Quite remarkably, the
subtle differences in the barcode distributions can be explained in term of
the different families of conjugacy classes of the Mathieu group
$\mathbb{M}_{24}$.

We consider this note as a first step to understand the role that topological
data analysis can play in string/field theory problems. Clearly much is left
to be understood and more extensive computations of enumerative invariants are
needed to put the results of this paper on firmer grounds. Also there are
several other problems which seem to be amenable to a topological analysis as
we have done in this paper. We list here a few which we are currently
investigating:

  * •

It would be interesting to generalize this formalism to the study of string
vacua. In particular a lot of compactifications of string theory are known,
with ${\cal N}=2$ or ${\cal N}=1$ supersymmetry. The vacuum selection problem
consists in the choice of one compactification, or a class thereof, over the
others due to some particular features. It is natural to wonder if persistent
homology has anything to say about this, if for examples the distributions of
vacua with certain features have distinctive topological properties with
respect to others. A first step in this direction appears in [20].

  * •

A similar problem concerns the distribution of attractor points in the Calabi-
Yau moduli space. These points corresponds to black holes via the attractor
mechanism, and as argued in [63] they are deeply related to certain arithmetic
aspects of string compactifications.

  * •

One of the original motivations of this paper was if there is any particular
distinctive feature of BPS states in quantum field theory in the presence of
defects. When a theory is modified by the presence of a defect, new BPS states
appear, those which can bound to the defect. For example in six dimensional
topological quantum field theories such a modification is related to the
conjectural enumerative problem of Donaldson-Thomas invariants for moduli
spaces of parabolic sheaves [64]. In four dimensions, it was shown in [65]
that for theories in which a line defect can be engineered via laminations on
a curve [66], line defects come in distinct families, which are generated by
the action of a cluster algebra. Each family is generated by the repeated
action of a certain sequence of cluster transformations and can contain an
infinite number or a finite number of elements (including just one element).
It is natural to wonder if the analysis we have performed in Section 4 can be
of any use in the classification problem for defects.

  * •

Persistent homology is very closed in spirit to Morse theory. It would be
interesting to give a more physical description of persistent homology classes
via the correspondence between Morse theory and supersymmetric quantum
mechanics [67]. Viceversa, generalization of this correspondence, for example
along the lines of Floer theory, are likely to provide interesting variants of
persistent homology.

  * •

The equivariant elliptic genera that we have discussed in Section 9 are just a
small set of functions which arise in the field of moonshine. It would be very
interesting to generalize our formalism to other similar problems. For
example, it is an open problem to understand the distribution of coefficients
of the McKay-Thompson series for the monster moonshine modules. Partial
results on the asymptotics are in [56] and it would be interesting to see what
are the persistent features of these distributions.

We hope to report on these matters in the near future.

## Acknowledgements

I wish to thank the Applied Topology group at Stanford University for making
javaplex available in [22]. I also thank the theory division at CERN and IHES
for hospitality and support during the preparation of this note. This work was
partially supported by FCT/Portugal and IST-ID through UID/MAT/04459/2013,
EXCL/MAT-GEO/0222/2012 and the program Investigador FCT IF2014, under contract
IF/01426/2014/CP1214/CT0001.

##  Appendix A Moonshine barcodes

In this appendix we collect all of the barcodes computed in Section 9. Every
point cloud ${\mathsf{X}}^{g}$ is constructed out of vectors of the form
$\mathsf{x}=(\log|c_{g}(n,l)|,n,l)$ computed via the equivariant elliptic
genus $\mathscr{E}_{g}(\tau,z)$. The labels of the conjugacy classes $[g]$ are
indicated on top of each Figure. Each point cloud consists of 259 states. For
each conjugacy class $[g]$ we use the Vietoris-Rips complex to compute the
persistence modules
$H_{i}\left(\mathsf{VR}_{\epsilon}({\mathsf{X}}^{g});{\mathbb{Z}}_{2}\right)$.
This Appendix contains all the associated barcodes distributions, each Figure
labelled by the conjugacy class.

## References

  * [1] N. Seiberg and E. Witten, “Electric - magnetic duality, monopole condensation, and confinement in N=2 supersymmetric Yang-Mills theory,” Nucl. Phys. B 426 (1994) 19 [Nucl. Phys. B 430 (1994) 485] [hep-th/9407087].
  * [2] A. Strominger and C. Vafa, “Microscopic origin of the Bekenstein-Hawking entropy,” Phys. Lett. B 379 (1996) 99 [hep-th/9601029].
  * [3] S. H. Katz, A. Klemm and C. Vafa, “M theory, topological strings and spinning black holes,” Adv. Theor. Math. Phys. 3 (1999) 1445 [hep-th/9910181].
  * [4] N. Seiberg and E. Witten, “Gauge dynamics and compactification to three-dimensions,” In *Saclay 1996, The mathematical beauty of physics* 333-366 [hep-th/9607163].
  * [5] D. Gaiotto, G. W. Moore and A. Neitzke, “Four-dimensional wall-crossing via three-dimensional field theory,” Commun. Math. Phys. 299 (2010) 163 [arXiv:0807.4723 [hep-th]].
  * [6] M. Kontsevich and Y. Soibelman, “Stability structures, motivic Donaldson-Thomas invariants and cluster transformations,” arXiv:0811.2435 [math.AG].
  * [7] D. Gaiotto, G. W. Moore and A. Neitzke, “Wall-crossing, Hitchin Systems, and the WKB Approximation,” arXiv:0907.3987 [hep-th].
  * [8] M. Alim, S. Cecotti, C. Cordova, S. Espahbodi, A. Rastogi and C. Vafa, “$\mathcal{N}=2$ quantum field theories and their BPS quivers,” Adv. Theor. Math. Phys. 18 (2014) 27 [arXiv:1112.3984 [hep-th]].
  * [9] S. Cecotti and M. Del Zotto, “$Y$ systems, $Q$ systems, and 4D $\mathcal{N}=2$ supersymmetric QFT,” J. Phys. A 47 (2014) 47, 474001 [arXiv:1403.7613 [hep-th]].
  * [10] H. Edelsbrunner, D. Letscher, and A. Zomorodian, “ Topological persistence and simplification”, Discrete $\&$ Computational Geometry, 28 (2002), pp. 511–533.
  * [11] A. Zomorodian and G. Carlsson, “Computing persistent homology”, Discrete Comput. Geom., 33 (2005), pp. 249–274.
  * [12] G. Carlsson, “Topology and Data”, Bulletin of the American Mathematical Society, 46 (2009), pp. 255–308.
  * [13] H. Edelsbrunner and J. Harer, “Persistent Homology — A Survey”, in Surveys on Discrete and Computational Geometry, vol. 453, Amer Mathematical Society, 2008, p. 257.
  * [14] R. Ghrist, “Elementary Applied Topology”, Createspace, 2014. ed. 1.0.
  * [15] G. Carlsson, T. Ishkhanov, V. de Silva, and A. Zomorodian, “On the local behavior of spaces of natural images”, International Journal of Computer Vision, 76 (2008), pp. 1–12.
  * [16] J. M. Chan, G. Carlsson, and R. Rabadan, “Topology of viral evolution”, Proceedings of the National Academy of Sciences, 110 (2013), pp. 18566–18571.
  * [17] J. Binchi, E. Merelli, M. Rucco, G. Petri, and F. Vaccarino, “jHoles: A tool for understanding biological complex networks via clique weight rank persistent homology”, Electronic Notes in Theoretical Computer Science, 306 (2014), pp. 5 – 18. Proceedings of the 5th International Workshop on Interactions between Computer Science and Biology (CS2Bio14).
  * [18] M. Nicolau, A. J. Levine, and G. Carlsson, “Topology based data analysis identifies a subgroup of breast cancers with a unique mutational profile and excellent survival”, Proceedings of the National Academy of Sciences of the United States of America, 108 (2011), pp. 7265–7270.
  * [19] A. Port, I. Gheorghita, D. Guth, J. M. Clark, C. Liang, S. Dasu, M. Marcolli, “Persistent topology of syntax”, arXiv:1507.05134 [cs.CL]
  * [20] M. Cirafici, “Persistent homology and string vacua”, to appear
  * [21] All the matlab programs used in this paper and the accompanying datasets can be found at http://www.math.tecnico.ulisboa.pt/~cirafici/BPSbarcodesFiles/
  * [22] A. Tausz, M. Vejdemo-Johansson and H. Adams, “JavaPlex: a research software package for persistent (co)homology”. In Han Hong and Chee Yap, editors, Proceedings of ICMS 2014, volume 8592 of Lecture Notes in Computer Science, pages 129 136, 2014. Software available at http://appliedtopology.github.io/javaplex/
  * [23] D. Joyce and Y. Song, “A theory of generalized Donaldson-Thomas invariants,” Mem. Amer. Math. Soc. 217 (2012) 1–199 [arXiv:0810.5645 [math.AG]].
  * [24] J. Manschot, B. Pioline and A. Sen, “Wall Crossing from Boltzmann Black Hole Halos,” JHEP 1107 (2011) 059 [arXiv:1011.1258 [hep-th]].
  * [25] G. W. Moore, “Applications of the six-dimensional (2,0) theory to physical mathematics”, Felix Klein lectures, available at http://www.physics.rutgers.edu/~gmoore/FelixKleinLectureNotes.pdf
  * [26] M. Kontsevich and Y. Soibelman, “Wall-crossing structures in Donaldson-Thomas invariants, integrable systems and Mirror Symmetry,” arXiv:1303.3253 [math.AG].
  * [27] D. Galakhov, P. Longhi, T. Mainiero, G. W. Moore and A. Neitzke, “Wild Wall Crossing and BPS Giants,” JHEP 1311 (2013) 046 [arXiv:1305.5454 [hep-th]].
  * [28] K. Behrend, “Donaldson–Thomas invariants via microlocal geometry,” Ann. Math. 170 (2009) 1307–1338 [arXiv:math.AG/0507523].
  * [29] D. L. Jafferis and G. W. Moore, “Wall crossing in local Calabi Yau manifolds,” arXiv:0810.4909 [hep-th].
  * [30] E. Andriyash, F. Denef, D. L. Jafferis and G. W. Moore, “Bound state transformation walls,” JHEP 1203 (2012) 007 [arXiv:1008.3555 [hep-th]].
  * [31] B. Szendrői, “Noncommutative Donaldson-Thomas theory and the conifold,” Geom. Topol. 12 (2008) 1171–1202 [arXiv:0705.3419 [math.AG]].
  * [32] K. Nagao and H. Nakajima, “Counting invariants of perverse coherent sheaves and its wall-crossing,” Int. Math. Res. Notices 17 (2011) 3885–3938 [arXiv:0809.2992 [math.AG]].
  * [33] F. Denef and G. W. Moore, “Split states, entropy enigmas, holes and halos,” JHEP 1111 (2011) 129 [hep-th/0702146].
  * [34] A. Iqbal, N. Nekrasov, A. Okounkov and C. Vafa, “Quantum foam and topological strings,” JHEP 0804 (2008) 011 [hep-th/0312022].
  * [35] M. Cirafici, A. Sinkovics and R. J. Szabo, “Cohomological gauge theory, quiver matrix models and Donaldson-Thomas theory,” Nucl. Phys. B 809 (2009) 452 [arXiv:0803.4188 [hep-th]].
  * [36] M. Cirafici, A. Sinkovics and R. J. Szabo, “Instantons, Quivers and Noncommutative Donaldson-Thomas Theory,” Nucl. Phys. B 853 (2011) 508 [arXiv:1012.2725 [hep-th]].
  * [37] M. Cirafici and R. J. Szabo, “Curve counting, instantons and McKay correspondences,” J. Geom. Phys. 72 (2013) 54 [arXiv:1209.1486 [hep-th]].
  * [38] M. Cirafici, “On framed quivers, BPS invariants and defects”, contribution to the special volume BPS states, Hitchin systems and quivers, Confluentes Math. To appear.
  * [39] R. Gopakumar and C. Vafa, “M theory and topological strings. 1.,” hep-th/9809187.
  * [40] R. Gopakumar and C. Vafa, “M theory and topological strings. 2.,” hep-th/9812127.
  * [41] M. x. Huang, A. Klemm and S. Quackenbush, “Topological string theory on compact Calabi-Yau: Modularity and boundary conditions,” Lect. Notes Phys. 757 (2009) 45 [hep-th/0612125].
  * [42] M. x. Huang, A. Klemm, M. Marino and A. Tavanfar, “Black holes and large order quantum geometry,” Phys. Rev. D 79 (2009) 066001 [arXiv:0704.2440 [hep-th]].
  * [43] G. Lopes Cardoso, B. de Wit and T. Mohaupt, “Corrections to macroscopic supersymmetric black hole entropy,” Phys. Lett. B 451 (1999) 309 [hep-th/9812082].
  * [44] H. Ooguri, A. Strominger and C. Vafa, “Black hole attractors and the topological string,” Phys. Rev. D 70 (2004) 106007 [hep-th/0405146].
  * [45] G. Lopes Cardoso, B. de Wit, J. Kappeli and T. Mohaupt, “Black hole partition functions and duality,” JHEP 0603 (2006) 074 [hep-th/0601108].
  * [46] G. L. Cardoso, B. de Wit and S. Mahapatra, “Deformations of special geometry: in search of the topological string,” JHEP 1409 (2014) 096 [arXiv:1406.5478 [hep-th]].
  * [47] R. Dijkgraaf, E. P. Verlinde and H. L. Verlinde, “Counting dyons in N=4 string theory,” Nucl. Phys. B 484 (1997) 543 [hep-th/9607026].
  * [48] G. Lopes Cardoso, B. de Wit, J. Kappeli and T. Mohaupt, “Asymptotic degeneracy of dyonic N = 4 string states and black hole entropy,” JHEP 0412 (2004) 075 [hep-th/0412287].
  * [49] D. Shih and X. Yin, “Exact black hole degeneracies and the topological string,” JHEP 0604 (2006) 034 [hep-th/0508174].
  * [50] A. Dabholkar, F. Denef, G. W. Moore and B. Pioline, “Precision counting of small black holes,” JHEP 0510 (2005) 096 [hep-th/0507014].
  * [51] A. Dabholkar, S. Murthy and D. Zagier, “Quantum Black Holes, Wall Crossing, and Mock Modular Forms,” arXiv:1208.4074 [hep-th].
  * [52] G. L. Cardoso, M. Cirafici and S. Nampuri, “Indefinite theta functions for counting attractor backgrounds,” JHEP 1410 (2014) 17 [arXiv:1407.0197 [hep-th]].
  * [53] G. Lopes Cardoso, M. Cirafici, R. Jorge and S. Nampuri, “Indefinite theta functions and black hole partition functions,” JHEP 1402 (2014) 019 [arXiv:1309.4428 [hep-th]].
  * [54] F. Denef and M. R. Douglas, “Distributions of flux vacua,” JHEP 0405 (2004) 072 [hep-th/0404116].
  * [55] T. Eguchi, H. Ooguri and Y. Tachikawa, “Notes on the K3 Surface and the Mathieu group $M_{24}$,” Exper. Math. 20 (2011) 91 [arXiv:1004.0956 [hep-th]].
  * [56] J. F. R. Duncan, M. J. Griffin and K. Ono, “Moonshine,” arXiv:1411.6571 [math.RT].
  * [57] M. C. N. Cheng and J. F. R. Duncan, “The Largest Mathieu Group and (Mock) Automorphic Forms,” arXiv:1201.4140 [math.RT].
  * [58] M. R. Gaberdiel and R. Volpato, “Mathieu Moonshine and Orbifold K3s,” Contrib. Math. Comput. Sci. 8 (2014) 109 [arXiv:1206.5143 [hep-th]].
  * [59] M. C. N. Cheng, “K3 Surfaces, N=4 Dyons, and the Mathieu Group M24,” Commun. Num. Theor. Phys. 4 (2010) 623 [arXiv:1005.5415 [hep-th]].
  * [60] M. R. Gaberdiel, S. Hohenegger and R. Volpato, “Mathieu twining characters for K3,” JHEP 1009 (2010) 058 [arXiv:1006.0221 [hep-th]].
  * [61] M. R. Gaberdiel, S. Hohenegger and R. Volpato, “Mathieu Moonshine in the elliptic genus of K3,” JHEP 1010 (2010) 062 [arXiv:1008.3778 [hep-th]].
  * [62] T. Eguchi and K. Hikami, “Note on twisted elliptic genus of $K3$ surface,” Phys. Lett. B 694 (2014) 446 [arXiv:1008.4924 [hep-th]].
  * [63] G. W. Moore, “Arithmetic and attractors,” hep-th/9807087.
  * [64] M. Cirafici, “Defects in Cohomological Gauge Theory and Donaldson-Thomas Invariants,” arXiv:1302.7297 [hep-th].
  * [65] M. Cirafici, “Line defects and (framed) BPS quivers,” JHEP 1311 (2013) 141 [arXiv:1307.7134 [hep-th]].
  * [66] D. Gaiotto, G. W. Moore and A. Neitzke, “Framed BPS States,” Adv. Theor. Math. Phys. 17 (2013) 241 [arXiv:1006.0146 [hep-th]].
  * [67] E. Witten, “Supersymmetry and Morse theory,” J. Diff. Geom. 17 (1982) 661.


# A GCICA Grant-Free Random Access Scheme for M2M Communications in Crowded
Massive MIMO Systems

Huimei Han1, Lushun Fang1, Weidang Lu1, Wenchao Zhai2, Ying Li3, Jun Zhao4
Huimei Han, Lushun Fang, and Weidang Lu are with College of Information
Engineering, Zhejiang University of Technology, Hangzhou, Zhejiang, 310032,
P.R. China, Wenchao Zhai is with College of Information Engineering, China
Jiliang University, Hangzhou, Zhejiang, 310018, P.R. China, Ying Li is with
the State Key Lab of Integrated Services Networks, Xidian University, Xi’an,
710071, P.R. China, and Jun Zhao is with School of Computer Science and
Engineering, Nanyang Technological University, Singapore. (Email: {hmhan1215,
2111903007<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A high success rate of grant-free random access scheme is proposed to support
massive access for machine-to-machine communications in massive multiple-input
multiple-output systems. This scheme allows active user equipments (UEs) to
transmit their modulated uplink messages along with super pilots consisting of
multiple sub-pilots to a base station (BS). Then, the BS performs channel
state information (CSI) estimation and uplink message decoding by utilizing a
proposed graph combined clustering independent component analysis (GCICA)
decoding algorithm, and then employs the estimated CSIs to detect active UEs
by utilizing the characteristic of asymptotic favorable propagation of massive
MIMO channel. We call this proposed scheme as GCICA based random access
(GCICA-RA) scheme. We analyze the successful access probability, missed
detection probability, and uplink throughput of the GCICA-RA scheme. Numerical
results show that, the GCICA-RA scheme significantly improves the successful
access probability and uplink throughput, decreases missed detection
probability, and provides low CSI estimation error at the same time.

###### Index Terms:

Grant-free random access, independent component analysis (ICA), massive MIMO,
M2M communications.

##  I Introduction

THE machine-to-machine (M2M) communications are centered on the intelligent
interaction of user equipments (UEs) without human intervention, which is the
enabler for the Internet of Things (IoT) to achieve the envision of the
“Internet of Everything” [1, 2, 3, 4]. In recent years, M2M communications
develop rapidly and have been applied to many scenarios, such as smart
medical, smart vehicle, smart logistics, etc. Cisco visual networking index
and forecast predicts that there will be around 28.5 billion connected UEs by
2022 [5]. The massive multiple-input multiple-output (MIMO) technology, which
achieves significant improvements in energy and spectral efficiency and can
serve massive UEs in the same time-frequency resource, is well suited for M2M
communications [6, 7, 8].

For M2M communications in massive MIMO systems, random access procedure is a
critical step to initiate a data transmission [9, 10]. Since payload data is
usually in small size and the number of UEs in a cell is envisioned in the
order of hundreds or thousands in M2M communications [11], random access
procedure utilized in the long term evolution (LTE) network may induce
excessive signaling overhead and cannot support massive access [9, 12].

Researchers are exploring new random access schemes for M2M communications in
massive MIMO systems, and the grant-based random access schemes have been
proposed in recent years. E. Björnson et al. proposed a strongest-user
collision resolution (SUCRe) scheme, which allocates pilots to active UEs with
largest channel gain among the contenders [13]. However, the number of
successful accessing UEs decreases with the increase of the number of
contenders [14]. To improve the pilot resource utilization of the SUCRe
scheme, SUCR combined idle pilots access (SUCR-IPA) scheme was proposed in
[15], where the weaker UEs randomly select idle pilots to increase the number
of successful accessing UEs. A user identity-aided pilot access scheme was
proposed for massive MIMO with interleave-division multiple-access systems,
where the interleaver of each UE is available at the BS according to the one-
to-one correspondence between UE’s identity number (ID) and its interleaver
[16]. However, since the grant-based random access schemes require two
handshake processes between the base station (BS) and UEs, considering the
small packet transmission in M2M communication, such kind of random access
scheme will introduce heavy signaling overhead and low data transmission
efficiency.

To address this problem, the grant-free random access schemes have attracted
much attention in recent years, which allow active UEs to transmit their
pilots and uplink messages to the BS directly and perform activity detection,
channel state information (CSI) estimation, and uplink message decoding in one
shot. J. Ahn et al. proposed a Bayesian based random access scheme to detect
the UE’s activity and estimate the CSI jointly by utilizing the expectation
propagation algorithm, considering the BS with one antenna [17]. L. Liu et al.
proposed a approximate message passing (AMP) based grant-free scheme to
achieve the joint activity detection and CSI estimation for massive MIMO
systems [18]. However, this AMP-based grant-free random access scheme requires
long pilot sequence to achieve better performance, resulting in heavy pilot
overhead. To support massive connectivity with low access delay and overhead,
a novel grant-free random access scheme was proposed to detect active UEs and
uplink message without CSI estimation in one shot [19]. However, the CSI
estimation is needed for some M2M traffics for downlink beamforming in the
case of time-division duplex (TDD) model, such as the telemedicine surgery
requiring the downlink control signal transmission. The above-mentioned grant-
free random access schemes consider the single pilot structure. Jiang et al.
proposed to concatenate multiple orthogonal sub-pilots into a super pilot
sequence, and the super pilot sequence is utilized for activity detection and
CSI estimation [20]. Simulation results show that, compared to the single
pilot structure, this super pilot structure can improve the number of
successful UEs by utilizing their proposed detection approach. However, with
the increase of the number of active UEs, the probability that more UEs have
the same super pilot sequence increases, resulting in low row rank pilot
selection matrix and thus degrading the number of successful UEs. This
motivates interesting researches on how to support massive access by utilizing
the super pilot structure.

To achieve this goal, we propose a graph combined clustering independent
component analysis based random access (GCICA-RA) scheme. This scheme utilizes
the super pilot structure and mainly consists of two steps. During the first
step, all the active UEs concatenate their randomly selected sub-pilots to
obtain their super pilot sequences, and send their super pilot sequences and
modulated uplink messages to the BS. During the second step, the BS performs
CSI estimation and uplink message decoding jointly by utilizing a proposed
GCICA decoding algorithm, employs the estimated CSIs to detect active UEs by
utilizing the characteristic of asymptotic favorable propagation of massive
MIMO channel, and sends random access response (RAR) to active UEs. We analyze
the successful access probability, missed detection probability, and uplink
throughput of the GCICA-RA scheme. Numerical results show that, compared to
the random access scheme proposed in [20], the GCICA-RA scheme significantly
improves the successful access probability and uplink throughput, decreases
missed detection probability, and provides low CSI estimation error at the
same time.

The main contributions of this paper are summarized as follows.

  * •

We propose a GCICA algorithm to perform CSIs estimation and uplink message
decoding jointly. Specifically, based on the received sub-pilot signal, the BS
first utilizes a successive interference cancellation (SIC) algorithm to
estimate CSIs of active UEs, and then decodes the uplink message by utilizing
the estimated CSIs. Then, the BS utilizes a proposed CICA algorithm to decode
uplink message and estimate CSIs of UEs whose CSIs cannot be estimated by the
SIC algorithm. Simulation results show that, compared to the detection method
proposed in [20], the number of successful UEs can be improved even multiple
UEs having the same super pilot sequence.

  * •

For the SIC algorithm, we propose to employ the estimated CSIs without other
overhead to find edges connecting to a recovered variable node, by utilizing
the asymptotic favorable propagation of massive MIMO systems. This is
different from the existing method where indexes of selected pilots should be
inserted into the uplink message to find edges connected to a recovered
variable node [14], which degrades data transmission efficiency.

  * •

We analyze the upper bound for the successful access probability and uplink
throughput of the proposed GCICA-RA scheme, and derive the lower bound for the
missed detection probability of the proposed GCICA-RA scheme.

The remainder of this paper is organized as follows. Section II introduces
system model and the proposed GCICA-RA scheme. Section III describes SIC and
CICA algorithms in the the proposed GCICA-RA scheme. We present the
performance analysis in Section IV. Simulation results and the conclusion are
given in Section V and VI, respectively.

Notations utilized throughout this paper are described in Table I.

TABLE I: notations. Notations | Description
---|---
Italic letters | Scalars
Boldface lower-case | Vectors
Boldface upper-case letters | Matrices
$(\cdot)^{T}$ and $(\cdot)^{H}$ |  | The transpose and conjugate transpose
---
of a vector or a matrix
‘*’ |  Complex conjugate of a vector or a matrix
$\boldsymbol{x}_{i}$ |  The $i^{th}$ element of a vector $\boldsymbol{x}$
$\mathbb{R}$ | The set of all real numbers
$\mathcal{C}\mathcal{N}(\mu,\sigma^{2})$ |  | A circularly-symmetric complex Gaussian
---
distribution with mean $\mu$ and variance $\sigma^{2}$
$||\cdot||$ | The Euclidean norm of a vector
$|\cdot|$ | The cardinality of a set
arg($d$) |  The phase of complex $d$
$\lfloor\cdot\rceil$ | The rounding operation
$[\bm{x}]_{n}$ |  The $n^{th}$ element of vector $\bm{x}$
$\bm{X}(i)$ |  The $i^{th}$ column of matrix $\bm{X}$

##  II SYSTEM MODEL and The proposed GCICA-RA Scheme

###  II-A SYSTEM MODEL

In this paper, we consider massive MIMO communication systems in time-division
duplex (TDD) mode. There is a BS with $M$ antennas and $K$ single-antenna UEs
in a cell, where the number of active UEs in a random access procedure is
$N_{a}$.

In the proposed GCICA-RA scheme, each UE transmits its modulated uplink
message and super pilot to the BS directly. The BS performs CSI estimation and
uplink message decoding by utilizing a proposed GCICA decoding algorithm, and
then employs the estimated CSIs to detect active UEs. The frame structure of
the GCICA-RA scheme is shown in Fig. 1. Specifically, the super pilot of each
UE consists of $L$ consecutive sub-pilots, and the symbol length of each sub-
pilot is $\tau_{p}$. Therefore, the symbol length of each super pilot is
$L\tau_{p}$. In addition, the proposed GCICA decoding algorithm utilizes
multiple ICA classifiers to separate source signals. However, the phase of the
separated signals is unpredictable.

To solve this problem, each active UE inserts a reference symbol (termed RS)
with symbol length 1 into its uplink message [19]. Let
$\bm{v_{k}}=[v_{k}(1),v_{k}(2),\ldots,v_{k}({N_{\text{m}}})]^{{T}}$ denote the
modulated uplink message of UE $k$, where $N_{\text{m}}$ represents the symbol
length of the modulated uplink message. Then, $v_{k}(1)$,
$[v_{k}(2),\ldots,v_{k}(N_{\text{m}})]$ represent the RS and the payload data
with symbol length $N_{\text{PD}}$, respectively. Furthermore, since any
complex symbol can be represented as a real symbol [21], we employ the BPSK
modulation scheme in the proposed GCICA-RA scheme, i.e.
$\upsilon_{k}(i)=\left\\{{\pm 1}\right\\}$.

Figure 1: The frame structure of the proposed GCICA-RA scheme.
Figure 2: Block diagram of the proposed GCICA-RA scheme. Each active UE sends
its super pilot and modulated uplink message to the BS. The BS performs CSI
estimation and uplink message decoding by utilizing a proposed GCICA decoding
algorithm, employs the estimated CSIs to detect active UEs, and sends RAR
information to active UEs. In addition, the proposed GCICA decoding algorithm
utilizes a SIC algorithm to estimate CSIs and thus decodes the uplink message
of active UEs, performs the number of active UEs estimation, and utilizes a
proposed CICA algorithm to decode uplink message and estimate CSIs of active
UEs whose CSIs cannot be recovered by the SIC algorithm.

###  II-B GCICA-RA scheme description

As shown in Fig. 2, the proposed GCICA-RA scheme allows each UE to transmit
its modulated uplink message and super pilot to the BS directly. Then, the BS
performs CSI estimation and uplink message decoding by utilizing a proposed
GCICA decoding algorithm, employs the estimated CSIs to detect active UEs, and
sends RAR information to active UEs. In addition, the proposed GCICA decoding
algorithm first utilizes a SIC algorithm to estimate CSIs and decode the
uplink message of active UEs, performs the number of active UEs estimation,
and utilizes a proposed CICA algorithm to decode uplink messages and estimate
CSIs of active UEs whose CSIs cannot be recovered by the SIC algorithm. The
details are described as follows.

####  II-B1  Step 1: Super pilot and modulated uplink message transmission

There are $L$ sub-pilot phases, and each active UE randomly selects its sub-
pilot from a set of mutually orthogonal normalized pilots
$\bm{{{P}}_{{o}}}=\left\\{\bm{S_{1}},\bm{S_{2}},\cdots,\bm{S_{\tau_{p}}}\right\\}\in{\mathbb{R}^{{\tau_{p}}\times\tau_{p}}}$
during each sub-pilot phase. Concatenating these randomly selected sub-pilot
sequences, each active UE obtains its super pilot. Following the frame
structure as shown in Fig. 1, each active UE sends it super pilot and
modulated uplink message to the BS. In addition, the value of RS is the same
for all active UEs, and we set the RS to be 1, i.e., $\upsilon_{k}(1)=1\
\forall k$.

Then, let $\bm{g_{k}}={(g_{k}^{1},g_{k}^{2},\cdots,g_{k}^{M})^{{T}}}$ denote
the independent and identically distributed (i.i.d.) small-scale fading
channel from UE $k$ to the BS, i.e., $\bm{g_{k}}\sim\text{
}\mathcal{N}(0,{\bm{{{I}}_{M}}})$. Let ${\rho_{k}}$ denote UE $k$’s
transmitting power. Let ${{{\beta_{k}}}}$ stand for the channel gain between
UE $k$ and the BS, which can be known for UE $k$ before the random access
procedure [13]. Then, the received sub-pilot signal during the $l^{th}\
(l=1,2,\ldots,L)$ sub-pilot phase at the BS is

$\bm{Y_{\text{p}}^{l}}=\sum\limits_{t=1}^{{\tau_{p}}}\sum\limits_{k\in\bm{A_{t}^{l}}}{\sqrt{{\beta_{k}}{\rho_{k}}}}{{{\bm{g_{k}}(\bm{S_{t}})^{{T}}}}}+\bm{Z_{\text{p}}^{l}},~{}l=1,2,\ldots,L,$
(1)

where $\bm{A_{t}^{l}}$ stands for the set of UEs that select pilot $S_{t}$
during the $l^{th}$ sub-pilot phase,
$\bm{Z_{\text{p}}^{l}}\in{{\mathbb{R}}^{\text{ }M\times{\tau_{p}}}}$ denotes
the additive white noise matrix with each element following from a
distribution of $\mathcal{N}(0,{{\sigma}^{2}})$. In this paper, $\bm{Z}$ with
different subscripts or superscripts will follow the same distribution. In
addition, we consider ${{\rho_{k}}{{\beta}_{k}}}=1$ to make the received
signal from all active UEs have the same power [22].

The received superimposed uplink message
$\bm{Y_{\text{m}}}\in{{\mathbb{R}}^{\text{ }M\times{N_{\text{m}}}}}$ at the BS
is

$\bm{Y_{\text{m}}}=\sum\limits_{k=1}^{{N_{a}}}{{{\bm{g_{k}}(\bm{v_{k}})^{{T}}}}}+\bm{Z_{\text{m}}}.$
(2)

####  II-B2  CSI estimation, uplink message decoding, active UE detection and
RAR transmission

Based on the received sub-pilot signal and the superimposed uplink message,
the BS performs ${①}$ CSI estimation and uplink message decoding jointly, and
${②}$ active UE detection and RAR information transmission. The details are
described as follows.

${①}$ CSI estimation and uplink message decoding

The BS utilizes the proposed GCICA decoding algorithm to perform CSI
estimation and uplink message decoding. The details of the GCICA decoding
algorithm are described as follows.

The BS first estimates the CSI corresponding to pilot $\bm{S_{t}}$ during the
$l^{th}$ sub-pilot phase by utilizing least squares (LS) method, denoted by
$\bm{h_{t}^{l}}$

$\bm{h_{t}^{l}}=\bm{Y_{\text{p}}^{l}}\bm{(S_{t})^{*}}=\sum\limits_{k\in\bm{A_{t}^{l}}}{\bm{{g_{k}}}}+\bm{Z_{t}^{l}}.$
(3)

It can be seen from (3) that, during the $l^{th}$ sub-pilot phase, the
estimated CSI corresponding to pilot $\bm{S_{t}}$ is the sum of CSIs of active
UEs selecting pilot $\bm{S_{t}}$. So, considering all $L$ sub-pilot phases, a
bipartite graph can be used to describe the proposed GCICA-RA scheme, and thus
the SIC algorithm can be utilized to estimate UEs’ CSIs on the bipartite
graph. The procedure of the SIC algorithm is described in Section III-A in
detail. We use
$\bm{\hat{H}_{\text{SIC}}}=[\bm{\hat{h}_{\text{SIC}}^{1},\cdots,\hat{h}_{\text{SIC}}}^{N_{\text{SIC}}^{\text{s}}}]\in{\mathbb{R}}^{\text{
}M\times{N_{\text{SIC}}^{\text{s}}}}$ to denote the estimated CSIs via the SIC
algorithm, where $N_{\text{SIC}}^{\text{s}}$ is the number of estimated CSIs
via the SIC algorithm. Thus, the detected uplink message can be obtained by
utilizing LS method

$\bm{\hat{x}_{\text{SIC}}^{i}}=\frac{(\bm{\hat{h}_{\text{SIC}}^{i}})^{H}\bm{{Y_{\text{m}}}}}{(\bm{\hat{h}_{\text{SIC}}^{i}})^{H}(\bm{\hat{h}_{\text{SIC}}^{i}})},1\leq
i\leq N_{\text{SIC}}^{\text{s}}.$ (4)

After demodulating and decoding the detected uplink message
$\bm{\hat{x}_{\text{SIC}}^{i}}$, we obtain the corresponding decoded uplink
message $\bm{\hat{v}_{\text{SIC}}^{i}}$.

Then, based on the estimated CSIs corresponding to all pilots during the
$l^{th}$ sub-pilot phase, the BS estimates the number of active UEs by
utilizing the estimation method proposed in [19], denoted by $\hat{{N_{a}}}$,
which is given by

$\begin{array}[]{l}\hat{{N_{a}}}={\sum\limits_{t=1}^{\tau_{p}}\bm{\frac{(h_{t}^{l})^{H}h_{t}^{l}}{M}}}\to{\sum\limits_{t=1}^{\tau_{p}}|\bm{A_{t}^{l}}|}={N_{a}},M\to\infty.\end{array}$
(5)

We can observe from (5) that $\hat{{N_{a}}}$ approaches $N_{a}$ with large
$M$. Therefore, the estimated number of active UEs whose CSIs cannot be
obtained via the SIC algorithm is

$\hat{N}_{\text{r}}=\hat{N}_{a}-N_{\text{SIC}}^{\text{s}}.$ (6)

Finally, by utilizing $\hat{N}_{\text{r}}$,
$\bm{\hat{v}_{\text{SIC}}^{i}}(1\leq i\leq N_{\text{SIC}}^{\text{s}})$, and
$\bm{{Y_{m}}}$, the BS employs a proposed CICA algorithm to obtain the CSIs
and decoded uplink message of active UEs whose CSIs cannot be obtained via the
SIC algorithm, denoted by
$\bm{\hat{H}_{\text{CICA}}}=[\bm{\hat{h}_{\text{CICA}}^{1}},\cdots,\bm{\hat{h}_{\text{CICA}}}^{{\hat{N}_{\text{r}}}}]\in{\mathbb{R}}^{M\times{\hat{N}_{\text{r}}}}$
and $\bm{\hat{v}_{\text{CICA}}^{i}}$, respectively. The details of the CICA
algorithm are descried in Section III-B.

${②}$ Active UE detection and RAR information transmission

The BS detects an active UE by judging whether its CSI can be estimated (i.e.,
whether the corresponding estimated CSI is a valid CSI estimation), and this
UE realizes whether it has been detected by the BS in a distributed manner.
The details are described as follows.

We utilize
$\hat{\bm{H}}=[\bm{\hat{H}_{\text{SIC}}},\bm{\hat{H}_{\text{CICA}}}]$ to
denote the estimated CSIs, and assume that $\hat{\bm{H}}(1)$ is a valid CSI
estimation. More specifically, $\hat{\bm{H}}(1)$ is the estimated CSI for UE
$m$ (i.e., $\hat{\bm{H}}(1)=\bm{g_{m}+\Delta g_{m}}$, where $\bm{\Delta
g_{m}}$ denotes the estimation error and takes small values). By correlating
$\hat{\bm{H}}(1)$ with the estimated CSI corresponding to pilot $\bm{S_{t}}$
during the $l^{th}$ sub-pilot phase (i.e.,
$\bm{h_{t}^{l}},~{}t=1,2,\ldots,\tau_{p}$), we have

$\displaystyle\frac{(\bm{h_{t}^{l}})^{H}\hat{\bm{H}}(1)}{M}$
$\displaystyle=\frac{\left(\sum\limits_{k\in\bm{A_{t}^{l}}}{\bm{{g_{k}}}}+\bm{Z_{t}^{l}}\right)^{H}(\bm{g_{m}+\Delta
g_{m}})}{M}$ (7)
$\displaystyle=\frac{(\bm{{g_{m}}})^{H}\bm{g_{m}}}{M}+\frac{\sum\limits_{k\in\bm{A_{t}^{l}},k\neq
m}({\bm{{g_{k}}})^{H}\bm{g_{m}}}}{M}+$
$\displaystyle\frac{\sum\limits_{k\in\bm{A_{t}^{l}}}\bm{{(\bm{{g_{k}}})^{H}}\Delta
g_{m}}}{M}+\frac{(\bm{Z_{t}^{l}})^{H}(\bm{g_{m}+\Delta g_{m}})}{M}$
$\displaystyle\xlongrightarrow[\quad]{(a)}1,M\to\infty,{\rm{\text{if}}}~{}t~{}{\rm{\text{is
the index of}}}$ $\displaystyle{\rm{\text{the pilot selected by UE}}}~{}m$
$\displaystyle\frac{(\bm{h_{t}^{l}})^{H}\hat{\bm{H}}(1)}{M}$
$\displaystyle=\frac{\sum\limits_{k\in\bm{A_{t}^{l}}}({\bm{{g_{k}}}})^{H}\hat{\bm{H}}(1)+(\bm{Z_{t}^{l}})^{H}\hat{\bm{H}}(1)}{M}$
$\displaystyle\xlongrightarrow[\quad]{(b)}0,M\to\infty,{\rm{otherwise}}$

where (a) follows from the fact that, when $M$ goes into infinity, terms
$\frac{(\bm{{g_{m}}})^{H}\bm{g_{m}}}{M}=1$,
$\frac{\sum\nolimits_{k\in\bm{A_{t}^{l}},k\neq
m}{\bm{{g_{k}}}^{H}}\bm{g_{m}}}{M}=0,\frac{(\bm{Z_{t}^{l}})^{H}(\bm{g_{m}+\Delta
g_{m}})}{M}=0$ by utilizing the asymptotic favorable propagation of massive
MIMO channel, term
$\frac{\sum\nolimits_{k\in\bm{A_{t}^{l}}}{(\bm{{g_{k}}})^{H}}\bm{\Delta
g_{m}}}{M}$ is close to $0$ since elements in $\bm{\Delta g_{m}}$ are small
values, and (b) is derived based on the random matrix theory. In addition, it
is easy to derive that, when $\hat{\bm{H}}(i)$ is not a valid CSI estimation,
$\frac{\left(\bm{h_{t}^{l}}\right)^{H}\hat{\bm{H}}(i)}{M}$ is close to 0 by
utilizing the random matrix theory. Therefore, the BS utilizes the following
rule to judge whether $\hat{\bm{H}}(i)$ is a valid CSI estimation: for the
$l^{th}\left(l\in\\{1,...,L\\}\right)$ sub-pilot phase, if there only exits an
element in set
$[\frac{(\bm{h_{l}^{1}})^{H}\hat{\bm{H}}(i)}{M},\cdots,\frac{(\bm{h_{l}^{t}})^{H}\hat{\bm{H}}(i)}{M},\cdots,\frac{(\bm{h_{l}^{\tau_{p}}})^{H}\hat{\bm{H}}(i)}{M}]$
close to 1, $\hat{\bm{H}}(i)$ is a valid CSI estimation because each active UE
only selects one pilot during each sub-pilot phase; otherwise,
$\hat{\bm{H}}(i)$ is not a valid CSI estimation.

After obtaining all valid CSIs, the BS generates and broadcasts the RAR
information $\bm{V}=\sum\nolimits_{u\in\bm{B}}(\bm{\hat{H}(u)})^{H}$ to all
active UEs where $\bm{B}$ stands for the index set of valid CSI estimations.

The received RAR at UE $m$ is

${R_{m}}=\sqrt{{\beta_{m}}}\bm{{V}{g_{m}}}+{Z_{d}}.$ (8)

Then, we have

$\begin{array}[]{l}\frac{{{R_{m}}}}{{\sqrt{{\beta_{m}}}M}}=\frac{{\sqrt{{\beta_{m}}}\bm{{{V}}{g_{m}}}}}{{\sqrt{{\beta_{m}}}M}}+\frac{{{Z_{d}}}}{{\sqrt{{\beta_{m}}}M}}\vspace{1ex}\\\
=\frac{{\sqrt{{\beta_{m}}}\left({\sum\limits_{u\in\bm{B}}{{(\hat{\bm{H}}(u))}^{H}}{\bm{g_{m}}}}\right)}}{{\sqrt{{\beta_{m}}}M}}+\frac{{{Z_{d}}}}{{\sqrt{{\beta_{m}}}M}}\vspace{1ex}\\\
\xlongrightarrow[\quad]{(a)M\to\infty}\left\\{\begin{array}[]{l}1{\quad{\text{\rm{if
UE}}}~{}m~{}{\text{{\rm{has been detected}}}}}\\\ \\\
0{\quad\rm{otherwise}}\end{array}\right.\end{array}$ (9)

where (a) is derived by utilizing the asymptotic favorable propagation of
massive MIMO channel, similar to (7). Based on the value of
$\frac{{{R_{m}}}}{{\sqrt{{\beta_{m}}}M}}$, UE $m$ can realize whether it has
been detected. If not, UE $m$ will reattempt the random access procedure in
the upcoming random access slot.

##  III  SIC and CICA algorithms

###  III-A SIC Algorithm

As we described in Section II-B, a bipartite graph can be used to describe the
proposed GCICA-RA scheme, and thus the SIC algorithm can be utilized to
estimate UEs’ CSIs on the bipartite graph. In this subsection, we first
present the bipartite graph representation, and then describe the SIC
algorithm in detail.

####  III-A1  Bipartite Graph Representation

A bipartite graph consists of variable nodes, factor nodes, and edges
connecting variable nodes and factor nodes. We utilize active UEs and
available pilots during each sub-pilot phase to represent variable nodes and
factor nodes, respectively. Furthermore, if pilot $\bm{S_{t}}$ is selected by
UE $k$ during sub-pilot phase $l$, we use $t^{l}$ to denote the corresponding
factor node, and there exists an edge connecting the variable node $k$ with
the factor node $t^{l}$. The CSI of UE $k$ is marked on edges connected to UE
$k$. Note that, since we use the power control mechanism to make
${{\rho_{k}}{{\beta}_{k}}}=1$, we only need to estimate the small-scale fading
channel. So, we only mark the small-scale fading coefficients between the BS
and active UE. Furthermore, the degree of a node is referred to the number of
edges connected to this node.

Example: Fig. 3(a) illustrates a bipartite graph based on 5 active UEs, 2 sub-
pilot phases, and 3 available pilots. During each sub-pilot phase, active UEs
(represented by circles) are variable nodes and all available pilots
(represented by squares) are factor nodes. Since UE 1 and UE 2 select pilot
$\bm{S_{1}}$ during sub-pilot phase $1$, UE 1 and UE 2 are connected to factor
node $1^{1}$ by a black solid line. Since UE 3 selects pilot $\bm{S_{2}}$
during sub-pilot phase $1$, UE 3 is connected to factor node $2^{1}$ by a
black solid line. Since UE 4 and UE 5 select pilot $\bm{S_{3}}$ during sub-
pilot phase $1$, UE 4 and UE 5 are connected to factor node $3^{1}$ by a black
solid line. Since UE 1 and UE 3 select pilot $\bm{S_{1}}$ during sub-pilot
phase $2$, UE 1 and UE 3 are connected to factor node $1^{2}$ by a black solid
line. Since UE 2, UE 4, and UE 5 all select pilot $\bm{S_{3}}$ during sub-
pilot phase $2$, UE 2, UE 4 and UE 5 are connected to factor node $3^{2}$ by a
black solid line.

Figure 3: An example of SIC process.

####  III-A2  SIC Algorithm

The BS employs SIC algorithm to estimate UE’s CSI on the bipartite graph, and
then utilizes the estimated CSIs to decode the uplink message. In the
following, we briefly discuss the SIC algorithm, present the pseudo-code of
the SIC algorithm in Algorithm 1, and give an example to show the SIC
algorithm intuitively.

By performing the autocorrelation operation on CSI corresponding to each pilot
during each sub-pilot phase, we have

$\displaystyle\frac{(\bm{h_{t}^{l}})^{H}\bm{h_{t}^{l}}}{M}$
$\displaystyle=\frac{\sum\limits_{k\in\bm{A_{t}^{l}}}(\bm{{g_{k}}})^{H}{\bm{{g_{k}}}}+\left({\sum\limits_{k\in\bm{A_{t}^{l}}}{\bm{{g_{k}}}}+\bm{Z_{t}^{l}}}\right)^{H}\bm{Z_{t}^{l}}}{M}$
(10) $\displaystyle\xlongrightarrow[\quad]{(a)}|\bm{A_{t}^{l}}|,M\to\infty,$

where (a) can be easily derived by utilizing the characteristic of asymptotic
favorable propagation of massive MIMO channel.

(10) shows that, with large $M$,
$\frac{(\bm{h_{t}^{l}})^{H}\bm{h_{t}^{l}}}{M}$ approaches the number of active
UEs selecting pilot $\bm{S_{t}}$ during the $l^{th}$ sub-pilot phase. In other
words, $\frac{(\bm{h_{t}^{l}})^{H}\bm{h_{t}^{l}}}{M}$ represents the degree of
factor node $t^{l}$.

During the first SIC iteration, the BS first finds factor nodes with degree 1
by utilizing (10). A factor node with degree 1 indicates that there is only
one UE selecting the corresponding sub-pilot and the corresponding CSI can be
recovered. Then, the BS finds other selected sub-pilots of UEs whose CSIs are
recovered. More specifically, for each CSI corresponding to factor node with
degree larger than 1, after subtracting each recovered CSI, the BS utilizes
(10) to compute the new degree of this factor node. If the new degree is less
than the previous degree, this means that the sub-pilot corresponding to this
factor node is selected by active UE corresponding to this recovered CSI. We
subtract the recovered CSI from the CSI corresponding to this factor node and
update the new degree as the previous degree of this factor node. This is
equivalent to delete the edges between this active UE and its selected sub-
pilots on the bipartite graph. Thus, a new bipartite graph can be obtained,
and SIC iteration continues until there is no factor node with degree 1.

Since each active UE selects its sub-pilots randomly, there exits the case of
UE $k$ having no pilot collision during multiple sub-pilot phases. This
results in the case of the recovered multiple CSIs corresponding to the same
UE. Luckily, by utilizing the asymptotic favorable propagation of massive MIMO
channel, we can find these CSIs from the recovered CSIs and select one of them
as the recovered CSI of this UE. For example, if estimated $\bm{\hat{h}_{1}}$
and $\bm{\hat{h}_{2}}$ correspond to the same active UE, by utilizing the
asymptotic favorable propagation of massive MIMO channel,
$\frac{(\bm{\hat{h}_{1}})^{H}\bm{\hat{h}_{2}}}{M}$ is close to 1 with large
$M$. Otherwise, $\frac{(\bm{\hat{h}_{1}})^{H}\bm{\hat{h}_{2}}}{M}$ is close to
0 with large $M$. Finally, the recovered CSIs via the SIC algorithm are
$\bm{\hat{H}_{\text{SIC}}}=[\bm{\hat{h}_{\text{SIC}}^{1},\cdots,\hat{h}_{\text{SIC}}}^{N_{\text{SIC}}^{\text{s}}}]\in{\mathbb{R}}^{\text{
}M\times{N_{\text{SIC}}^{\text{s}}}}$.

Input : $\bm{h_{t}^{l}},(t=1,2,\ldots,\tau_{p};l=1,2,\ldots,L)$

Output : $\bm{\hat{H}_{\text{SIC}}}$: a set of recovered CSIs

1 $\bm{h}_{1,t}^{l}=\bm{h}_{t}^{l},flag=1,\bm{\hat{H}_{\text{SIC}}}=\phi,i=1$;

2 while  _$(flag\neq 0)$_ do

3 $flag=0$;

4 for  _all
$\bm{h}_{i,t}^{l},l\in\\{1,2,\ldots,L\\},t\in\\{1,2,\ldots,\tau_{p}\\}$ _ do

5 Compute the degree of the corresponding factor node by utilizing (10);

6 if  _the degree is 1_ then

7 Set $flag=1$;

8 $\bm{\hat{H}_{\text{SIC}}}=\bm{\hat{H}_{\text{SIC}}}\cup\bm{h}_{i,t}^{l}$;

9 for  _all
$\bm{h}_{i,m}^{n},m\in\\{1,2,\ldots,L\\},n\in\\{1,2,\ldots,\tau_{p}\\}$ _ do

10 if  _ $\frac{(\bm{h}_{i,m}^{n})^{H}\bm{h}_{i,t}^{l}}{M}$ is less than the
previous degree _ then

11 update the channel response:
$\bm{h_{i+1,m}^{n}}=\bm{h_{i,m}^{n}}-\bm{h}_{i,t}^{l}$;

12

13 end if

14 else

15 update the channel response: $\bm{h_{i+1,m}^{n}}=\bm{h_{i,m}^{n}}$;

16

17 end if

18

19 end for

20

21 end if

22

23 end for

24 $i=i+1$;

25

26 end while

Return : After deleting the recovered CSIs corresponding to the same UEs in
set $\bm{\bm{H}_{\text{SIC}}}$ and selecting one of them as the recovered CSI,
we have
$\bm{\hat{H}_{\text{SIC}}}=[\bm{\hat{h}_{\text{SIC}}^{1},\cdots,\hat{h}_{\text{SIC}}^{N_{\text{SIC}}^{s}}}]$

Algorithm 1 The SIC Algorithm

Algorithm 1 gives the detailed iteration process of SIC algorithm, where
$\bm{h_{i,j}^{k}}$ to represent the CSI corresponding to factor node $j^{k}$
during the $i^{th}$ iteration.

To better understand the SIC algorithm, an example is given to illustrate the
SIC algorithm. The details are described as follows:

During the first iteration (Fig. 3(a)), the BS finds that the degree of factor
node $2^{1}$ is 1 based on (10). So the BS can directly recover the CSI of UE
3, and incorporates this CSI (i.e., $\bm{h_{1,2}^{1}}$) into the set
$\bm{\hat{H}_{\text{SIC}}}$. Then, after subtracting $\bm{h_{\text{SIC}}}^{1}$
from CSIs corresponding to factor node with degree larger than 1, the BS finds
that the new degree of factor node $1^{2}$ is less than the previous degree
according to (10), and thus realizes that the other index of sub-pilot
selected by this UE is 1 in the second sub-pilot phase. Subtracting the
recovered CSI from the CSI corresponding to factor node $1^{2}$ (i.e, deleting
the lines connected to variable node UE 3), a new bipartite graph can be
obtained. Same as above operation, the BS can recover the CSI of UE 1 during
the second iteration (i.e., $\bm{h_{2,1}^{2}}$ Fig.3(b)) and the CSI of UE 2
during the third iteration (i.e., $\bm{h_{3,1}^{1}}$ Fig.3(c)). In the end,
since there is no factor node with degree 1 (Fig.3(d)), SIC ends. Finally, by
utilizing the asymptotic favorable propagation of massive MIMO channel, we
find that the recovered CSIs (i.e., $\bm{h_{1,2}^{1}}$, $\bm{h_{2,1}^{2}}$,
and $\bm{h_{3,1}^{1}}$) belong to different UEs. Hence, we have
$\bm{\hat{h}_{\text{SIC}}^{1}}=\bm{h_{1,2}^{1}}$,
$\bm{\hat{h}_{\text{SIC}}^{2}}=\bm{h_{2,1}^{2}}$, and
$\bm{\hat{h}_{\text{SIC}}^{3}}=\bm{h_{3,1}^{1}}$.

###  III-B CICA algorithm

We can see from Fig. 3 that, for these 5 active UEs, we can only recover CSIs
of UE 1, UE 2, and UE 3 by utilizing the SIC algorithm, and UE 4 and UE 5 fail
to access the network. We propose a CICA algorithm to further obtain the CSIs
and decoded uplink messages of UE 4 and UE 5. We briefly describe the CICA
algorithm in Algorithm 2. In the following, we use the term “the remaining
active UEs” to represent active UEs whose CSIs cannot be recovered by the SIC
algorithm for simplify representation. The proposed CICA algorithm first ${①}$
utilizes multiple ICA classifiers to separate the source signals and obtain
the decoded signals, and then ${②}$ employs a proposed clustering algorithm to
obtain estimated CSIs and decoded uplink messages of the remaining active UEs.
The details are described as follows.

Input :
$\bm{Y_{\text{m}}^{{}^{\prime}}}=[\bm{Y_{\text{m}}}(1),\bm{Y_{\text{m}}}(2),\cdots,\bm{Y_{\text{m}}}(M)]^{{T}}$,
$\hat{N}_{r}$

Output : $\bm{\hat{v}_{\text{CICA}}}$: a set of decoded uplink message,
$\bm{\hat{H}_{\text{CICA}}}$: a set of estimated CSIs;

1
${r}=N_{\text{I}},\bm{\hat{v}_{\text{CICA}}}=\phi,\bm{\hat{H}_{\text{CICA}}}=\phi$;

2 while  _$(r\neq 0)$_ do

3 $\bm{f_{i}}\leftarrow$ Select $\hat{N}_{r}$ row vectors from
$\bm{Y_{m}^{{}^{\prime}}}$ randomly to obtain $\bm{Y_{\text{ICA}}^{i}}$, and
then feed it to the $i^{th}$ ICA classifier to obtain separated source signal
$\bm{f_{i}}=[\bm{f_{i}^{1}},\bm{f_{i}^{2}},...,\bm{f_{i}^{\hat{N}_{r}}}]$;

4 $r=r-1$;

5

6 end while

7for  _all $i\in\\{1,2,...,N_{\text{I}}\\},t\in\\{1,2,...,{\hat{N}_{r}}\\}$ _
do

8 $\bm{\hat{f_{i}^{t}}}\leftarrow$ Use the first decoding symbol in
$\bm{f_{i}^{t}}$ to solve the ICA phase ambiguity problem via (12);

9 $\text{DEC}(\bm{\hat{f_{i}^{t}}})\leftarrow$ After the first symbol in
$\bm{\hat{f_{i}^{t}}}$ is deleted, $\bm{\hat{f_{i}^{t}}}$ is fed into the
decoder to obtain the decoded signal $\text{DEC}(\bm{\hat{f_{i}^{t}}})$;

10 $j\leftarrow$ Judge which class $\text{DEC}(\bm{\hat{f_{i}^{t}}})$ belongs
to by using the proposed clustering algorithm (assume it belongs to class $j$)
;

11 $\bm{C_{j}}\leftarrow$ Incorporate this signal to the decoded signal set of
the corresponding class ;

12 $\bm{\hat{d}_{j}}\leftarrow$ Add the remodulated signal via (13);

13 end for

14for  _$j\in\\{1,2,...,{\hat{N}_{r}}\\}$ _ do

15 $\bm{\hat{v}_{\text{CICA}}^{j}}=\text{sign}(\bm{\hat{d}_{j}})$,
$\bm{\hat{v}_{\text{CICA}}}=\bm{\hat{v}_{\text{CICA}}}\cup\bm{\hat{v}_{\text{CICA}}^{j}}$;

16
$\bm{\hat{h}_{\text{CICA}}^{j}}=\frac{\bm{{Y_{m}^{{}^{\prime}}}({{\hat{\bm{v}}_{\text{CICA}}^{j}}})^{*}}}{({{\hat{\bm{v}}_{\text{CICA}}^{j}}})^{H}({{\hat{\bm{v}}_{\text{CICA}}^{j}}})}$
;

17
$\bm{\hat{H}_{\text{CICA}}}=\bm{\hat{H}_{\text{CICA}}}\cup\bm{\hat{h}_{\text{CICA}}^{j}}$
;

18 end for

Return :
$\bm{\hat{v}_{\text{CICA}}=[\hat{v}_{\text{CICA}}^{1},\hat{v}_{\text{CICA}}^{2},\cdots,\hat{v}_{\text{CICA}}^{\hat{N}_{r}}]}$;

$\bm{\hat{H}_{\text{CICA}}}=[\bm{\hat{h}_{\text{CICA}}^{1}},\cdots,\bm{\hat{h}_{\text{CICA}}}^{{\hat{N}_{\text{r}}}}]$.

Algorithm 2 CICA decoding algorithm

${①}$ ICA classifiers

The BS first obtains the sum of received superimposed uplink message of the
remaining active UEs, denoted by
$\bm{Y_{\text{m}}^{{}^{\prime}}}\in\mathbb{R}^{M\times N_{\text{m}}}$.

$\begin{array}[]{l}\bm{Y_{\text{m}}^{{}^{\prime}}}=\bm{Y_{\text{m}}}-\sum\limits_{i=1}^{N_{\text{SIC}}^{\text{s}}}{\bm{h_{\text{SIC}}^{i}}(\bm{\hat{v}_{\text{SIC}}^{i}})^{T}}\\\
={\sum\limits_{i=1}^{{N_{\text{r}}}}{{{\bm{g_{i}}(\bm{v_{i}})^{\text{T}}}}}}+\bm{Z_{\text{m}}^{{}^{\prime}}}\\\
=\bm{G}_{\text{r}}\bm{V}_{\text{r}}+\bm{Z_{\text{m}}^{{}^{\prime}}}.\end{array}$
(11)

where ${N_{\text{r}}}$ denotes the number of the remaining active UEs. Recall
that we can obtain the estimated value of ${N_{\text{r}}}$ (i,e.,
$\hat{N}_{\text{r}}$) by utilizing (6).

Then, by utilizing the ICA method in [19], the BS utilizes $N_{\text{I}}$ ICA
classifiers to separate the source signals $\bm{V_{\text{r}}}$ in
$\bm{Y_{\text{m}}^{{}^{\prime}}}$. More specifically, the BS randomly selects
$\hat{N}_{\text{r}}$ different row vectors from
$\bm{Y_{\text{m}}^{{}^{\prime}}}=[\bm{Y_{\text{m}}^{{}^{\prime}}}(1),\cdots,\bm{Y_{\text{m}}^{{}^{\prime}}}(M)]^{T}$
as the input signal of the $i^{th}$ ICA classifier, and the $i^{th}$ ICA
classifier outputs the separated source signal, denoted by
$\bm{f_{i}}=[\bm{f_{i}^{1}},\cdots,\bm{f_{i}^{\hat{N}_{\text{r}}}}],~{}i=1,\ldots,N_{\text{I}}$.
We further use the first symbol in each separated source signal to solve the
phase ambiguity problem introduced by the ICA classifier as follows

$\bm{\hat{f_{i}^{t}}}={f_{i}^{1}}(t)\times\bm{f_{i}^{t}},(i\in\\{1,2,...,N_{\text{I}}\\},t\in\\{1,2,...,{\hat{N}_{\text{r}}}\\}).$
(12)

Finally, after deleting the first symbol in $\bm{\hat{f_{i}^{t}}}$,
$\bm{\hat{f_{i}^{t}}}$ is fed into the channel decoder to obtain the decoded
signal $\text{DEC}(\bm{\hat{f_{i}^{t}}})$. Thus, we can obtain the decoded
signal from the $i^{th}$ ICA classifier, i.e.,
$\text{DEC}(\bm{f_{i}})=[\text{DEC}(\bm{f_{i}^{1}}),\cdots,\text{DEC}(\bm{f_{i}^{t}}),\cdots,\text{DEC}(\bm{f_{i}^{\hat{N}_{\text{r}}}})]$.

${②}$ Clustering algorithm

The order of the ICA classifier output is undetermined, and thus we cannot
sure which UE each decoded signal from the $i^{th}$ ICA classifier belongs to.
To solve this problem, we propose a clustering algorithm to cluster the
decoded signals of the same UE from these $N_{\text{I}}$ ICA classifiers, and
further obtain the estimated CSIs and decoded uplink messages of the remaining
active UEs.

We regard each remaining active UE as a class, and the decoded signals of one
UE belong to the same class. Therefore, there are ${\hat{N}_{\text{r}}}$
classes, indexed by $1,2,\ldots,{\hat{N}_{\text{r}}}$. Actually, since each
ICA classifier outputs decoded signals of ${\hat{N}_{\text{r}}}$ UEs, we can
utilize decoded signals from any ICA classifiers as the first sample in each
class. Here, we assume that ${\hat{N}_{\text{r}}}$ decoded signals from the
first ICA classifier (i.e., $\text{DEC}(\bm{f_{1}})$) are the first sample in
each class accordingly.

The BS first judges which class each decoded signal from the
$i^{th}~{}(i=2,\ldots,{\hat{N}_{\text{I}}})$ ICA classifier belongs to. More
specifically, considering that the fewer the number of different bits between
two decoded signals, the more likely these two decoded signals come from the
same UE. The BS compares the number of different bits between the decoded
signal $\text{DEC}(\bm{f_{i}^{t}}),~{}(t=1,2,\ldots,{\hat{N}_{\text{r}}})$ and
each element in $\text{DEC}(\bm{f_{1}})$, and judges
$\text{DEC}(\bm{f_{i}^{t}})$ belonging to the class which has the minimum
number of different bits with $\text{DEC}(\bm{f_{i}^{t}})$. We use
$\bm{C_{j}}$ to denote a set of decoded signals in the $j^{th}$ class.

Then, for the $j^{th}$ class, the BS remodulates each element in set
$\bm{C_{j}}$ by utilizing the BPSK modulator and performs the sum operation as
follows

$\bm{\hat{d}_{j}}=\sum\limits_{i=1}^{|\bm{C_{j}}|}2\times{\bm{C_{j}}(i)}-1,$
(13)

where $2\times{\bm{C_{j}}(i)}-1$ refers to the remodulated signal.

Finally, the decoded uplink message of the remaining active UE can be obtained
by

${{\hat{\bm{v}}_{\text{CICA}}^{j}}}=\text{sign}(\bm{\hat{d}_{j}}),~{}j=1,\ldots,{\hat{N}_{\text{r}}},$
(14)

and the estimated CSIs of the remaining active UEs can be obtained by

$\bm{\hat{h}_{\text{CICA}}^{j}}=\frac{\bm{{Y_{m}^{{}^{\prime}}}({{\hat{\bm{v}}_{\text{CICA}}^{j}}})^{*}}}{({{\hat{\bm{v}}_{\text{CICA}}^{j}}})^{H}({{\hat{\bm{v}}_{\text{CICA}}^{j}}})},~{}j=1,\ldots,{\hat{N}_{\text{r}}}.$
(15)

##  IV  Performance Analysis

In this section, we analyze the performance of the proposed GCICA-RA scheme,
including the successful access probability, missed detection probability, and
the uplink throughput.

###  IV-A Successful access probability analysis

Based on the procedure of the proposed GCICA-RA scheme, if the CSI of one UE
can be obtained, this UE is a detected UE and its uplink message can also be
decoded. This means that this UE can access the network successfully. The
GCICA-RA scheme first utilizes the SIC algorithm to obtain the CSIs of active
UEs, and then employs the proposed CICA algorithm to obtain the CSIs of active
UEs whose CSI cannot be recovered by the SIC algorithm. Let $S_{\text{GCICA-
RA}}$, $S_{\text{SIC}}$ and $S_{\text{CICA}}$ denote the number of successful
UEs obtained by the proposed GCICA-RA scheme, the SIC algorithm and the CICA
algorithm, respectively. Then, the successful access probability of the
proposed GCICA-RA scheme is given by

$P_{s}=\frac{S_{\text{GCICA-
RA}}}{N_{a}}=\frac{S_{\text{SIC}}+S_{\text{CICA}}}{N_{a}}.$ (16)

Next, we discuss how to compute terms $S_{\text{SIC}}$ and $S_{\text{CICA}}$
in (16).

####  IV-A1 $S_{\rm{SIC}}$ computation

To compute $S_{\text{SIC}}$, we follow the two assumptions below [14]:

_Assumption 1_ : If an active UE has no pilot collision with other active UEs,
its CSI can be recovered successfully.

_Assumption 2_ : When an active UE selects the same pilot with other $l$
active UEs, if CSIs of these $l$ UEs can be recovered, its CSI can be
recovered successfully.

Let $P_{\text{fail}}$ denote the probability that UE’s CSI cannot be recovered
by the SIC algorithm. Then, we have

$S_{\text{SIC}}=N_{a}\left(1-P_{\text{fail}}\right).$ (17)

The and-or tree principle is used to study random process, and we employ this
principle to compute $P_{\text{fail}}$. In the following, we first derive the
degree distributions of the proposed GCICA-RA scheme, and then discuss how to
utilize the and-or tree principle to compute $P_{\text{fail}}$ based on the
derived degree distributions.

${①}$ Degree distributions of the proposed GCICA-RA scheme derivation

Let $\Lambda_{l}$ and $\Psi_{l}$ denote the probability that the degree of a
variable node is $l$ and the probability that the degree of a factor node is
$l$, respectively. Then, the corresponding polynomial representations are
given by

$\Lambda(x)\triangleq\sum\nolimits_{l}{\Lambda_{l}x^{l}},\quad\Psi(x)\triangleq\sum\nolimits_{l}{\Psi_{l}x^{l}}.$
(18)

The proposed GCICA-RA scheme allows each active UE to select its sub-pilot
during each sub-pilot phase, and there are $L$ sub-pilot phases. Therefore,
the degree of each variable node is $L$, i.e.,

$\Lambda_{l}=\left\\{\begin{aligned} 1,&&l=L,\\\ 0,&&{\rm{otherwise}}.\\\
\end{aligned}\right.$ (19)

In addition, since each active UE randomly selects its sub-pilot during each
sub-pilot phase, the event that $l$ active UEs select the same sub-pilot
follows a Bernoulli distribution with parameters $N_{a}$ and ${1}/{\tau_{p}}$.
Thus, the probability of a factor node with degree $l$ (i.e., $\Psi_{l}$) is

$\Psi_{l}=\binom{N_{a}}{l}\left(\dfrac{1}{\tau_{p}}\right)^{l}\left(1-\dfrac{1}{\tau_{p}}\right)^{\left({N_{a}}-{l}\right)}.$
(20)

Let $\lambda_{l}$ denote the probability of an edge connecting to a variable
node with degree $l$, and $\rho_{l}$ denote the probability of an edge
connecting to a factor node with degree $l$. Then, we transfer the node-
perspective presentations in (19) and (20) into the edge perspective
presentations, respectively

$\lambda_{l}=\dfrac{\Lambda_{l}l}{\sum\nolimits_{l}{\Lambda_{l}l}},\quad\rho_{l}=\dfrac{\Psi_{l}l}{\sum\nolimits_{l}{\Psi_{l}l}},$
(21)

and the polynomial representations are

$\lambda(x)\triangleq\sum\nolimits_{l}{\lambda_{l}x^{l-1}},\rho(x)\triangleq\sum\nolimits{\rho_{l}x^{l-1}}.$
(22)

${②}$ $P_{\text{fail}}$ computation

Figure 4: Probability updates

Considering a factor node with degree $l$ (i.e., there are $l$ UEs selecting
the corresponding pilot), an edge connecting to this factor node can be
recovered when the other $l-1$ edges connecting to this factor node have been
recovered, as shown in Fig. 4(a). Similarly, as shown in Fig. 4(b),
considering a variable node with degree $l$ (i.e., UE corresponding to this
variable node sends $l$ sub-pilots to the BS), an edge connecting to a
variable node with degree $l$ can be recovered when at least one of the other
$l-1$ edges have been recovered. Let $n(m)$ denote the probability that an
edge connecting to a factor (variable) node with degree $l$ cannot be
recovered. Then, we have

$n=1-\left(1-m\right)^{l-1},$ (23)

and

$m=n^{l-1}.$ (24)

By averaging $m$ and $n$ over the edge distributions, the evolution of the
average erasure probability during the $i^{th}~{}(1\leq i\leq I)$ iteration
where $I$ is the maximum number of iterations, can be written as

$\begin{array}[]{rcl}n_{i}&=&\sum\nolimits_{l}{\rho_{l}\left(1-\left(1-m_{i-1}\right)^{l-1}\right)\hfill}\\\
\\\ &=&1-\rho\left(1-m_{i-1}\right),\hfill\\\ \\\
m_{i}&=&\sum\nolimits_{l}{\lambda_{l}n_{i}^{l-1}}=\lambda\left(n_{i}\right).\end{array}$
(25)

The update process of $m_{i}$ and $n_{i}$ is shown in Fig. 4(c). Specifically,
given initial values $m_{0}$ and $n_{0}$, $m_{i}$ and $n_{i}$ can be derived
by (25).

It can be seen from Fig. 4(b) that, a variable node with degree $l$ is unknown
when all of these $l$ edges cannot be recovered in the $I^{th}$ iteration. So,
the probability of the variable node being unknown is $(n_{I})^{l}$. By
averaging $(n_{I})^{l}$ over the degree distributions of variable nodes,
$P_{\text{fail}}$ can be calculated by [23]

$P_{\text{fail}}=\sum\nolimits_{l}{\Lambda_{l}(n_{I})^{l}}=\Lambda(n_{I}).$
(26)

Substituting (26) into (17), we have

$S_{\text{SIC}}=N_{a}\left(1-P_{\text{fail}}\right)=N_{a}(1-\Lambda(n_{I})).$
(27)

####  IV-A2 $S_{\rm{CICA}}$ computation

The number of successful UEs obtained by the CICA algorithm can be written as

$S_{\text{CICA}}=(N_{a}-S_{\text{SIC}})P_{\text{CICA}}^{s},$ (28)

where $P_{\text{CICA}}^{s}$ denotes the probability that the CSI of one UE can
be obtained successfully.

Based on the procedure of the CICA algorithm described in Subsection III-B, we
can see that, the CICA algorithm first obtains the decoded uplink message via
(14), and then obtains UE’s CSI via (15). Thus, if the uplink message of one
UE can be decoded successfully, the CSI of this UE can be recovered
successfully. Then, $P_{\text{CICA}}^{s}$ can be computed by

$P_{\text{CICA}}^{s}=\left(1-P_{e}\right)^{N_{\text{m}}}.$ (29)

where $P_{e}$ is the bit error rate (BER) of the proposed CICA algorithm,
which is larger than that using maximum likelihood (ML) decoding method [24],
and can be written as

$\displaystyle{P_{e}}$
$\displaystyle\geq\dfrac{1}{{\sqrt{\pi}}}{\int}_{\sqrt{\rho_{1}\beta_{1}\sum\nolimits_{m=1}^{M}{{{({g_{1,m}})}^{2}}}}}^{\infty}{{e^{-{x^{2}}}}dx}\hfill$
$\displaystyle\mathop{=}\limits^{(a)}\dfrac{1}{{\sqrt{\pi}}}{\int}_{0}^{\infty}{{\int}_{\sqrt{y}}^{\infty}{\dfrac{{{y^{M/2-1}}{e^{-y/2}}}}{{{2^{M/2}}\Gamma(\frac{M}{2})}}}{e^{-{x^{2}}}}}dxdy,$
(30)

where step (a) follows from the fact that
$\sum\nolimits_{m=1}^{M}{{{({g_{1,m}})}^{2}}}$ is a central chi-squared
distribution with degree of freedom $M$ and that ${\rho_{1}\beta_{1}}=1$, and
$\Gamma(z)$ is the gamma function with parameter $z$ which is given by

$\Gamma(z)=\int_{0}^{\infty}{{t^{z-1}}{e^{-t}}dt}.$ (31)

Substituting (30) into (29), we can derive the upper bound for
$P_{\text{CICA}}^{s}$,

$P_{\text{CICA}}^{s}\leq
P_{\text{CICA}}^{\text{U}}=\leavevmode\resizebox{305.33788pt}{}{$\left(1-\dfrac{1}{{\sqrt{\pi}}}{\int}_{0}^{\infty}{{\int}_{\sqrt{y}}^{\infty}{\dfrac{{{y^{M/2-1}}{e^{-y/2}}}}{{{2^{M/2}}\Gamma(\frac{M}{2})}}}{e^{-{x^{2}}}}}dxdy\right)^{N_{\text{m}}}$},$
(32)

and thus derive the upper bound for $S_{\text{CICA}}$

$\displaystyle S_{\text{CICA}}$ $\displaystyle\leq S_{\text{CICA}}^{\text{U}}$
$\displaystyle=P_{\text{CICA}}^{\text{U}}(N_{a}-S_{\text{SIC}})\hfill$
$\displaystyle=\leavevmode\resizebox{352.31624pt}{}{$\left(1-\dfrac{1}{{\sqrt{\pi}}}{\int}_{0}^{\infty}{{\int}_{\sqrt{y}}^{\infty}{\dfrac{{{y^{M/2-1}}{e^{-y/2}}}}{{{2^{M/2}}\Gamma(\frac{M}{2})}}}{e^{-{x^{2}}}}}dxdy\right)^{N_{\text{m}}}$}$
$\displaystyle\quad\times(N_{a}-S_{\text{SIC}}).$ (33)

Substituting (33), (27) into (16), we can get the upper bound for the
successful access probability of the proposed GCICA-RA scheme

$\displaystyle P_{s}\leq P_{s}^{\text{U}}$
$\displaystyle=\frac{\leavevmode\resizebox{352.31624pt}{}{$(N_{a}\Lambda(n_{I}))\left(1-\dfrac{1}{{\sqrt{\pi}}}{\int}_{0}^{\infty}{{\int}_{\sqrt{y}}^{\infty}{\dfrac{{{y^{M/2-1}}{e^{-y/2}}}}{{{2^{M/2}}\Gamma(\frac{M}{2})}}}{e^{-{x^{2}}}}}dxdy\right)^{N_{\text{m}}}$}}{N_{a}}$
$\displaystyle\qquad+\frac{N_{a}(1-\Lambda(n_{I}))}{N_{a}}.$ (34)

###  IV-B Missed detection probability analysis

The missed detection probability refers to the ratio of the number of
undetected active UEs to the total number of active UEs, which can be
calculated as

$P_{\text{md}}=\dfrac{N_{a}-S_{\text{GCICA-RA}}}{N_{a}}=1-P_{s}.$ (35)

Substituting (34) into (35), we can get the lower bound for the missed
detection probability

$\begin{array}[]{rcl}P_{\text{md}}&\geq&P_{\text{md}}^{\text{L}}\hfill\\\
&=&1-P_{s}^{\text{U}}\hfill\\\
&=&1-[\frac{N_{a}(1-\Lambda(n_{I}))}{N_{a}}+(N_{a}\Lambda(n_{I}))\times\\\
&&\frac{\leavevmode\resizebox{352.31624pt}{}{$\left(1-\dfrac{1}{{\sqrt{\pi}}}{\int}_{0}^{\infty}{{\int}_{\sqrt{y}}^{\infty}{\dfrac{{{y^{M/2-1}}{e^{-y/2}}}}{{{2^{M/2}}\Gamma(\frac{M}{2})}}}{e^{-{x^{2}}}}}dxdy\right)^{N_{\text{m}}}$}}{N_{a}}].\end{array}$
(36)

###  IV-C Uplink throughput analysis

We define the uplink throughput as follows

$\begin{array}[]{l}\gamma=\dfrac{{{S_{\text{SchemeName}}}{N_{\text{PD}}}R}}{{N_{\text{PD}}}+O_{\text{SchemeName}}}=\dfrac{{{N_{a}}{P_{s}}{N_{\text{PD}}}R}}{{N_{\text{PD}}}+O_{\text{SchemeName}}},\end{array}$
(37)

where ${S_{\text{SchemeName}}}$ denotes the number of successful UEs where
subscript “SchemeName” stands for the name of the scheme,
$O_{\text{SchemeName}}$ denotes the overhead of the random access scheme, and
${N_{\text{PD}}+O_{\text{SchemeName}}}$ is the frame length. Fig. 1 indicates
that the amount of total overhead of the proposed GCICA-RA scheme is
$O_{\text{GCICA-RA}}=\tau_{p}\times L+1$.

Substituting (34) into (37), we can derive the upper bound for the uplink
throughput of the proposed GCICA-RA scheme

$\gamma\leq\gamma^{\text{U}}=\dfrac{{{N_{a}P_{s}^{\text{U}}}{N_{\text{PD}}}R}}{{N_{\text{PD}}}+\tau_{p}\times
L+1}.$ (38)

##  V Numerical Results

In this section, we first present the performance of the proposed GCICA-RA
scheme, including the successful access probability, missed detection
probability, and MSE performance of the CSI estimation method. Then, we make a
comparison between the proposed GCICA-RA scheme, multiple-preamble grant-free
RA scheme proposed in [20] and the traditional random access scheme in terms
of the uplink throughput. For the traditional random access scheme, active UEs
sends their randomly selected pilots along with their uplink messages to the
BS, and the BS can only decode the uplink message of UE free from pilot
collision.

In the simulation, we consider the urban micro scenario, where the path loss
exponent is 3.8 [25]. The radius of the cell is 200 meters and all UEs are
uniformly distributed among the location farther than 25 meters from the BS.
Simulation parameters are shown in Table II. In addition, we utilize a LDPC
code with 1/2 rate to code the payload data, and then employ the BPSK scheme
to module the coded payload data. In the simulation, since the length of
payload data is a few hundred bytes in M2M communications [26], the payload
data size is set to 128 bytes to meet the communication requirements of M2M
communications, and further to verify the proposed GCICA-RA scheme.

TABLE II: Simulation parameters Parameter | Value
---|---
$R$ (Code rate) | 1/2
$N_{\text{PD}}$ | (Symbol length of the coded
---
payload data)
2048
$M$ (The number of antennas) | 400
$\text{SNR}_{e}$ | (Median signal to noise ratio (SNR)
---
of the UEs at the corner of cell)
5, 10, 15, 20 dB
$\tau_{p}$ (The number of pilots during each sub-pilot phase) | 6,9,10
$L$ (The number of sub-pilot phases) | 2,3
$N_{\rm{I}}$ (The number of ICA classifiers) |  5$\sim$30
$N_{a}$ | (The number of active UEs during
---
one random access procedure)
10 $\sim$ 40

Figure 5:  Successful access probability versus the number of ICA classifiers
$N_{\text{I}}$.
Figure 6: Missed detection probability versus the number of ICA classifiers
$N_{\text{I}}$.
Figure 7: MSE performance versus the number of ICA classifiers $N_{\text{I}}$.
Figure 8: Uplink throughput comparison between GCICA-RA and multiple-preamble
grant-free RA scheme.

Fig. 5 shows how the successful access probability $P_{s}$ changes with the
number of ICA classifiers $N_{\text{I}}$, where $\text{SNR}_{e}=5,10,15,20$
dB, $M=400$, $L=2$, $N_{a}=20$ and $\tau_{p}=10$. The upper bound for the
successful access probability can be obtained according to (34). We can see
from Fig. 5 that, with the increase of $N_{\text{I}}$, the successful access
probability increases dramatically first and then increases at a slower pace
for different $\text{SNR}_{e}$. The reason is that, the more the number of ICA
classifiers, the more signals of one UE we can obtain, which further improves
the decoding performance. We can also observe that, with the increase of
$\text{SNR}_{e}$, the successful access probability increases, and is close to
the upper bound.

Fig. 6 shows how the missed detection probability changes with the number of
ICA classifiers $N_{\text{I}}$, where $\text{SNR}_{e}=5,10,15,20$ dB, $M=400$,
$L=2$, $N_{a}=20$ and $\tau_{p}=10$. The lower bound for the missed detection
probability can be obtained according to (36). Fig. 6 illustrates that, with
the increase of $N_{\text{I}}$, the missed detection probability decreases
dramatically first and then decreases at a slower pace for different
$\text{SNR}_{e}$. We can also note that, with the increase of
$\text{SNR}_{e}$, the missed detection probability takes small values and
decreases, and is close to the lower bound.

Fig. 7 illustrates MSE performance of the proposed GCICA-RA, where
$\text{SNR}_{e}=5,10,15,20$ dB, $M=400$, $L=2$, $N_{a}=20$ and $\tau_{p}=10$.
We can see from Fig. 7 that, with the increase of $N_{\text{I}}$, MSE
decreases dramatically first and then decreases at a slower pace for different
$\text{SNR}_{e}$. We can also note that MSE takes small values and decreases
with the increase of $\text{SNR}_{e}$.

Fig. 8 compares the uplink throughput between the proposed GCICA-RA, multiple-
preamble grant-free RA schemes in [20], and the traditional random access
scheme with different $L$. We set $M=400$, $\text{SNR}_{e}=10$ dB and
$N_{\text{I}}=30$. To see how the uplink throughput changes with the number of
sub-pilot phases, we fix the length of the super pilot to 18, i.e.,
$L\times\tau_{p}=18$, and set $L$ to 2 and 3. To compare fairly, we ensure all
random access scheme have the same pilot length and the same symbol length of
the modulated uplink message. The upper bound for the uplink throughput of the
proposed GCICA-RA is obtained via (38). Note that, since upper bounds for the
uplink throughput of the proposed GCICA-RA are same for cases of $L=2$ and
$L=3$ based on (38), we not distinguish these two upper bounds in this figure.
It is observed that, the uplink throughput of $L=3$ is slightly higher than
that of $L=2$ for the proposed GCICA-RA scheme. The reason is that the number
of iterations of the SIC for $L=2$ is larger than that for $L=3$ on average,
which results in the higher CSI estimation error [14] and further increase the
equivalent noise introduced by the SIC algorithm. We also see that, the uplink
throughput of the proposed GCICA-RA scheme is close to the upper bound, and
significantly higher than those of baselines. With the increase of the number
of active UEs, the gap between the uplink throughput of the proposed GCICA-RA
scheme and those of baselines increases. This means that our proposed GCICA-RA
is more suitable to crowded scenarios.

##  VI Conclusion

An GCICA-RA grant-free scheme was proposed to support massive connectivity for
M2M communications in massive MIMO systems. The GCICA-RA scheme allows each
active UE to transmit their super pilots and modulated uplink messages to the
BS. Based on the received pilot signal and modulated uplink message, the BS
utilizes a proposed GCICA algorithm to estimate the CSI and decode the uplink
message of each active UE jointly, and then employs the estimated CSIs to
detect active UEs. Simulation results demonstrated that, our proposed GCICA-RA
scheme achieves performance gain even in the very crowded scenarios.

## References

  * [1] J. A. Stankovic, “Research directions for the internet of things,” _IEEE Internet of Things Journal_ , vol. 1, no. 1, pp. 3–9, Feb. 2014.
  * [2] P. Schulz, M. Matthe, H. Klessig, M. Simsek, G. Fettweis, J. Ansari, S. A. Ashraf, B. Almeroth, J. Voigt, I. Riedel, A. Puschmann, A. Mitschele-Thiel, M. Muller, T. Elste, and M. Windisch, “Latency critical IoT applications in 5G: Perspective on the design of radio interface and network architecture,” _IEEE Communications Magazine_ , vol. 55, no. 2, pp. 70–78, Feb. 2017.
  * [3] K. Ashton, “That ’internet of things’ thing,” _RFiD Journal_ , vol. 22, pp. 97–114, Jan. 2009.
  * [4] G. Wu, S. Talwar, K. Johnsson, N. Himayat, and K. Johnson, “M2m: From mobile to embedded internet,” _Communications Magazine, IEEE_ , vol. 49, pp. 36 – 43, May. 2011.
  * [5] C. VNI, “Cisco visual networking index: Forecast and trends, 2017Â¨c2022,” _White Paper_ , 2018.
  * [6] H. Q. Ngo, E. G. Larsson, and T. L. Marzetta, “Energy and spectral efficiency of very large multiuser MIMO systems,” _IEEE Transactions on Communications_ , vol. 61, no. 4, pp. 1436–1449, Apr. 2013.
  * [7] C. Xu, Y. Hu, C. Liang, J. Ma, and L. Ping, “Massive MIMO, non-orthogonal multiple access and interleave division multiple access,” _IEEE Access_ , vol. 5, pp. 14 728–14 748, 2017.
  * [8] C. Anton and M. Dohler, _Machine-to-machine (M2M) Communications Architecture, Performance and Applications_. Elsevier, Jan. 2015.
  * [9] J. Yuan, H. Shan, A. Huang, T. Q. S. Quek, and Y. Yao, “Massive machine-to-machine communications in cellular network: Distributed queueing random access meets MIMO,” _IEEE Access_ , vol. 5, pp. 2981–2993, 2017\.
  * [10] L. Sanguinetti, A. A. D’Amico, M. Morelli, and M. Debbah, “Random access in uplink massive mimo systems: How to exploit asynchronicity and excess antennas,” in _2016 IEEE Global Communications Conference (GLOBECOM)_ , Dec. 2016, pp. 1–5.
  * [11] Z. Dawy, W. Saad, A. Ghosh, J. G. Andrews, and E. Yaacoub, “Toward massive machine type cellular communications,” _IEEE Wireless Communications_ , vol. 24, no. 1, pp. 120–128, Feb. 2017.
  * [12] J. H. Sørensen, E. de Carvalho, and P. Popovski, “Massive MIMO for crowd scenarios: A solution based on random access,” in _IEEE Globecom Workshops_ , Dec. 2014, pp. 352–357.
  * [13] E. Björnson, E. de Carvalho, J. H. Sørensen, E. G. Larsson, and P. Popovski, “A random access protocol for pilot allocation in crowded massive MIMO systems,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 4, pp. 2220–2234, Apr. 2017.
  * [14] H. Han, Y. Li, and X. Guo, “A graph-based random access protocol for crowded massive MIMO systems,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 11, pp. 7348–7361, Nov. 2017.
  * [15] H. Han, X. Guo, and Y. Li, “A high throughput pilot allocation for M2M communication in crowded massive MIMO systems,” _IEEE Transactions on Vehicular Technology_ , vol. 66, no. 10, pp. 9572–9576, Oct. 2017.
  * [16] H. Han, Y. Li, and X. Guo, “User identity-aided pilot access scheme for massive MIMO-IDMA system,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 6, pp. 6197–6201, Jun. 2019.
  * [17] J. Ahn, B. Shim, and K. B. Lee, “Ep-based joint active user detection and channel estimation for massive machine-type communications,” _IEEE Transactions on Communications_ , vol. 67, no. 7, pp. 5178–5189, Jul. 2019.
  * [18] L. Liu and W. Yu, “Massive connectivity with massive MIMO-Part I: Device activity detection and channel estimation,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 11, pp. 2933–2946, Jun. 2018.
  * [19] H. Han, Y. Li, W. Zhai, and L. Qian, “A grant-free random access scheme for M2M communication in massive MIMO systems,” _IEEE Internet of Things Journal_ , vol. 7, no. 4, pp. 3602–3613, Apri. 2020.
  * [20] H. Jiang, D. Qu, J. Ding, and T. Jiang, “Multiple preambles for high success rate of grant-free random access with massive MIMO,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 10, pp. 4779–4789, Oct. 2019.
  * [21] L. Shen, Y. Yao, and H. Wang, “ICA based semi-blind decoding method for a multicell multiuser massive MIMO uplink system in Rician/Rayleigh fading channels,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 11, pp. 7501–7511, Nov. 2017.
  * [22] J. Shen, J. Zhang, and K. B. Letaief, “Downlink user capacity of massive MIMO under pilot contamination,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 6, pp. 3183–3193, Jun. 2015.
  * [23] G. Liva, “Graph-based analysis and optimization of contention resolution diversity slotted aloha,” _IEEE Transactions on Communications_ , vol. 59, no. 2, pp. 477–487, Feb. 2011.
  * [24] G. P. John, _Digital Communications_. Columbus, OH, USA: McGraw-Hill, 2001.
  * [25] M. Fallgren and B. Timus, “Spatial channel model for multiple input multiple output (MIMO) simulations (Release 13),” document TR 25.996, 3GPP, Tech. Rep., Dec. 2015.
  * [26] H. Lee, “The seven charicteristics of M2M traffic: Observations and implications,” _IEEE Transactions on Vehicular Technology_ , vol. 2, pp. 341–349, Dec. 2013.