* Van Aertryck and Jensen [2003] Lionel Van Aertryck and Thomas Jensen.  Uml-casting: Test synthesis from uml models using constraint resolution.  Proc. AFADL’2003 (Approches Formelles dans l’Assistance au Dévelopment de Logiciel). Springer, 2003.
  * Vargha and Delaney [2000] András Vargha and Harold D Delaney.  A critique and improvement of the cl common language effect size statistics of mcgraw and wong.  _Journal of Educational and Behavioral Statistics_ , 25(2):101–132, 2000.  doi: 10.3102/10769986025002101.
  * Warmer and Kleppe [2003] Jos B Warmer and Anneke G Kleppe.  _The object constraint language: getting your models ready for MDA_.  Addison-Wesley Professional, 2003.
  * Weißleder and Schlingloff [2007] Stephan Weißleder and Bernd-Holger Schlingloff.  Deriving input partitions from uml models for automatic test generation.  International Conference on Model Driven Engineering Languages and Systems, pages 151–163. Springer, 2007.
  * Weyuker et al. [1994] Elaine Weyuker, Tarak Goradia, and Ashutosh Singh.  Automatically generating test data from a boolean specification.  _IEEE Transactions on Software Engineering_ , 20(5):353–363, 1994.  doi: 10.1109/32.286420.
  * Williams et al. [2003] Ryan Williams, Carla P Gomes, and Bart Selman.  Backdoors to typical case complexity.  IJCAI, pages 1173–1178. Citeseer, 2003.
  * Woodward and Hennell [2006] Martin R Woodward and Michael A Hennell.  On the relationship between two control-flow coverage criteria: all jj-paths and mcdc.  _Information and Software Technology_ , 48(7):433–440, 2006.  doi: 10.1016/j.infsof.2005.05.003.
  * Wotawa et al. [2010] Franz Wotawa, Mihai Nica, and Bernhard K Aichernig.  Generating distinguishing tests using the minion constraint solver.  2010 Third International Conference on Software Testing, Verification, and Validation Workshops, pages 325–330. IEEE, 2010.
  * Wu [2016] Hao Wu.  Generating metamodel instances satisfying coverage criteria via smt solving.  2016 4th International Conference on Model-Driven Engineering and Software Development (MODELSWARD), pages 40–51. IEEE, 2016.
  * Wu [2022] Hao Wu.  A query-based approach for verifying uml class diagrams with ocl invariants.  18th European Conference on Modelling Foundations and Applications. Springer, 2022.
  * Wu [2023] Hao Wu.  Qmaxuse: A new tool for verifying uml class diagrams and ocl invariants.  _Science of Computer Programming_ , 228:102955, 2023.  doi: 10.1016/j.scico.2023.102955.
  * Wu and Timoney [2020] Hao Wu and Joseph Timoney.  Verifying ocl operational contracts via smt-based synthesising.  MODELSWARD, pages 249–259. Springer, 2020.
  * Yu and Lau [2006] Yuen Tak Yu and Man Fai Lau.  A comparison of mc/dc, mumcut and several other coverage criteria for logical decisions.  _Journal of Systems and Software_ , 79(5):577–590, 2006.  doi: 10.1016/j.jss.2005.05.030.
  * Zafar et al. [2021] Muhammad Nouman Zafar, Wasif Afzal, Eduard Enoiu, Athanasios Stratis, Aitor Arrieta, and Goiuria Sagardui.  Model-based testing in practice: An industrial case study using graphwalker.  14th Innovations in Software Engineering Conference (formerly known as India Software Engineering Conference), pages 1–11. ACM, 2021.


11institutetext: Centre for mathematical Plasma Astrophysics, Department of
Mathematics, KU Leuven, Celestijnenlaan 200B, B-3001, Leuven, Belgium
22institutetext: Centre for Computational Helio Studies, Ilia State
University, G. Tsereteli street 3, 0162 Tbilisi, Georgia
33institutetext: Evgeni Kharadze Georgian National Astrophysical Observatory,
M. Kostava street 47/57, 0179 Tbilisi, Georgia
44institutetext: Space Research Institute, Austrian Academy of Sciences,
Schmiedlstrasse 6, 8042 Graz, Austria
55institutetext: Combinatorial Optimization and Decision Support, KU Leuven
campus Kortrijk, E. Sabbelaan 53, 8500 Kortrijk, Belgium
66institutetext: Institute of Physics, IGAM, University of Graz,
Universitätsplatz 5, 8010 Graz, Austria
77institutetext: Ilia State University, Cholokashvili Ave 3/5, 0162 Tbilisi,
Georgia
88institutetext: Skobeltsyn Institute of Nuclear Physics, Moscow State
University, Moscow 119992, Russia
99institutetext: Institute of Physics, University of Maria Curie-Skłodowska,
Pl. M. Curie-Skłodowska 5, PL-20-031 Lublin, Poland
1010institutetext: Institute of Astronomy, Russian Academy of Sciences, Moscow
119017, Russia

# Eigenspectra of solar active region long-period oscillations

G. Dumbadze 112233 B.M. Shergelashvili 22334455 S. Poedts 1199 T.V.
Zaqarashvili 336677 M. Khodachenko 44881010 and P. De Causmaecker 55

###### Abstract

Context. We studied the low-frequency $\lesssim 0.5\;$h-1 (long-period
$\gtrsim 2\;$h) oscillations of active regions (ARs). The investigation is
based on an analysis of a time series built from Solar Dynamics
Observatory/Helioseismic and Magnetic Imager (SDO/HMI) photospheric
magnetograms and comprises case studies of several types of AR structures.

Aims. The main goals are to investigate whether ARs can be engaged in long-
period oscillations as unified oscillatory entities and, if so, to determine
the spectral pattern of such oscillations.

Methods. Time series of characteristic parameters of the ARs, such as, the
total area, total unsigned radial magnetic flux, and tilt angle, were measured
and recorded using the image moment method. The power spectra were built out
of Gaussian-apodised and zero-padded datasets.

Results. There are long-period oscillations ranging from 2 to 20 h, similarly
to the characteristic lifetimes of super-granulation, determined from the
datasets of the AR total area and radial magnetic flux, respectively. However,
no periodicity in tilt angle data was found.

Conclusions. Whatever nature these oscillations have, they must be
energetically supported by convective motions beneath the solar surface. The
possible interpretations can be related to different types of
magnetohydrodynamic (MHD) oscillations of the multi-scale structure of the AR
magnetic field, which is probably linked with the characteristic turnover
timescales of the super-granulation cells. The presence of oscillations in the
radial magnetic flux data may be connected to periodic flux emergence or
cancellation processes.

###### Key Words.:

Sun: magnetic fields; Active Regions; Methods: data analysis; Sun:
oscillations

##  1 Introduction

Active regions (ARs) are areas of strong magnetic field concentration visible
on the solar surface and in its lower atmosphere. These are areas of intense
and often complex magnetic activity in which sunspots frequently form. ARs are
often related to solar activity; for example, solar flares and coronal mass
ejections. The number, the locations, and the sizes of these ARs vary with
time, making ARs valuable indicators or tracers of the magnetic activity of
the Sun. The available high-resolution (both in time and space) observations
vividly show the complexity of the AR morphology and their intriguing
dynamics. One important aspect of the observable dynamical properties of ARs
and of the sunspots they contain concerns the huge variety of waves and
oscillations they exhibit. The various oscillatory phenomena can be classified
in different types such as umbral chromospheric oscillations with a typical
period of three minutes (Chorley et al., 2010; Centeno et al., 2006; Fleck &
Schmitz, 1991; Kuridze et al., 2009), umbral photospheric oscillations with a
typical period of five minutes (Thomas et al., 1984; Shergelashvili & Poedts,
2005), long-period oscillations with typical periods of several hours
(Dumbadze et al., 2017; Efremov et al., 2007; Goldvarg et al., 2005), and
ultra-long-period oscillations with typical periods of several days
(Khutsishvili et al., 1998; Gopasyuk, 2004). The long-period oscillations in
sunspots were discovered both in radio emission measurements and in magnetic
field data (Gelfreikh et al., 2006; Efremov et al., 2007; Solov’ev & Kirichek,
2008; Smirnova et al., 2013; Abramov-Maximov et al., 2013). The ground-based
and space-based observations yielded very similar long-period oscillations of
sunspots that confirm the solar nature of these oscillations (Abramov-Maximov
et al., 2013; Nagovitsyn & Nagovitsyna, 2011). State-of-the-art, space-based
observational facilities provide data with unprecedented accuracy, which
encourages both observational and theoretical studies of the physics behind AR
dynamics (Moradi, 2012).

The image moment method represents a tool for image processing to recognise
visual patterns of the image objects and to characterise their position, size,
and orientation (Flusser, 2007; Hu, 1962). The invariants are computed using
their boundary shape and interior area, and provide characteristics of an
object to represent its properties (Hu, 1962; Prokop & Reeves, 1992). Hu
(1962) used the mathematical foundation for two-dimensional moment invariants
to derive seven known invariants with respect to the rotation of 2D objects.
Furthermore, Hu’s moment invariant method has been used in template matching
and the registration of satellite images (e.g. Goshtasby, 1985; Flusser & Suk,
1994) for character recognition of objects (Belkasim et al., 1991; Tsirikolias
& Mertzios, 1993) and in many other applications.

In Dumbadze et al. (2017, hereafter Paper I), we reported and discussed long-
period oscillations in ARs, with periods of up to several hours. In that
paper, we also focused on oscillations of the tilt angle of the AR’s cross-
sectional area. The tilt angle oscillations were interpreted as the second
harmonic of standing kink modes, assuming the ARs are related to magnetic flux
tubes with feet embedded in the dense plasma of the solar interior and
connected to the lighter (less dense) part of the loop observed above the
surface. We assumed that the related kink modes have nodes in the apex above
the solar surface and that the ARs oscillate as a whole. This simplifying
assumption enabled us to roughly determine the possible distribution of their
phase speed values along the corresponding magnetic tubes. Our estimations
showed that the characteristic depth of the ARs is about $40\,000\;$km, and
that below this depth the sunspots are presumably fragmented into smaller
magnetic tube structures as was suggested by Parker (1979).

The analysis presented here is in some sense a continuation of the previous
phenomenological study of ARs reported in Paper I. The objective of the
present paper is to more systematically examine the existence of the
oscillations in different types of ARs. For this purpose, we performed a set
of case studies for typical ARs along their transit across the solar disk,
which were selected according to their morphological structure. The remainder
of this paper is organised as follows: the observations of the ARs and the
data processing methods are described in Sect. 2. Section 3 is dedicated to
the analysis of the significant periods discovered. The discussion and
conclusions are presented in Sect. 4.

##  2 Observations and data processing methodology

###  2.1 Observational data

Figure 1: Illustration of Helioseismic and Magnetic Imager (HMI) magnetogram
snapshots of the five studied ARs.

The aim of the study reported here is to reveal temporal oscillatory patterns
in the datasets related to a number of ARs. We carried out case studies for
the five ARs that are shown in Fig. 1. The list of these five studied ARs with
some of their properties is also given in Table 1. During the selection of
these specific ARs, we paid attention to their sizes, shapes, and locations on
the solar disk. The observational time span is chosen so that we observe the
transition of the ARs from roughly $-70^{\circ}$ to $+70^{\circ}$ longitude.

Table 1: The selected ARs.

1 warning AR | Obs. time |  Longitude (1) |  Latitude (1) |  Area (2) |  Threshold (3)
---|---|---|---|---|---
no. | (h) |  (∘) |  (∘) |  (km2) | (G)
12381 | $240$ |  [$-69.9\;$ $65.2$] |  [$14.7\;$ $14.5$] | $1.6\cdot 10^{9}$ | $400$
12378 | $240$ |  [$-64.5\;$ $60.2$] |  [$-15.7\;$ $-14.9$] | $1.5\cdot 10^{8}$ | $300$
12435 | $225$ |  [$-62.5\;$ $54.4$] |  [$-15.3\;$ $-15.4$] | $6.0\cdot 10^{7}$ | $350$
12437 | $225$ |  [$-57.9\;$ $59.2$] |  [$-21.6\;$ $-21.7$] | $6.0\cdot 10^{7}$ | $350$
12524 | $240$ |  [$-69.2\;$ $69.6$] |  [$14.9\;$ $14.3$] | $6.1\cdot 10^{8}$ | $400$

  * •

Notes. (1)The location of the AR on the solar surface at the beginning and the
end moment of the observation period. (2)The area values are taken at the
moment when the centre of mass of the AR is at or very close to the visible
central meridian. (3)The thresholds of unsigned magnetograms for the
corresponding unsigned binary images are used to isolate magnetically active
pixels (Hoeksema et al., 2014) for each AR.

