{'fraction_non_alphanumeric': 0.045504644937805074, 'fraction_numerical': 0.00999842544481184, 'mean_word_length': 4.353704796425862, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sensing RISs: Enabling Dimension-Independent CSI Acquisition for Beamforming

Reconfigurable intelligent surfaces (RISs) are envisioned as a potentially transformative technology for future wireless communications. However, RISs' inability to process signals and the attendant increased channel dimension have brought new challenges to RIS-assisted systems, including significantly increased pilot overhead required for channel estimation. To address these problems, several prior contributions that enhance the hardware architecture of RISs or develop algorithms to exploit the channels' mathematical properties have been made, where the required pilot overhead is reduced to be proportional to the number of RIS elements. In this paper, we propose a dimension-independent channel state information (CSI) acquisition approach in which the required pilot overhead is independent of the number of RIS elements. Specifically, in contrast to traditional signal transmission methods, where signals from the base station (BS) and the users are transmitted in different time slots, we propose a novel method in which signals are transmitted from the BS and the user simultaneously during CSI acquisition. With this method, an electromagnetic interference random field (IRF) will be induced on the RIS, and we propose the structure of sensing RIS to capture its features. Moreover, we develop three algorithms for parameter estimation in this system, in which one of the proposed vM-EM algorithm is analyzed with the fixed-point perturbation method to obtain an asymptotic achievable bound. In addition, we also derive the Cram\'er-Rao lower bound (CRLB) and an asymptotic expression for characterizing the best possible performance of the proposed algorithms. Simulation results verify that our proposed signal transmission method and the corresponding algorithms can achieve dimension-independent CSI acquisition for beamforming.

Abstract-Reconfigurable intelligent surfaces (RISs) are envisioned as a potentially transformative technology for future wireless communications. However, RISs' inability to process signals and the attendant increased channel dimension have brought new challenges to RIS-assisted systems, including significantly increased pilot overhead required for channel estimation. To address these problems, several prior contributions that enhance the hardware architecture of RISs or develop algorithms to exploit the channels' mathematical properties have been made, where the required pilot overhead is reduced to be proportional to the number of RIS elements. In this paper, we propose a dimension-independent channel state information (CSI) acquisition approach in which the required pilot overhead is independent of the number of RIS elements. Specifically, in contrast to traditional signal transmission methods, where signals from the base station (BS) and the users are transmitted in different time slots, we propose a novel method in which signals are transmitted from the BS and the user simultaneously during CSI acquisition. With this method, an electromagnetic interference random field (IRF) will be induced on the RIS, and we propose the structure of sensing RIS to capture its features. Moreover, we develop three algorithms for parameter estimation in this system, in which one of the proposed vM-EM algorithm is analyzed with the fixedpoint perturbation method to obtain an asymptotic achievable bound. In addition, we also derive the Cramér-Rao lower bound (CRLB) and an asymptotic expression for characterizing the best possible performance of the proposed algorithms. Simulation results verify that our proposed signal transmission method and the corresponding algorithms can achieve dimension-independent CSI acquisition for beamforming.

I. INTRODUCTION
Reconfigurable intelligent surfaces (RISs) are considered to be a potentially important technology for future wireless communications. The characteristics of low cost and power consumption make RISs a promising solution for overcoming blockages, improving capacity, and reducing transmit power [1]- [3]. Specifically, an RIS is a large-scale array composed of passive elements, which can achieve significant beamforming gain by appropriately imposing phase shifts on the incident electromagnetic waves [4]. To achieve this beamforming gain, accurate channel state information (CSI) should be acquired beforehand, which makes channel estimation an essential problem for RIS-assisted communications [5].
Although channel estimation has been well investigated in conventional communication systems, the additional employment of RISs brings about two challenges [6]. Firstly, in contrast to traditional antenna array capable of transmitting, receiving, and processing the signals, RISs can only passively reflect the incident signals. Secondly, since the number of RIS elements is usually large, the dimensions of the channels increase sharply compared with conventional communication systems, which results in unacceptably high pilot overhead for channel estimation. These two main challenges result in the need for new channel estimation techniques for beamforming in RIS-assisted communications.

A. Prior Works
Generally, for RIS-assisted systems, channel estimation and beamforming are two separate procedures. Channel estimation is performed first, and then the obtained CSI is utilized for beamforming. The beamforming gain relies heavily on the channel estimation accuracy.
To tackle the challenges mentioned above in channel estimation for RIS-assisted communications, some solutions have been proposed, which can be generally divided into two categories. The first category modifies the hardware architecture of the RIS, which enables some signal processing capability. For example, by sparsely replacing some of the RIS elements with active sensors capable of baseband processing, the authors of [7] proposed a compressive sensing and deep learning-based channel estimation scheme with negligible pilot overhead. To further reduce the number of active sensors, in [8] the authors proposed an alternating direction method-based channel estimation procedure with a single radio-frequency (RF) chain, with the help of an extended analog combiner [9]. Since RF chains have been attached to RISs for channel sensing in these works, they can also be modified to perform signal relaying. In [10], a hybrid relay-reflecting architecture is considered, where a few passive RIS elements are replaced by active amplify-and-forward relay modules. These active modules are then capable of channel estimation.
The second category preserves the original hardware architecture of RISs, while conceiving new algorithms to exploit the specific features that RISs bring about. Exploiting the two-timescale channel property in the RIS-assisted system, the authors of [11] proposed a two-timescale channel estimation algorithm that reduces the pilot overhead in the time-averaged sense. In other works [5], [12], by exploiting the sparsity of channels in the angular domain, compressive sensing-based algorithms are developed with reduced pilot overhead. In [13], based on the shared reflected channels among multiple users, the authors proposed a three-phase channel estimation framework to further reduce the pilot overhead.
Note that, in order to achieve high beamforming gain, all the above channel estimation methods are designed to estimate the full channel matrix as accurately as possible. However, since the channel matrix is at least of size O(N ), where N denotes the number of RIS elements, the required pilot overhead for channel estimation is usually proportional to the number of RIS elements [7]- [14] in practical systems, which makes most of the existing channel estimation schemes dimensiondependent. This is unacceptable especially when the RIS is fabricated with a large number of elements (e.g. 1100 elements in [15], 2304 elements in [16], and 10240 elements in [17]). Therefore, the following question naturally arises: Does there exist a dimension-independent approach where the required pilot overhead is independent of the number of RIS elements?

B. Our Contributions
We point out that the main drawback of conventional channel estimation methods is that they only extract the mathematical features of channels, while neglecting their electromagnetic nature. This limitation has caused the problem of high pilot overhead required for channel estimation in RIS-assisted systems. Thus, we propose an interference random field (IRF)based approach, where the IRF induced on an RIS is utilized for channel estimation 1 . Specifically, the contributions of this paper can be summarized as follows.
• Inspired by optical interference where phase information can be obtained from the interference fringes, we reveal that phase information about the channel can also be gathered from the phenomenon of electromagnetic interference that occurs on RISs, which we name the IRF. To induce the IRF on an RIS, we propose a novel pilot transmission method called simultaneous rotational signaling, where signals are transmitted from the base station (BS) and the user to the RIS simultaneously, and the two signals bear a slight frequency difference during CSI acquisition. • To exploit this IRF for CSI acquisition, we employ a sensing RIS that integrates power detectors into the RIS elements to capture the features of the IRF. Each of the power detectors can acquire its phase information independently from the IRF, so that the required pilot overhead is independent of the number of RIS elements. We then develop discrete Fourier transform (DFT), maximum likelihood (ML), and von Mises-expectation maximization (vM-EM) phase estimation algorithms to extract the phase information from the power detector signals in order to perform beamforming. • By the fixed-point perturbation method, we analytically prove that the proposed vM-EM algorithm achieves an error decay of O(γ −1 ), which is the best expectable asymptotic precision that a phase estimator could attain. Furthermore, we derive the Cramér-Rao lower bound (CRLB) of the IRF phase estimation problem as well as its approximated asymptotic expression. Our numerical results verify that the developed vM-EM algorithm achieves the theoretical error decay rate, and that it is close to the CRLB. Our simulation results further demonstrate that the proposed signal transmission method and the corresponding algorithms can realize dimensionindependent CSI acquisition for beamforming, and can achieve near-optimal system spectral efficiency.

C. Organization and Notation
Organization: The rest of the paper is organized as follows. In Section II, we introduce the system model of RIS-aided communications, and review the existing separate procedures of beamforming and channel estimation. In Section III, we propose a novel signal transmission method for CSI acquisition in RIS-aided system, and reveal the IRF phenomenon under this method. In Section IV, we introduce the hardware architecture of a sensing RIS to exploit the IRF. Based on the sensing RIS, we propose three algorithms to realize channel estimation. In Section VI, we derive the CRLB and its asymptotic expression for channel estimation in the sensing RIS-assisted system. In Section VII, simulation results are provided for quantifying the performance of our proposed sensing RIS-based channel estimation as a dimension-independent approach. Finally, in Section VIII, we provide our conclusions followed by promising future research ideas.
Notation: C and R denote the set of complex and real numbers, respectively; i denotes the imaginary unit; {L} represents the set of integers {0, 1, · · · , L−1}; A −1 , A * , A T , and A H denote the inverse, conjugate, transpose, and conjugate transpose of matrix A, respectively; · 2 is the L 2 -norm of its argument function; · F denotes the Frobenius norm of its argument matrix; x, y := x H y denotes the inner product of complex vectors x, y; arg(x) and exp(x) denote the phase angle and exponential of the complex scalar x, respectively; |x| denotes the amplitude of a complex scalar x; diag(·) is the diagonal operation; CN µ, σ 2 represents the complex univariate Gaussian distribution with the mean µ and the variance σ 2 ; VM(µ, κ) denotes the von Mises distribution with circular mean µ ∈ [0, 2π] and centrality κ; N C χ 2 k (λ, σ 2 ) and N C χ k (λ) are the degree-k non-central chi-squared distribution with noncentrality parameter λ and variance parameter σ 2 , and the non-central chi distribution with non-centrality parameter λ, respectively; I L denotes the L × L identity matrix; x BB (t) denotes the baseband representation of a passband signal x(t); I ν (z) is the ν-th order modified Bessel function of the first kind.

II. SYSTEM MODEL
In this section, we will first specify the system model of the RIS-assisted multiple-input single-output (MISO) system. Then, we will clarify the power allocation between CSI acquisition and data transmission in II-A. Finally, traditional approaches for the corresponding beamforming design and channel estimation will be introduced in Subsection II-B and II-C, respectively.
Let us consider an RIS-aided MISO system, where an N -element RIS is employed for enhancing the transmission from an M -antenna BS to a single-antenna user. Assume furthermore that the phase-shift of each element of the RIS can be continuously and independently controlled with unit reflective gain [3]. Then, the precoding matrix of the RIS can be represented by where θ n (n ∈ {N }) denotes the phase-shift of the n-th RIS element, satisfying |θ n | = 1. Therefore, the signal received by the user can be written as where f ∈ C N ×1 and G ∈ C N ×M denote the channel spanning from the RIS to the user and the channel spanning from the BS to the RIS, respectively; w ∈ C M ×1 denotes the beamformer at the BS transmitter, with power constraint w 2 2 ≤ P max ; s denotes the normalized BS transmitted symbol satisfying E[ss * ] = 1; z ∼ CN 0, σ 2 z denotes the additive white Gaussian noise (AWGN) introduced at the receiver user. To focus on RIS beamforming, other possible links are neglected in this paper 2 .

A. Power Allocation
We assume a block-fading Rayleigh channel with B symbols in each block, i.e., the i.i.d. Rayleigh random channel matrices G and f are generated once for every consecutive B symbols. Within these B symbols, there are N p pilot symbols (i.e., the number of time slots for pilot signals) for channel estimation, and N d data symbols for data transmission, where B = N p + N d . The pilot symbols and the data symbols are subject to a total energy of E p and E d , respectively. Suppose the average transmit power of each block is P , then the energy values E p and E d should satisfy Thus, a smaller number of pilot symbols N p will lead to a higher average pilot signal-to-noise ratio (SNR) γ p = E p /(N p σ 2 z ) during CSI acquisition. This property ensures fair comparison among different channel estimation schemes, since equal energy E p is injected into the channel for the purpose of CSI acquisition.
By defining the total pilot energy E p and the total data enery E d , the actual transmit power values during the CSI 2 In fact, our proposed CSI acquisition procedure automatically works when there exist direct BS-user links, which will be explained in Sec. IV-D. acquisition phase and the data transmission phase are given by where P max is the average BS transmit power budget, P max is the average user transmit power budget.

B. Beamforming Design
Based on the received signal (1), we can formulate the signal-to-noise ratio (SNR) maximization problem for designing the BS beamformer w and precoding matrix Θ, i.e., A near-optimal solution to this joint active-passive beamforming problem can be solved by the following alternating optimization (AO) loop [3]: where norm(z), z ∈ C M ×1 denotes z/ z , and t is the iteration index. Note that this AO-based beamforming method depends on the knowledge of the full channel matrices G and f , thus the active beamformer w and the passive beamformer Θ have to be frequently re-optimized once the channel changes. However, a more often assumption in real-world RIS-aided systems is that the locations of the BS and the RIS are usually fixed, implying that the BS-RIS channel G enjoys a much longer channel coherence time compared with the RIS-user channel f [11]. Thus, in order to mitigate the computational cost, a suboptimal solution can be immediately obtained by designing w and Θ in a one-shot manner, i.e., where G = [g 1 , · · · , g N ] T . Note that this method is equivalent to optimizing the upper bound of the target function (4a):

C. Channel Estimation
Observe from (5) that, it suffices to know the cascaded channel H = diag(f * )G for beamforming design. To accurately acquire this cascaded channel, the user sends pilot signal x ∈ C to the BS with N p different RIS configurations. These different RIS configurations are usually designed to be the first N p columns of the DFT matrix F N , denoted as F N,Np . Due to the channel reciprocity, the received signal at the BS with the p-th RIS configuration Θ p = diag(θ p ) can be written as [18] y BS,p = H T θ p w s + n, p = 1, 2, · · · , N p , where w is the precoding scalar of the user satisfying |w | 2 ≤ P max , s is the normalized user transmitted symbol with E[s (s ) * ] = 1, and n ∼ CN 0 M , σ 2 z I M is the AWGN at the BS receiver. Equivalently, (8) can be expressed in the matrix form where Y BS = y BS,1 , · · · , y BS,Np and N = n 1 , · · · , n Np . Given the received signal Y BS and, without loss of generality, assuming w s = P max , the channel estimation problem in RIS-assisted system can be solved by the MMSE estimator [14] asĤ However, the evaluation of the MMSE estimator requires exact knowledge of the prior p.d.f. of the channel p(H), which is difficult to be written in an explicit form. This is because the entries of the cascaded channel matrix H ij is the product of two Gaussian-distributed random variables f * i ∼ CN (0, σ 2 f ) and G ij ∼ CN (0, σ 2 g ), which is generally not a Gaussian random variable.
An alternative way to circumvent this difficulty is to constrain the estimatorĤ in (10) to be a linear transform of Y BS . This linear constraint produces a linearly approximated MMSE estimator, which is called the linear MMSE (LMMSE) estimator. By assuming σ 2 h = σ 2 f σ 2 g , the LMMSE estimator can be obtained as which provides a feasible solution to the channel estimation problem in the RIS-assisted MISO system. In this paper, we name it the linear minimum mean squared error for H (LMMSE-H) channel estimation method.
Furthermore, by assuming known BS-RIS link G, the channel estimation problem can be significantly simplified to the estimation problem of f * . By exchanging the position of θ and f * in (8), we obtain the equivalent uplink channel estimation model as which immediately leads to the following MMSE 3 estimator for f : whereỹ = vec(Y BS ) ∈ C NpM ×1 , and Since this channel estimation scheme assumes known G and estimates the unknown f by an MMSE estimator, we refer to it as the MMSE-f scheme.
It is worth noting that, both the LMMSE-G and the MMSEf methods consume a linear number of pilot slots in N . For the LMMSE-G scheme, In order to determine the RIS-user channel f ∈ C N ×1 , at least N p = N/M = Θ(N ) pilots are required, even with the assumption that the slow-varying [11] BS-RIS channel G is known in advance. In fact, without exploiting the inherent structure of f , it seems impossible to obtain the N channel coefficients with pilot overhead smaller than Θ(N ). However, in the following section, we will show that it is possible to realize dimension-independent CSI acquisition by exploiting a common physical phenomenon in electrodynamics.

III. INTERFERENCE RANDOM FIELD
In this section, the signal model for the IRF will be first introduced in Subsection III-A. Then, in order to acquire the CSI, the simultaneous rotational signaling method will be proposed in Subsection III-B.

A. Signal Model for IRF
Physical intuition is vital in designing RIS-aided systems [19], and interference is a fundamental physical phenomenon that appears in all kinds of waves. The most wellknown example is the double-slit optical interference [20], as is shown in Fig. 1 (a), where the interference fringes created during the wave superposition reveal the phase difference of the two optical paths. Similarly, the same interference phenomenon occurs when two RF electromagnetic waves meet together on the RIS, as is shown in Fig. 1 (b), and the IRF created by this interference reveals the CSI.
Suppose two signals impinge upon the RIS simultaneously, creating an interference field at the n-th RIS element. Denote the symbol transmitted from the BS by s, and the symbol transmitted from the user by s . In order to perform channel estimation by exploiting the IRF, we need to probe the interference fringes, i.e., the power of the IRF. To clarify the power problem associated with an electromagnetic signal, the relationship between physical signal power and baseband equivalent signal power is stressed in the following formulas. If we represent both of the BS-RIS and user-RIS signals that appear on each RIS element by their baseband equivalent signals E BB (t), then the corresponding physical passband electric field induced on each RIS element is characterized by where ω c = 2πf c denotes the carrier frequency, and the coefficient √ 2 ensures that the passband signal power is equal to the baseband signal power, i.e., E BB 2 2 = E 2 2 . Thus, taking the square of the baseband signal is equivalent to calculating the power of the physical electromagnetic signal.  Now we consider the IRF case, where a superposition of the BS-RIS signal and the user-RIS signal is considered. Following the notations in Section II and due to the linearity of (15), we obtain the noisy IRF signal by adding up the two impinging baseband signals where v(t) ∼ CN (0, σ 2 v ) is the electromagnetic noise signal in its baseband representation, w and w are the beamformer at the BS side and the user side satisfying w 2 2 ≤ P max and |w | 2 ≤ P max , respectively, s is the symbol transmitted from the BS to RIS, and s = e iψ(t) is the time-varying transmitted symbol from the user to the RIS. Since the interference fringes are not sensitive to a global phase change, we can safely assume that s = 1, and then the relative phase between the user and the BS can be fully characterized by a time-varying phase function ψ(t).
Furthermore, by defining α = g T n w , β = |f * n w |, and the phase difference between the BS-RIS link 4 and the RIS-user channel as ϕ = arg (f * n w ) − arg g T n w , then the power of 4 Since it is usually difficult to obtain the optimal BS beamformer w, during the IRF-based CSI acquisition procedure, we fix w to a vector that guides most of the signal energy to the RIS aperture according to (6). As a result, g T n w can be treated as the equivalent BS-link from the BS to the n-th RIS element.
the IRF can be written as where ζ is the noise introduced by digital signal processing after measuring the IRF power, A is the amplification factor of the power sensor, and the equivalent electromagnetic noise v (t) = e i·arg(g T n w) v * (t) follows the same distribution as v(t) for any time t. The name IRF comes from the randomness of the unknown ϕ, the unknown noise realization, and the interferential nature of the field.
The phase difference ϕ defined here will play an important role in our following channel estimation, since it carries enough CSI for beamforming, which will be justified in the next section. To estimate the phase difference ϕ, we need enough observations from the sensor-detected power signal P (t). For simplicity and without loss of generality, we can assume ψ(t) = 2π Ts t, where T s is the symbol period. Furthermore, assume that L observations P [l] = P (t l ) located at equally-spaced instants t l = l L T s for all l ∈ {L} are used for phase estimation ϕ.

B. Simultaneous Rotational Signaling and Interference Detection
However, at first glance it seems impractical to use only power signals to obtain enough information for RIS beamforming, in the absence of phase information. Generally, signal processing with only amplitude signals are called noncoherent detection. The target of RIS beamforming is to strengthen the signal at the user, hence it requires coherent signal combining at the user antenna. As a result, for noncoherent detection, it seems that the lack of phase information makes beamforming impossible. Fortunately, the physical phenomenon interference makes it possible to convert phase difference into intensity difference, allowing noncoherent devices to perform coherent detection.
To clarify this idea, we take only one RIS element into consideration. As is introduced in the previous subsection, if we allow the BS and the user to transmit electromagnetic waves simultaneously, then the interference phenomenon will occur on each RIS element. However, equal carrier frequencies of the BS and the user create stable interference fringes [20] on the RIS, which do not carry information about the channel. In contrast, a rotational symbol s = e iψ(t) at the user, which is equivalent to a slightly higher carrier frequency, enables the IRF power to vary over time. Thus, the desired CSI can be drawn from the varying IRF power signal received by sensing RIS. Fig. 2 shows the simultaneous rotational signaling procedure and the IRF waveform.
To further analyze the interference, we focus on the received IRF signal P (t) of the n-th RIS element. Note that the power of interference field measured at instant t, expressed by (17) Fig. 2. Simultaneous rotational signaling scheme. The IRF power is the instantaneous total power of the BS-RIS signal and the user-RIS signal. The composite power waveform that appears on RIS enables our algorithm to obtain the desired CSI. The parameters α and β are estimated from the signals Pα(t), P β (t) measured during the first two time slots, and ϕ is estimated by IRF signal P (t) that occurs during the third time slot.
H-MIMO noise initial phase initial phase Fig. 3. Phasor representation of the IRF on each RIS element. Red vector represents the complex BS-RIS signal g T n w; Blue vector represents the complex user-RIS signal f * n w ; α and β denote the amplitude of these signals respectively. In our signaling method, the BS transmits a fixed symbol s = 1, while the user transmits a rotating symbol s = e iψ(t) , causing the output of the power sensor P (t) to vary in a waveform that is similar to a sine curve. exhibits a sinusoidal waveform in the time domain, as is shown in Fig. 3. The initial phase (t = 0) of this sinusoidal waveform is uniquely determined by ϕ. Thus, if we have access to the power received by the n-th element at successive instants t l , it is possible to retrieve the phase difference ϕ. However, according to (6), it is the phase sum of the RIS-user channel and the BS-RIS channel, i.e., arg(g T n w) + arg(f * n w ), that determines the optimal phase-shift of the n-th element. Given the phase difference ϕ = arg (f * n w ) − arg g T n w , we assume the phase arg(g T n w) to be known in our algorithms in order to acquire the phase sum. The underlying reason for this assumption is that, in contrast to the fast time-varying RIS-user link, the BS-RIS link often exhibits a quasi-static property, thus two-timescale methods in [11] can be applied, in which the quasi-static BS-RIS channel G is estimated in a longer timescale, while the frequently-varying RIS-user channel f is estimated in a shorter timescale. Since the BS-RIS link G is usually stable, it can be estimated only once for many following data frames. Thus, in our simulations for the IRFbased methods, the phase arg(g T n w) is assumed to be known.

IV. SENSING RIS-BASED CHANNEL ESTIMATION
In this section, we will first briefly introduce the hardware architecture required for detecting the IRF in Subsection IV-A. Then, the core concept of IRF channel estimation and beamforming procedures will be introduced in Subsection IV-B. After that, Subsection IV-C will include the analysis on pilot overhead and computational complexity of the proposed IRFbased CSI acquisition methods. Finally, the extension to the co-existence case of direct-reflective links are presented in Subsection IV-D.

A. Hardware Architecture of Sensing RIS
To obtain the amplitude of the induced IRF, inspired by the sensing metasurface [21], we propose a hardware architecture named sensing RIS. In contrast to traditional RIS architecture shown in Fig. 4 (a), where each RIS element includes a phaseshift circuit and a patch antenna, each sensing RIS element additionally integrates a power sensor which is responsible for detecting the amplitude of the IRF [22], as shown in Fig. 4 (b). Since all of the RIS elements can sense and adjust the phase independently from one another, low-cost microcontroller units (MCUs) can be attached locally to each RIS element to allow parallel computation of the optimal phases. However, for extremely large-scale RIS systems, sparse sensing and controlling may be preferred to reduce cost and hardware complexity.

B. IRF Channel Estimation and Beamforming
In this subsection, we will thoroughly introduce the IRFbased algorithms for adjusting the phases of the RIS elements. The IRF channel estimation and beamforming procedure can be divided into three steps: 1) Simultaneous rotational signaling and power sensing; 2) Phase estimation based on power data; 3) Integrate phase information and other CSI to perform beamforming. In step 1), three power signals are recorded: P α (t), P β (t), and P (t) (see Fig. 2 for details). The IRF appears during the third signal P (t). In order to create an IRF, simultaneous rotational signaling must be performed. In fact, the signaling requirements can be realized by transmitting a symbol s = 1 on the zeroth subcarrier at the BS, while transmitting a symbol s = 1 on the k-th subcarrier at the user. As is mentioned above, if we denote the OFDM symbol period by T s , then the equivalent IRF angular frequency ω is given by ω = 2πk/T s , which means the noiseless power signal exhibits exactly k sinusoidal periods during an OFDM symbol. Without loss of generality, in this paper, we always assume k = 1. In fact, for multi-antenna MIMO systems, rotating pilot symbols at different frequencies k can be simultaneously transmitted by different antennas to enable CSI acquisition. Then, multiantenna CSI acquisition can be fulfilled by performing Fourier analysis on the composite IRF power signal P (t), and extracting the phase angles for each antenna at distinct frequencies.
Similarly, this pseudo-frequency division idea can be extended to multi-user MIMO systems, where different users can also be distinguished by different IRF frequencies. Fortunately, since Algorithm 1 Near-optimal RIS Beamforming by IRF Input: Number of RIS elements N , IRF power signals detected on each RIS element P α (t), P β (t) and P (t).

