Investigation of Phase-Locked Loop Statistics via Numerical Implementation of the Fokker–Planck Equation

Abstract: The goal of this paper is to explore the effect of various parameters on the information geometric structure of the phase-locked loop (PLL) statistics, both transient and stationary. Comprehensive treatment on the behavior of PLL statistics will be given. The behavior of the phase-error statistics of the first-order PLL, in the presence of additive white Gaussian noise (WGN) is investigated through solving the differential equations known as the Fokker–Planck (FP) equation using the implicit Crank–Nicolson finite-difference method. The PLL is one of the most commonly used circuits in electrical engineering. A full knowledge of probability density functions (PDFs) of the phase-error statistics becomes essential in understanding the PLLs. Several illustrative examples are presented to provide profound insights on understanding the PLL statistics both qualitatively and quantitatively. Results covered include the transient and stationary statistics for the nonmodulo-2π probability density function, modulo-2π probability density function, and cycle slipping density function, of the phase error. Various numerical settings of PLL parameters are involved, including the detuning factor and signal-to-noise ratio (SNR). The results presented in this paper elucidate the link between various parameters and the information geometry of the phase-error statistics and form a basis for future investigation on PLL designs.


Introduction
The phase-locked loop or phase lock loop (PLL) [1][2][3][4] is basically a nonlinear feedback control system that tracks the phase of the input signal by detecting the phase error and then adjusting the phase of the output. There are several different types; the simplest is an electronic circuit consisting of a variable frequency oscillator and a phase detector in a feedback loop. The oscillator generates a periodic signal, and the phase detector compares the phase of that signal with the phase of the input periodic signal, adjusting the oscillator to keep the phases matched.
The GPS receiver [5][6][7][8][9] contains a code tracking loop and a carrier tracking loop for tracking the Doppler-shifted carrier simultaneously. A carrier-tracking loop tracks the carrier phase and the code-tracking loop tracks the signal code to within a small fraction of the chip duration. Most GPS receiver designs have a mode of operation that employs a non-coherent delay-locked loop (DLL) for code tracking and a Costas PLL for tracking the Doppler-shifted carrier (the frequency change due to vehicle dynamics). The pseudoranges obtained from the code tracking provide a position fix, while the pseudorange rate (or delta range) estimates obtained from the Costas loop provide a velocity fix, which is accomplished by counting the number of Doppler frequency shifts of carrier cycles that occur over a finite time interval.
The statistics of a first-order PLL output with a white noise input have been shown to be a first-order Markov process. Furthermore, the time varying probability density for any Markov process can be described by the Fokker-Planck (FP) equation [10][11][12][13][14]. The FP equation is a second-order nonlinear, partial differential equation (PDE) that describes the time propagation of probability density functions and the transition probability for these stochastic processes. Viterbi employed the FP equation to describe the variation of probability distribution of the phase error over time. He attempted to solve that equation for time approaching infinity to obtain the steady-state phase-error probability distribution [1,10]. Since the analytical solution to the FP equation is relatively difficult, the numerical approach provides a good alternative. In fact, the transport (advection-diffusion) equation that frequently appears in the field of fluid dynamics [15][16][17][18] bears a very good analogy with the modified FP equation. Therefore, the transport PDE is selected as the model equation for the modified FP equation. Among the available numerical methods [19,20], the finite difference technique is the most attractive one to solve this type of equation.
Explicit and implicit methods are approaches used in the numerical analysis for obtaining numerical approximations to the solutions of differential equations of physical processes. Explicit methods calculate the state of a system at a later time from the state of the system at the current time, while implicit methods find a solution by solving an equation involving both the current state of the system and the later one. The problem of choosing appropriate discretization is present in all finite difference methods-both explicit and implicit. The requirement based on the explicit methods places a serve restriction on the grid spacing size in the time direction with a resultant increase in computation. On the other hand, implicit finite-difference methods are unconditionally stable and need not satisfy this constrain. Therefore, the implicit finite difference methods are computationally advantageous methods with respect to the explicit ones for solving the PDE's. Implicit methods employed in the work include the backward-time/central space (BTCS) implicit method and the Crank-Nicolson (CN) method, which is second-order accurate both in time and space domains. The CN method offers a considerable improvement over the FTCS and BTCS implicit schemes, which are only first-order accurate in time.
Previous numerical solutions of various FP equations in both one-dimensional and two-dimensional have been implemented based on the CN approach [21][22][23][24]. Hollerbach, Dimanche, and Kim [21] investigated the effect of different orders of nonlinear interaction on the geometric structure by considering the case when the equilibrium is a stable point. Jacquet, Kim, and Hollerbach [22] presented the time-evolution of probability density functions (PDFs) in a toy model of self-organized shear flows. Using theoretical analyses of limiting cases, as well as numerical solutions of the full FP equation, a thorough parameter study of PDFs for different values of the correlation time and amplitude of stochastic forcing was presented. Hollerbach, Kim, and Mahi [23] reported the PDFs using direct numerical solutions of the FP equation for a classical double-well stochastic resonance, showing clear signature in information length. Kim, Jacquet, and Hollerbach [24] investigated the case where a stationary PDF that is unimodal for δ-correlated stochastic noise becomes bimodal when the noise has a finite correlation time. As for the periodic-in-space systems [25][26][27], the Fourier series approach might be potentially better than finite-differencing in space with the CN method. For example, in the work by Hollerbach and Kim [25], the effect of different spatially periodic, deterministic forces on the information geometry of stochastic processes was explored. Kim and Hollerbach [26] presented time-dependent PDFs for a nonlinear stochastic process with a cubic force using analytical and computational studies. Numerical simulation of the FP equation revealed detailed evolution of the time-dependent PDF. Guérin and Dean [27] studied the time-dependent dispersion properties of overdamped tracer particles diffusing in a near critically tilted one-dimensional periodic potential under the influence of an additional constant tilting force.
Stochastic noise is ubiquitous and plays a crucial role in the evolution of many different systems such as the PLL. Given uncertainty inherent in the systems, a probabilistic description is essential for understanding the dynamics of stochastic systems. A direct numerical solution to the transient behavior of the PLL phase-error process was carried out by Dominiak and Pickholtz [28] by comparing the FP equation to the one-dimensional heat-flow equation. The numerical method for solving the transient behavior of the PLL phase-error process was carried out by Dominiak and Pickholtz [28] by comparing the FP equation to the one-dimensional heat-flow equation. The numerical method for solving the modified FP equation employed in references [29][30][31] is the Euler's forward time/central space (FTCS) explicit method. To elucidate the link between various parameters and the information geometry of the phase-error statistics, in the presence of additive white Gaussian noise, it is useful to quantify the difference among different probability density functions (PDFs) and compare their evolution. Based on the preliminary study presented in reference [29], extension of the work for more comprehensive treatment and profound insights on the behavior of PLL statistics will be presented in this paper.
The remainder of this paper is organized as follows. In Section 2, the phase-locked loop and normalized FP equation are reviewed; applications of the FP equation are also addressed. In Section 3, the numerical methods using finite difference technique for discretization of the model equation, known as the transport (advection-diffusion) PDE are discussed; numerical formulations are then developed for the FP equation. Numerical implementation and analysis on the PLL statistics are given in Section 4. Illustration and discussion for both the transient and stationary statistics, including the nonmodulo-2π probability density function (PDF), modulo-2π PDF, and cycle slipping density function, with various parameters, such as the detuning factor and signal-to-noise ratio (SNR), are presented. Section 5 summarizes the main results of this paper.