We used data from the Solar Dynamics Observatory (SDO)/Helioseismic and
Magnetic Imager (HMI) (Schou et al., 2012; Scherrer et al., 2012). The line-
of-sight magnetogram ($B_{LOS}$) dataset of 45 s cadence was used for the
systematic investigation of active region long-period oscillatory patterns.
The projection effects were corrected using the ’mtrack’ algorithm in the
SDO/HMI ring-diagram pipeline (Bogart et al., 2011a, b), which employs the
azimuthal equidistant (Postel) projection. This involves a mapped rectangular
frame following the AR at selected heliographic coordinates and tracks it
during the time interval with the Carrington rotation rate. Also, we
considered all three components of the magnetic field (i.e. radial $B_{r}$,
meridional $B_{\theta,}$ and azimuthal $B_{\varphi}$), which were obtained
from Space-weather HMI Active Region Patches (SHARP, Bobra et al. (2014)),
available as sharp_cea_720s data with a cadence rate of 720 s. For these data,
the projection effects were corrected using cylindrical equal-area mapping.

###  2.2 Generation of time series

As mentioned before, we used the image moment method to process the obtained
series of rectangular frames described in the previous sub-section. The
standard procedure implies a transformation of the original $B_{LOS}$
magnetograms (or $B_{r}$, $B_{\theta}$ and $B_{\varphi}$) into unsigned binary
images of the ARs using certain thresholds for each of them (cf. last column
of Table 1). For each of the different ARs, the value of the threshold is
chosen individually from the $\pm\leavevmode\nobreak\ 300-400\;$G range (see
Table 1) so that the spurious small magnetic features around the AR are
filtered out, while the principal components of the AR body remain visible,
and more or less unchanged throughout the observational time span. These
thresholds are kept the same within the duration of observation for a given
AR. After the implementation of thresholding and binary imaging procedures,
the binary snapshots were used for the calculation of the image moments. The
moments of invariants are characterised with their two-dimensional $(p+q)$th
order defined by Hu (1962):

$M_{pq}=\iint{x^{p}y^{q}I(x,y)\;dx\,dy},$ (1)

where $p,q=0,1,2,...$; $x$ and $y$ are the coordinates of the pixel in the
respective scaling; $I(x,y)$ denotes the binary rate (set of active pixels) of
the unsigned magnetic field. For clarity with regard to the analysis, we
introduce the notation for total areas as $D_{i}$ ($i=$ $B_{LOS}$, $B_{r}$,
$B_{\theta}$ or $B_{\varphi}$), which are defined by the corresponding
unsigned magnetic field binary rate $I$.

For the analysis, we use the discretised counterpart of Eq. (1):

$M_{pq}=\sum_{\hbox{\tiny active pixels}}{{x^{p}y^{q}I(x,y)}},$ (2)

where $I(x,y)$ represents the mentioned discrete binary rate of $B_{LOS}$(or
$B_{r}$, $B_{\theta}$, $B_{\varphi}$). In this case, the zeroth-order moment
invariant $M_{00}$ determines the area of the object measured in the pixel
areas. We can also calculate the principal axes of the object (Hu, 1962) using
the first- and second-order moment invariants. Here, the principal axes of the
object correspond to the major and minor axes of the AR. The inclination of
the AR major axis with respect to the equatorial line of the solar disk, which
is identified as the tilt angle of the major axis (the largest principal axis
of the object), can be derived using the first- and second-order moments of
the image object:

$\theta=\frac{1}{2}\arctan\left(\frac{2\left(M_{11}/M_{00}-\bar{x}\cdot\bar{y}\right)}{\left(M_{20}-M_{02}\right)/M_{00}-\left(\bar{x}^{2}-\bar{y}^{2}\right)}\right).$
(3)

Using the sequence of snapshots within the observational time span, we thus
produced time series for the tilt angles and the total areas of the studied
ARs, and this for all components of the magnetic field. In addition, we
produced time series of the total unsigned radial magnetic flux, determined as

$|\Phi|=\sum_{\hbox{\tiny active pixels}}|B_{r}(x,y)|S(x,y),$ (4)

where $B_{r}(x,y)$ is the radial magnetic field component in each active pixel
identified by the corresponding binary image, and $S(x,y)=1.33\times
10^{5}\;$km2 is the area covered by these pixels on the solar surface
corresponding to the cylindrical equal area remapping (Bobra et al., 2014).

###  2.3 Data pre-processing

We focussed our study on long-period ($>2\;$h) (low-frequency) oscillation
signals with a satisfactory level of confidence. Further discussion reveals
that such periodicities are seen in the Fourier spectra of the considered time
series. However, they are ’drenched’ in the background power-law noise, making
it difficult to prove the significance of the detected periods. Therefore, in
order to elaborate the characteristic periods with a satisfactory level of
confidence and to distinguish the real quasi-oscillations in the time series
from artefacts, we used three different methods for the period detection:

####  2.3.1 Method 1 (M1) – Power spectra of the datasets

Vaughan (2005) suggested a method that finds significant peaks on top of the
power-law noise and calculates their uncertainties. We applied this method to
our datasets to find significant peaks in the power spectra, implying the
following steps. (i) De-trending and apodisating the initial data (as
described in Appendix A, items 1 and 2). (ii) Computing of the periodograms of
the detrended and apodised datasets; and along with this, we omit the power of
the Nyquist frequency following the recipe given in Vaughan (2005) and convert
the frequency ($f$) and power ($W$) into log-space (see example in Fig. 4).
(iii) Making a linear model fit ($\lg W_{model}=m\lg f+b$) of the obtained
power spectrum in log-scale (cf. solid blue lines in Fig. 4). (iv) Determining
the residuals of the computed power from the model line assuming that they
follow the second order two-dimensional chi-square ($\chi_{2}^{2}$)
distribution as a null hypothesis. In particular, we test the goodness of the
model with the Kolmogorov-Smirnov (KS) test (Press et al., 2007). When the KS
test does not satisfy the null hypothesis, we adjust the slope $m$ and the
offset $b$ until the test is fulfilled. (v) Estimating the 95% probability
limit $\gamma_{e}$ by assuming that the noise is two-dimensional $\chi^{2}$
distributed following Eq. (16) from Vaughan (2005) and using the computed $m$
and $b$ coefficients (we also plot this confidence limit into log-space; cf.
dashed red lines in Fig. 4). (vi) Selecting the significant peak of the lowest
frequency (longest period) in the spectrum and comparing it with similar
period of the sinusoidal fitting function $A_{1}\sin(A_{2}t+A_{3})$ applied to
the considered initial dataset (for further details, we refer the reader to
Appendix A). The corresponding peak is then removed from the power spectrum by
directly subtracting the found model sine function from the dataset. (vii)
Steps (ii)-(v) are iteratively repeated until all outstanding peaks are
removed from the spectrum, and a pure power-law noise remains.

####  2.3.2 Method 2 (M2) – Re-binned power spectra of the datasets

In order to enhance significant peaks in the power spectrum, we used the
spectral re-binning method (Appourchaux, 2003). Pugh et al. (2017) showed how
the Vaughan (2005) method can be updated for use with re-binned power spectra.
They used this method for solar and stellar flare datasets to determine quasi-
periodic pulsations. We used the same method and proceeded with the following
steps: (i) performing steps (i) and (ii) of M1; (ii) compute the sum of
consecutive $p$ frequency bins using the original power spectrum and then
divide it by $p$ (in our calculations we use $p=2$); (iii) repeat steps
(iii)-(iv) of M1 (however, for the description of the noise distribution,
instead of a two-dimensional $\chi_{2}^{2}$ distribution function, we used a
$2p$-dimensional distribution function $\chi_{2p}^{2}$ (Pugh et al., 2017);
then, the corresponding 95% confidence line is estimated and plotted
accordingly); (iv) repeat steps (vi)-(vii) of M1 (see example in Fig. 4).

####  2.3.3 Method 3 (M3) – Average power spectra for four consecutive equal
parts (time windows) of entire observational time

In order to provide additional proof of the presence of the discovered
significant periods in the signal, we used a third method. It is based on the
division of the total observational time span in four equivalent mutually non-
overlapping intervals of time, and we derived the power spectra for each of
these intervals. Considering the shorter datasets, we sacrifice the spectral
resolution (the uncertainties increase by a factor four compared to the first
two methods) and, in this case, the noise no longer follows a chi-squared
distribution. However, we are able to generate the ’averaged’ Fourier power
spectrum out of these four observational windows.

The analysis steps of M3 are as follows. (i) Divide the analysed time series
into four equivalent non-overlapping parts. (ii) Compute the power spectrum of
each individual part, as described in step (ii) of M1. (iii) Sum the four
power spectra obtained. The peaks revealed by M1 and M2 are ’real’ if they
still remain in the summed spectrum, while the spurious ones disappear (see
example in Fig. 4). Then, we apply the same procedures as those outlined in
steps (iii)-(vii) of M1 again.

Both methods M1 and M2 enable us to isolate the individual periods of
oscillations with a spectral accuracy of approximately the same magnitude. The
purpose of using these two different methods is to cross-check the validity of
the discovered significant peaks. Thus, we are able to compare the periods
obtained from these two main methods and evaluate their average values.
Strictly speaking, however, these two methods of measurement are not
completely independent of each other. Therefore, when we calculate the mean
values of the periods, including measurements of areas and fluxes related to
particular components of the magnetic field obtained by one of the two methods
and cross confirmed by them where possible, the aim is not to demonstrate the
statistical mean of the mutually independent ensemble of the measurements
(according to its classical definition). Instead, we use the average values
just as a method of approximation to evaluate the globally averaged values we
discovered. With regard to the errors, the simple averaging of uncertainties,
when we take the arithmetic mean of measurements, is not consistent. In
particular, it is not consistent in the case when the methods of measurements
are not completely independent of each other. Therefore, we choose the
uncertainty approximation method in order to be prudent with the assumption
that the uncertainty of the mean of two measurements cannot be larger than the
maximum error of each individual method; hence, whenever we use the maximum
uncertainties. Additionally, M3 enables the verification and confirmation
(with a reduced but still reasonable accuracy) of the characteristic periods
inferred by the main methods, M1 and M2, where possible. This yields the
maximal robustness of the analysis that can be achieved.

##  3 Analysis and comparison of spectra for different cases

###  3.1 Revealing of major periodicities

Figure 2: Characteristic mean oscillation periods (measured in hours)
calculated as a global average of those found in different time series
separately (for $D_{B_{LOS}}$, $D_{B_{r}}$, $D_{B_{\theta}}$,
$D_{B_{\varphi,}}$ and total radial magnetic flux). On the horizontal axes,
the M1, M2 and M3 labels, respectively, denote the global mean of the periods
obtained by M1 (coloured ’$*$’ with error bars), M2 (coloured ’$\circ$’ with
error bars), and M3 (coloured ’$\times$’ with error bars). Moreover, we also
indicate the total mean periods by taking the average of the mentioned global
means of M1 and M2 (i.e. the mean periods obtained by M3 are excluded from the
total mean calculation) via the solid-coloured horizontal lines (with the
uncertainties indicated by the dashed coloured lines). The colour scheme
follows the sequence of quasi-harmonics described in Sect. 3.2 and also
corresponds to the colouring used in the tables (see Appendix B). Figure 3:
Ratios of global mean periods ($P_{i}/P_{i+1}$, $i=1,2,\ldots$) of M1, M2 and
the total mean periods shown in the panels of Fig. 2. The formatting of the
data points and horizontal lines is the same as in Fig. 2. The colouring in
each panel coincides with the period standing in the denominator in each
ratio. On the vertical axes, the levels of the period ratios in the reference
spectrum are labelled and the corresponding horizontal solid black lines are
plotted in all panels.

We study the long-period ($>2$ h) oscillatory patterns in ARs with four
different typical morphological structures (e.g. following the classification
given in Jaeggli & Norton, 2016). The considered periods are obtained from the
Fourier spectra of the time series, which have certain (limited) observational
time spans and sampling rates. The characteristic parameters of these datasets
determine the period resolution and the overall accuracy of the obtained power
spectra. As a result, some significant periods in the spectra are embedded in
strong background noise showing a power-law nature (see Appendix C). We do not
attempt to provide a physical interpretation of the mechanism(s) that might be
responsible for the appearance of these oscillatory patterns. Instead, we give
a few plausible indications concerning the possible scenarios. The final
interpretation of the detected oscillations requires the construction of
rigorous mathematical models, and this is beyond the scope of the present,
purely observational, study. The results of the performed analysis for each
considered AR are summarised in Tables 2-5 (see Appendix B).

Most of the periodic features are revealed for $D_{B_{LOS}}$. However, the
analysis of the estimated $D_{B_{r}}$, $D_{B_{\theta}}$ and $D_{B_{\varphi}}$
and the total unsigned radial magnetic flux yield a smaller number of detected
significant periods. In most cases, the values of these periods are similar to
those in $D_{B_{LOS}}$.