4:
Estimate ψ n = arg(g T n w) from known locations of BS and RIS 5: different users are usually separated in the angular domain from the RIS's perspective, these users can be further identified by performing joint detection and estimation across all the RIS sensors. The joint CSI estimation is beyond the scope of this paper, and is left for our future work.
The second step 2) is called phase estimation. This step is designed to extract the phase difference ϕ between the BS-RIS and RIS-user channel. This is the key step of dimensionindependent IRF channel estimation and beamforming, since this step can be done independently among all the RIS elements.
The third step 3) is to utilize the phase information provided by step 2) and calculate the near-optimal phase shifts for each RIS element in order to perform beamforming. The additional CSI refers to the phase information arg(g T n w).
Step 2) and step 3) are executed by processors, so we collect these two steps together into pseudo codes shown in Algorithm 1: Note that in each iteration n of Algorithm 1, the data P (t) is independent of other iterations. Thus, we can perform the calculations in parallel for each phase-shift θ n . For example, we can install an MCU for each of the elements on the RIS. Each MCU is only responsible for gathering the data from its own power sensor and adjusting the phase-shift of its own RIS element. Since all MCUs can work in parallel, the computational time is independent of the number of RIS elements, resulting in an O(1) time complexity.

C. Pilot Overhead and Computational Complexity
In our proposed IRF-based CSI acquisition method, the pilot overhead is fixed to O(1), which is independent of the RIS dimension N . The reason is that, no matter how many RIS elements are employed, the IRF appears on them simultaneously. Thus, both the channel estimation and beamforming can be fulfilled within only three pilot symbols, as is depicted in Fig. 2. Thus, the pilot overhead is independent of the RIS dimension N . To the best of our knowledge, this dimension-independent property is unprecedented, if we cannot endure the cost of attaching a dedicated RF chain to every RIS element. Some hybrid solutions do exist, such as connecting all the RIS elements to several RF chains with analog combiners [8], [23] to enable explicit channel estimation at RIS. Unfortunately, the cost of analog combiners and RF chains are usually much higher than power sensors. Table I compares the pilot overhead of our proposed IRF method with other different CSI acquisition methods. All the CSI acquisition methods that appear in Table I assume a narrowband system, so the pilot overhead is equal to the number of pilot time slots. K denotes the number of users, S is the sparsity of the channel assumption, M is the number of antennas at BS, and α is the ratio of the large-timescale channel coherence time and the small-timescale channel coherence time [11]. Except for CS-based methods [5], all the channel estimation methods [11], [13], [24] require the pilot overhead to be linearly dependent on N , while our IRF method needs a constant number of exactly 3 pilot slots per user, regardless of the number of RIS elements.
The computational complexity of processing the obtained power signals depends on the hardware implementation. The O(1) time complexity only holds when one MCU is installed for each RIS element. If all the RIS elements are collectively controlled by one processor, the computational time would be O(N ), but the pilot overhead still remains O(1). These conclusions of complexity can be easily extended to multi-RIS schemes, where multiple RISs are employed to serve a single user at the same time. The independent nature of IRF methods allows the RISs to work without the need to exchange data with the BS or other RISs. This property makes it much  [24] N + 1 RIS General Multi-user [13] 1 + N/K + (K−1)N M /K RIS General CS [5] O(S log N ) RIS Sparse Two-timescale [11] 2(N +1) αK 3 Sensing RIS Quasi-static BS-RIS easier to integrate a new RIS into an existing communication system, which greatly enhances the extendibility of the system.