The Phase-Locked Loop and Normalized Fokker-Planck Equation
In this section, a brief review on the FP equation and the related application issues to PLLs is provided. Detailed discussion can be referred to references [1][2][3][30][31][32][33]. The procedure is computationally efficient, and readily provides numerical results describing the phase-error process, the modulo-2π phase-error process, and the cycle-slipping statistics.
The basic PLL configuration contains a phase detector (PD), for detecting the phase error between the output and the input through feedback system; a loop filter; and a voltage-controlled oscillator (VCO), which is the oscillator to be synchronized by adjusting the phase difference. The mathematical block diagram of the phase-locked loop is shown as in Figure 1. The transfer function ) (s F represents the loop filter of a tracking loop. The loop filter is usually a low-pass filter used to suppress the noise and high-frequency signal components from the phase-detector/delay discriminator and provide a dc-controlled signal for the VCO. A first-order tracking loop is obtained when Figure 1. Mathematical block diagram of the phase-locked loop (PLL). The basic PLL configuration contains a phase detector (PD), for detecting the phase error between the output and the input through feedback system; a loop filter; and a voltage-controlled oscillator (VCO). The loop filter ) (s F is usually a low-pass filter used to suppress the noise and high-frequency signal components from the phase-detector/delay discriminator and provide a dc-controlled signal for the VCO [1][2][3]. The equation describing the loop operation is the first-order differential equation given as [1][2][3][4]10,30]  . The basic PLL configuration contains a phase detector (PD), for detecting the phase error between the output and the input through feedback system; a loop filter; and a voltage-controlled oscillator (VCO). The loop filter F(s) is usually a low-pass filter used to suppress the noise and high-frequency signal components from the phase-detector/delay discriminator and provide a dc-controlled signal for the VCO [1][2][3].
The equation describing the loop operation is the first-order differential equation given as [1][2][3][4]10,30] where φ = Φ i − Φ is the error between the phase angle of the incoming signal and the output of the VCO, K is the total loop gain, A is the input signal voltage amplitude, ω o is the quiescent frequency of the VCO, and n (t) is a stationary white Gaussian noise (WGN) process. Assume that the input signal to the PLL is a pure sinusoid with the average power of A 2 Watt, a constant frequency ω rad/s, and that the signal is corrupted by a stationary white Gaussian noise with a single-sided power spectral density N 0 Watt/Hz. For a continuous Markov process, the instantaneous probability density function (PDF) of the phase error φ at time t, p(φ, t), satisfies the partial differential equation with the initial condition where φ o is the initial value of φ; A(φ) and B(φ) are normalized moments of the infinitesimal increment [10]. The equation that describes the stochastic behavior of the phase error φ(t) of the first order PLL is the Fokker-Planck equation After substituting Equation (5) into Equation (4), the Fokker-Planck equation that describes the PDF p(φ, t) of the Markov process phase error φ(t) can be written as Although the equation is linear in p, the complete solution for p(φ, t) is somewhat complicated by the nonlinear behavior of the variable coefficients [10,14].