Tables 2-5 show that each AR yields a different number of modes, and further
investigations with more ARs are needed to verify the dependence of the number
of oscillation modes on the morphology. Another finding is that long periods
($>10$ h) are detected less frequently than the shorter periods ($<10$ h) in
all analysed datasets. The apparent reason for this is the gradual reduction
of the spectral resolution with the growth of the periods. We know also that
the longest period parts of the spectra are contaminated by two artificial
instrumental peaks at $12$ and $24$ h (Liu et al., 2012), which might cause
spectral power leakage and mixture between real, long periods and artificial
spectral peaks. Presumably, these artificial signals have less impact on the
part of the spectrum below $10$ h. Altogether, we show that in all datasets
there are long-period oscillations with characteristic periods in the range
from $2$ to $20$ h (beyond $20$ h, the derived spectra are rather unreliable
because of a very large uncertainty), which are similar to the characteristic
lifetimes of super-granulation (Rincon & Rieutord, 2018). So, it is intuitive
to connect the observed oscillation periods with the characteristic turnover
times of super-granulation cells.

None of the applied analysis methods described in Sect. 2.3 gave us any
significant periodicities in the AR major axis tilt angle data. Therefore, we
assume that the power spectrum for this quantity does not contain significant
periods and it remains of a pure power-law type. At the same time, we would
like to note that this non-detection does not mean that tilt angle
oscillations do not exist, because the uncertainty tolerance of the image
moment method used here, does not allow for sufficient confidence of such
oscillations. The difference between this result and that discussed in Paper I
can be related to the methodology of the analysis. In particular, the method
discussed in Paper I, identifies the principal axis of the AR and its tilt on
the basis of the least-squares fitting of its elliptical shape. While with the
method of image moments used in the present paper, the tilt angle $\theta$
(Eq. 3) is defined in terms of zeroth-, first-, and second-order moments of
each image, which probably produces more noise in the tilt angle datasets.
Thus, the noise might overwhelm the oscillation amplitude itself. Of course, a
more precise analysis of this difference requires a thorough cross-comparison
of the two methods of measurements.

The radial unsigned magnetic flux shows some significant periodicities (see
Table 5) for all ARs, comparably to those of the area oscillations (see Tables
2-4). For example, the oscillations with periods of about $6$ and $4$ h are
the most frequently detected ones. They are confirmed, with some exceptions,
by the two main methods M1 and M2, as well as being additionally proven by M3
(with some discrepancy tolerated by the error of M3). These periods are found
in the radial unsigned magnetic flux oscillations in practically all ARs. The
frequent presence of the six- and four-hour periods both in area and radial
magnetic flux datasets may be connected with periodic flux emergence or
cancellation processes, which can be another possible component of the
physical scenario behind the observed dynamic picture.

###  3.2 Possible interpretation as a signature of quasi-standing oscillation
patterns

Our analysis enables the detection of a discrete spectrum of significant
periods. Tables 2-5 (in Appendix B) show that the particular and averaged
periods measured by various methods are grouped in the following ranges: (i)
’harmonic’ $P_{1}$ – $17.1\pm 0.71-17.4\pm 0.74$ h (marked in cyan), (ii)
’harmonic’ $P_{2}$ – $7.66\pm 0.14-9.54\pm 0.90$ h (marked in magenta), (iii)
’harmonic’ $P_{3}$ – $5.69\pm 0.08-6.37\pm 0.10$ h (marked in red), (iv)
’harmonic’ $P_{4}$ – $4.45\pm 0.19-4.88\pm 0.23$ h (marked in green), and (v)
’harmonic’ $P_{5}$ – $3.66\pm 0.03-3.88\pm 0.05$ h (marked in blue). In Fig.
2, we plot all these period groups together. In Fig. 3, we show the ratios of
the average periods from each group in the same sequence as in Fig. 2. For the
calculation of the uncertainties found in the period ratios, we use the
following relation:

$\vartriangle\frac{P_{i}}{P_{i+1}}=\left|\frac{1}{P_{i+1}}\vartriangle
P_{i}-\frac{P_{i}}{P_{i+1}^{2}}\vartriangle
P_{i+1}\right|\leq\left|\frac{\vartriangle
P_{i}}{P_{i+1}}\right|+\left|-\frac{P_{i}}{P_{i+1}^{2}}\vartriangle
P_{i+1}\right|,$

where $P_{i}>P_{i+1}$ ($i=1,2,\ldots$ is the number of ’harmonics’) are the
global and total mean periods (see definition of these terms in the caption of
Fig. 2) evaluated for each AR. They sequentially follow each other in
descending order. One can see that the ratios of the averaged periods resemble
the sequence of oscillation harmonics (Fig. 3) typical for a standing
oscillations. The latter are shown with black solid horizontal lines as the
reference values of the ratios. We note that the observed ratios do not
exactly coincide with the sequence of ratios of periods corresponding to pure
standing oscillations ($2/1,3/2,4/3,\ldots$). Nevertheless, the observed
ratios with the entire uncertainty interval are well separated and
approximately follow the reference spectrum ratios. Hence, with the present
level of accuracy we can say that the observed discrete spectrum follows the
sequence of quasi-standing oscillations still showing some shift of the
periods and corresponding ratios from the reference values. This shift of
periods might be due to the incomplete line-tying of the magnetic loop system,
Doppler shifts due to the internal flows in ARs, and/or the inhomogeneity of
the magnetic loop system constituting the AR below and above the solar
surface. The presence of the structures in the emerged magnetic fields with
power-law-distributed spatial scales is well known and actively used in sand-
pile models of flares (see Baiesi et al. (2008) and references therein). The
formation of such transfer structures can be ascribed to non-equilibrium
statistical process of driven diffusion (Maes et al., 2009) of the emerged
elements of the magnetic flux.

##  4 Discussion and conclusions

One can link the discovered periods with the characteristic turnover timescale
of the super-granulation cells, which is known to be of the same order of
magnitude. Another possible interpretation can be related to the different
types of MHD oscillations of the multi-scale structure of the AR magnetic
field. It is well known that the MHD oscillations and waves sustained in
structures with non-uniform magnetic field can be engaged in the processes of
resonant (Terradas et al., 2008) or swing (Shergelashvili et al., 2005)
absorption. It is evident that, whatever nature these oscillations have, they
must be energetically supported by convective motions beneath the solar
surface, which represents the source of entropy production in the system
through macroscopic stochastic super-granular motions of the plasma. The
entropy variations can be due to microscopic transport processes as well
(Shergelashvili et al., 2007). The presence of oscillations in the radial
magnetic flux data may be connected with a periodic flux
emergence/cancellation process. One could interpret some parts of the
discovered spectra in terms of sequences of quasi-standing eigen-oscillations
with a full or partial line-tying configuration of the magnetic structures.
The MHD oscillations can be also driven by convective turbulence. To confirm
which of these scenarios prevails, further dedicated studies are needed.

In the very simple set-up we considered in Paper I, we interpreted the
oscillations of the major and minor axis lengths of the ellipses we used to
frame the ARs as higher (than sausage mode) degree flute modes. However, this
interpretation was rather intuitive and speculative. The goal of the current
work is to reveal the basic spectral characteristic of the oscillations that
will be used in the next work for a more rigorous modelling. The comparison of
the theoretical model with the observed spectra discussed here will enable us
to draw definite conclusions about the nature of these oscillations. According
to our preliminary results in Paper I, one can intuitively expect that the
large period oscillations should be related to the modes with phase speeds
comparable to the Alfvén velocity in the high-beta plasma in sunspots. Another
general conclusion can be that, whatever nature these modes have, if for some
reason there is a line-tying at some depth enabling the presence of nodes for
velocity or pressure perturbations, then such discrete spectra with fixed
sequences of period ratios can be formed. Since we plan further theoretical
investigations of the issue, we avoided any speculations about the physical
nature of the observed periodicities in the present paper, with an ultimate
goal of simply addressing the observational facts; that is to say, the
oscillation spectra observed in the range of $2-20$ h periods. Concerning the
pure standing oscillations, these would occur when the ratios of the periods
of all overtones strictly following the sequence $2/1,3/2,4/3,\ldots$. The
reason why the ratios differ form the reference spectrum requires
investigation. Most likely, the inhomogeneity along the flux tubes causes this
effect, as it does in coronal loops (Andries et al., 2009). Nevertheless, our
results indicate that the ratios of the detected periods roughly follow the
reference spectrum (shown on the ordinate axis of Fig. 3); however, they still
deviate from it somewhat. This might be caused by the uncertainty tolerance we
admit in our measurements and the foot point leakage of modes (not complete
line-tying) or their combination.

###### Acknowledgements.

The work was supported by Shota Rustaveli National Science Foundation grants
DI-2016-52, FR17_609. Work of G.D. was supported under Shota Rustaveli
National Science Foundation grants for doctoral students – PhDF2016_177 and
grant for young scientists for scientific research internships abroad
IG/50/1/16. Work of B.M.S. was supported by the Austrian Fonds zur Forderung
der wissenschaftlichen Forschung (FWF) under the project P25640-N27. S.P. was
supported by the projects C14/19/089 (C1 project Internal Funds KU Leuven),
G.0D07.19N (FWO-Vlaanderen), SIDC Data Exploitation (ESA Prodex-12). The work
of TVZ was supported by the Austrian Science Fund (FWF, project P30695-N27).
M.L.K. additionally acknowledges the projects I2939-N27 and S11606-N16 of the
Austrian Science Fund (FWF), grant No.075-15-2020-780 (GA No.13.1902.21.0039)
of the Russian Ministry of Education and Science and the project ”Study of
stars with exoplanets”, supported by grant No.075-15-2019-1875 from the
government of Russian Federation. We are thankful to anonymous referee for
constructive remarks on our manuscript that led to the significant improvement
of content.

## References

  * Abramov-Maximov et al. (2013) Abramov-Maximov, V. E., Efremov, V. I., Parfinenko, L. D., Solov’ev, A. A., & Shibasaki, K. 2013, PASJ, 65, S12
  * Andries et al. (2009) Andries, J., van Doorsselaere, T., Roberts, B., et al. 2009, Space Sci. Rev., 149, 3
  * Appourchaux (2003) Appourchaux, T. 2003, A&A, 412, 903
  * Baiesi et al. (2008) Baiesi, M., Maes, C., & Shergelashvili, B. M. 2008, Physica A Statistical Mechanics and its Applications, 387, 167
  * Belkasim et al. (1991) Belkasim, S., Shridhar, M., & Ahmadi, M. 1991, Pattern Recognition, 24, 1117
  * Bobra et al. (2014) Bobra, M. G., Sun, X., Hoeksema, J. T., et al. 2014, Sol. Phys., 289, 3549
  * Bogart et al. (2011a) Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & Rabello-Soares, M. C. 2011a, Journal of Physics: Conference Series, 271, 012008
  * Bogart et al. (2011b) Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & Rabello-Soares, M. C. 2011b, Journal of Physics: Conference Series, 271, 012009
  * Centeno et al. (2006) Centeno, R., Collados, M., & Trujillo Bueno, J. 2006, ApJ, 640, 1153
  * Chorley et al. (2010) Chorley, N., Hnat, B., Nakariakov, V. M., Inglis, A. R., & Bakunina, I. A. 2010, A&A, 513, A27
  * Dumbadze et al. (2017) Dumbadze, G., Shergelashvili, B. M., Kukhianidze, V., et al. 2017, A&A, 597, A93
  * Efremov et al. (2007) Efremov, V. I., Parfinenko, L. D., & Solov’ev, A. A. 2007, Astronomy Reports, 51, 401
  * Fleck & Schmitz (1991) Fleck, B. & Schmitz, F. 1991, A&A, 250, 235
  * Flusser (2007) Flusser, J. 2007, International Journal of Computer, Electrical, Automation, Control and Information Engineering, 1, 3708
  * Flusser & Suk (1994) Flusser, J. & Suk, T. 1994, IEEE Transactions on Geoscience and Remote Sensing, 32, 382
  * Gelfreikh et al. (2006) Gelfreikh, G. B., Nagovitsyn, Y. A., & Nagovitsyna, E. Y. 2006, PASJ, 58, 29
  * Goldvarg et al. (2005) Goldvarg, T. B., Nagovitsyn, Y. A., & Solov’Ev, A. A. 2005, Astronomy Letters, 31, 414
  * Gopasyuk (2004) Gopasyuk, O. S. 2004, in IAU Symposium, Vol. 223, Multi-Wavelength Investigations of Solar Activity, ed. A. V. Stepanov, E. E. Benevolenskaya, & A. G. Kosovichev, 249–250
  * Goshtasby (1985) Goshtasby, A. 1985, IEEE Transactions on Geoscience and Remote Sensing, 7, 338
  * Hoeksema et al. (2014) Hoeksema, J. T., Liu, Y., Hayashi, K., et al. 2014, Sol. Phys., 289, 3483
  * Hu (1962) Hu, M.-K. 1962, IRE Transactions on Information Theory, 8
  * Jaeggli & Norton (2016) Jaeggli, S. A. & Norton, A. A. 2016, ApJ, 820, L11
  * Khutsishvili et al. (1998) Khutsishvili, E., Kvernadze, T., & Sikharulidze, M. 1998, Sol. Phys., 178, 271
  * Kuridze et al. (2009) Kuridze, D., Zaqarashvili, T. V., Shergelashvili, B. M., & Poedts, S. 2009, A&A, 505, 763
  * Liu et al. (2012) Liu, Y., Hoeksema, J. T., Scherrer, P. H., et al. 2012, Sol. Phys., 279, 295
  * Maes et al. (2009) Maes, C., Netočný, K., & Shergelashvili, B. M. 2009, Phys. Rev. E, 80, 011121
  * Moradi (2012) Moradi, H. 2012, Astronomische Nachrichten, 333, 1003
  * Nagovitsyn & Nagovitsyna (2011) Nagovitsyn, Y. A. & Nagovitsyna, E. Y. 2011, Geomagnetism and Aeronomy, 51, 1049
  * Parker (1979) Parker, E. N. 1979, ApJ, 234, 333
  * Press et al. (2007) Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing (Cambridge University Press)
  * Prokop & Reeves (1992) Prokop, R. J. & Reeves, A. P. 1992, CVGIP: Graph. Models Image Process., 54, 438
  * Pugh et al. (2017) Pugh, C. E., Broomhall, A. M., & Nakariakov, V. M. 2017, A&A, 602, A47
  * Rincon & Rieutord (2018) Rincon, F. & Rieutord, M. 2018, Living Reviews in Solar Physics, 15, 6
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229
  * Shergelashvili et al. (2007) Shergelashvili, B. M., Maes, C., Poedts, S., & Zaqarashvili, T. V. 2007, Phys. Rev. E, 76, 046404
  * Shergelashvili & Poedts (2005) Shergelashvili, B. M. & Poedts, S. 2005, A&A, 438, 1083
  * Shergelashvili et al. (2005) Shergelashvili, B. M., Zaqarashvili, T. V., Poedts, S., & Roberts, B. 2005, A&A, 429, 767
  * Smirnova et al. (2013) Smirnova, V., Riehokainen, A., Solov’ev, A., et al. 2013, A&A, 552, A23
  * Solov’ev & Kirichek (2008) Solov’ev, A. A. & Kirichek, E. A. 2008, Astrophysical Bulletin, 63, 169
  * Sun (2013) Sun, X. 2013, arXiv e-prints, arXiv:1309.2392
  * Terradas et al. (2008) Terradas, J., Arregui, I., Oliver, R., et al. 2008, ApJ, 679, 1611
  * Thomas et al. (1984) Thomas, J. H., Cram, L. E., & Nye, A. H. 1984, ApJ, 285, 368
  * Tsirikolias & Mertzios (1993) Tsirikolias, K. & Mertzios, B. 1993, Pattern Recognition, 26, 877
  * Vaughan (2005) Vaughan, S. 2005, A&A, 431, 391