D. Extension to Systems with Direct BS-User Link
Since the direct link does not affect the creation of the IRF, the RIS phase-shift matrix Θ can be obtained as if there were no such a direct link, up to an undetermined global phase e iφ . Note that this IRF-based procedure for obtaining Θ consumes only 3 time slots. After this procedure, two additional uplink training symbols can be transmitted from the user. The uplink channel model is where h d ∈ C M ×1 is the direct BS-user link. In the first symbol period φ is set to 0, and in the second φ = π. Thus, the direct BS-user link h d and the reflective link h RIS can be recovered byĥ After obtaining the estimators of these two links, the global phase φ should be tuned to maximize the total channel energy h d + e iφ h RIS 2 , and this is done by settinĝ Though the final solution exp(iφ)Θ for the RIS phase-shift matrix is suboptimal in general, it only consumes two additional time slots, resulting in 5 time slots in total for configuring N reflective elements of the RIS. This pilot overhead is still dimension-independent.

V. PHASE ESTIMATION ALGORITHMS
In this section, three phase estimation algorithms will be proposed in Subsection V-A, V-B, and V-C respectively, which constitute the core of the IRF channel estimation and beamforming algorithm.

A. DFT method
The key challenge of the IRF channel estimation and beamforming is the phase estimation step, i.e., how to obtain the phase difference ϕ. Since the interferential power P (t) exhibits a sinusoidal waveform, Fourier transforms can be applied to extract its phase. Apply L-point discrete Fourier transform (DFT) to the discrete-time observed sensor detection signals P [0], · · · , P [L − 1], and we have Specifically, we have the complex amplitude of the first harmonic p [1] as Then, the phase ϕ can be estimated aŝ Note that this DFT method simply ignores the non-Gaussian noise. For Gaussian noise, the DFT method is optimal. However, in fact, the noise in (17) contains a squared term of Gaussian noise, resulting in the noise being non-Gaussian. Thus, we further conceive an ML method, as described in the following Subsection V-B.