Nonmodulo-2π Phase-Error PDF
An initial distribution of the dependent variable and two sets of boundary conditions are required for a complete description of the problem. Equation (7) can be solved for the PDF of phase error with the initial condition lim and normalizing condition where δ(φ − φ o ) is the unit impulse function, also known as the Dirac delta function possessing the following property The PDF p(φ, τ) can be written as which can be substituted into Equation (7) and lead to the following equation The function k(φ) can be found to be [31] and consequently, Equation (12) becomes where the density function g(φ, τ) satisfies the following one-dimensional advection-diffusion parabolic PDE, called the modified FP equation Solving for g(φ, τ) is computationally simpler than solving for p(φ, τ), mainly due to the fact that coarser grid spacing is allowed for g(φ, τ). Additionally, p(φ, τ) can be reconstructed readily using Equation (15) without difficulty when g(φ, τ) is obtained.

Modulo-2π Phase-Error PDF
A common form of partial solution to Equation (7) is the modulo-2π probability [1][2][3][4][10][11][12][13] where the capital letter P is used to denote the modulo-2π PDF while little letter p denotes the unbounded nonmodulo-2π PDF, respectively. In the steady state, p(φ, τ) is spread over all φ since cyclic slipping event (lock on at k2π) have occurred. The relation between these two PDF's can be expressed as P(−π < φ < π, τ) = p(φ, τ)modulo 2π (18) which shows that P(φ, τ) has a period of 2π rad. This is like the original function except that each 2π-wide region contains an identical distribution and P begins with impulses at k2π for all k. As P diffuses with time, the probability that the phase will jump out of the primary interval −π < φ < π is the same as the probability that it will jump into that interval. As time increases, p(φ, τ) diffuses so the ultimate value approaches zero everywhere. Consider only the interval −π < φ < π and set the boundary and normalized conditions, respectively: Appl. Sci. 2020, 10, 2625 6 of 25 Hence the modulo-2π PDF describing the phase error is obtained.