##  Appendix A Data processing details

For the application of original, re-binned and windowed power spectra to the
analysed datasets, we use data de-trending, apodisation, significant peak
identification, and removal by sine function fitting, which are specified
below.

###  A.1 Data de-trending and preparation for different spectral methods

Data processing begins with our obtaining time series using the method
described in Sect. 2.2 of image moments dealing with the selected values of
magnetic field thresholds for individual ARs. However, apart from the
stochastic and quasi-periodic parts, the initial data include aperiodic trends
with a characteristic timescale of variation equal to the total observation
time. We remove these aperiodic variations by subtracting a third-order
polynomial function fitted to the original data. Furthermore, as mentioned in
Sect. 2.1, for 45-second cadence data of $B_{LOS}$ we use Postel’s equidistant
projection mapping. Additionally, we used 720-second cadence data from the
Joint Science Operations Center (JSOC, Sun (2013); Hoeksema et al. (2014)).
The latter data are de-projected by using cylindrical equal area mapping for
radial $B_{r}$, meridional $B_{\theta}$ and azimuthal $B_{\varphi}$ components
of the magnetic field. Therefore, the projection effect related to the
spherical shape of the Sun is removed from the magnetogram image data.

###  A.2 Data apodisation

We perform a Gaussian apodisation of the de-trended data to reduce the effect
of finiteness of the total sampling time. The procedure implies multiplication
of the original datasets by a Gaussian function of the following form:

$w(n)=e^{-\frac{1}{2}\left(\alpha\frac{2n}{(N-1)}\right)^{2}},$ (5)

where $N$ is total number of sample points, while the integer $n$ is varying
in the range $-(N-1)/2\leq n\leq(N-1)/2,$ and the parameter $\alpha$ controls
the shape of the Gaussian. For all area datasets of ARs 12381, 12435, 12437,
and 12524, for $D_{B_{LOS}}$ of AR 12378, and for the radial flux datasets of
ARs 12381, 12435, and 12524, we use $\alpha=3$. For the remaining datasets, we
use $\alpha=2$. These $\alpha$ values were chosen for each AR separately by
gradually increasing it until we obtain the value for which any further
increase does not show any further significant affect on the obtained
spectrum.

###  A.3 Significant peak identification

To identify the significant peaks using this method, we take the following
steps:

  1. 1.

Perform a regression of the dataset on the sum of sine functions. We choose
the number of terms in this sum depending on the best regression goodness and
detect the one with the longest period (with its amplitude, angular frequency,
and phase).

  2. 2.

Test the relevance of the periods found, analysing an auto-correlation
function of the de-trended datasets, by taking its power spectra and making
sure that the period deduced directly in step (1) is indeed present in the
spectrum. Then this term is subtracted from the original dataset which almost
completely removes the longest oscillation period from the signal.

  3. 3.

Take the remaining dataset, after the subtraction of the largest periodic
signal found, as a new dataset and apply steps (1) and (2) to identify and
extract the next longest period. These iterations are repeated until the
shortest period of interest is detected and extracted.

##  Appendix B Tables

Here, we present the tables that summarise the significant periods obtained
from the total area variation in time corresponding to AR 12378, AR 12381, AR
12435, AR 12437, and AR 12524 (see Tables 2-4). The first columns of these
tables show the periods obtained by the power spectra of the datasets (M1).
The second columns indicate the periods detected by the method using the re-
binned power spectra of the datasets (M2). The third columns indicate the
averaged periods calculated based on the results of M1 and M2, and the last
columns denote the outcome of average power spectra for four consecutive equal
parts of entire observational time (M3). All periods and uncertainties shown
in the tables are expressed in hours. The colours correspond to the colour-
coding of the period groups described in Sect. 3.2.

In addition, we provide Table 5, which contains similar periodicities inferred
from the total radial magnetic flux of all considered ARs. The average
characteristic parameters of the oscillations, derived from the content of
these tables, are visualised in Fig. 2, while the ratios of the derived global
and total mean periods are plotted in Fig. 3. When multiple peaks occur in one
particular method only, while a single isolated one is detected in other
methods, we calculate the weighted mean of the frequency and define them as
’multiplets’. Accordingly, we assume that such multiplets do not carry real
physical information, but they simply represent spurious artefacts related to
analysis accuracy.

Table 2: Detected periods (with corresponding errors) for the area of the ARs
12378 and 12381.

AR 12378 | AR 12381
---|---
M1 | M2 | Averaged periods | M3 | M1 | M2 | Averaged periods | M3
$D_{B_{LOS}}$ | $D_{B_{LOS}}$
– | – | – | – | $17.1\pm 0.71$ | $17.4\pm 0.74$ | $17.3\pm 0.74$ | –
– |  $8.92\pm 0.20$* | $8.92\pm 0.20$ |  $9.54\pm 0.90$* | – |  $7.91\pm 0.15$* | $7.91\pm 0.15$ | –
$5.94\pm 0.09$ | $5.89\pm 0.08$ | $5.92\pm 0.09$ | $5.69\pm 0.32$ |  $6.37\pm 0.10$* |  $6.14\pm 0.09$* | $6.26\pm 0.10$ | –
– | $4.68\pm 0.05$ | $4.68\pm 0.05$ | – | $4.60\pm 0.05$ | $4.58\pm 0.05$ | $4.59\pm 0.05$ | $4.45\pm 0.19$
$3.79\pm 0.04$ | $3.81\pm 0.04$ | $3.80\pm 0.04$ | – | $3.69\pm 0.03$ | $3.67\pm 0.03$ | $3.68\pm 0.03$ | $3.79\pm 0.14$
$3.28\pm 0.03$ | – | $3.28\pm 0.03$ | $3.10\pm 0.10$ | – | – | – | –
– | $2.82\pm 0.04$ | $2.82\pm 0.04$ | $2.84\pm 0.08$ | – | – | – | –
$D_{B_{r}}$ | $D_{B_{r}}$
– | $7.66\pm 0.14$ | $7.66\pm 0.14$ | – | – | – | – | –
$6.11\pm 0.09$ | $5.89\pm 0.08$ | $6.00\pm 0.09$ | $6.02\pm 0.36$ | – | – | – | –
$4.65\pm 0.05$ | $4.68\pm 0.06$ | $4.67\pm 0.06$ | $4.65\pm 0.21$ | – | – | – | –
$D_{B_{\theta}}$ | $D_{B_{\theta}}$
– |  $8.14\pm 0.17$* | $8.14\pm 0.17$ | – | – | – | – | –
$5.77\pm 0.08$ |  $5.95\pm 0.09$* | $5.86\pm 0.09$ | $6.02\pm 0.36$ | $5.94\pm 0.08$ | $5.89\pm 0.08$ | $5.92\pm 0.08$ | $5.69\pm 0.63$
$4.72\pm 0.05$ |  $4.78\pm 0.06$* | $4.75\pm 0.06$ | $4.65\pm 0.21$ | – | – | – | –
– | – | – | – | – | $3.88\pm 0.05$ | $3.88\pm 0.05$ | –
$D_{B_{\varphi}}$ | $D_{B_{\varphi}}$
– | $7.66\pm 0.14$ | $7.66\pm 0.14$ | – | – | – | – | –
$6.02\pm 0.09$ |  $5.95\pm 0.09$* | $5.99\pm 0.09$ | $6.02\pm 0.36$ | – | – | – | –
– | – | – | – | – | $4.68\pm 0.03$ | $4.68\pm 0.03$ | –

All values marked by asterisks indicate the weighted mean periods calculated
from multiplets.

Table 3: Same as in Table 2, but for AR 12435 and AR 12437.

AR 12435 | AR 12437
---|---
M1 | M2 | Averaged periods | M3 | M1 | M2 | Averaged periods | M3
$D_{B_{LOS}}$ | $D_{B_{LOS}}$
– | $16.1\pm 0.63$ | $16.1\pm 0.63$ | – | $17.1\pm 0.71$ | $17.4\pm 0.74$ | $17.3\pm 0.73$ | –
$8.90\pm 0.19$ | $9.00\pm 0.20$ | $8.95\pm 0.20$ | $8.53\pm 0.71$ | – | $7.95\pm 0.15$ | $7.95\pm 0.15$ | $7.88\pm 0.61$
$6.40\pm 0.10$ | $6.45\pm 0.10$ | $6.43\pm 0.10$ | $6.02\pm 0.35$ | – | – | – | –
$D_{B_{r}}$ | $D_{B_{r}}$
– | $7.95\pm 0.15$ | $7.95\pm 0.15$ | – | – | – | – | –
– | $5.89\pm 0.08$ | $5.89\pm 0.08$ | $6.02\pm 0.35$ | $5.85\pm 0.08$ | $5.89\pm 0.08$ | $5.87\pm 0.08$ | $6.02\pm 0.36$
$4.71\pm 0.06$ |  $4.77\pm 0.06$* | $4.74\pm 0.06$ | $4.88\pm 0.23$ | – | – | – | $4.65\pm 0.21$
$D_{B_{\theta}}$ | $D_{B_{\theta}}$
– |  $7.79\pm 0.15$* | $7.79\pm 0.15$ | $7.88\pm 0.61$ | – | $7.95\pm 0.15$ | $7.95\pm 0.15$ | $7.88\pm 0.61$
$5.95\pm 0.09$* |  $6.01\pm 0.09$* | $5.98\pm 0.09$ |  $6.13\pm 0.37$* | $6.02\pm 0.09$ | $6.07\pm 0.09$ | $6.04\pm 0.09$ | $6.02\pm 0.36$
$5.00\pm 0.06$ | $5.03\pm 0.06$ | $5.02\pm 0.06$ | – | – | – | – | –
$4.75\pm 0.05$* |  $4.77\pm 0.06$* | $4.76\pm 0.06$ | $4.88\pm 0.23$ | $4.82\pm 0.06$ | $4.79\pm 0.06$ | $4.81\pm 0.06$ | $4.88\pm 0.23$
– | $3.40\pm 0.06$ | $3.40\pm 0.06$ | $3.41\pm 0.11$ | – | – | – | –
$D_{B_{\varphi}}$ | $D_{B_{\varphi}}$
– | $7.95\pm 0.15$ | $7.95\pm 0.15$ | $7.88\pm 0.61$ | $8.53\pm 0.18$ | $8.62\pm 0.18$ | $8.58\pm 0.18$ | –
$5.85\pm 0.08$ | $5.89\pm 0.08$ | $5.87\pm 0.08$ | $5.69\pm 0.32$ |  $5.89\pm 0.09$* |  $5.96\pm 0.09$* | $5.93\pm 0.09$ |  $5.86\pm 0.34$*
$4.82\pm 0.06$ | $4.79\pm 0.06$ | $4.81\pm 0.06$ | $4.65\pm 0.21$ | – | $4.68\pm 0.05$ | $4.68\pm 0.05$ | $4.65\pm 0.21$