B. Newton-ML method
Suppose the noise field v (t) = v R (t)+iv I (t) ∼ CN (0, σ 2 v ), and the noise of the power sensor is ζ ∼ N (0, σ 2 ζ ). Without loss of generality, we can assume that the sensor noise power σ 2 ζ is much weaker than the electromagnetic noise field v (t). Thus, we assume σ 2 ζ = 0 in the following discussion. As a result, the distribution of P (t) = A |E BB,IRF (t)| 2 is a noncentral chi-squared distribution N C χ 2 2 (A(µ 2 R + µ 2 I ), Aσ 2 v /2) with degrees of freedom k = 2, and mean values µ R , µ I given by µ R = α + β cos(ψ(t) + ϕ), µ I = β sin(ψ(t) + ϕ). (24) Thus, according to (17), the output signal of the power sensor is given by Let us define the noncentral parameter λ(t) as then, according to the definition of the N C χ 2 2 , the p.d.f. of P (t) is given by the zeroth-order modified Bessel function of the first kind (27) Then, the log likelihood function of ϕ based on the observations P [l] can be represented by where λ l := λ(t l ), and the derivative of (28) is where the function R(z) is defined as R(z) = I 1 (z)/I 0 (z). Since the derivative of the function R(z) satisfies the property [25] R the second derivative of the likelihood function L can be expressed as where z l = P [l]λ(t l )/(Aσ 2 v /2). Then, we can perform the Newton iteration to obtainφ, by iteratively using the updating formulaφ where L (φ (k) ) and L (φ (k) ) are given by (29) and (31) respectively. Note that during the calculation of (32), A, σ v , α, β and the received signal P [l] are all assumed to be known. Thus, in fact, (29) and (31) are functions of a single variable ϕ.