Cyclic Slipping Density Function
A cycle slip is defined as an increase (positive cycle slip) or decrease (negative cycle slip) of the phase error by 2π. Cycle slipping in a phase-locked oscillator occurs where the dynamic phase error is larger than 2π. If the error should become too large that the VCO skips cycles, the loop is considered to have lost lock. The integral of this density becomes 2π −2π As long as the loop phase error remains within the limits of a specified subinterval, the density function q(φ, τ), accompanied with the initial condition q(φ, 0) = δ(φ), will satisfy the FP equation and thus The use of q(φ, τ) for the bounded case is to distinguish it from the unbounded nonmodulo-2π function p(φ, τ). The numerical algorithm used to generate q(φ, τ) is the same as for generating p(φ, τ), except that the end points q(φ l , τ) and q(−φ l , τ) are forced to be zero, whereas the end points for p(φ, τ) were found by extrapolation, which simulates an infinite interval.

Numerical Formulations for the Normalized Fokker-Planck Equation
The solution procedure of a PDE depends on the type of the equation. Imposition of initial and/or boundary conditions also depends on the type of PDE. To make this paper self-contained, the main steps for discretization of the typical model equations are introduced first. After that, numerical formulations using the implicit method for the normalized Fokker-Planck equation are presented.

Discretization of the Model Equation
The heat conduction equation is a parabolic PDE with the form which is typically characterized by a diffusivity (or dispersion) coefficient ρ. The convective (or advection) equation such as the wave equation is categorized as a hyperbolic type given by Consider the one-dimensional transport equation with both convection and diffusion terms, where u is a scalar, being convected with a known velocity v and being diffused. Equation (25) is a one-dimensional non-linear convective diffusion equation, or a hyperbolic parabolic equation; it is hyperbolic due to the convective term, v∂u/∂x, and it is parabolic due to the diffusive term, ρ∂ 2 u/∂x 2 . Like the diffusion equation the transport equation is strictly parabolic. Closed-form solution u(x, t) of Equation (25) is usually not available except for very simplified cases. It can, however, be obtained numerically in most applications, for example, by the finite difference and the finite element methods. The transport equation, as the model equation in this paper, will be implemented. Explicit methods are easier to solve than implicit methods in terms of time and effort spent on a computer. However, there exist several disadvantages in explicit methods. For the explicit finite difference method to be stable, the grid spacing in the time domain and grid spacing in the space domain need to satisfy certain requirements. The choice of grid spacing is made in accordance with the condition of stability. Implicit methods offer great advantage on the stability of the finite difference equations, since most are unconditionally stable. Therefore, a larger grid size in time is permitted. It should be noted that the selection of a larger time step is limited due to accuracy consideration because an increase in time step will increase the truncation error of the finite difference equation.
In the finite-difference methods, partial derivatives are approximated using Taylor's series; hence, solutions can be obtained using the first-, second-, or higher-order accurate methods. In this paper, the following finite difference approximations are selected for the partial derivatives: An equation is defined as being implicit when more than one unknown appears in the finite difference equation. For the implicit methods, a set of simultaneous equations needs to be solved, which requires more computation time per time step for the implicit methods. Implicit methods generally possess better stability properties and higher accuracy than the explicit methods. For implicit formulation, most are unconditionally stable. Implicit methods include (1) the backward time/central space (BTCS) implicit method, and (2) the Crank-Nicolson (CN) method. For the model equation, the BTCS implicit scheme, also known as the simple (Laasonen) implicit scheme, leads to the following formulation which can be rearranged as By considering all spatial nodes Equation (31) produces a tridiagonal system of equations that can be solved efficiently using the Thomas algorithm [17].
The CN scheme evaluates the spatial derivative at the average of the nth and (n + 1)-th time levels, i.e., at the (n + 1/2)-th time level. For the model equation, the CN scheme gives which can be rewritten as Equation (33) is consistent with Equation (25) with a truncation error of O[(∆t) 2 , (∆x) 2 ] since the Taylor expression is made about ( j, n + 1/2), hence second-order accurate spatial differences are inserted into this algorithm. The CN method offers a considerable improvement over the FTCS and BTCS implicit schemes whose local truncation errors are only first-order accurate in time. A von Neumann stability analysis indicates that the CN scheme is unconditionally stable.