All values marked by asterisks indicate the weighted mean periods calculated
from multiplets.

Table 4: Same as in Table 2, but for AR 12524.

M1 | M2 | Averaged periods | M3
---|---|---|---
$D_{B_{LOS}}$
– |  $3.73\pm 0.04$* | $3.73\pm 0.04$ | –
– | $2.84\pm 0.04$ | $2.84\pm 0.04$ | $2.84\pm 0.08$
$D_{B_{r}}$
$6.02\pm 0.09$ | $6.06\pm 0.09$ | $6.04\pm 0.09$ | –
$4.60\pm 0.05$ | $4.58\pm 0.05$ | $4.59\pm 0.05$ | $4.65\pm 0.21$
$3.98\pm 0.04$ | $3.96\pm 0.04$ | $3.97\pm 0.04$ | –
$D_{B_{\theta}}$
– | $5.73\pm 0.08$ | $5.73\pm 0.08$ | $5.69\pm 0.32$
– | $4.58\pm 0.06$ | $4.58\pm 0.06$ | $4.88\pm 0.23$
$D_{B_{\varphi}}$
– | – | – | –

All values marked by asterisks indicate the weighted mean periods calculated
from multiplets.

Table 5: Detected periods (with corresponding errors) for the radial unsigned
magnetic flux for each active region. The formatting is the same as in
previous tables.

M1 | M2 | Averaged periods | M3
---|---|---|---
AR 12378
$6.21\pm 0.09$ | $6.25\pm 0.09$ | $6.23\pm 0.09$ | $6.02\pm 0.36$
AR 12381
– | – | – | –
AR 12435
– | $7.38\pm 0.14$ | $7.38\pm 0.14$ | $7.31\pm 0.52$
AR 12437
$6.03\pm 0.09$* | – | $6.03\pm 0.09$ | $6.02\pm 0.36$
$4.80\pm 0.06$ | $4.91\pm 0.06$ | $4.86\pm 0.06$ | $4.65\pm 0.21$
AR 12524
– | $6.07\pm 0.09$ | $6.07\pm 0.09$ | –
– | $4.58\pm 0.05$ | $4.58\pm 0.05$ | $4.65\pm 0.21$
– | $3.96\pm 0.04$ | $3.96\pm 0.04$ | –

All values marked by asterisks indicate the weighted mean periods calculated
from multiplets.

##  Appendix C Example eigenspectra

Here, we provide an example of the eigenspectra of the AR 12381 area. Fig. 4
shows the obtained frequency spectra for AR 12381 with all three methods that
are shown with the detected peaks, as well as the 12-hour artefact frequencies
(Liu et al. 2012).

Figure 4: Example of power spectra obtained for the area of AR 12381 using
methods M1, M2, and M3 for $B_{LOS}$ (the first row of panels from top),
$B_{r}$ (the second row), $B_{\theta}$ (the third row), $B_{\varphi}$ (the
fourth row) data. The solid blue line presents the model line ($\lg
W_{model}=m\lg f+b$) and the red dashed line indicates the 95% confidence
level. The arrows indicate the significant peaks for each case, accordingly.
The colouring of all arrows follows that used in Table 2. The thick black
arrows denote the artefact periods of 12 hours (Liu et al. 2012).

## References

  * Abramov-Maximov et al. (2013) Abramov-Maximov, V. E., Efremov, V. I., Parfinenko, L. D., Solov’ev, A. A., & Shibasaki, K. 2013, PASJ, 65, S12
  * Andries et al. (2009) Andries, J., van Doorsselaere, T., Roberts, B., et al. 2009, Space Sci. Rev., 149, 3
  * Appourchaux (2003) Appourchaux, T. 2003, A&A, 412, 903
  * Baiesi et al. (2008) Baiesi, M., Maes, C., & Shergelashvili, B. M. 2008, Physica A Statistical Mechanics and its Applications, 387, 167
  * Belkasim et al. (1991) Belkasim, S., Shridhar, M., & Ahmadi, M. 1991, Pattern Recognition, 24, 1117
  * Bobra et al. (2014) Bobra, M. G., Sun, X., Hoeksema, J. T., et al. 2014, Sol. Phys., 289, 3549
  * Bogart et al. (2011a) Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & Rabello-Soares, M. C. 2011a, Journal of Physics: Conference Series, 271, 012008
  * Bogart et al. (2011b) Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & Rabello-Soares, M. C. 2011b, Journal of Physics: Conference Series, 271, 012009
  * Centeno et al. (2006) Centeno, R., Collados, M., & Trujillo Bueno, J. 2006, ApJ, 640, 1153
  * Chorley et al. (2010) Chorley, N., Hnat, B., Nakariakov, V. M., Inglis, A. R., & Bakunina, I. A. 2010, A&A, 513, A27
  * Dumbadze et al. (2017) Dumbadze, G., Shergelashvili, B. M., Kukhianidze, V., et al. 2017, A&A, 597, A93
  * Efremov et al. (2007) Efremov, V. I., Parfinenko, L. D., & Solov’ev, A. A. 2007, Astronomy Reports, 51, 401
  * Fleck & Schmitz (1991) Fleck, B. & Schmitz, F. 1991, A&A, 250, 235
  * Flusser (2007) Flusser, J. 2007, International Journal of Computer, Electrical, Automation, Control and Information Engineering, 1, 3708
  * Flusser & Suk (1994) Flusser, J. & Suk, T. 1994, IEEE Transactions on Geoscience and Remote Sensing, 32, 382
  * Gelfreikh et al. (2006) Gelfreikh, G. B., Nagovitsyn, Y. A., & Nagovitsyna, E. Y. 2006, PASJ, 58, 29
  * Goldvarg et al. (2005) Goldvarg, T. B., Nagovitsyn, Y. A., & Solov’Ev, A. A. 2005, Astronomy Letters, 31, 414
  * Gopasyuk (2004) Gopasyuk, O. S. 2004, in IAU Symposium, Vol. 223, Multi-Wavelength Investigations of Solar Activity, ed. A. V. Stepanov, E. E. Benevolenskaya, & A. G. Kosovichev, 249–250
  * Goshtasby (1985) Goshtasby, A. 1985, IEEE Transactions on Geoscience and Remote Sensing, 7, 338
  * Hoeksema et al. (2014) Hoeksema, J. T., Liu, Y., Hayashi, K., et al. 2014, Sol. Phys., 289, 3483
  * Hu (1962) Hu, M.-K. 1962, IRE Transactions on Information Theory, 8
  * Jaeggli & Norton (2016) Jaeggli, S. A. & Norton, A. A. 2016, ApJ, 820, L11
  * Khutsishvili et al. (1998) Khutsishvili, E., Kvernadze, T., & Sikharulidze, M. 1998, Sol. Phys., 178, 271
  * Kuridze et al. (2009) Kuridze, D., Zaqarashvili, T. V., Shergelashvili, B. M., & Poedts, S. 2009, A&A, 505, 763
  * Liu et al. (2012) Liu, Y., Hoeksema, J. T., Scherrer, P. H., et al. 2012, Sol. Phys., 279, 295
  * Maes et al. (2009) Maes, C., Netočný, K., & Shergelashvili, B. M. 2009, Phys. Rev. E, 80, 011121
  * Moradi (2012) Moradi, H. 2012, Astronomische Nachrichten, 333, 1003
  * Nagovitsyn & Nagovitsyna (2011) Nagovitsyn, Y. A. & Nagovitsyna, E. Y. 2011, Geomagnetism and Aeronomy, 51, 1049
  * Parker (1979) Parker, E. N. 1979, ApJ, 234, 333
  * Press et al. (2007) Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing (Cambridge University Press)
  * Prokop & Reeves (1992) Prokop, R. J. & Reeves, A. P. 1992, CVGIP: Graph. Models Image Process., 54, 438
  * Pugh et al. (2017) Pugh, C. E., Broomhall, A. M., & Nakariakov, V. M. 2017, A&A, 602, A47
  * Rincon & Rieutord (2018) Rincon, F. & Rieutord, M. 2018, Living Reviews in Solar Physics, 15, 6
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229
  * Shergelashvili et al. (2007) Shergelashvili, B. M., Maes, C., Poedts, S., & Zaqarashvili, T. V. 2007, Phys. Rev. E, 76, 046404
  * Shergelashvili & Poedts (2005) Shergelashvili, B. M. & Poedts, S. 2005, A&A, 438, 1083
  * Shergelashvili et al. (2005) Shergelashvili, B. M., Zaqarashvili, T. V., Poedts, S., & Roberts, B. 2005, A&A, 429, 767
  * Smirnova et al. (2013) Smirnova, V., Riehokainen, A., Solov’ev, A., et al. 2013, A&A, 552, A23
  * Solov’ev & Kirichek (2008) Solov’ev, A. A. & Kirichek, E. A. 2008, Astrophysical Bulletin, 63, 169
  * Sun (2013) Sun, X. 2013, arXiv e-prints, arXiv:1309.2392
  * Terradas et al. (2008) Terradas, J., Arregui, I., Oliver, R., et al. 2008, ApJ, 679, 1611
  * Thomas et al. (1984) Thomas, J. H., Cram, L. E., & Nye, A. H. 1984, ApJ, 285, 368
  * Tsirikolias & Mertzios (1993) Tsirikolias, K. & Mertzios, B. 1993, Pattern Recognition, 26, 877
  * Vaughan (2005) Vaughan, S. 2005, A&A, 431, 391


# Blind Channel Estimation for MIMO Systems via Variational Inference

Jiancheng Tang, Qianqian Yang2, Zhaoyang Zhang 0*This work is partly supported
by the SUTD-ZJU IDEA Grant (SUTD-ZJU (VP) 202102), and partly by the
Fundamental Research Funds for the Central Universities under Grant
2021FZZX001-20. College of Information Science and Electronic Engineering,
Zhejiang University, Hangzhou 310007, China
Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we investigate the blind channel estimation problem for MIMO
systems under Rayleigh fading channel. Conventional MIMO communication
techniques require transmitting a considerable amount of training symbols as
pilots in each data block to obtain the channel state information (CSI) such
that the transmitted signals can be successfully recovered. However, the pilot
overhead and contamination become a bottleneck for the practical application
of MIMO systems with the increase of the number of antennas. To overcome this
obstacle, we propose a blind channel estimation framework, where we introduce
an auxiliary posterior distribution of CSI and the transmitted signals given
the received signals to derive a lower bound to the intractable likelihood
function of the received signal. Meanwhile, we generate this auxiliary
distribution by a neural network based variational inference framework, which
is trained by maximizing the lower bound. The optimal auxiliary distribution
which approaches real prior distribution is then leveraged to obtain the
maximum a posterior (MAP) estimation of channel matrix and transmitted data.
The simulation results demonstrate that the performance of the proposed blind
channel estimation method closely approaches that of the conventional pilot-
aided methods in terms of the channel estimation error and symbol error rate
(SER) of the detected signals even without the help of pilots.

##  I Introduction

With the capability of substantially improving the channel capacity and
spectrum efficiency, multi-input multi-output (MIMO) approach has been
recognized as one of the key enabling techniques for reliable high data rate
transmission[1],[2]. And the orthogonal frequency division multiplexing (OFDM)
approach is deployed in MIMO systems to eliminate the intersymbol interference
(ISI) resulted form the frequency-selective fading[3] by breaking the
broadband system into many different narrowband subchannels. However, the
transmitted signals in subchannels still experience power attenuation caused
by flat fading which impedes the reliable detection of received signals.
Therefore the accurate estimation of channel state information (CSI) is
essential to compensate the channel responses and recover the transmitted
information[4].

The channel estimation techniques for MIMO systems have been studied
extensively in the past [5],[6], and the existing techniques can be summarized
in two categories: pilot-aided estimation and blind estimation. The former
approach requires to transmit orthogonal pilot sequences in order to estimate
the channel matrices, which occupies a considerable portion of spectrum
resources thus reduces spectrum efficiency [7]. Furthermore, the overlap of
frequency by adjacent cells in cellular system impedes the orthogonality of
pilot, which is referred to as pilot contamination[8]. On the other hand, the
blind estimation methods are studied to overcome the pilot contamination
problem and reduce the spectrum overhead. The second-order statistics of
received signals is utilized to derive the blind channel estimation algorithm
in [9], where the covariance matrix of the received signal is divided into an
orthogonal signal subspace and a noise subspace by singular value
decomposition (SVD). However the SVD operation of large-scale matrix works in
an iterative manner, which result in prohibitive computational cost [10].

Recently, deep learning (DL) based channel estimation methods have achieved
promising results and attracted a lot of research interests [11]. In
[12],[13], the time-frequency response of pilots going through fading channels
is regarded as a low-resolution image, which is then enhanced by a DL-based
image super-resolution (SR) network to estimate the CSI. These works are
designed for single user scenarios and still require a certain amount of
pilots to initialize the time-frequency response images. In [14], a blind
channel estimator based on a denoising convolutional neural network (DnCNN) is
proposed for massive MIMO systems, where the DnCNN is employed to remove the
residual noise effects which cannot be averaged out through channel asymptotic
orthogonality. However, the performance of this method degrades in limited-
scale MIMO systems as the asymptotic orthogonality of the channel vectors only
exists in large scale MIMO systems.