C. von Mises-EM method
The Newton-ML algorithm, if convergent, is asymptotically optimal [26]. However, the computation of the Newton-ML estimator is quite complicated due to the intensive calculation of modified Bessel functions. Now we introduce an iterative method for estimating ϕ without any computation of such special functions. Our method is based on the von Mises distributions [27].
The von Mises-EM algorithm is based on the Bayesian inference of von Mises distributions [28]. The von Mises distribution VM(µ, κ) is a two-parameter distribution on [0, 2π], with the probability density function given by where µ ∈ [0, 2π] and κ > 0 being the cyclic location parameter and the concentration parameter. Note that the von Mises distribution is a distribution on a circle, thus it acts as a perfect prior distribution of a phase estimation problem. More fortunately, the von Mises distribution is also closely related to the complex Gaussian distribution, thus implying the possibility of designing an iterative EM algorithm [26] based on the interactions between the von Mises distribution and the complex Gaussian noise distribution. The following two lemmas: Lemma 1 and Lemma 2 reveal these interactions.
Proof: See Appendix A.
In this paper, we also use VM(κe iµ ) to denote the von Mises distribution VM(µ, κ). This representation provides convenience for the calculation of the posterior distribution of the von Mises distribution in Bayesian inference.

Lemma 2 (Circular CN posterior is VM)
Suppose z ∼ CN (z 0 , σ 2 ), where z 0 ∈ C, and a positive radius r > 0. Then the posterior distribution of angle θ = arg(z), constrained on a circle |z| = r obeys the von Mises distribution p(θ| |z| = r) ∼ VM arg(z 0 ), Proof: Replacing z by re iθ give rise to the conclusion immediately.
After inferring the latent variablesθ l,ML , we can update the estimation of ϕ using Bayesian rule in Lemma 1 where the coefficient β/(σ 2 v /2) comes from scaling the phasor in Fig. 3 by a factor β −1 and applying Lemma 1. Performing E-step with (35) and M-step with (36) alternately, then the Algorithm 2 von Mises-EM phase estimation (vM-EM algorithm) Input: Incident wave intensity α, β; sensor data P [l]; amplification factor A and noise variance σ 2 v ; predefined phase shifts ψ l = ωt l .