Discretization of the Normalized Fokker-Planck Equation
The modified FP equation is a one-dimensional non-linear hyperbolic parabolic PDE in which both convection and diffusion terms are involved. The numerical algorithms for the FP equation can be obtained based on the procedures discussed in Section 3.1. The equation to be solved is Equation (16), which is analogous to the model equation, Equation (25), with the relations: The implicit finite difference methods are computationally advantageous methods with respect to the explicit ones since they are unconditionally stable. In addition to the Euler's BTCS implicit method, the CN method is second-order accurate both in time and space domains.
Applying the Euler's BTCS implicit scheme for the FP equation, we have the following formulation for numerical implementation and which can be rearranged Equation (37) can be represented as where ξ = (sin ϕ − γ) ∆τ ∆ϕ , representing the Courant number, and η = ∆τ α(∆φ) 2 , representing the diffusion number.
Employing the CN scheme, we have With the finite-difference operators defined as Appl. Sci. 2020, 10, 2625 9 of 25 the CN implicit formulation, after rearrangement, is given by

Numerical Implementation and Analysis of the Phase-Locked Loop Statistics
In this section, numerical implementation and analysis of the first-order phase-locked loop statistics are presented by solving the normalized FP equation. Computer codes were developed using the Matlab ® software. According to Equation (9), the initial condition is supposed to be a Dirac delta function. Since the delta function is a continuous-time signal and therefore, for numerically resolved purpose, the value 1/∆φ is used for approximating the delta function. In the following illustrative examples, the initial value of phase error φ 0 = 0 is utilized for discussion.
First of all, the time-evolution of PDFs is investigated to reveal detailed evolution of the time-dependent PDFs of the PLL. Typical results for the time-evolution of nonmodulo-2π transient statistics for SNR α = 1 and detuning factor γ = 0 are shown in Figures 2-4. The grid spacing for the phase-error domain is chosen as ∆φ = π/50 rad; the grid spacing in the time marching direction is set as 0.01. Figure 2 shows the transition phase-error density function g(φ, τ) and the PDF p(φ, τ). The numbers beside individual curves indicate the normalized times, switching τ from 0.625 to 20. The corresponding cumulative probability distribution of the phase-error is depicted in Figure 3. Typical time-evolution of p(φ, τ), as shown in Figure 2b, revealed that as P diffused with time, the probability that the phase would jump out of the primary interval −π < φ < π was the same as the probability that it would jump into that interval, implying that the possibility of larger phase error increased. In the case of γ = 0, both of the error density functions g(φ, τ) and p(φ, τ) diffuse with time while remaining symmetric with respect to the phase-error axis, and thus the initial global peak values of the density functions g and p at zero phase error (φ = 0) decreased as time progresses. As functions of φ, the local peak positions of p appeared at a period of 2π rad. The amplitudes of local peaks were monotonically nonincreasing when moving outwards from their mean values with respect to the phase-error axis. The figures and later ones as well only show part of the results. In Figures 2  and 3, the results at a phase error between −20 and 20 rad were the area of interest. Figure 4 provides the three-dimensional plots for nonmodulo-2π transition g and p, respectively, for more profound insights into the statistics information, especially the relation between g/p and phase error, as time progresses. Both the front view and back view are shown such that better understanding on the behavior of PLL statistics is accessible. The result enabled us to mathematically quantify the difference among different PDFs, providing a meaningful conceptual link between the stochastic process and information geometry, and could be utilized for optimizing various desired outcomes.
Figures 5-7 provide information basically similar to that of Figures 2-4, and thus the parameter values were identical except that the detuning factor was now changed to γ = sin(π/4). The identical numerical values including the normalized time (τ), the grid size for the phase-error domain (∆φ) and the grid size in the time marching direction (∆τ) were utilized for comparison. It can be seen that in this case the density functions diffused to one direction as time progresses, given a certain SNR (e.g.,α = 1 as demonstrated here). The error probabilities g and p diffused virtually in opposite directions. In this example, g virtually diffused leftwards while p diffused rightwards. As can be seen, p had a right-skewed (or positively skewed) distribution, and most data fell to the right (positive direction) of the graph's peak. The histogram skewed in such a way that its right side (or "tail") was longer than its left side. Again, for the density function p, there were local peaks appearing with a period of 2π rad. It should be noticed that the position of global peak of p, which initially appeared at the phase error near by zero, moved rightwards as time progressed. As can be seen, the peak position appeared at 2π rad away from its original position when τ exceeded a certain value, e.g., τ = 10 in this illustration. According to Equation (10), if the solution is properly resolved, the total probability remains constant (=1) for all time. The cumulative probability distribution can be used to check the result. As compared to the actual total probability, the numerically calculated cumulative probabilities for various normalized time are summarized in Table 1. The errors using the numerical approach was of the order of 0.02% to 0.5% depending on the various τ, for γ = 0 and γ = sin(π/4), respectively. For better resolution that is required, a smaller grid can be utilized. seen, the peak position appeared at 2π rad away from its original position when τ exceeded a certain value, e.g., τ = 10 in this illustration. According to Equation (10), if the solution is properly resolved, the total probability remains constant (=1) for all time. The cumulative probability distribution can be used to check the result. As compared to the actual total probability, the numerically calculated cumulative probabilities for various normalized time are summarized in Table 1. The errors using the numerical approach was of the order of 0.02% to 0.5% depending on the various τ, for 0 γ = and ) 4 / sin(π γ = , respectively. For better resolution that is required, a smaller grid can be utilized.            Figure 7. Three-dimensional plots for the nonmodulo-2π transition g and p, respectively, for SNR α = 1 as time progresses: front view (left) and back view (right; γ = sin(π/4)). In the next example, an illustration of how the SNR α influences the information geometry of the phase-error statistics is presented. Figures 8-10 provide the results of the nonmodulo-2π statistics at normalized time τ = 10 for SNR α in the range 0.3-1 where the detuning factor γ = 0 was utilized. The results are shown to capture the effect of different SNR values and understand the information changes on the error-phase statistics of PLL. Figure 8 presents the phase-error density function g(φ, τ) and the PDF p(φ, τ) and Figure 9 shows the corresponding cumulative probability distribution of p(φ, τ). The numbers beside individual curves indicate the various SNR values. Again, three-dimensional plots for nonmodulo-2π g and p, respectively, are given in Figure 10. As in the presented case where γ = 0 was utilized, both the density functions g and p remained symmetric with respect to the phase-error axis as α varied. As α decreased, the global peak values of the density functions at φ = 0 decreased and thus the initial narrow density functions became 'flatter', implying that the variance the phase error increased, and the possibility of a larger phase error increased. As functions of φ, the local peak positions appeared at a period of 2π rad when moving outwards, with the peak amplitudes monotonically nonincreasing. Figures 11-13 elucidated information similar to that of Figures 8-10, and the parameter values were identical except that the detuning factor γ = sin(π/4) was utilized. It can be seen that the density functions g and p virtually stretched (or expanded) to one direction where they stretched to opposite directions. As shown in this example, g had a left-skewed distribution, while p had a right-skewed distribution. The density p shows that the peak values appeared at the phase error near by 2π rad. As α decreased, p tended to be flatter, and the variance of the phase error increased. Three-dimensional plots of information geometry were provided such that more comprehensive understanding could be gained on the relationship among the parameters, p/q, φ, and α.
To elucidate the effect of detuning factor γ on the error-phase statistics of PLL, the following illustrations are shown. Figures 14-16 provide the results for the nonmodulo-2π statistics at normalized time τ = 10 for various detuning factors γ, where the signal was assumed to be sufficient and a high SNR value α = 1 was utilized for illustration. The numbers beside individual curves indicate the detuning factor γ, from γ = 0 to sin(0.5π). The density functions stretched to one direction as the detuning factor γ increased, where g and p stretched to opposite directions. As can be seen in this example, g moved leftwards while p went rightwards. Furthermore, the increase of the γ value caused the density functions to flatten, revealing the possibility of the larger phase error increasing. The global peak of the density function p, which initially appeared at the phase error near by zero, moved rightwards with the increase of γ. As can be seen from the graphs, the peak position appeared at 2π rad away from its original position when γ was sufficiently large, e.g., sin(0.3π) in this example. The three-dimensional plots provided elucidated more comprehensive geometry information on the relationship among the parameters, p/q, φ, and γ.         In the following discussion, the modulo-2π PDF, ( , ) P φ τ (denoted by a capital letter P) as well as the cumulative probability distribution, and the cycle slipping density function ) , In the following discussion, the modulo-2π PDF, P(φ, τ) (denoted by a capital letter P) as well as the cumulative probability distribution, and the cycle slipping density function q(φ, τ) of the phase error for the first-order PLL are presented. The time-evolution of modulo-2π phase-error PDF P(φ, τ) and cumulative probability distribution is shown in Figure 17. In this illustrative example, the SNR α = 1 and detuning factor γ = 0 were utilized. The numbers beside individual curves indicate the normalized times, switching τ from 0.625 to 5. The steady state was reached when the normalized time τ ≥ 5. Since the total time steps for simulation were set as 1000, the grid spacing in the time marching direction was ∆τ = τ/1000, which is dependent on the normalized time τ. As time τ progressed, the peak value of P(φ, τ) at zero phase error (φ = 0) was decreasing, while remaining symmetric with respect to the phase-error axis.
that the phase was completely unknown.
The function ) , for the primary cycle slip was essentially identical to the ) , ( τ φ p function until a significant amount of probability mass was lost at the absorbing boundaries at = φ π ± . The inclusion of absorbing boundaries around some subinterval of φ , e.g., [−π, π], enabled one to determine the probability that the phase error had not yet exceeded the limits of that subinterval at any prior time. The end points were held at zero to simulate the effect of the slip boundary. Figure 20 shows the time-evolution of primary cyclic slipping density function ) , The curves were plotted for the conditions SNR 1 = α , where three values of detuning factor, 0 = γ , sin( / 4) π and sin( / 2) π , respectively, are depicted.
The results provide us a meaningful conceptual link between the stochastic process and information geometry. Through the above examples of physical realizations, the results underscore the natural, intuitive interpretation of information profiles and changes associated with the stochastic process in the phase-locked loop, and can be utilized for optimizing various desired outcomes.   Figure 18 provides the comparison on variation of modulo-2π phase-error steady state PDF P(φ) and cumulative probability distribution function for various detuning factor γ where SNR α = 1 was utilized. The numbers beside individual curves indicate the detuning factors in the range 0-sin(π/2). It can be seen that, initially with a symmetric shape, P(φ) became the left-skewed distribution. The peak moved rightwards with its amplitude decreasing as γ increased. Figure 19 provides the variation of modulo-2π phase-error steady state PDF and cumulative probability distribution function, respectively, of the phase error at τ = 0.625 for the SNR values α in the range 0-1 in the case of γ = 0. The numbers beside individual curves indicate the SNR values α. It can be seen that, the peak value of P(φ) appeared at φ = 0 and was monotonically decreasing, while remaining symmetric with respect to the phase-error axis as α was increased. For a very high SNR (α = 1),P(φ) resembled a unity Gaussian distribution, while for very low SNR ( α → 0 ), it became flat and approached a uniform distribution from -π to π, with the PDF P(φ) = 1/(2π rad), implying that the phase was completely unknown.
The function q(φ, τ) for the primary cycle slip was essentially identical to the p(φ, τ) function until a significant amount of probability mass was lost at the absorbing boundaries at φ = ±π. The inclusion of absorbing boundaries around some subinterval of φ, e.g., [−π, π], enabled one to determine the probability that the phase error had not yet exceeded the limits of that subinterval at any prior time. The end points were held at zero to simulate the effect of the slip boundary. Figure 20 shows the time-evolution of primary cyclic slipping density function q(φ, τ). The curves were plotted for the conditions SNR α = 1, where three values of detuning factor, γ = 0, sin(π/4) and sin(π/2), respectively, are depicted.
The results provide us a meaningful conceptual link between the stochastic process and information geometry. Through the above examples of physical realizations, the results underscore the natural, intuitive interpretation of information profiles and changes associated with the stochastic process in the phase-locked loop, and can be utilized for optimizing various desired outcomes. Appl. Sci. 2020, 10, x FOR PEER REVIEW 22 of 25   For very low SNR (α→0), P approaches a uniform distribution sin(π/4) sin(π/2) γ = 0 sin(π/8) sin(π/4) sin(π/2) γ = 0 sin(π/8) Figure 18. Variation of modulo-2π phase-error steady state: (a) PDF and (b) cumulative probability distribution function, for SNR α = 1 with various detuning factors γ (τ ≥ 5 or steady-state).