In this paper, inspired by the widely applied approximate inference
algorithms, i.e., variational inference (VI) which could provide a good
approximation to the complex distribution[15], we propose a VI based blind
channel estimation technique for MIMO systems. The conventional signal
detection method exploits the maximum a posterior (MAP) receiver to perform
the exact inference of transmitted signals, which requires the posterior
probability of the CSI. However, this may be intractable due to the
dimensionality of the channel matrix. To deal with this challenge, we utilize
a neural network based VI framework to generate an auxiliary distribution,
which is optimized by maximizing a lower bound to the log-likelihood of the
received signal. The derived auxiliary distribution approximates the real
distribution, from which, the estimation of the CSI can be sampled, and thus
the transmitted signals can be detected using maximum likelihood estimation.
The numerical results validate the effectiveness of the proposed blind channel
estimation approach, compared to the conventional pilot-aided methods, for
MIMO systems of different scales.

The rest of this paper is organized as follows. In Section II, we present the
system model for the uplink transmission of MIMO system. In Section III, we
introduce the proposed blind channel estimation framework based on VI. Then,
the performance of the proposed approach is numerically evaluated in Section
IV. Finally, we conclude the paper in Section V.

Notions: matrices and vectors are denoted by boldface symbols, $\mathbb{C}$
denotes the set of complex values and ${\rm N}_{+}$ represents positive
integer. ${\left(\cdot\right)^{T}}$ and ${\left(\cdot\right)^{H}}$ denote
transpose, Hermitian respectively. The elements $i$ and $\left({i,j}\right)$
of the vector $a$ and matrix ${\bf{A}}$ are represented by ${a_{i}}$ and
${{\bf{A}}_{ij}}$, ${\bf{I}}_{K}$ is $K\times K$ identity matrix, and
${tr\left(\cdot\right)}$ represents the trace of matrix. ${\rho^{2}}$is the
transmitting power. $E\left(\cdot\right)$ and $Var\left(\cdot\right)$ denote
the expectation and variance of a random variable,
$\mathcal{CN}\left({{\bf{\mu}},{\bf{\Sigma}}}\right)$ denotes circularly
symmetric complex Gaussian random vectors with mean $\bf{\mu}$ and covariance
matrix $\bf{\Sigma}$.

##  II SYSTEM MODEL AND PRELIMINARIES

Consider the uplink transmission of a MIMO system under the Rayleigh fading
channel, where $K$ single-antenna mobile stations (MS) are uploading data to a
base station (BS) with $N$ antennas. Typical transmission operates in two
phases; (i) _the channel estimation phase_ where known pilot sequences are
transmitted to the BS to obtain the CSI, and (ii) _the signal detection phase_
when the estimated CSI is used to detect transmitted signals [14]. However,
different from this pilot-aided approach which can result in the loss of
bandwidth efficiency and pilot contamination, the proposed method in this
paper aims to recover the sent signal without the help of pilots. In the
channel estimation phase, instead of sending pilots, the MSs send data in an
one by one order such that while one is sending during its allocated time
slot, the others remain silent to maintain orthogonality. The received signal
by the BS at the end of the _t_ th time slot can be represented as

${{\bf{y}}_{t}}={\bf{H}}{{\bf{x}}_{t}}+{{\bf{n}}_{t}},$ (1)

where the channel matrix ${\bf{H}}\in{\mathbb{C}^{N\times K}}$ is considered
to remain constant over a number of transmissions and then changes to a new
state according to block fading distribution, and the elements of $\bf{H}$ are
independent random variables that follow zero-mean complex Gaussian
distribution with unit variance, i.e., ${h_{ij}}\sim\mathcal{CN}(0,1)$.
${{\mathbf{x}}_{t}}={\left[{{x_{1}}\left(t\right),...,{x_{K}}\left(t\right)}\right]^{T}}\in{\mathbb{C}^{K\times
1}}$ is the transmitted signal vector during the _t_ th time slot of the
channel estimation phase with only one of the elements is non-zero
corresponding to the user allocated to transmit in this time slot. We also
note that the independent bit streams at each MB are mapped using the same
constellation map.
${{\mathbf{y}}_{t}}={\left[{{y_{1}}\left(t\right),...,{y_{N}}\left(t\right)}\right]^{T}}\in{\mathbb{C}^{N\times
1}}$is the received signals vector by BS, and ${{\mathbf{n}}_{t}}$ is the
additive noises vector following zero-mean complex Gaussian distribution.

The goal of receiver design in channel estimation phase is to recover the
transmitted signal ${{\mathbf{x}}_{t}}$ and then estimate the channel matrix
$\bf{H}$ from the received signals ${{\mathbf{y}}_{t}}$, and then the obtained
channel matrix $\bf{H}$ is utilized in the signal detection phase to detect
the transmitted signals. The MAP estimation of the sent signal $\bf{x}_{t}$
can be expressed as

$\displaystyle{{{\bf{\hat{x}}}}_{t}}^{MAP}$
$\displaystyle=\mathop{\arg\max}\limits_{\bf{x}}\int\limits_{\bf{H}}{p\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right){\rm{}}}d{\bf{H}}$
(2)
$\displaystyle=\mathop{\arg\max}\limits_{\bf{x}}\int\limits_{\bf{H}}{p\left({{{\bf{x}}_{t}}|{\bf{H}},{{\bf{y}}_{t}}}\right)}p\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)d{\bf{H}},$

where the second equation can be derived by applying Bayesian criterion.
Unfortunately, due to the high dimensionality of channel matrix $\bf{H}$, it
is infeasible to directly calculate ${{{\bf{\hat{x}}}}_{t}}^{MAP}$ by (2).
Hence, we propose an approximate inference method instead that approximates
${{{\bf{\hat{x}}}}_{t}}^{MAP}$ and estimates $\bf{H}$ simultaneously via a
neural network based variance inference framework which will be introduced in
details in next section.

##  III OPTIMIZATION FRAMEWORK FOR BLIND ESTIMATION

In this section, we present the proposed method for simultaneously estimating
the transmitted symbols and channel matrix without the assistance of pilots.
Note that the log-likelihood function of the received signal vector
${{\mathbf{y}}_{t}}$ during the $t$th time slot of the channel estimation
phase can be written as

$\displaystyle\log p\left({{{\bf{y}}_{t}}}\right)$
$\displaystyle=\int_{{{\bf{x}}_{t}},{\bf{H}}}\log
p\left({{{\bf{y}}_{t}}}\right)\cdot{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}d{{\bf{x}}_{t}}d{\bf{H}}$
(3a)
$\displaystyle=\int_{{{\bf{x}}_{t}},{\bf{H}}}{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}$
$\displaystyle\cdot\log\frac{{p\left({{{\bf{x}}_{t}},{\bf{H}},{{\bf{y}}_{t}}}\right)}}{{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\frac{{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}{{p\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}d{{\bf{x}}_{t}}d{\bf{H}}$
(3b) $\displaystyle={E_{{{\bf{x}}_{t}}{\bf{,H}}\sim
q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{\log\frac{{p\left({{{\bf{x}}_{t}},{\bf{H}},{{\bf{y}}_{t}}}\right)}}{{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}}\right]$
$\displaystyle+\underbrace{{E_{{{\bf{x}}_{t}}{\bf{,H}}\sim
q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{\log\frac{{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}{{p\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}}\right]}_{Kullback-
Leibler\enspace divergence}$ (3c)

where an auxiliary distribution
$q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)$ is introduced in
(3a), and (3b) is derived by expanding ${p\left({{{\bf{y}}_{t}}}\right)}$. We
note that the _Kullback-Leibler divergence_ (KL-D)[16] between the posterior
distribution $p\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)$ and the
auxiliary distribution
$q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)$ is non-negative,
therefore a lower bound of the intractable likelihood function of received
signal ${{\mathbf{y}}_{t}}$ is derived through VI framework. We have

$\displaystyle\log p\left({{{\bf{y}}_{t}}}\right)$
$\displaystyle\geq{E_{{{\bf{x}}_{t}}{\bf{,H}}\sim
q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{\log\frac{{p\left({{{\bf{y}}_{t}},{\bf{H}},{{\bf{x}}_{t}}}\right)}}{{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}}\right]$
(4) $\displaystyle\buildrel\Delta\over{=}-L\left(q\right).$

Since $\log p\left({{{\bf{y}}_{t}}}\right)$ is a unknown constant by (3c), we
can minimize the KL-divergence
${D_{KL}}\left({p\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)||q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}\right)$
by maximizing the lower bound $-L\left(q\right)$ over the parameters
${{\bf{x}}_{t}}$ and ${\bf{H}}$, by which the distribution
$q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)$ approximates the
posterior distribution
$p\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)$. Thus we can obtain
the maximum likelihood estimation of channel matrix ${\bf{H}}$ and transmitted
signals ${{\bf{x}}_{t}}$ with
$q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)$. We exploit the mean-
field approximation[17] to further simplify the maximum likelihood estimation
problem by assuming
$q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)=q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)$.
We also assume $q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)$ and
$q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)$ follow complex Gaussian
distribution as follows

$\displaystyle
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)\sim\mathcal{CN}\left({{m_{{{\bf{x}}_{t}}}},{S_{{{\bf{x}}_{t}}}}}\right),$
(5) $\displaystyle
q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)\sim\mathcal{CN}\left({{m_{{{\bf{H}}}}},{S_{{{\bf{H}}}}}}\right),$

where ${m_{{{\bf{x}}_{t}}}}$, ${m_{{{\bf{H}}}}}$, $S_{{{\bf{x}}_{t}}}$,
$S_{{{\bf{H}}}}$ denote the means and variances of the Gaussian distribution
respectively. We obtain these parameters by two trainable neural networks
denoted by _g_ and _f_ respectively as shown in Fig. 1

$\displaystyle\left[{{m_{{{\bf{x}}_{t}}}},{S_{{{\bf{x}}_{t}}}}}\right]{\rm{=}}g\left({{{\bf{y}}_{t}}\left|\phi\right.}\right),$
(6)
$\displaystyle{\rm{}}\left[{{m_{{{\bf{H}}}}},{S_{{{\bf{H}}}}}}\right]=f\left({{{\bf{y}}_{t}}\left|\varphi\right.}\right),$

referred to as $Encoder1$ and $Encoder2$ of the proposed blind channel
estimation framework. The inputs to these two encoders are the received
signals ${{\bf{y}}_{t}}$, where the two channels corresponding to the real and
imaginary elements. Both $Encoder1$ and $Encoder2$ have two fully connected
layer with $\phi$ and $\varphi$ being the trainable parameters of the
networks. The estimated ${m_{{{\bf{x}}_{t}}}}$, ${m_{{{\bf{H}}}}}$,
$S_{{{\bf{x}}_{t}}}$, $S_{{{\bf{H}}}}$ are used to generate the complex
Gaussian distribution
$\mathcal{CN}\left({{m_{{{\bf{x}}_{t}}}},{S_{{{\bf{x}}_{t}}}}}\right)$ and
$\mathcal{CN}\left({{m_{{{\bf{H}}_{t}}}},{S_{{{\bf{H}}_{t}}}}}\right)$
respectively, from which we sample the estimation of ${\bf{x}}_{k}$ and
${\bf{H}}$, denoted by $\widehat{{{\bf{x}}_{t}}}$, ${\widehat{{\bf{H}}}}$.
Finally, we generate $\widehat{{{\bf{y}}_{t}}}$ by multiplying
$\widehat{{{\bf{x}}_{t}}}$ and ${\widehat{{\bf{H}}}}$, which is referred to as
$Decoder$ of the proposed framework. We train the whole network structure by
minimizing the objective function $L\left(q\right)$ where
$q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)=q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)$,
that is the product of the two normal distribution generate by $Encoder1$ and
$Encoder2$, $L\left(q\right)$ can be written as