VI. PERFORMANCE ANALYSIS
In this section, we first provide the achievability proof the proposed vM-EM phase estimation algorithm. Then, we derive various CRLB expressions for the phase estimation problems.
A. Asymptotic Achievability Bound of the vM-EM algorithm Theorem 1 (Fixed-point perturbation bound) Assume α > β > 0 are fixed. Then, there exists a sufficiently large integer L 0 , such that for any L ≥ L 0 , if the proposed vM-EM algorithm converges, the returned estimatorφ achieves an MSE performance of E[(φ − ϕ) 2 ] = O(γ −1 ), where ϕ is the true parameter, andγ is the interferential SNR defined asγ = (α 2 + β 2 )/σ 2 v . Before proving this theorem, we first introduce the following Lemma 3 that characterizes how fast the vM-EM estimator ϕ in Algorithm 2 varies with the random noise vector v.
Lemma 3 (The speed of the fixed-point) Suppose that for the input data sequence s = P /A = (s 0 , · · · , s L−1 ) T , the vM-EM algorithm converges to the estimatorφ. Then, the estimatorφ =φ(v) can be viewed as a function of the noise v, and the Wirtinger derivative [29] ofφ w.r.t. v satisfies where x = (x 0 , · · · , x L−1 ) T is defined as and Proof: The proof is provided in Appendix E.
The aim of Lemma 3 is to establish the connection between the vM-EM estimatorφ and the noise vector v. In the following, we denote the estimation error as Specifically, it is justified in Lemma 3 that, the squared error |∆ϕ| 2 of the vM-EM estimator is intrinsically bounded by the noise energy v 2 , since the differential d∆ϕ of a general complex-valued function ∆ϕ : C → C can be written as where for real-valued functionφ(v) : C L → [0, 2π] ⊂ R this differential relation reduces to Thus, an integration inequality will hold to upper-bound the MSE of the estimatorφ, i.e., Before explaining this idea in detail, we first introduce some interesting results about some intermediate variables.

(46)
Proof Sketch. The integral expressions for H ∞ and G ∞ as well as the asymptotic residuals O(L −1 ) can be obtained by applying the numerical trapezoidal integration formula [30] to the definition of the sequence H L and G L . The integrals on [0, 2π] can be evaluated by the following replacements [31]: and the applying the Residue Theorem to all poles inside the closed curve |z| = 1.
Remark 1 These limiting expressions will play an important rule in obtaining the upper bound of the estimation error |∆ϕ|, since this invariance can asymptotically eliminate the dependency of the Wirtinger derivative (Lemma 3) on the unknown estimatorφ.
Remark 2 In fact, we can also prove that by the same numerical integration approximation technique and the Residue Theorem [31].
Lemma 5 (ODE bound) Suppose for the input data sequence s = P /A = (s 0 , · · · , s L−1 ) T , the vM-EM algorithm converges, and the estimation error is denoted by ∆ϕ =φ − ϕ. Then, for any sufficiently large L ≥ L 0 , there exists some positive δ = δ(L) > 0 and C = C(L) > 0, such that for any noise vector v satisfying v / √ L ≤ δ, the estimation error is upper-bounded by Proof: The proof is provided in Appendix F. From the above Lemma 5, we can directly prove Theorem 1. Choose L ≥ L 0 , and then let us compute the MMSE of the vM-EM estimator, i.e., where (a) comes from applying the Markov inequality, and r > 0 can be arbitrarily chosen. Particularly, by setting r = 2, we obtain which completes the proof of Theorem 1.

Remark 3
The conclusion of Theorem 1 guarantees that as γ → ∞, the MSE of the vM-EM estimator decays at a rate of at least (γ) −1 . This conclusion is verified in the following numerical simulation, where it is shown that the MSE curve of the proposed vM-EM algorithm has an asymptotic slope of −1 in the logarithmic coordinate.

B. Expressions for the CRLB
In Section IV, we have introduced three phase estimation methods to solve the probabilistic parameter estimation problem. To analyze the performance limit of the proposed schemes, we derive the CRLB [26] of the estimation in Lemma 6.
Proof: See Appendix B. From Lemma 6, we note that the CRLB of the phase estimation problem is completely determined by the two parameters K andγ, which are the interferential contrast and the average interferential SNR, respectively. However, this exact CRLB expression is difficult to calculate due to the sophisticated evaluation of the integral in (53). Thus, we will provide a simpler approximated CRLB in the following Theorem 2.

Theorem 2 (Asymptotic CRLB)
Suppose the probabilistic model is specified by (25). Then the approximated CRLB of this L-point phase estimation problem is given by where (x) + represents x1 {x≥0} , and the definition of parameters are the same as in Lemma 6. The asymptotic approximation functionĝ(γ) is defined aŝ π γ e −γ/2 ((1 + 1/γ)I 0 (γ/2) + I 1 (γ/2)) . (55) Proof: See Appendix C. The asymptotic expression of the CRLB still relies only on two parameters K andγ. But the calculation is much simpler compared to the exact expression. The discussions of the properties of the exact and asymptotic expressions of CRLB are continued in Subsection VI-C.

C. Properties of the CRLB
From the above analysis of CRLB, we can observe two properties of the phase estimation problem: 1) The CRLB relies almost only on the physical parameters K andγ.
2) The CRLB is insensitive to the value ϕ. The first point can be concluded from the expressions, and the second point comes from the symmetry of the sin 2 (·) function when L is large enough. We have plotted 1/CRLB(ϕ) as a two-variable function in Fig. 5. The best prediction accuracy occurs when |K| = 1, i.e., the RIS-received signal from BS is as strong as that from the user. Also, the larger value of the average interferential SNRγ also contributes to a more accurate phase estimation.
We also perform error analysis for the use of asymptotic expansion technique, which is stated in the following Theorem 3.

Theorem 3 (Asymptotic Optimality of CRLB)
The relative error r of 1/γ −ĝ(γ), as an approximation of 1/γ − g(γ), decreases at a rate that is inverse proportional to the interferential SNR γ, i.e., Proof: See Appendix D. From Theorem 3, we can conclude that our derived expression of approximated CRLB is asymptotically optimal.

VII. SIMULATION RESULTS
In this section, we will present our simulation results. In Subsection VII-A, we will show the performance comparison of the phase estimation algorithms. In Subsection VII-B, we will show the achievable rate comparison of our IRF method and other RIS-aided channel estimation and beamforming algorithms.