Conclusions
Information geometry of the phase-error statistics for the first-order phase-locked loop was presented. A thorough parameter study of PDFs for different values of the PLL parameters was presented by solving the Fokker-Planck equation numerically. Due to the fact that the modified Fokker-Planck equation bears a very good analogy with the transport equation, the numerical algorithms were derived for the transport equation and then applied to solve the Fokker-Planck equation. The implicit Crank-Nicolson finite-difference method was employed in that it is computationally advantageous in terms of both accuracy and stability.
Illustrative examples include the transient and stationary statistics for the nonmodulo-2π probability density functions and cumulative probability distribution functions, modulo-2π probability density functions and probability distribution functions, and cycle slipping density function of the phase error, with various settings of parameters involved, including the detuning factor γ and SNR α. The property of skewness on the probability density/histogram due to the non-zero detuning factor γ was discussed. Additionally, shown for the nonmodulo-2π PLL statistics were the three-dimensional plots for g and p, respectively, of the phase error for profound insights into the PLL behavior. As for the modulo-2π statistics, the time-evolutions of phase-error PDF and cumulative probability distribution function, respectively, were shown. Comparison on the variation of steady-state modulo-2π phase-error statistics, for a specific SNR with various detuning factors was presented. Furthermore, variation of modulo-2π statistics of the phase error for various SNR values was also involved. Furthermore, the time-evolution of cyclic slipping density function for varying detuning factor was addressed.
The results presented in this paper provided comprehensive treatment and elucidated the statistics of the PLL phase error from different perspectives, qualitatively and quantitatively in some extent, and are useful for future PLL investigation and designs.