$\displaystyle L\left(q\right)$
$\displaystyle=\int_{{{\bf{x}}_{t}},{\bf{H}}\sim
q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}\log\frac{{q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}{{p\left({{{\bf{x}}_{t}},{\bf{H}},{{\bf{y}}_{t}}}\right)}}d{{\bf{x}}_{t}}d{\bf{H}}$
(7) $\displaystyle=\underbrace{{E_{{\bf{H}},{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{\log\frac{{q\left({{{\bf{x}}_{t}}{\bf{|}}{{\bf{y}}_{t}}}\right)}}{{p\left({{{\bf{x}}_{t}}}\right)}}}\right]}_{los{s_{1}}}$
$\displaystyle+\underbrace{{E_{{\bf{H}},{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{\log\frac{{q\left({{\bf{H|}}{{\bf{y}}_{t}}}\right)}}{{p\left({\bf{H}}\right)}}}\right]}_{los{s_{2}}}$
$\displaystyle-\underbrace{{E_{{\bf{H}},{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}},{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
p\left({{{\bf{y}}_{t}}{\bf{|H,}}{{\bf{x}}_{t}}}\right)}\right]}_{los{s_{3}}}.$

We have

$\displaystyle los{s_{1}}$ (8) $\displaystyle={E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log\frac{{q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}{{p\left({{{\bf{x}}_{t}}}\right)}}}\right]$
$\displaystyle={E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\right]-{E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
p\left({{{\bf{x}}_{k}}}\right)}\right].$

Note that the term ${E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log{q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\right]$
is the entropy of multivariate normal distribution
${q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\sim\mathcal{CN}\left({{m_{{{\bf{x}}_{t}}}},{S_{{{\bf{x}}_{t}}}}}\right)$.
We have

$\displaystyle{E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\right]=-\frac{1}{2}\log\left|{{S_{{{\bf{x}}_{t}}}}}\right|+C_{1},$
(9)

where $C$ is a constant term and we use $C_{i}$, $i\in{{\rm N}_{+}}$, to
denote different constants in the sequel. And the term ${E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log{p\left({{{\bf{x}}_{t}}}\right)}}\right]$
in (8) can be written as

$\displaystyle-{E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log{p\left({{{\bf{x}}_{t}}}\right)}}\right]$
$\displaystyle=\frac{1}{{2{\rho^{2}}}}{E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{{{\bf{x}}_{t}}^{H}{{\bf{x}}_{t}}}\right]+C_{2}$
(10)
$\displaystyle=\frac{1}{{2{\rho^{2}}}}\left({tr({S_{{{\bf{x}}_{t}}}})+{m_{{{\bf{x}}_{t}}}}^{T}{m_{{{\bf{x}}_{t}}}}}\right)+C_{2},$

where the discrete sample space constrain of ${{\bf{x}}_{t}}$ is relaxed to be
continuous via the assumption
$p\left({{{\bf{x}}_{t}}}\right)\sim\mathcal{CN}\left({0,2{\rho^{2}}{\bf{I}}_{K}}\right)$
for the convenience of the calculation. And hence the $los{s_{1}}$ can be
rewritten as

$\displaystyle
los{s_{1}}=\frac{1}{{2{\rho^{2}}}}\left({tr({S_{{{\bf{x}}_{t}}}})+{m_{{{\bf{x}}_{t}}}}^{T}{m_{{{\bf{x}}_{t}}}}}\right)-\frac{1}{2}\log\left|{{S_{{{\bf{x}}_{t}}}}}\right|+C_{3}.$
(11)

The detailed proof of equations (9)-(11) can be found in Appendix I.
Similarly, we have

$\displaystyle los{s_{2}}$ (12) $\displaystyle={E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{\log\frac{{q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}}{{p\left({\bf{H}}\right)}}}\right]$
$\displaystyle={E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{k}}}\right)}}\left[{\log
q\left({{\bf{H}}|{{\bf{y}}_{k}}}\right)}\right]-{E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{k}}}\right)}}\left[{\log
p\left({\bf{H}}\right)}\right].$

Following the same procedure of deriving (9)-(11), we have

$\displaystyle{E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{k}}}\right)}}\left[{\log
q\left({{\bf{H}}|{{\bf{y}}_{k}}}\right)}\right]=-\frac{1}{2}\log\left|{{S_{\bf{H}}}}\right|+C_{4},$
(13)

and the term ${E_{{{\bf{H}}}\sim
q\left({{{\bf{H}}}|{{\bf{y}}_{t}}}\right)}}\left[{\log{p\left({{{\bf{H}}}}\right)}}\right]$
in (12) can be written as

$\displaystyle-{E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{k}}}\right)}}\left[{\log
p\left({\bf{H}}\right)}\right]$ $\displaystyle=\frac{1}{2}{E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{{{\bf{H}}^{H}}{\bf{H}}}\right]+C_{5}$
(14)
$\displaystyle=\frac{1}{2}tr({S_{\bf{H}}})+\frac{1}{2}{m_{\bf{H}}}^{T}{m_{\bf{H}}}+C_{5}.$

Hence we rewritten $los{s_{2}}$ as

$\displaystyle
los{s_{2}}=\frac{1}{2}_{t}r({S_{\bf{H}}})+\frac{1}{2}{m_{\bf{H}}}^{T}{m_{\bf{H}}}-\frac{1}{2}\log\left|{{S_{\bf{H}}}}\right|+C_{6}.$
(15)

We employ Monte Carlo method to compute the $los{s_{3}}$ in (7), we have

$\displaystyle los{s_{3}}$ $\displaystyle=-{E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right),{\bf{x}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
p\left({{{\bf{y}}_{t}}|{\bf{H}},{{\bf{x}}_{t}}}\right)}\right]$ (16)
$\displaystyle\approx\frac{1}{L}\sum\limits_{l=1}^{L}{}{tr({\widehat{{\bf{H}}}_{\rm{l}}}{S_{{{\bf{x}}_{t}}}}{\widehat{{\bf{H}}}_{\rm{l}}}^{H})}$
$\displaystyle+\frac{1}{L}\sum\limits_{l=1}^{L}{}{{{({\widehat{{\bf{H}}}_{\rm{l}}}{m_{{{\bf{x}}_{t}}}}-{{\bf{y}}_{t}})}^{H}}({\widehat{{\bf{H}}}_{\rm{l}}}{m_{{{\bf{x}}_{t}}}}-{{\bf{y}}_{t}})},$

Figure 1: The blind estimation framework structure.

where ${\widehat{{\bf{H}}}_{\rm{l}}}$ are sampled from
${q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}\sim\mathcal{CN}\left({{m_{{{\bf{H}}_{t}}}},{S_{{{\bf{H}}_{t}}}}}\right)$
and $L$ is the number of sample points. The proof of (16) can also be found in
Appendix.

We note that (8) and (12) represent the KL-D between the posterior
distributions generated by the proposed neural network and the actual priors
distribution of ${{\bf{x}}_{t}}$ and ${{\bf{H}}}$ respectively. $los{s_{3}}$
represents the reconstruction error $\widehat{{{\bf{y}}_{t}}}$ with the
variational distributions ${q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}$
and ${q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}$. Hence minimizing the
objective function $L\left(q\right)=los{s_{1}}+los{s_{2}}+los{s_{3}}$ pushes
the generated posterior distribution approach the prior distribution and the
reconstructed signal $\widehat{{{\bf{y}}_{t}}}$ to approach the actual
received signal ${{{\bf{y}}_{t}}}$. Thus when the loss function converges,
reasonable estimation results about ${\bf{H}}$ and decision results
${{\bf{x}}_{t}}$ about can be obtained.

##  IV NUMERICAL ANALYSIS

(a) signals before equalization with QPSK

(b) signals after equalization with QPSK

(c) signals before equalization with 16QAM

(d) signals after equalization with 16QAM

Figure 2: Blind equalization result of receiver setup with 4 antennas in
medium SNR scenario: (a) QPSK modulation before blind equalization; (b) QPSK
modulation after blind equalization; (c) 16QAM modulation before blind
equalization; (d) 16QAM modulation after blind equalization.

In this section, we numerically evaluate the performance of our proposed blind
channel estimation framework. A MIMO system in which $K=4$ users are
communicating with a BS through QPSK/16QAM modulation is considered here. Both
$Encoder1$ and $Encoder2$ of the proposed framework are fully connected neural
networks, each of which consists of an input layer, a 16 node hidden layer
with tanh activation and an output layer with tanh activation. The Adam
optimizer[18] with an initial learning rate of 0.05 is used to train the whole
network. The pilot-aided channel estimation methods [7] is used as the
benchmark.

We first demonstrate in Fig. 2 the equalization results by the proposed blind
estimation method through the constellation graph. It is seen from Fig. 2(a)
that the transmitted signals during the signal detection phase interfered by
multi-path fading are overlapped with each other, and the proposed method is
able to separate the overlapped signal points as shown in Fig. 2(b) which
implies the effectiveness of our blind estimation algorithm. When the
modulation order increased from QPSK to 16QAM, our algorithm still works
effectively as shown in Fig. 2(c) and Fig. 2(d), which demonstrates the
generalization ability of our method among different order modulation schemes.

Figure 3: The MSE with respect to SNR for $K=4$ users and $N=4$ receive
antennas with QPSK modulation. Figure 4: The SER performance with respect to
SNR for $K=4$ users and $N=4$ antennas with QPSK modulation.  Figure 5: The
SER performance with respect to SNR for $K=4$ users and $N=40$ antennas with
QPSK modulation.

Then we compare the channel estimation performance of conventional pilot-aided
approach with minimum mean squared error algorithm (Aided-MMSE), pilot-aided
with least square algorithm (Aided-LS), and the proposed method in terms of
through mean square error (MSE) depicted in Fig. 3. We can observe that in the
low-to-medium SNR region, knowing the exact pilot symbols and utilizing the
statistics of channel to eliminate the AWGN, Aided-MMSE achieves the best
performance, while the estimation error of our blind method is slightly higher
due to the strong AWGN interference. However, the proposed method still
outperforms Aided-LS which does not use the channel statistics even without
the assistance of pilots. It can also be observed that in the medium-to-high
SNR region, the estimation performance of the three methods is very close
which validates the effectiveness of the proposed channel estimation method.

We also simulate the symbol error rate (SER) in signal detection phase, where
the maximum-likelihood detector (MLD) with CSI estimated in channel estimation
phase is employed, we also use the results by the MLD with perfect CSI
estimation as a benchmark. As shown in Fig. 4, it can be seen that the signal
detection performance of proposed blind method in terms of SER approaches the
Aided-MMSE detector while slightly outperforms the Aided-LS detector in low-
to-medium SNR scenario due to the channel estimation performance as shown in
Fig. 3. With the increase of SNR, the performance gap between the three
methods becomes neglectable. Meanwhile, the performance gap between the
detector with perfect CSI and the three detectors with estimated CSI decreases
as well due to the increase of estimation accuracy. We emphasize that our
method could save spectrum resources and mitigate the pilot contamination
since there is no need for pilot signals. We also evaluate signal detection
performance when there are more receiver antennas ($K=40$) in Fig. 5 and
observe the same performance trend, which proves the compatibility of our
blind method with MIMO systems of different scales.

##  V CONCLUSIONS

In this paper, we proposed a novel blind channel estimation approach for MIMO
systems experiencing Rayleigh fading by exploiting variational inference and
neural network. We derived a lower bound to the intractable log-likelihood of
received signal by introducing an auxiliary which is generated by a neural
network based framework. By training the neural network to maximize the lower
bound, the auxiliary posterior distribution closely approaches the real
distribution, by sampling from which the estimation of CSI can be obtained. We
then numerically compared the proposed blind estimation method with the
conventional pilot-aided methods in terms of the channel estimation error and
SER of the detected signals which demonstrated that the proposed method
outperforms the pilot-aided scheme with LS algorithm, and closely approaches
the performance of the pilot-aided scheme with MMSE algorithm while saving
spectrum resource and mitigating the pilot contamination problem.

## APPENDIX I

In this section, the detailed derivation of each term in the blind estimation
framework (7) is given. For ${E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\right]$, the probability
density function (PDF) of the multivariate normal distribution
$q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)\sim\mathcal{CN}\left({{m_{{{\bf{H}}_{t}}}},{S_{{{\bf{H}}_{t}}}}}\right)$
can be written as

$\displaystyle q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)$
$\displaystyle={\left({2\pi}\right)^{{\rm{-}}\frac{N\times
K}{2}}}{\left|{{S_{\bf{H}}}}\right|^{{\rm{-}}\frac{1}{2}}}$ (17)
$\displaystyle\times\exp\left[{{\rm{-}}\frac{1}{2}{{\left({{\bf{H}}-{m_{\bf{H}}}}\right)}^{H}}{S_{\bf{H}}}^{{\rm{-1}}}\left({{\bf{H}}-{m_{\bf{H}}}}\right)}\right]d{{\bf{x}}_{t}}d{\bf{H}},$

and substituting the PDF (17) into the entropy of normal distribution
${E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\right]$, we have

$\displaystyle\int_{x,{\bf{H}}}{q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}\log
q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)d{{\bf{x}}_{t}}d{\bf{H}}$ (18)
$\displaystyle=\int_{x,{\bf{H}}}{q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}\log{\left({2\pi}\right)^{{\rm{-}}\frac{N\times
K}{2}}}{\left|{{S_{\bf{H}}}}\right|^{{\rm{-}}\frac{1}{2}}}$
$\displaystyle\times\exp\left[{{\rm{-}}\frac{1}{2}{{\left({{\bf{H}}-{m_{\bf{H}}}}\right)}^{H}}{S_{\bf{H}}}^{{\rm{-1}}}\left({{\bf{H}}-{m_{\bf{H}}}}\right)}\right]d{{\bf{x}}_{t}}d{\bf{H}}$
$\displaystyle{\rm{=}}-\frac{N\times K}{2}\left({\log
2\pi+1}\right)-\frac{1}{2}\log\left|{{S_{\bf{H}}}}\right|$
$\displaystyle=-\frac{1}{2}\log\left|{{S_{\bf{H}}}}\right|+C_{4},$