A. Phase Estimation Algorithms
In Section IV, we have already introduced three phase estimation algorithms: DFT, Newton-ML and vM-EM. Now we compare the performance of these algorithms together with the CRLB which has been derived in Section VI. In all the simulations for phase estimation algorithms, the amplification factor A is fixed to be 1, the number of power samples 5 is fixed to L = 2 6 = 64, and σ ζ = 0.05; The interferential SNRγ = 20 in Fig. 6. The CRLB and CRLB-approx curves are calculated from (52) and (54), respectively. For Newton- ML algorithm, the Newton iteration is performed 4 times according to (32). For vM-EM algorithm, the iteration number is also fixed to 4. The true value of random variable ϕ under estimation is drawn from a uniform distribution on [0, 2π]. In the simulations in Fig. 7, all the simulation parameters, except for K andγ, are the same with that of Fig. 6. It can be concluded from Fig. 6 and Fig. 7 that, the vM-EM algorithm has comparable performance with the Newton-ML algorithm, but the computational cost is significantly lower, since Algorithm 2 does not require the evaluation of the complicated modified Bessel functions.
It is shown in Fig. 7 that the vM-EM algorithm achieves an asymptotic error decay of O(γ −1 ), which coincides with the statement in Theorem 1. Furthermore, both the Newton-ML and vM-EM algorithms are close to the CRLB, and both of them outperform the simple DFT algorithm. The CRLB approximation (54) is also satisfactory under a wide range of K andγ. However, the DFT estimator in Fig. 7 exhibits a near-constant performance gap toward the CRLB in the high-SNR region, while the Newton-ML and vM-EM estimators are able to bridge this gap and finally approaches the CRLB asγ → ∞.

B. Spectral Efficiency with IRF-based CSI Acquisition
All our simulation data are acquired under f c = 3.5 GHz, P max = 300 mW, n 0 = −174 dBm/Hz, subcarrier bandwidth BW = 180 kHz, thermal noise at the receiver σ 2 = BW × n 0 , σ 2 v = 100 MHz × n 0 F p , where F p = 10 dB is the noise factor of the power sensor. Both the BS and the RIS are equipped with λ/2-spaced uniform planar arrays (UPAs), while the user has a single antenna. The BS is located at a distance from the RIS that is drawn from a uniform distribution between 20 m and 100 m, and the user appears uniformly within a distance ranging from 10 m to 100 m. The size of the RIS is set to be N = 20×10, and that of the BS antenna is 4×2. The detailed simulation settings about the pilot energy and data energy (3) are listed in Table II.
Moreover, for the IRF methods, the BS beamformer w is chosen to be the right singular vector that corresponds to the most significant singular value of the BS-RIS channel G, which maximizes the total signal energy received by the RIS from the BS, i.e., and the user's transmit scaler is set to the maximum allowed power w = P max for all the CSI acquisition schemes to ensure a high SNR. For the MMSE/LMMSE algorithms, after the channel estimation procedure, the beamforming step is fulfilled by alternately updating the BS beamformer w and the RIS phase-shift matrix Θ = diag(θ) according to (5).
The initial values w (0) and θ (0) are generated randomly, and are subject to their corresponding power constraint and unit-modulus constraint, respectively. In our simulations, the iteration is performed for at most T = 10 times. The proposed IRF channel estimation and beamforming method are compared with the conventional MMSE-f , LMMSE-H, and other idealized settings in Fig. 8. In the simulation of the IRF method, the vM-EM phase estimation algorithm is utilized with the assumption that the phases of the BS-RIS line-of-sight (LoS) path are known. The "Oracle" (black dashed line) assumes perfectly-known CSI with iterative beamforming (5). In the LMMSE-H method (pink dashed line), we first utilize (11) to obtain estimates for the cascaded channel [5], [14], and then perform iterative RIS beamforming (5) based on the estimated channel. The MMSEf method (blue line) is calculated by (13) before applying the beamforming algorithm (5). The matched filter (MF) scheme (cyan line) is obtained by assuming σ 2 z → ∞ in (11). In addition, we also consider the randomly-phased RIS scheme as a benchmark. Fig. 8 shows the performance curve of the ergodic spectral efficiency C PSI against the BS transmitted power. The ergodic spectral efficiency is defined as which is normalized by the pilot overhead, and the subscript PSI represents perfect side information [32]. The achievable user SNR is defined as From Fig. 8, we can conclude that the proposed IRF method with vM-EM phase estimation is nearly optimal even with extremely low pilot overhead. The proposed IRF method is even close to the oracle scheme, and achieves satisfactory performance throughout a wide range of the BS transmitted power.

VIII. CONCLUSION
In this paper, we have introduced a dimension-independent CSI acquisition method for sensing RIS-assisted MISO wireless communication systems. Combined with our proposed vM-EM phase estimation algorithm, the pilot overhead of our CSI acquisition method is made independent of the number of RIS elements with low computational cost, enabling the implementation of extremely large-scale RISs to achieve significant beamforming gain. Theoretical analysis have demonstrated the asymptotic optimality of the proposed vM-EM algorithm, which is further supported by the CRLB analysis. Simulation results have also verified the near-optimality of our vM-EM algorithm. Furthermore, due to the elementwise independent property of our IRF-based CSI acquisition method, the nearfield effect cannot corrupt the precision of the CSI. Also, due to the simultaneous signaling protocol, the "multiplicative fading" effect of RIS [33], [34] during channel estimation is automatically avoided. Thus, the IRF method has promising applications to high-frequency large-scale systems.
For future work, the spatial interferential fringes on the sensing RIS may be exploited to recover the CSI at higher precision, and the data obtained by the power sensors may be utilized to perform joint channel estimation and beamforming with sparse assumptions on the channel. In addition, different interferential frequencies can be assigned to different users to perform multi-user CSI acquisition simultaneously, but the waveforms should be re-designed to avoid interference among users. Furthermore, when equipped with a sensing RIS, traditional MIMO systems can also benefit from the additional CSI provided by the phase estimation methods based on the IRF.

APPENDIX A PROOF OF LEMMA 1
Denote z = z r + iz i ; then the posterior density p(θ|z) ∝ p(θ)p(z|θ) can be expressed as Since the density of the von Mises distribution VM(µ, κ) can also be expressed as p(θ) ∝ exp Re e iθ (κe iµ ) * , we can also parameterize the von Mises distribution by a single complex parameter κe iµ . Thus, the above p(θ|z) is a von Mises density with parameter κ e iµ , satisfying κ e iµ = κe iµ + 2z/σ 2 . This completes the proof.

APPENDIX B PROOF OF LEMMA 6
According to the definition of the CRLB, taking the negative expectation of (31) yields the reciprocal CRLB of the estimators for ϕ. Note that in (31) The expectation E R(z l ) P [l] can be evaluated by calculating the derivative w.r.t λ l on both sides of the identity According to the definition of z l , the expectation E [R(z l )P [l]/z l ] is of the same form as the expectation E R(z l ) P [l] . Thus, the expectation result is As for the third kind of expectation E (1 − R 2 (z l )) P [l] λ l , we first define Let x = λ l t, and write down the integral expression of (64), we obtain Combining the above equations of expectations (61), (62), (63) and (65), we obtain where step (a) comes from substituting these three types of expectation into (31), and step (b) comes from replacing the trickiest expectation by the definition of the function g(·) in (53).
Note that the precise value of the CRLB can be determined by the exact single-variable g function (53), whose variable γ l = λ l /a is the interferential SNR. After further derivations, it can be observed from the expressions that the CRLB only relies on two intrinsic physical parameters: the interferential SNR γ l , and the interferential contrast K [20], which coincides with the physical intuition. Imitating the notations in optics, we define the parameter K to be and K automatically satisfies −1 ≤ K ≤ 1. Define the average interferential SNRγ to be the arithmetic average of γ l , 0 ≤ l < L:γ Thus, the CRLB can be expressed as where the values γ l =γ (1 + K cos(ψ l + ϕ)) are jointly determined by both the average interferential SNRγ and the interferential contrast K. This completes the proof.