Conclusions
Information geometry of the phase-error statistics for the first-order phase-locked loop was presented. A thorough parameter study of PDFs for different values of the PLL parameters was presented by solving the Fokker-Planck equation numerically. Due to the fact that the modified Fokker-Planck equation bears a very good analogy with the transport equation, the numerical algorithms were derived for the transport equation and then applied to solve the Fokker-Planck equation. The implicit Crank-Nicolson finite-difference method was employed in that it is computationally advantageous in terms of both accuracy and stability.
Illustrative examples include the transient and stationary statistics for the nonmodulo-2π probability density functions and cumulative probability distribution functions, modulo-2π probability density functions and probability distribution functions, and cycle slipping density function of the phase error, with various settings of parameters involved, including the detuning factor γ and SNR α. The property of skewness on the probability density/histogram due to the non-zero detuning factor γ was discussed. Additionally, shown for the nonmodulo-2π PLL statistics were the three-dimensional plots for g and p, respectively, of the phase error for profound insights into the PLL behavior. As for the modulo-2π statistics, the time-evolutions of phase-error PDF and cumulative probability distribution function, respectively, were shown. Comparison on the variation of steady-state modulo-2π phase-error statistics, for a specific SNR with various detuning factors was presented. Furthermore, variation of modulo-2π statistics of the phase error for various SNR values was also involved. Furthermore, the time-evolution of cyclic slipping density function for varying detuning factor was addressed.
The results presented in this paper provided comprehensive treatment and elucidated the statistics of the PLL phase error from different perspectives, qualitatively and quantitatively in some extent, and are useful for future PLL investigation and designs.