where the constant term $-\frac{N\times K}{2}\left({\log 2\pi+1}\right)$ is a
constant and denoted by $C_{4}$, and ${E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}\right]$
can be derived directly in the same way

$\displaystyle\int_{{{\bf{x}}_{t}},{\bf{H}}}{q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\log
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)d{{\bf{x}}_{t}}d{\bf{H}}=-\frac{1}{2}\log\left|{{S_{{{\bf{x}}_{t}}}}}\right|+C_{1}.$
(19)

For $-{E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
p\left({{{\bf{x}}_{t}}}\right)}\right]$, the PDF of
$p\left({{{\bf{x}}_{t}}}\right)\sim\mathcal{CN}\left({0,2{\rho^{2}}{{\bf{I}}_{K}}}\right)$
can be written as

$\displaystyle p\left({{{\bf{x}}_{t}}}\right)$
$\displaystyle={\left({2\pi}\right)^{{\rm{-}}\frac{P}{2}}}{\left|{2{\rho^{2}}{{\bf{I}}_{K}}}\right|^{{\rm{-}}\frac{1}{2}}}$
(20)
$\displaystyle\times\exp\left[{{\rm{-}}{{\bf{x}}_{t}}^{H}{{\left({2{\rho^{2}}{{\bf{I}}_{K}}}\right)}^{{\rm{-1}}}}{{\bf{x}}_{t}}}\right]d{{\bf{x}}_{t}}d{\bf{H}},$

and substituting the PDF (20) into $-{E_{{{\bf{x}}_{t}}}}\left[{\log
p\left({{{\bf{x}}_{t}}}\right)}\right]$, we have

$\displaystyle-\int_{{{\bf{x}}_{t}},{\bf{H}}}{q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\log
p\left({{{\bf{x}}_{t}}}\right)d{{\bf{x}}_{t}}d{\bf{H}}$ (21)
$\displaystyle{\rm{=}}-\int_{{{\bf{x}}_{t}},{\bf{H}}}{q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\log{\left({2\pi}\right)^{{\rm{-}}\frac{K}{2}}}{\left|{2{\rho^{2}}{{\bf{I}}_{K}}}\right|^{{\rm{-}}\frac{1}{2}}}$
$\displaystyle\times\exp\left[{{\rm{-}}{{\bf{x}}_{t}}^{H}{{\left({2{\rho^{2}}{{\bf{I}}_{K}}}\right)}^{{\rm{-1}}}}{{\bf{x}}_{t}}}\right]d{{\bf{x}}_{t}}d{\bf{H}}$
$\displaystyle=\int_{{{\bf{x}}_{t}},{\bf{H}}}{q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\left({\log{{\left({2\pi}\right)}^{\frac{K}{2}}}{{\left|{2{\rho^{2}}{{\bf{I}}_{K}}}\right|}^{\frac{1}{2}}}+\frac{1}{{2{\rho^{2}}}}{{\bf{x}}_{t}}^{H}{{\bf{x}}_{t}}}\right)d{{\bf{x}}_{t}}d{\bf{H}}$
$\displaystyle=\frac{1}{{2{\rho^{2}}}}{E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{{{\bf{x}}_{t}}^{H}{{\bf{x}}_{t}}}\right]+C_{2},$

where the expectation of
${\log{{\left({2\pi}\right)}^{\frac{N}{2}}}{{\left|{2{\rho^{2}}{{\bf{I}}_{K}}}\right|}^{\frac{1}{2}}}}$
is constant and detonated by $C_{2}$, and since
$E\left[{{x^{2}}}\right]=Var\left(x\right)+{E^{2}}\left[x\right]$, (21) can be
represented as

$\displaystyle\frac{1}{{2{\rho^{2}}}}{E_{{{\bf{x}}_{t}}\sim
q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}}\left[{{{\bf{x}}_{t}}^{H}{{\bf{x}}_{t}}}\right]$
$\displaystyle{\rm{}}=\frac{1}{{2{\rho^{2}}}}\left({Var({{\bf{x}}_{t}})+{E^{2}}\left[{{{\bf{x}}_{t}}}\right]}\right)$
(22)
$\displaystyle{\rm{}}=\frac{1}{{2{\rho^{2}}}}\left({tr({S_{{{\bf{x}}_{t}}}})+{m_{{{\bf{x}}_{t}}}}^{T}{m_{{{\bf{x}}_{t}}}}}\right).$

Similarly, for $-{E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{\log
p\left({\bf{H}}\right)}\right]$

$\displaystyle-\int_{\bf{H}}{q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}\log
p\left({\bf{H}}\right)d{\bf{H}}$ $\displaystyle=\frac{1}{2}{E_{{\bf{H}}\sim
q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)}}\left[{{{\bf{H}}^{H}}{\bf{H}}}\right]+C_{5}$
(23)
$\displaystyle=\frac{1}{2}tr({S_{\bf{H}}})+\frac{1}{2}{m_{\bf{H}}}^{T}{m_{\bf{H}}}+C_{5}.$

For the decoder term in (16), the Monte Carlo method is used to approximate
the expectation

$\displaystyle-\int_{{{\bf{x}}_{t}},{\bf{H}}}{q\left({{\bf{H}}|{{\bf{y}}_{t}}}\right)q\left({{{\bf{x}}_{t}}|{{\bf{y}}_{t}}}\right)}\log
p\left({{{\bf{y}}_{t}}|{\bf{H}},{{\bf{x}}_{t}}}\right)d{{\bf{x}}_{t}}d{\bf{H}}$
(24)
$\displaystyle={E_{{\bf{H}},{{{\bf{x}}_{t}}}}}\left[{{{({{\bf{y}}_{t}}-{\bf{H}}{{\bf{x}}_{t}})}^{H}}({{\bf{y}}_{t}}-{\bf{H}}{{\bf{x}}_{t}})}\right]$
$\displaystyle{\rm{=}}{E_{{\bf{H}},{{{\bf{x}}_{t}}}}}\left[{{{\bf{y}}_{t}}^{H}{{\bf{y}}_{t}}-{{\bf{y}}_{t}}^{H}{\bf{H}}{{\bf{x}}_{t}}-{{\bf{x}}_{t}}^{H}{{\bf{H}}^{H}}{{\bf{y}}_{t}}-{{({\bf{H}}{{\bf{x}}_{t}})}^{H}}({\bf{H}}{{\bf{x}}_{t}})}\right]$
$\displaystyle={E_{\bf{H}}}\left[{tr({\bf{H}}{S_{{{\bf{x}}_{t}}}}{{\bf{H}}^{H}})+{{({\bf{H}}{m_{{{\bf{x}}_{t}}}}-{{\bf{y}}_{t}})}^{H}}({\bf{H}}{m_{{{\bf{x}}_{t}}}}-{{\bf{y}}_{t}})}\right]$
$\displaystyle\approx\frac{1}{L}\sum\limits_{{\widehat{{\bf{H}}}_{\rm{l}}}}{tr({\widehat{{\bf{H}}}_{\rm{l}}}{S_{{{\bf{x}}_{t}}}}{\widehat{{\bf{H}}}_{\rm{l}}}^{H})+{{({\widehat{{\bf{H}}}_{\rm{l}}}{m_{{{\bf{x}}_{t}}}}-{{\bf{y}}_{t}})}^{H}}({\widehat{{\bf{H}}}_{\rm{l}}}{m_{{{\bf{x}}_{t}}}}-{{\bf{y}}_{t}})},$

where the _reparameterization trick_ is used in sampling operation as the
Monte Carlo method is not differentiable which impediments towards BP
operation, sampling a ${\widehat{{\bf{H}}}_{\rm{l}}}$ from distribution
$\mathcal{CN}\left({{m_{{{\bf{H}}_{t}}}},{S_{{{\bf{H}}_{t}}}}}\right)$ is
equivalent to sampling a $\bf{h}$ from distribution $\mathcal{CN}(0,I)$ and
let
${\widehat{{\bf{H}}}_{\rm{l}}}={{m_{{{\bf{H}}_{t}}}}}+{{S_{{{\bf{H}}_{t}}}}}\times\bf{h}$,
and then the encode network can be trained as the sampling operation does not
need to participate in the gradient descent process.

## References

  * [1] F. Rusek et al., “Scaling Up MIMO: Opportunities and Challenges with Very Large Arrays,” in IEEE Signal Processing Magazine, vol. 30, no. 1, pp. 40-60, Jan. 2013, doi: 10.1109/MSP.2011.2178495.
  * [2] E. Basar, “Reconfigurable Intelligent Surface-Based Index Modulation: A New Beyond MIMO Paradigm for 6G,” in IEEE Transactions on Communications, vol. 68, no. 5, pp. 3187-3196, May 2020, doi: 10.1109/TCOMM.2020.2971486.
  * [3] H. Ye, G. Y. Li and B. Juang, “Power of Deep Learning for Channel Estimation and Signal Detection in OFDM Systems,” in IEEE Wireless Communications Letters, vol. 7, no. 1, pp. 114-117, Feb. 2018, doi: 10.1109/LWC.2017.2757490.
  * [4] H. He, C. Wen, S. Jin and G. Y. Li, “Deep Learning-Based Channel Estimation for Beamspace mmWave Massive MIMO Systems,” in IEEE Wireless Communications Letters, vol. 7, no. 5, pp. 852-855, Oct. 2018, doi: 10.1109/LWC.2018.2832128.
  * [5] B. A. Kenney, A. J. Majid, H. Moradi and B. Farhang-Boroujeny, “Frequency Domain Detection and Precoding for Massive MIMO with Single Carrier Modulation,” in IEEE Transactions on Wireless Communications, doi: 10.1109/TWC.2021.3119486.
  * [6] X. Yuan, C. Fan and Y. J. Zhang, “Fundamental Limits of Training-Based Uplink Multiuser MIMO Systems,” in IEEE Transactions on Wireless Communications, vol. 17, no. 11, pp. 7544-7558, Nov. 2018, doi: 10.1109/TWC.2018.2868100.
  * [7] E. Björnson, E. G. Larsson and M. Debbah, “Massive MIMO for Maximal Spectral Efficiency: How Many Users and Pilots Should Be Allocated?,” in IEEE Transactions on Wireless Communications, vol. 15, no. 2, pp. 1293-1308, Feb. 2016, doi: 10.1109/TWC.2015.2488634.
  * [8] B. Wang, F. Gao, S. Jin, H. Lin and G. Y. Li, “Spatial- and Frequency-Wideband Effects in Millimeter-Wave Massive MIMO Systems,” in IEEE Transactions on Signal Processing, vol. 66, no. 13, pp. 3393-3406, 1 July1, 2018, doi: 10.1109/TSP.2018.2831628.
  * [9] L. Chen and X. Yuan, “Blind Multiuser Detection in Massive MIMO Channels With Clustered Sparsity,” in IEEE Wireless Communications Letters, vol. 8, no. 4, pp. 1052-1055, Aug. 2019, doi: 10.1109/LWC.2019.2906168.
  * [10] C. Shin, R. W. Heath and E. J. Powers, “Blind Channel Estimation for MIMO-OFDM Systems,” in IEEE Transactions on Vehicular Technology, vol. 56, no. 2, pp. 670-685, March 2007, doi: 10.1109/TVT.2007.891429.
  * [11] J. Guo, C. -K. Wen and S. Jin, “CAnet: Uplink-aided Downlink Channel Acquisition in FDD Massive MIMO Using Deep Learning,” in IEEE Transactions on Communications, doi: 10.1109/TCOMM.2021.3120294.
  * [12] M. Soltani, V. Pourahmadi, A. Mirzaei, and H. Sheikhzadeh, “Deep learning-based channel estimation,” IEEE Communications Letters, 2019.
  * [13] H. He, C. Wen, S. Jin, and G. Y. Li, “Deep learning-based channel estimation for beamspace mmwave massive mimo systems,” IEEE Wireless Communications Letters, vol. 7, pp. 852–855, Oct 2018.
  * [14] P. Sabeti, A. Farhang, I. Macaluso, N. Marchetti and L. Doyle, “Blind Channel Estimation for Massive MIMO: A Deep Learning Assisted Approach,” ICC 2020 - 2020 IEEE International Conference on Communications (ICC), 2020, pp. 1-6, doi: 10.1109/ICC40277.2020.9149291.
  * [15] C. M. Bishop, “Pattern Recognition and Machine Learning (Information Science and Statistics).” Springer-Verlag New York, Inc., 2006.
  * [16] Z. Jiang, Y. Zheng, H. Tan, B. Tang, and H. Zhou, “Variational deep embedding: An unsupervised and generative approach to clustering,” 2016.
  * [17] D. D. Lin and T. J. Lim,“The Variational Inference Approach to Joint Data Detection and Phase Noise Estimation in OFDM,” in IEEE Transactions on Signal Processing, vol. 55, no. 5, pp. 1862-1874, May 2007, doi: 10.1109/TSP.2006.890915.
  * [18] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” Computer Ence, 2014.