APPENDIX C PROOF OF THEOREM 2
Since an exact expression of g(γ) in (53) is difficult to calculate, we evaluate it approximately by utilizing the asymp- [25]. According to the definition, evaluating g(γ) is the same as evaluating In order to evaluate the expectation E[(1 − R 2 (z l ))P [l]/λ l ], we first introduce some preliminaries about the noncentral chi distribution N C χ k (λ) with noncentrality parameter λ > 0.
The distribution N C χ k (λ) is the law of the length (2-norm) of a k-dimensional standard normal distribution N (µ, I k ), with λ = µ 2 . Specifically, we are interested in the case where k = 2, since this is the case of the 2-dimensional complex plane. For k = 2, let Y ∼ N C χ2 (m), then we have [35] where L 1/2 denotes the generalized Laguerre function of order 1/2. The function L 1/2 (x) has explicit expression Recall that the asymptotic expansion x(1−R 2 (2 √ x)) ∼ √ x/2 holds for large x, and the random variable P [l] obeys a non-central chi-squared distribution which can be equivalently expressed as Thus, we obtain Finally, substituting the approximation (74) into the exact CRLB (52) yields the conclusion (54), which completes the proof. Note that the operation (x) + in (54) ensures that each term of the CRLB is non-negative.

APPENDIX D PROOF OF THEOREM 3
Since the only imprecise step of the above derivation is the asymptotic expansion, the approximation error of the expectation can be upper-bounded by the asymptotic expansion error. Assume that the asymptotic expansion error does not exceed δ, i.e., Then, by error analysis on (73), we have It can be numerically confirmed from Fig. 9 that δ ≤ 0.07, and that as x → ∞, the approximation error tends to zero. Thus, (74) is a nearly perfect approximation when the interferential We can see from the curve that the approximation error δ ≤ 0.07. SNR γ is large. To be more precise, at γ → +∞, the function g(γ) can be asymptotically expanded aŝ Thus, 1/γ −ĝ(γ) ∼ 1/2 Θ(1/γ), and the relative error r of a single term in the CRLB expression is upper bounded by introducing a positive parameter , which is where (a) comes from applying the triangle inequality, and (b) comes from assuming sufficiently large γ such that the denominator is positive (note that this can be done due to (77)). In order to let this upper bound hold, it should be satisfied that γ > (2 + )δ, and the parameter should satisfy which is equivalent to Such > 0 exists. Since the function 1/(1 − γĝ(γ)) is decreasing for sufficiently large γ, the inequality (80) is equivalent to In fact, choosing = 4 will satisfy all the conditions above when δ = 0.07. Thus, from (78), the relative error is upper bounded by We can easily see from the above inequality that this approximation becomes arbitrarily good when γ → ∞, and the decreasing rate is of order O(1/γ). This completes the proof.

APPENDIX E PROOF OF LEMMA 3
Proof. Let z = exp(iφ) be the complex representation of the estimatorφ, and be the scaled intermediate result in the proposed vM-EM algorithm (Algorithm 2, line 7). Then, on convergence,φ is the fixed point of the algorithmic iteration, i.e., the complex number w is parallel to z, which is equivalent to w * z ∈ R.
Define p(z, w) = (zw * − z * w)/(2i) : C 2 → R, then the convergence of this algorithm is equivalent to w * z ∈ R, which is further equivalent to p(z, w) = 0. Since both z = z(φ) and w = w(φ, s) are functions of the estimatorφ, the output of the vM-EM algorithm is a root of the equation p = 0, i.e., p(z(φ), w(φ, s)) = 0. (84) Notice that if the observed signal is noiseless, i.e., s = µ (0) := α+βexp(i(ψ +ϕ)), then the true value ϕ is a solution to the equation (84). This can be easily seen from the fact that ifφ = ϕ, then µ = µ (0) . However, generally the input of the algorithm s is noisy, which is modeled by s = |µ (0) + v |, where v ∼ i.i.d. CN (0, σ 2 v ) are the thermal noise at theth power sensor. Thus, the estimatorφ can be viewed as a perturbed version of the true value ϕ, i.e.,φ =φ(v), which is a function of the noise v. Thus, we aim to find the Wirtinger derivative [29] ∇ v (φ(v)), where the gradient operator ∇ v acting on f : C L → C is defined through its components: ∂f ∂v 0,R − i ∂f ∂v 0,I , · · · , ∂f ∂v L−1,R − i ∂f ∂v L−1,I T .
(85) Choose ∈ {L}, we first evaluate ∂φ/∂v . Taking the derivative of both sides of (84) w.r.t. v ∈ C and using the derivative formula for implicit functions, we obtain .
Finally, the squared norm of the gradient ∇ vφ is given by which completes the proof.
APPENDIX F PROOF OF LEMMA 5 From the differential representation of the estimation error (42), for all t ∈ [0, 1], we obtain Since H L = H ∞ + O(L −1 ), we can conclude that the upper limit and the lower limit of the sequence H L converges to the same limit H ∞ as L → ∞, and this convergence is uniform in ϕ ∈ [0, 2π]. Specifically, if we denote Thus, by Lemma 3, we can further bound the increasing rate of |∆ϕ| from above, i.e., where s(t) = |µ (0) + tv|, and the components of µ (0) ∈ C L is defined as µ (0) := α + βexp(i(ψ + ϕ)). In order to obtain an upper bound for the increasing rate of |∆ϕ| w.r.t. t, we first try to lower-bound the inner product s(t), x . ∀t ∈ [0, 1]: 99) where (a) comes from applying the triangle inequality to the definition of s(t) and µ (0) . Combining the above inequality with the differential inequality (98), we obtain where the initial condition is |∆ϕ|(t = 0) = |∆ϕ|(v = 0) = 0. In order to solve this differential inequality, we construct an ordinary differential equation (ODE) with zero initial condition as dy du = g L (y, u) := 2 H L G L − X L (u + βy) , y(0) = 0. (101) Since u = 0 is a non-singular point of the ODE, there exists some small positive number δ = δ(L) > 0 such that the solution to (101) uniquely exists in some closed interval I δ = [0, δ] [36]. Note that this δ is chosen for some given L, but the existence of the ODE solution, as well as all the following properties, hold true for arbitrary larger values of L in the same interval I δ , because of the monotone property of the sequences H L , G L , and X L . Furthermore, there exists some positive constant C > 0 such that ∀u ∈ I δ , y(u) ≤ Cu. Particularly, since g L (y, u) is an increasing function of y and u, C can be chosen to be C(δ, L) = 2 H L G L − X L (δ + βy(δ)) .
Since the solution y(u) of this ODE (101) characterizes the maximum value attainable for arbitrary smooth function subject to the increasing rate constraint (100), the solution y(u) in I δ serves as an upper bound to the unknown function |∆ϕ|(tv) as a function of t, i.e., Thus, by further applying y(u) ≤ Cu, we can infer that ∀t : If we apply the constraint v / √ L ≤ δ, then t can be chosen to t = 1, which completes the proof.