\section{Introduction}
The phase noise (PN) remains one of the main limiting factors for wireless communications. Phase noise due to instabilities of local oscillators, both at the transmitter and at the receiver, can in fact cause a severe performance degradation.
Ideally, a local oscillator would produce a sinusoidal signal, whose power spectrum is a delta function at the carrier frequency, but in reality its output power appears also in a band around the desired frequency.

Since the PN can strongly limit the performance, the study of the PN effects in communication systems has attracted a lot of interest in the literature of the last decades, see, e.g., ~\cite{PoBlMo95,MeMo97,To98,ArCa98, Ar01,WuBa04,CoBaCa05,BaCo07,MuErSv08,NiPa09,Co14,KrKhErSv13,KhKuPa14,KrCoGrEr15b}. In this context, the assumed PN statistical model plays a crucial role. For example, having a good statistical model for the PN can support the design of estimation/detection algorithms, or allowing accurate performance prediction.
Hence, in the communication community, the studies about PN are normally based on quite simple discrete-time models.

One of the simplest discrete-time models for PN assumes that the samples are independent random phases, each distributed uniformly in the interval $[-\pi, \pi)$. As an example, this model is adopted in~\cite{KaSh04}, with the aim of obtaining capacity bounds for the considered noncoherent channels in the case of ideal interleaving. In~\cite{PeSh98,CoRa01,NuAn05},  a constant channel phase error over a block of symbols is considered for computing the noncoherent capacity of the MPSK channel.
Discrete-time models that are used to describe the residual phase error after phase tracking~\cite{FoGiWe73,KrGrErCo13b,CoYa18,UgPiEr19} assume that the distribution of the phase samples is Gaussian or Tikhonov, which is similar to a wrapped $[-\pi,\pi]$ Gaussian. Further, oscillator PN is strongly correlated with a non-flat spectrum, as we will study in this paper. A widely studied discrete-time PN model, which considers the correlation between PN samples, is the Wiener model~\cite{CoBaCa05,De06,DeMeRo00 ,GhKr13,GhKr17}. The strength of this model resides in the simplicity, since it is defined by a single parameter, i.e., the variance of the phase increments. Blockwise Wiener channel models are also studied~\cite{GaMaNi20}.

In the literature, several models for the power spectral density (PSD) of the PN can be found. The analysis of real oscillators is based on direct measurements of the PSD of the PN $\theta(t)$ or of the phasor $e^{j\theta(t)}$. These measurements are related to the carrier frequency and are given in dBc/Hz~\cite{De06}. PSD models that try to match the measurements can be found in~\cite{TuCo05,HaDrHa98,Mu95,RoKa95}. These works use piece-wise linear PSD models for phase-locked (PLL-locked) oscillators. The phase noise model of 3GPP, described in~\cite{3GPPTR38803,R4-1701165}, is a multi-pole/zero model extended to fractional orders that matches well the measurements. In some contexts, the PN is described by \textit{masks}, that specify the frequencies and the maximum acceptable PN level at those frequencies, given in  dBc/Hz, that guarantee the target performance. Several examples of PN masks adopted in satellite communications can be found in~\cite{DVB-S2X_CM}, for the Digital Video Broadcasting DVB-S2~\cite{DVB-S2X}.

Despite the large number of publications on the subject, there are a lot of questions that are not concisely answered yet. For example, how can we relate the parameters of the PN model to the oscillator measurements?
In fact, common time-domain PN models are not explicitly linked to the oscillator measurements, for example, there is no correlation with oscillator data sheets, such as~\cite{datasheet}. Moreover, the above mentioned PSD models, which try to match measurements, can be used for computer simulations but do not lead to an analytical time domain model for the PN. Therefore, a good time-domain PN model that is related to the measurements is still missing in the literature.
Another question that is worth to consider is what is a good discrete-time model for a free-running versus a PLL-locked oscillator? The largely used Wiener model can describe a free-running oscillator, while less is said in the literature about discrete-time models for PLL-locked oscillator, one example is~\cite{Ar01}.
Considering that the PN can be described either by the phase process or by the phasor process, we can ask what is the relationship between the PSDs of these two processes.
An important point regards a largely used discrete-time channel model, that ignores the power loss and the intersymbol interference (ISI) affecting the received symbols due to the PN. The signal with
PN will not be exactly matched to the matched filter in the receiver. In~\cite{GhKr17}, the effect of filtering on the phase noise is considered by using a multi-sample receiver, in the case of Wiener PN and in the absence of ISI.
In~\cite{MaTrRa13}, the ISI effect is discussed in terms information rate by using an oversampled model for the received signal, in the case of Wiener PN. However, to the authors' knowledge, the quantification of these neglected effects  is missing in the literature.

\subsection{Contribution}
In this work, we give an answer to these questions starting from a general analytical model that can describe the PSD of the PN of real oscillators.

\textbf{Phase noise spectrum model:} We propose a novel phase noise spectrum model
\begin{equation}\nonumber
S_\theta(f)=\frac{10^{10}\ell^2_{100}}{f^2_{\mathrm{3dB}}+f^2}+\ell^2_\infty\, ,
\end{equation}
where $f_{\mathrm{3dB}}$ is the 3dB bandwidth, and $\ell^2_{100}$ and $\ell^2_\infty$ are the spectrum levels for $f\!=\!100$~kHz and for high frequency, respectively. The proposed model can represent a PLL-locked oscillator including white phase noise, and is very suitable for theoretical calculations. The proposed PSD is characterized by a $-20$~dB/decade slope and by a flat part at low frequencies, representing the attenuation of the PN stabilized by means of a PLL, and by another flat part at high frequencies.

\textbf{Phasor spectrum model:}
Starting from the PN PSD, we derive the PSD of the phasor $e^{j\theta(t)}$ and express it in terms of the parameters of the PN PSD, that are related to the oscillator measurements. The phasor process   has the spectrum
\begin{align}
&S_{h} (f)=e^{-\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}}  \delta(f) + \nonumber\\ \nonumber
&\;\;\;e^{-\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}}  \! \! \int_{-\infty}^{\infty}  \!\!\!\!  (  e^{     \frac{  \pi 10^{10} \ell^2_{100} e^{-2 \pi f_{\mathrm{3dB}} |\tau | } }{f_{\mathrm{3dB}}}   } \!\!\! -\!\!1 \!   ) \mathrm{cos}(2\pi f \tau) \mathrm{d} \tau + \ell^2_\infty\, ,
\end{align}
where $\delta(f) $ denotes the Dirac delta function.
Two important special cases are the \textit{free-running oscillator} case when \mbox{$f_{\mathrm{3dB}} \rightarrow 0$},
\begin{equation}\nonumber
S_h(f)_{f_{\mathrm{3dB}} \to 0}  =  \frac{ 10^{10}\ell^2_{100}}{\pi^2 10^{20} \ell_{100}^4  + f^2} +\ell^2_\infty\, ,
\end{equation}
and the \textit{PLL-locked oscillator} case when $f_{\mathrm{3dB}}\gg \pi 10^{10} \ell^2_{100}$,
\begin{equation}\nonumber
S_h(f)_{\mathrm{high}f_{\mathrm{3dB}} }
=\Big( 1- \frac{\pi 10^{10} \ell^2_{100}  }{f_{\mathrm{3dB}}} \Big) \delta(f)   + \frac{  10^{10} \ell^2_{100} }{f^2_{\mathrm{3dB}}+f^2}+\ell^2_\infty\, .
\end{equation}
It is notable that the phasor spectrum overlaps with the phase noise spectrum for high frequencies.

\textbf{Time domain model:}
We show that a discrete-time model for phase noise is
\begin{align}\nonumber
z_k=&x_k e^{j \theta_k}  + w_k \\\nonumber
\theta_k=&\theta^c_k+\theta^f_k\, ,
\end{align}
where $\theta^c_k$ and $\theta^f_k$ are independent processes, the former inducing the $-20$~dB/decade slope in the PN PSD, which dominates close-to-carrier, the latter inducing the flat part of the PSD, that dominates at high frequencies. In particular, $\theta^f_k$ are independent and identically distributed (i.i.d.) zero-mean Gaussian variables with variance $\sigma^2_f=\ell^2_\infty/T_s$, while $\theta^c_k$ can be described as
\begin{equation}\nonumber
\theta^c_k=a\theta^c_{k-1}+u_k \, ,
\end{equation}
where $u_k$ are i.i.d. zero-mean Gaussian variables with variance $\sigma^2_u$ and the parameters are computed from
\begin{align}\nonumber
a&=e^{-2\pi f_{\mathrm{3dB}} T_\mathrm{s}}\\ \nonumber
\sigma^2_u&=\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}} (1-e^{-4\pi f_{\mathrm{3dB}} T_\mathrm{s}})\, .
\end{align}
In the case of a free-running oscillator, the proposed time domain model reduces to the Wiener model, with increment variance
\begin{equation} \nonumber
\sigma^2_u= 4\pi\rho \, ,
\end{equation}
where $\rho\triangleq \pi 10^{10} \ell^2_{100} T_\mathrm{s}$ represents the ratio between the bandwidth of the phasor and the signal bandwidth.

\textbf{Error analysis:}
The PN process has an infinite bandwidth, and some approximation errors will occur when using a discrete-time model. Further, the communication signal with PN will not be exactly matched to the matched filter in the receiver. We bound the approximation errors that the discrete-time representation suffers from, and show that for most realistic cases those approximation errors are limited. For example, in the case of a free-running oscillator, to quantify the ISI due to the presence of PN, we derive the signal-to-interference (SIR) ratio in closed form as
\begin{equation}\nonumber
\text{SIR}=\frac{   \text{tan}^{-1}   \Big( \frac{4 \pi}{\sigma_u^2} \Big) \Big[ 1 -\frac{\sigma_u^4}{16\pi^2} \Big] -\frac{\sigma_u^2}{4\pi} \text{log}\Big(1+\frac{16\pi^2}{\sigma_u^4}\Big) + \frac{\sigma_u^2}{4\pi} }{ \frac{\sigma_u^4}{16\pi^2} \mathrm{tan}^{-1}   \Big( \frac{4 \pi}{\sigma_u^2}   \Big)  +  \frac{\sigma_u^2}{8 \pi}  \mathrm{log}\Big(1+\frac{16\pi^2}{\sigma_u^4}\Big)       -\frac{\sigma_u^2}{4 \pi} }\, .
\end{equation}
Interestingly, from the above equation, we can find a limit on $\sigma_u$ to have a SIR higher than a given value, i.e., the SIR is higher that $25$~dB if $\sigma_u<0.1$~rad.

\subsection{Outline}
The remainder of this paper is structured as follows. In Section~\ref{s:model}, we introduce the system model. In Section~\ref{s:PNmodel}, the proposed statistical model of the PN is described, in particular, we propose an analytical form for the PN PSD and from that we derive the PSD of the phasor. In Section~\ref{s:DT channel}, we analyse the baseband discrete-time PN channel. We describe a new  discrete-time model for the PN, which is given in Theorem~\ref{thm:AR model}. Then we analyze the error due to the adoption of the discrete-time channel model and quantify the power loss and the ISI due to the presence of PN. Section~\ref{s:results} collects some simulation results and Section~\ref{s:conclusions} concludes the paper.

\section{Communication system model}\label{s:model}
We consider a wireless communication system, where linearly modulated symbols $x_n$ are transmitted through the channel. At the receiver, after down-conversion, the received signal is
\begin{equation}\label{e:r1}
r(t)= \sum_n x_n p(t-n T_\mathrm{s}) e^{j \theta_\text{T}(t)} e^{j \theta_\text{R}(t)}+ v(t)e^{j \theta_\text{R}(t)}\, ,
\end{equation}
where $T_\mathrm{s}$ is the symbol time, $p(t)$ is the shaping pulse with unitary energy, $v(t)$ is complex additive white Gaussian noise (AWGN) with PSD $2N_0$, and $\theta_\text{T}(t)$ and $\theta_\text{R}(t)$ are the PN processes that arise from local oscillator instabilities at the transmitter and at the receiver, respectively (see Section~\ref{s:PNmodel}).
The model~(\ref{e:r1}) can be rewritten as
\begin{equation}\label{e:r}
r(t)= \sum_n x_n p(t-n T_\mathrm{s}) e^{j \theta(t)} + w(t)\, ,
\end{equation}
where $w(t)$ is the AWGN with same statistics of $v(t)$ and $\theta(t)$ is the sum of the transmitter and receiver PN processes. The PSD of $\theta(t)$ is the sum of the PSDs of $\theta_\text{T}(t)$ and $\theta_\text{R}(t)$, where the two PN processes are independent. In many practical cases, receiver and transmitter oscillators have quite different performance, and hence only the dominant PN can be considered, i.e., the PN of the oscillator adopted in consumer grade equipment. In other cases, they are similar and $\theta(t)$ has the sum level of both PNs. Hence, our PSD model, described in the following section, is applicable in both cases.

At the receiver side, a filter matched to the shaping pulse $p(t)$ is employed, followed by a sampler at symbol time. The received samples are approximate sufficient statistics and are given by
\begin{equation}\label{e:rec}
y_k=r(t)\otimes p^*(-t)_{|_{t=k T_\mathrm{s}}}\, .
\end{equation}
The considered system model is shown in Figure~\ref{f:blockdiagram}.
The PN can be tracked/estimated at the receiver following several approaches, as discussed in~\cite{Co14}.

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/blockdiagram-paper}
\caption{Block diagram of the system model.} \label{f:blockdiagram}
\end{figure}

\section{Proposed statistical model of the phase noise}\label{s:PNmodel}
This section introduces the proposed statistical model of the PN. We first describe an analytical form for the PN PSD in the case of free-running and PLL oscillators, then we derive the corresponding PSD of the phasor.
\subsection{Power spectral density of the PN}
The single-side-band PN spectrum of many practical oscillators that can be found from measurements is characterized by a $-20$~dB/decade slope~\cite{TuCo05,HaDrHa98,Mu95} due to integration of white noise inside the oscillator circuitry, and by two flat parts: one at low frequencies, representing the attenuation of the PN stabilized by means of a PLL, the other at high frequencies, modelling the thermal noise at the oscillator output~\cite{De06}. Therefore, we propose to model the PN PSD as
\begin{equation}\label{e:psd}
S_\theta(f)=\frac{f^2_{\mathrm{3dB}}\ell^2_{0}}{f^2_{\mathrm{3dB}}+f^2}+\ell^2_\infty\, ,
\end{equation}
where $f_{\mathrm{3dB}}$ is the 3dB bandwidth, and $\ell^2_{0}$ and $\ell^2_\infty$ are the spectrum levels for $f\!=\!0$ and for high frequency, respectively.
In real oscillators, $\ell_\infty$ is very low and we can safely assume that $\ell_{0}\gg\ell_\infty$.
An example of such PSD is shown in Figure~\ref{f:psd}, where $\ell^2_{100}$ is the spectrum level at 100~kHz. In Appendix~\ref{a:mesures}, we report some experimental measurements of some free running and PLL-locked oscillators.

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/PSD}
\caption{Power spectral density of the phase noise, corresponding to equation~(\ref{e:psd}).  $\ell^2_{0}$ and $\ell^2_{100}$ are  the spectrum levels at 0 Hz and at 100 kHz, respectively. $\ell_\infty$ is the spectrum at high frequencies.} \label{f:psd}
\end{figure}

Oscillators are usually quantified in terms of the PN spectrum level in the $-20$~dB/decade region. For example, the spectrum level at $100$~kHz, $\ell^2_{100}$, is a common calibration point of oscillators, and
represents a more interesting quantity for the designers than the spectrum level at $0$~Hz~\cite{datasheet}. When $f_{\mathrm{3dB}}$ is much smaller than the frequency chosen as reference, in this case $100$~kHz, starting from (\ref{e:psd}), we have the following relationships
\begin{align}\label{e:app1}
\ell^2_{0}&=S_\theta(f)\frac{f^2+f^2_{\mathrm{3dB}}}{f^2_{\mathrm{3dB}}}-\ell^2_\infty\simeq\\\label{e:app2}
&\simeq\ell^2_{100}\frac{10^{10}}{f^2_{\mathrm{3dB}}}-\ell^2_\infty\simeq\\\label{e:approx1_l100}
&\simeq\ell^2_{100}\frac{10^{10}}{f^2_{\mathrm{3dB}}}\, ,
\end{align}
where~(\ref{e:app2}) is obtained by writing~(\ref{e:app1}) for $f=100$~kHz and using that $f_{\mathrm{3dB}}\ll100$~kHz,\footnote{When the condition $f_{\mathrm{3dB}} \ll 100$~kHz does not hold, we can simply choose to express the PSD as a function of another spectrum level in the $-20$~dB/decade region.} while~(\ref{e:approx1_l100}) is obtained by assuming that $\ell^2_\infty$ is negligible with respect to the first term in~(\ref{e:app2}). By using~(\ref{e:approx1_l100}), we can express the PSD of the PN as a function of $\ell_{100}$, that is
\begin{equation}\label{e:psd_l100}
S_\theta(f)=\frac{10^{10}\ell^2_{100}}{f^2_{\mathrm{3dB}}+f^2}+\ell^2_\infty\, .
\end{equation}
The flat part of the PSD at high frequencies dominates when $f> \sqrt{ \frac{10^{10}\ell^2_{100}}{\ell^2_\infty} -f_{\mathrm{3dB}}^2}$.
In order to simplify the exposure, in the following we assume $\ell_\infty=0$, unless otherwise specified. The extension of our results to the case where we consider that $\ell_\infty\neq0$ is reported in Appendices~\ref{a:linfinity} and~\ref{a:linfinity2}.

By taking the inverse Fourier transform of~(\ref{e:psd_l100}) with $\ell_\infty=0$, we obtain the autocorrelation function of the PN, that is
\begin{equation}\label{e:corr}
R_\theta(\tau)=\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}  e^{-2\pi f_{\mathrm{3dB}}  |\tau|}\, .
\end{equation}

In this work, we consider two limiting cases, i.e., free-running oscillator and PLL-locked oscillator. In particular, we study the case \mbox{$f_{\mathrm{3dB}} \rightarrow 0$}, which corresponds to a free-running oscillator. In this case, the PN process is a nonstationary Wiener process with a variance that increases linearly with time; on the other hand, the PN has stationary increments and can be described through the variance of the phase increments. The other case is the case of high $f_{\mathrm{3dB}}$, i.e., \mbox{$f_{\mathrm{3dB}}\gg \pi 10^{10} \ell^2_{100}$},  which can represent a PLL-locked oscillator with 3dB-bandwidth $f_{\mathrm{3dB}}$.
This condition is easily met in practical PLL-locked oscillators. In fact, the value of the spectrum at 100 kHz offset, $\ell^2_{100}$,  is most often between -110 and -80 dB, therefore the quantity $\pi 10^{10} \ell^2_{100}$ is in the range $[ 0.3, 300]$. Since the PLL bandwidth is typically in the order of $1$~kHz or 1-2 decades higher than this value, we can conclude that the condition \mbox{$f_{\mathrm{3dB}}\gg \pi 10^{10} \ell^2_{100}$} is normally satisfied for PLL-locked oscillators. We remark that this condition and the condition $f_{\mathrm{3dB}}\ll100$~kHz that we used in~(\ref{e:app2}) are usually met at the same time by common values of $f_{\mathrm{3dB}}$.\footnote{See footnote 1.} The range of values of parameters of typical oscillators are reported in Table~\ref{t:parameters}.
\begin{table}[ht]
\centering
\begin{tabular}{|c|c|}
\hline
parameter& range  \\ \hline
\hline
$f_{\mathrm{3dB}}$   & $[10^3,10^5]$ Hz  \\ \hline
$\ell^2_{100}$     & $[-110,-80]$ dB  \\ \hline
$\ell^2_\infty$     & $<-110$ dB \\ \hline
$\pi 10^{10} \ell^2_{100}$ & $[0.3,300]$ Hz  \\ \hline
\end{tabular}
\vspace{0.2cm}
\caption{Parameters of typical oscillators. Parameter $f_{\mathrm{3dB}}$ is given for PLL-locked oscillators.}
\label{t:parameters}
\end{table}

Our model can be easily extended to match more involved PSDs by assuming that the PN is given by the sum of independent processes, each described by~(\ref{e:psd_l100}). For example, in the case of PN given by the independent contribution of two processes, we get
\begin{equation}\label{e:extension}
    S_\theta(f)=S'_\theta(f)+S''_\theta(f)\, ,
\end{equation}
where both $S'_\theta(f)$ and $S''_\theta(f)$ are given by~(\ref{e:psd_l100}), with suitable different parameters. We use this extension to model a more involved spectrum in Section~\ref{s:results}, where we model also the PN of the reference oscillator in a PLL system (see Appendix~\ref{a:mesures}).

As a final remark, we note that our proposed model is similar to the piece-wise-linear spectrum models in~\cite{TuCo05,HaDrHa98,Mu95,RoKa95}, which strengthens the proposal; also other authors conclude that a similar model is useful for modeling of the PN. The main advantage with our proposal is the possibility to use it to mathematically derive properties of the PN and oscillator, and to derive a time-domain model with a spectrum as in~(\ref{e:psd_l100}), which we do in the following.
\subsection{Power spectral density of the phasor}\label{ss:phasorPSD}
Since both PN and phasor PSDs are used by industries and studied in the literature~\cite{PiMa02,BaLe90,GhKr17}, we now derive the PSD of the random process of the phasor $h(t)\triangleq e^{j\theta(t)}$.
Through simulations, it can be observed that the PSD of the phasor follows the one of the PN at high frequencies. In the literature, see for example~\cite{PiMa02}, this is motivated by using the low-PN approximation, i.e.,
\begin{equation}
h(t)\simeq1+j \theta(t)
\end{equation}
when $\theta(t)$ is small. However, this approximation is not always met (quite often it is not met in practice, actually), and we derive more exact expressions below.

Here, we derive the PSD of the phasor starting from the more general  model of the PN in~(\ref{e:psd_l100}), without resorting to the low PN approximation.
We first assume that the flat part of the PN PSD at high frequency is negligible and hence we set $\ell^2_\infty=0$.  The extension to the case where we consider that $\ell_\infty\neq0$ is reported in Appendix~\ref{a:linfinity}.
We start by expressing the autocorrelation of the phasor in Lemma~\ref{thm:phasorautocor}, and then get the spectrum for the two relevant cases of free-running and PLL-locked oscillators by Theorem~\ref{thm:phasorPSD}.
\begin{lemma}\label{thm:phasorautocor}
Consider the phase process $\theta(t)$ with PSD
\begin{equation}\nonumber
S_\theta(f)=\frac{10^{10}\ell^2_{100}}{f^2_{\mathrm{3dB}}+f^2}\, .
\end{equation}
The phasor process $e^{j\theta(t)}$ has autocorrelation
\begin{equation}\label{e:corr_phasor}
R_h(\tau)=e^{- \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}   (1- e^{-2\pi f_{\mathrm{3dB}}  |\tau|} )}\,.
\end{equation}
\end{lemma}
\begin{IEEEproof}
See Appendix~\ref{a:thmphasorautocor}
\end{IEEEproof}
\begin{theorem}\label{thm:phasorPSD}
Consider the phase process $\theta(t)$ with PSD
\begin{equation}\nonumber
S_\theta(f)=\frac{10^{10}\ell^2_{100}}{f^2_{\mathrm{3dB}}+f^2}\, .
\end{equation}
The phasor process $e^{j\theta(t)}$ has power spectral density
\begin{align}
&S_{h} (f)=e^{-\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}}  \delta(f) + \nonumber\\
&\;\;\;e^{-\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}}  \! \! \int_{-\infty}^{\infty}  \!\!\!\!  (  e^{     \frac{  \pi 10^{10} \ell^2_{100} e^{-2 \pi f_{\mathrm{3dB}} |\tau | } }{f_{\mathrm{3dB}}}   } \!\!\! -\!\!1 \!   ) \mathrm{cos}(2\pi f \tau) \mathrm{d} \tau\, ,
\end{align}
where $\delta(f) $ denotes the Dirac delta function.
Moreover, for \mbox{$f_{\mathrm{3dB}} \rightarrow 0$} we have
\begin{equation}\label{e:free_phasor_spectrum}
S_h(f)_{f_{\mathrm{3dB}} \to 0}  =  \frac{ 10^{10}\ell^2_{100}}{\pi^2 10^{20} \ell_{100}^4  + f^2}\, ,
\end{equation}
and for $f_{\mathrm{3dB}}\gg \pi 10^{10} \ell^2_{100}$
\begin{equation}\label{e:pll_phasor_spectrum}
S_h(f)_{\mathrm{high}f_{\mathrm{3dB}} }
=\Big( 1- \frac{\pi 10^{10} \ell^2_{100}  }{f_{\mathrm{3dB}}} \Big) \delta(f)   + \frac{  10^{10} \ell^2_{100} }{f^2_{\mathrm{3dB}}+f^2}\, .
\end{equation}
\end{theorem}
\begin{IEEEproof}
The PSD of the phasor process is obtained by taking the Fourier transform of the autocorrelation function of the phasor, given in Lemma~\ref{thm:phasorautocor}, that is
\begin{align}\nonumber
&S_{h} (f)=\mathcal{F} \{ e^{- \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}   (1- e^{-2\pi f_{\mathrm{3dB}}  |\tau|} )} \}=\\\nonumber
&=\mathcal{F} \{ e^{- \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}   (1- e^{-2\pi f_{\mathrm{3dB}}  |\tau|} )} + e^{- \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}} -e^{- \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}} } \} \\\nonumber
&=e^{-\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}}  \delta(f) + \\ \nonumber
&\;\;\;e^{-\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}}  \! \! \int_{-\infty}^{\infty}  \!\!\!\!  (  e^{     \frac{  \pi 10^{10} \ell^2_{100} e^{-2 \pi f_{\mathrm{3dB}} |\tau | } }{f_{\mathrm{3dB}}}   } \!\!\! -\!\!1 \!   ) \text{cos}(2\pi f \tau) d \tau\, .
\end{align}
We now consider the two limiting cases \mbox{$f_{\mathrm{3dB}} \rightarrow 0$} and \mbox{$f_{\mathrm{3dB}}\gg \pi 10^{10} \ell^2_{100}$}.

In the first case, using the approximation $e^{-2\pi f_{\mathrm{3dB}} |t|} \simeq  1 -2\pi f_{\mathrm{3dB}} |t|$, which is valid when $f_{\mathrm{3dB}}\rightarrow 0$, we can write
\begin{equation}\nonumber
S_h(f)_{f_{\mathrm{3dB}} \to 0}  \!= \! \mathcal{F} \{ e^{-2\pi^2 10^{10}\ell^2_{100} |\tau| }  \} \!\! =\!\!  \frac{ 10^{10}\ell^2_{100}}{\pi^2 10^{20} \ell_{100}^4  + f^2}\, .
\end{equation}

When $f_{\mathrm{3dB}}\gg \pi 10^{10} \ell^2_{100}$, the autocorrelation of the phasor can be approximated as
\begin{equation}
R_h(\tau)\simeq 1- \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}   (1- e^{-2\pi f_{\mathrm{3dB}}  |\tau|} )
\end{equation}
and the PSD can be computed as
\begin{align}\nonumber
S_h(f)_{\text{high}f_{\mathrm{3dB}} } & = \mathcal{F} \Big\{ 1- \frac{\pi 10^{10} \ell^2_{100}  }{f_{\mathrm{3dB}}}  + \frac{\pi 10^{10}  \ell^2_{100}  e^{-2\pi f_{\mathrm{3dB}} |\tau|} }{f_{\mathrm{3dB}}}  \Big\}\\ \nonumber
&=\Big( 1- \frac{\pi 10^{10} \ell^2_{100}  }{f_{\mathrm{3dB}}} \Big) \delta(f)   + \frac{  10^{10} \ell^2_{100} }{f^2_{\mathrm{3dB}}+f^2}\, .
\end{align}
\end{IEEEproof}

Theorem~\ref{thm:phasorPSD} gives the PSD in the case of free-running oscillator in~(\ref{e:free_phasor_spectrum}) and in the case of PLL-locked oscillator in~(\ref{e:pll_phasor_spectrum}). The PSD for free-running oscillator is a Lorentzian spectrum, with the value at 100~kHz given by $\ell^2_{100}$ and \mbox{3-dB} bandwidth equal to $\pi10^{10}\ell^2_{100}$, and overlaps with~(\ref{e:psd_l100}) for $f\gg \pi 10^{10} \ell^2_{100}$. In the case of PLL-locked oscillator, the spectrum has a delta function at $f\!=\!0$, and then follows exactly the PN spectrum. As expected, the power of the phasor computed by integrating the spectrum in~(\ref{e:free_phasor_spectrum}) or in~(\ref{e:pll_phasor_spectrum}) is unitary.

\section{Baseband discrete-time phase noise channel}\label{s:DT channel}
In this section, we consider the discrete-time channel affected by PN and AWGN, representative of the sample in~(\ref{e:rec}), obtained after matched filtering and sampling of the received signal. A commonly studied discrete-time channel model is~\cite{CoBaCa05,BaCo07b}
\begin{equation}\label{e:discretetime}
z_k=x_k e^{j \theta_k}  + w_k \, ,
\end{equation}
where $\theta_k$ is the PN sample obtained as $\theta_k =\theta(k T_\mathrm{s})$ and $w_k$ is the sample obtained by filtering and sampling the AWGN signal $w(t)$. Signals $\theta(t)$ and $w(t)$ are introduced in~(\ref{e:r}).
Model~(\ref{e:discretetime}) neglects the ISI and the power loss due to the presence of PN. In fact, the PN causes an enlargement of the signal bandwidth, therefore the filtering with the pulse matched to the shaping pulse $p(t)$ cuts some spectral components of the signal, resulting in a small power loss. Moreover, the received samples are affected by ISI due to mismatched filters. We quantify these effects in subsection~\ref{ss:error}. In the following subsection, we propose a model for the discrete-time PN process $\{\theta_k\}$. Our model is valid for common oscillators and symbol rates since it requires that $f_{\mathrm{3dB}} T_\mathrm{s} \ll 1$. Considering that $T_\mathrm{s}$ of practical systems is usually lower than $10^{-5}$ and the typical range of $f_{\mathrm{3dB}}$ reported in Table~\ref{t:parameters}, we can conclude that the above condition is easily met.

\subsection{Discrete-time phase noise model}\label{ss:DT model}
Here, we propose a discrete-time model for the PN. More precisely, we propose a model representing the PN sampled at symbol rate, whose PSD follows the  model~(\ref{e:psd_l100}) when $\ell^2_\infty=0$ in the signal band $ [-1/2T_\mathrm{s}, 1/2T_\mathrm{s}]$. The proposed model can describe the process $\theta_k=\theta(k T_s)$ to be used to simulate the simplified PN channel~(\ref{e:discretetime}). Again, we consider the case $\ell^2_\infty=0$, while we account for the case $\ell^2_\infty\neq0$ in Appendix~\ref{a:linfinity2}.

In the following theorem, we show that the process $\theta_k$ can be described by the autoregressive (AR) process.
\begin{theorem}\label{thm:AR model}
The discrete-time AR random process
\begin{equation}\label{e:dt model}
\theta_k=a\theta_{k-1}+u_k\, ,
\end{equation}
where $u_k$ are i.i.d. zero-mean Gaussian variables with variance $\sigma^2_u$, and parameters $a$ and $\sigma^2_u$ are
\begin{align}\label{e:a}
a&=e^{-2\pi f_{\mathrm{3dB}} T_\mathrm{s}}\\\label{e:s}
\sigma^2_u&=\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}} (1-e^{-4\pi f_{\mathrm{3dB}} T_\mathrm{s}})\, ,
\end{align}
has the following PSD
\begin{equation}\label{e:psdphi}
S_\theta(f)=\frac{10^{10}\ell^2_{100}}{f^2_{\mathrm{3dB}}+f^2}\,
\end{equation}
in the band $ [-1/2T_\mathrm{s}, 1/2T_\mathrm{s}]$, when $f_{\mathrm{3dB}} T_s\ll 1$.
\end{theorem}
\begin{IEEEproof}
See Appendix~\ref{a:thmARmodel}.
\end{IEEEproof}
The poof of the theorem uses the following Lemma, which quantifies the variance of the aliasing error due to the infinite bandwidth of the PN.
\begin{lemma}\label{l:aliasing}
Consider the discrete-time process $\theta_k$ obtained by sampling with sampling time $T_\mathrm{s}$ the process $\theta(t)$ with PSD
\begin{equation}\label{e:psdphiL}
S_\theta(f)=\frac{10^{10}\ell^2_{100}}{f^2_{\mathrm{3dB}}+f^2}\, .
\end{equation}
The aliasing affecting the PSD of the process $\theta_k$ in the band $[-1/2T_\mathrm{s}, 1/2T_\mathrm{s}]$ has variance
\begin{equation}\label{e:alias}
    \sigma^2_\mathrm{alias}= \frac{\pi 10^{10}\ell^2_{100}}{f_{\mathrm{3dB}}}\Big(1 - \frac{2}{\pi} \mathrm{tan}^{-1} \Big( \frac{1}{2 \Tsf_{\mathrm{3dB}}}\Big)\Big)\, .
\end{equation}
\end{lemma}
\begin{IEEEproof}
 The variance of the aliasing can be computed by integrating the PSD~(\ref{e:psdphiL}) of the continuous-time process $\theta(t)$ outside the band $[-1/2T_\mathrm{s}, 1/2T_\mathrm{s}]$, that is
\begin{align}\nonumber
    \sigma^2_\mathrm{alias}&=  2\int_{\frac{1}{2T_\mathrm{s}}}^\infty  \frac{10^{10}\ell^2_{100}}{f^2_{\mathrm{3dB}}+f^2} \, df\\ \nonumber
&= \frac{\pi 10^{10}\ell^2_{100}}{f_{\mathrm{3dB}}}\Big(1 - \frac{2}{\pi} \text{tan}^{-1} \Big( \frac{1}{2 \Tsf_{\mathrm{3dB}}}\Big)\Big)\, .
\end{align}
\end{IEEEproof}

The variance of the error due to the alias, given by~(\ref{e:alias}), normalized to the PN power inside the considered  bandwidth $P$, obtained by integrating~(\ref{e:psdphi}) in $[-1/2T_\mathrm{s}, 1/2T_\mathrm{s}]$ is given by
\begin{equation}
    \sigma^2_{\text{alias}}/P=\frac{\pi}{2}\frac{1}{\text{tan}^{-1}( \frac{1}{2 \Tsf_{\mathrm{3dB}}})} -1\, .
\end{equation}
In Figure~\ref{f:alias}, we report $\sigma^2_{\text{alias}}/P$ as a function of the term $f_{\mathrm{3dB}} T_s$. In the figure, we observe that  the error increases with $f_{\mathrm{3dB}} T_s$ and approaches $0$ for $f_{\mathrm{3dB}} T_s\to 0$.

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/ealias}
\caption{Normalized variance of the error due to the alias $\sigma^2_{\text{alias}}/P$ as a function of $f_{\mathrm{3dB}} T_s$.}\label{f:alias}
\end{figure}

\begin{comment}
A discrete time-domain model for the PN whose  PSD follows the one in equation~(\ref{e:psd_l100}) in the signal bandwidth is
\begin{equation}\label{e:dt model1}
\theta_k=\psi_k+v_k\, ,
\end{equation}
where $\psi_k$ is the random process representing the close-to-carrier PN, and $v_k$ is a white Gaussian process, independent of $\psi_k$, with variance $\sigma^2_v$.

Hence, the autocorrelation of the PN samples is
\begin{align}
R_\theta(n)&=\text{E} \{ \theta_k \theta_{k-n} \} \\
&=  \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}} \, e^{-2\pi f_{\mathrm{3dB}} T_\mathrm{s} |n|}   +\frac{\ell^2_\infty}{T_\mathrm{s}} \delta(n)\, ,
\end{align}
where $\delta(n)$ is the Kronecker delta.
\end{comment}

The case of a free-running oscillator implies that $f_{\mathrm{3dB}}\to 0$. In this case, the parameter $a$ is equal to 1 and the model~(\ref{e:dt model}) reduces to the Wiener (random-walk) model. The Wiener process is nonstationary, but the process of the phase increment $u_k$ is stationary and the PN can be described through the variance $\sigma^2_u$. We can compute the variance of the phase increments by taking the limit for $f_{\mathrm{3dB}}\rightarrow 0$ and we obtain
\begin{equation} \label{e:innovation variance}
\sigma^2_u= 4\pi^2 10^{10}\ell^2_{100} T_\mathrm{s} \, .
\end{equation}
As observed in Subection~\ref{ss:phasorPSD}, the \mbox{3-dB} bandwidth of the oscillator PSD is $\pi 10^{10} \ell^2_{100}$, therefore the variance $\sigma^2_u$ can be expressed as
\begin{equation} \label{e:innovation variance rho}
\sigma^2_u= 4\pi\rho \, ,
\end{equation}
where $\rho \triangleq \pi 10^{10} \ell^2_{100} T_\mathrm{s}$ represents the ratio between the phasor bandwidth and the signal bandwidth.

This analysis tells us that if we want to model the PN whose spectrum is available, we can generate the PN samples according to the model (\ref{e:dt model}), with $\ell_{100}$ and $f_{\mathrm{3dB}}$ taken from measurements.
For example, in the case of a free-running oscillator, we generate the PN as a Wiener model with the variance of the phase increments given by~(\ref{e:innovation variance}), where $\ell_{100}$ is simply the value of the measured spectrum at 100~kHz.

\subsection{Error analysis}\label{ss:error}
Here, we evaluate the accuracy of the communication model~(\ref{e:discretetime}). Thanks to the proposed PN model described in Section~\ref{s:PNmodel}, we can connect results of this analysis to the main PN parameters of practical oscillators. In particular we relate the following mean square error to the parameters of the PN
\begin{equation}\label{e:eta1}
\eta=E_{x,\theta,w}\Big \{ |z_k-y_k|^2\Big \}\, ,
\end{equation}
where $z_k$ is the received sample of the simplified model~(\ref{e:discretetime}) ignoring ISI and power loss, while $y_k$ is the received sample of the more accurate model~(\ref{e:rec}).
\begin{lemma}\label{l:eta}
Consider the system model described in Section~\ref{s:model}, and the phase noise $\theta(t)$ described in Section~\ref{s:PNmodel} in the case of a free-running oscillator. Let $e^{j \theta_k}$ be the sample $e^{j \theta(k T_s)}$ and $g_{\ell,k}=p(t)e^{j\theta(t+(k-\ell) T_s)}\otimes p^*(-t)|_{{t=\ell T_s}}$. The mean square error in~(\ref{e:eta1}) between $z_k$ defined in~(\ref{e:discretetime}) and $y_k$ defined in~(\ref{e:rec}) is given by
\begin{equation}\label{e:eta}
\eta=  \eta_\mathrm{D} + \eta_\mathrm{ISI}  \, ,
\end{equation}
where
\begin{align}\label{e:etad}
\eta_\mathrm{D}=&    E_{\theta}\Big \{ |  e^{j \theta_k} -  g_{0,k}|^2\Big \}    \, , \\ \label{e:etaisi}
\eta_\mathrm{ISI}=&   \sum_{\ell\neq 0}  E_{\theta}\Big \{ | g_{\ell,k-\ell}  |^2\Big \}    \, .  &
\end{align}
\end{lemma}
\begin{IEEEproof}
We expand~(\ref{e:rec}) and obtain
\begin{align}\nonumber
y_k&=\sum_n x_n \int_{-\infty}^{+\infty}   p(z-nT_s) e^{j\theta(z)} p^*(z-kT_s ) dz  +w_k\\ \nonumber
&=\sum_n x_n \!\! \int_{-\infty}^{+\infty}  \!\!\!  p(z) e^{j\theta(z+nT_s)} p^*(z-(k-n)T_s) dz  +w_k\\
&=\sum_\ell x_{k-\ell} g_{\ell,k-\ell}    +w_k\, ,\label{e:rec_k}
\end{align}
where we have defined
\begin{align}
g_{\ell,k-\ell}&\triangleq \int_{-\infty}^{+\infty}  p(z) e^{j\theta(z+(k-\ell)T_s)} p^*(z-\ell T_s ) dz \label{e:gkl}  \\
&= p(t)e^{j\theta(t+(k-\ell )T_s)}\otimes p^*(-t)|_{t=\ell T_s}\, ,
\end{align}
that represents a linear time-varying filter that describes the power loss and ISI.
Rewriting~(\ref{e:rec_k}), we get
\begin{equation}\label{e:y_k}
y_k= x_k g_{0,k} + \sum_{\ell\neq 0} x_{k-\ell} g_{\ell,k-\ell}   +w_k \, ,
\end{equation}
where $g_{0,k}$ is connected to the power loss and the summation represents the ISI.
Using (\ref{e:discretetime}) and (\ref{e:y_k}), we obtain
\begin{equation}
\eta=E_{x,\theta}\Big \{ |  x_k e^{j \theta_k} -  x_k g_{0,k} - \sum_{\ell\neq 0} x_{k-\ell} g_{\ell,k-\ell}      |^2\Big \}  \, .
\end{equation}
Assuming that the transmitted symbols are zero mean and i.i.d., with unitary energy, we obtain that $\eta$ is the sum of two errors, one due to the approximation of $g_{0,k}$ with the PN sampled in $kT_s$, the other due to the ISI, that is
\begin{equation}
\eta=  E_{\theta}\Big \{ |  e^{j \theta_k} -  g_{0,k}|^2\Big \}  + \sum_{\ell\neq 0}  E_{\theta}\Big \{ | g_{\ell,k-\ell}  |^2\Big \} \, .
\end{equation}
\begin{comment}
\begin{align}
\eta_\text{D}=&    E_{\theta}\Big \{ |  e^{j \theta_k} -  g_{0,k}|^2\Big \}    \, , \\
\eta_\text{ISI}=&   \sum_{\ell\neq 0}  E_{\theta}\Big \{ | g_{\ell,k-\ell}  |^2\Big \}    \, .  &
\end{align}
\end{comment}
\end{IEEEproof}
As expected, in the absence of PN and for Nyquist shaping pulse, we have that $g_{0,k}=1$ and $g_{\ell,k}=0$ for $\ell\neq0$. For convenience, we define the terms
\begin{equation}\label{e:gammal}
    \gamma_\ell= E_{\theta}  \{  | g_{\ell,k-\ell} |^2  \}\, ,
\end{equation}
that using~(\ref{e:gkl}) can be written as
\begin{align}\nonumber
\gamma_\ell& =E_{\theta}\Big \{  \!\! \int_{-\infty}^{+\infty}   \int_{-\infty}^{+\infty}   p(z)p^*(z_1)p^*(z \!-\! \ell T_s) p( z_1 \!\! - \!\!   \ell T_s)\cdot \\ \nonumber
& \cdot e^{j \theta(z+(k-\ell) T_s )}  e^{-j \theta(z_1+(k-\ell) T_s )}  d z d z_1 \Big\}\\
&\!= \!\! \int_{-\infty}^{+\infty}   \! \!\! \int_{-\infty}^{+\infty}   \!\!\! \! p(z)p^*\!(z_1)p^*\!(z \!-\! \ell T_s) p( z_1 \!\! - \!\!   \ell T_s) \!R_h(z_1\!\!-\!\!z) d z d z_1, \label{e:gamma_l}
\end{align}
where we recall that $R_h(\tau)$ is the autocorrelation of the phasor $h(t)=e^{j \theta(t)}$. Equation~(\ref{e:gamma_l}) shows that the terms $\gamma_\ell $ are independent of $k$. We can express the SIR when a detector that is not able to account for the ISI is employed, as a function of the terms $\gamma_\ell$.
Using the assumption of zero mean symbols, we have
\begin{equation}
\text{SIR}=  E_{x,\theta} \Big \{ |x_k g_{0,k}|^2  \Big\}/ E_{x,\theta} \Big\{  \Big| \sum_{\ell\neq 0} x_{k-\ell} g_{\ell,k-\ell} \Big|^2    \Big\}\, ,
\end{equation}
and considering that the symbols are independent, with unitary variance, we obtain
\begin{equation}\label{e:sir}
\text{SIR}=\frac{\gamma_0}{ \sum_{\ell\neq 0} \gamma_\ell}\, .
\end{equation}

In the following, to simplify the exposure, we assume that the filter $p(t)$ is given by
\begin{equation}\label{e:sinc}
p(t)=\frac{1}{\sqrt{T_s}} \text{sinc}\Big(\frac{t}{T_s}  \Big)
\end{equation}
with Fourier transform
\begin{equation}\label{e:rect}
P(f)=\sqrt{T_s} \Pi (f T_s )\, .
\end{equation}
In the above equations, sinc is the normalized sinc function, and $\Pi(x)$ is the rectangular function equal to one for \mbox{$-1/2<x<1/2$} and zero otherwise.
Moreover, we consider the case of a free-running oscillator, for which the PSD of the phasor is given by~(\ref{e:free_phasor_spectrum}). The extension of the following results to the case of the PLL-locked oscillator is straightforward.

The following theorem gives the expression of the mean square errors in Lemma~\ref{l:eta} as a function of the PN parameters.
\begin{theorem}\label{thm:eta}
Consider the system model described in Section~\ref{s:model}, where $p(t)$ is the normalized sinc function with unitary energy, and the phase noise $\theta(t)$ is described in Section~\ref{s:PNmodel} in the case of a free-running oscillator. The mean square errors $\eta$, defined in~(\ref{e:eta1}), $\eta_\mathrm{D}$ and $\eta_\mathrm{ISI}$, defined in Lemma~\ref{l:eta}, are given by
\begin{equation}\label{e:etarho}
    \eta= 1+\frac{2}{\pi}\Big[ \frac{\rho }{2} \mathrm{log}\Big(1+\frac{1}{\rho ^2}\Big) -\mathrm{tan}^{-1}   \Big( \frac{1}{\rho } \Big)  \Big]\, ,
\end{equation}
\begin{equation}\label{e:etaDfinal}
\eta_\mathrm{D}=1+\frac{2}{\pi}  \Big[  \rho  -(1+\rho ^2)     \mathrm{tan}^{-1}   \Big( \frac{1}{\rho } \Big) \Big]\,,
\end{equation}
\begin{equation}\label{e:etaisifinal}
\eta_\mathrm{ISI}= \frac{2}{\pi}\Big[ \rho ^2 \mathrm{tan}^{-1}   \Big( \frac{1}{\rho } \Big)  +  \frac{\rho }{2} \mathrm{log}\Big(1+\frac{1}{\rho ^2}\Big)       -\rho \Big]\, ,
\end{equation}
where $\rho=\pi 10^{10} \ell^2_{100} T_s$.
\end{theorem}
\begin{IEEEproof}
See Appendix~\ref{a:thmeta}.
\end{IEEEproof}

\begin{comment}
\begin{theorem}\label{thm:powerloss}
Consider the phase noise $\theta(t)$ described in Section~\ref{s:PNmodel} in the case of a free-running oscillator. Let $e^{j \theta_k}$ be the sample $e^{j \theta(k T_s)}$ and $g_{0,k}=p(t)e^{j\theta(t+k T_s)}\otimes p^*(-t)_{|_{t=0}}$, where $p(t)$ is the normalized sinc function with unitary energy.
The mean square error between the phasor $e^{j \theta_k}$ and $g_{0,k}$ is given by
\begin{equation}\label{e:etaDfinal}
\eta_\text{D}=1+\frac{2}{\pi}  \Big[  \rho  -(1+\rho ^2)     \mathrm{tan}^{-1}   \Big( \frac{1}{\rho } \Big) \Big]\,,
\end{equation}
where $\rho=\pi 10^{10} \ell^2_{100} T_s$.
\end{theorem}
\begin{IEEEproof}
See Appendix~\ref{a:thmpowerloss}.
\end{IEEEproof}

\begin{theorem}\label{thm:isi}
Consider the system model described in Section~\ref{s:model}, and the phase noise $\theta(t)$ described in Section~\ref{s:PNmodel} in the case of a free-running oscillator. Let $g_{\ell,k}=p(t)e^{j\theta(t+(k-\ell) T_s)}\otimes p^*(-t)_{|_{\ell T_s}}$, where $p(t)$ is the normalized sinc function with unitary energy.
The sum for $\ell \neq 0$ of the terms $|g_{\ell,k}|^2$ is given by
\begin{equation}\label{e:etaisifinal}
\eta_\text{ISI}= \frac{2}{\pi}\Big[ \rho ^2 \mathrm{tan}^{-1}   \Big( \frac{1}{\rho } \Big)  +  \frac{\rho }{2} \text{log}\Big(1+\frac{1}{\rho ^2}\Big)       -\rho \Big]\, ,
\end{equation}
where $\rho=\pi 10^{10} \ell^2_{100} T_s$.
\end{theorem}
\begin{IEEEproof}
See appendix~\ref{a:thmisi}.
\end{IEEEproof}
\end{comment}
From the expressions above, we observe that the error made when using the channel model~(\ref{e:discretetime}) depends only on the parameter $\rho$, that is the ratio between the oscillator bandwidth and the signal bandwidth. Considering the limits, we have
\begin{align}
\lim_{\rho \to 0} \eta&=0\\
\lim_{\rho\to \infty} \eta&=1\, .
\end{align}
In the case of sinc shaping pulse, we can express the SIR as a function of the phase increment standard deviation $\sigma_u$. Starting from~(\ref{e:sir}), we use the expression of $\gamma_0$ in~(\ref{e:gamma0}) in the numerator, while the denominator is $\eta_\text{ISI}$. Then, by using equation~(\ref{e:innovation variance rho}) which connects $\rho $ with the standard deviation of the phase increments, we obtain
\begin{equation}\label{e:sir_sigma}
\text{SIR}=\frac{   \text{tan}^{-1}   \Big( \frac{4 \pi}{\sigma_u^2} \Big) \Big[ 1 -\frac{\sigma_u^4}{16\pi^2} \Big] -\frac{\sigma_u^2}{4\pi} \text{log}\Big(1+\frac{16\pi^2}{\sigma_u^4}\Big) + \frac{\sigma_u^2}{4\pi} }{ \frac{\sigma_u^4}{16\pi^2} \mathrm{tan}^{-1}   \Big( \frac{4 \pi}{\sigma_u^2}   \Big)  +  \frac{\sigma_u^2}{8 \pi}  \mathrm{log}\Big(1+\frac{16\pi^2}{\sigma_u^4}\Big)       -\frac{\sigma_u^2}{4 \pi} }\, .
\end{equation}
From the above equation, we can compute a limit on $\sigma_u$ to have a SIR higher than a given value. For example, a SIR higher than $25$~dB requires $\sigma_u<0.1$ rad.

\begin{comment}
\begin{figure}
\includegraphics[width=\columnwidth]{PNoverleaf/Figures/Eta_vs_ro.pdf}
\caption{Mean square error $\eta$ related to the use of the discrete-time phase noise channel~(\ref{e:discretetime}) as a function of $\rho $, the ratio between oscillator and signal bandwidth.}\label{f:1}
\end{figure}

\begin{figure}
\includegraphics[width=\columnwidth]{PNoverleaf/Figures/Eta_vs_ro_zoom.pdf}
\caption{Mean square error $\eta$ as a function of $\rho $, for $\rho \in[5\cdot 10^{-4},5\cdot 10^{-3} ]$.}\label{f:2}
\end{figure}
\end{comment}

\section{Results}\label{s:results}
In this section, we show numerical results obtained through computer simulations to demonstrate the usefulness of the proposed analytical PN models. First, we show some results to validate the analysis of the discrete-time PN channel, reported in Section~\ref{s:DT channel}.
Then, as a case study, we consider a DVB-S2 system. We show how the proposed model for the PN PSD described in Section~\ref{s:PNmodel} can be used to represent the PN typical of such systems and compare the proposed discrete-time model of Section~\ref{s:DT channel} with the continuous-time counterpart in terms of bit error rate (BER). Finally, we use the extension of our model given by~(\ref{e:extension}) to describe the PN of next-generation cellular systems operating in millimeter-wave frequency bands.

\subsection{Validation of the error analysis}
We consider the transmission of linearly modulated QPSK symbols, with a root-raised cosine (RRC) filter $p(t)$ with roll-off factor $0.05$, $0.1$ and $0.5$. We simulate a continuous-time communication system by using an oversampling factor $5$. We add the PN of a free-running oscillator with $\ell_\infty=0$ in the continuous-time domain, then we compute the SIR after the matched filter and the downsampling operation. The simulated SIRs are reported in Figure~\ref{f:SIR} as a function of $\rho$. We remark that, under the assumption of free-running oscillator and $\ell_\infty=0$, the parameter $\rho$ fully describes the PN in our model.
The lowest curve in Figure~\ref{f:SIR} corresponds to the closed-form expression given in~(\ref{e:sir}), that is obtained under the assumption that the filter $p(t)$ is the ideal sinc function (roll-off zero). We observe that the numerical results are in line with the closed-form analysis, since all curves have the same behaviour of the one corresponding to the closed-form, and approach the closed-form as the roll-off decreases.

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/SIR2}
\caption{Signal-to-interference ratio after the matched filter at the receiver as a function of $\rho $, the ratio between oscillator and signal
bandwidth. The different curves are for different roll-off factors of the RRC filter.}\label{f:SIR}
\end{figure}

\subsection{Case study: DVB-S2 sysyem}
We consider the PN typical of satellite communications. In particular, we consider the PN with PSD given in~\cite{DVB-S2X}, table H. The PSD is given by points and the complete PSD is obtained by linear interpolation. Clearly, it is not possible to derive a closed-form expression of the time-domain PN signal based on this information. In Figures~\ref{f:mask1}, we show the DVB-S2 PSD in red and the PSD of the proposed model~(\ref{e:psd_l100}) in blue.

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/PSDcase1}
\caption{Power spectral density of the phase noise of DVB-S2 systems. Comparison with the proposed model.}\label{f:mask1}
\end{figure}
\begin{comment}
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Figures/PSDcase2}
\caption{Power spectral density of the phase noise of DVB-S2 systems. Case~2, DTH scenario and \textit{typical mask}. Comparison with the proposed model.}\label{f:mask2}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Figures/PSDcase3}
\caption{Power spectral density of the phase noise of DVB-S2 systems. Case~3, non DTH scenario. Comparison with the proposed model.}\label{f:mask3}
\end{figure}
\end{comment}

We now want to evaluate our discrete-time PN model. In particular, we compute the BER performance as a function of $E_s/N_0$, being $E_s$ the mean symbol energy. We consider coded \mbox{$16$-QAM} modulations affected by the PN generated according to our model in the scenario described above. We use RRC shaping pulse  with roll-off factor $0.3$. The information bits are encoded through a LDPC code of length $64800$ with rate $3/4$ from~\cite{DVB-S2X}. We consider two different symbol rates, that is $1/T_\mathrm{s}=100$~MBaud and $1/T_\mathrm{s}=10$~MBaud.
At the receiver, the PN tracking is obtained by means of linear interpolation between estimates on pilot fields of length $36$ symbols, inserted every $1476$ information symbols.
In Figures~\ref{f:BERmask1}, we show the BER performance obtained by using the continuous-time model and the discrete-time model for the PN. As in the previous paragraph, continuous-time signals are simulated in the oversampled domain, by using an oversampling factor equal to $5$. In the continuous-time model, the PN is added in the continuous-time domain, before the matched filter at the receiver. On the contrary, the discrete-time system is simulated according to equation~(\ref{e:discretetime}) and the PN samples are generated according to Theorem~\ref{thm:AR model}, extended to the case $\ell_\infty\neq 0$ in Appendix~\ref{a:linfinity2}. In the figure, we also show the curve of the coherent case as reference. The figure shows an excellent match between continuous-time and discrete-time model, revealing that the error due to the ISI, the power loss and the aliasing are small in the considered cases. In Tables~\ref{t:DTvsCTparam}, we report the parameters of the PN model, while in Table~\ref{t:DTvsCTerrors} we report the related errors $\eta$ and $\sigma_{\text{alias}}$ computed with the equations~(\ref{e:etaDfinal}) and (\ref{e:alias}), for the case $1/T_\mathrm{s}=100$~MBaud and $1/T_\mathrm{s}=10$~MBaud. We recall that $\eta$ and $\sigma_{\text{alias}}$ are computed in closed form under the simplified assumption of $\ell_\infty=0$, roll-off factor zero and, only for $\eta$, under the assumption of free-running oscillator. Table~\ref{t:DTvsCTerrors} show that these errors are small for the considered parameters.

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/BER_case1}
\caption{BER performance as a function of $E_s/N_0$. Comparison between the continuous-time (CT) model and the discrete-time (DT) model.}\label{f:BERmask1}
\end{figure}
\begin{comment}
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Figures/BER_case2}
\caption{BER performance as a function of $E_s/N_0$ for case 2. Comparison between the continuous-time (CT) model and the discrete-time (DT) model.}\label{f:BERmask2}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Figures/BER_case3}
\caption{BER performance as a function of $E_s/N_0$ for case 3. Comparison between the continuous-time (CT) model and the discrete-time (DT) model.}\label{f:BERmask3}
\end{figure}
\end{comment}
\begin{table}[ht]
\centering
\begin{tabular}{|c|c|}
\hline
parameter& value\\ \hline
\hline
 $f_{\mathrm{3dB}}$ & $10$~Hz\\ \hline
 $\ell^2_{100}$ & $-88$ dB   \\ \hline
 $\ell^2_\infty$ & $-114$ dB\\ \hline
\end{tabular}
\vspace{0.2cm}
\caption{Parameters of the proposed model to match the DVB-S2 PN.}
\label{t:DTvsCTparam}
\end{table}

\begin{table}[ht]
\centering
\begin{tabular}{|c|c|c|}
\hline
symbol rate & $\eta$  &$\sigma_{\text{alias}}$ \\ \hline
\hline
$10$ MBaud & $4.2\cdot 10^{-5}$ & $1.3\cdot 10^{-6}$ \\ \hline
$100$ MBaud & $4.9\cdot 10^{-6}$ & $1.3\cdot 10^{-7}$ \\ \hline
\end{tabular}
\vspace{0.2cm}
\caption{Errors related to the discrete-time model in the case of DVB-S2 PN.}
\label{t:DTvsCTerrors}
\end{table}

\subsection{Extended PSD model}
We saw in Section~\ref{s:PNmodel} that the proposed PSD model can be extended to describe more involved spectra.
\begin{comment}
\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/PSDcase3extended}
\caption{Power spectral density of the phase noise typical of non DTH services in DVB-S2 systems, case 3. The proposed extended model is obtained as the sum of 2 independent processes.}\label{f:mask3extended}
\end{figure}
\end{comment}
We consider the PN that can be found in the 3GPP documents~\cite{3GPPTR38803,R4-1701165}, whose PSD is given by the following expression
\begin{equation}\label{e:pds3gpp}
	S(f)=\mathrm{PSD0}\frac{\prod_{n=1}^N 1+\left(\frac{f}{f_{z,n}}\right)^{\alpha_{z,n}}}{\prod_{m=1}^M 1+\left(\frac{f}{f_{p,m}}\right)^{\alpha_{p,m}}}\,.
\end{equation}
When the carrier frequency is 45 GHz, the parameters of the above PSD are given in Table~\ref{t:param45}.
The PSD is shown in Figure~\ref{f:5g}. The spectrum is typical of a PLL-locked oscillator, where the PN at very low frequencies is due to the reference oscillator (see Appendix~\ref{a:mesures} for the description of PN of PLL systems). In Figure~\ref{f:5g}, we also show the proposed approximated model, obtained as the sum of two independent processes, each one with PSD of the form~(\ref{e:psd_l100}), whose parameters are given in Table~\ref{t:AR2}, where we have highlighted the frequencies (low and high) where each process dominates.

\begin{table}[ht]
\centering
\begin{tabular}{|c|c|c|c|c|}
\hline
$\mathrm{PSD0}$ & \multicolumn{4}{c|}{$3675$ ($35.65$ dB)}                                 \\ \hline
$n,m$           & $f_{z,n}$        & $\alpha_{z,n}$ & $f_{p,m}$        & $\alpha_{p,m}$ \\ \hline
$1$             & $3\mathrm{e}3$   & $2.37$         & $1$              & $3.3$          \\ \hline
$2$             & $451\mathrm{e}3$ & $2.7$          & $1.54\mathrm{e}6$ & $3.3$          \\ \hline
$3$             & $458\mathrm{e}6$ & $2.53$         & $30\mathrm{e}6$  & $1$            \\ \hline
\end{tabular}
\vspace{0.2cm}
\caption{Parameters for the 3GPP PN model given in~(\ref{e:pds3gpp}).}
\label{t:param45}
\end{table}

\begin{table}[ht]
\centering
\begin{tabular}{|c|c|| c|c|}
\hline
\multicolumn{2}{|c||}{low frequency}&\multicolumn{2}{c|}{high frequency}\\ \hline
$f_{\mathrm{3dB}}$     & $7\mathrm{e}2$  &$f_{\mathrm{3dB}}$        & $2\mathrm{e}6$\\ \hline
$\ell^2_{100}$     & $-105$ dB& $\ell^2_{100}$     & $-65$ dB\\ \hline
$\ell^2_\infty$     & $-200$ dB&$\ell^2_\infty$      & $-140$ dB\\ \hline
\end{tabular}
\vspace{0.2cm}
\caption{Proposed model: parameters for the two processes for the 3GPP PN.}
\label{t:AR2}
\end{table}

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/PSDcase4}
\caption{Phase noise for a millimeter-wave frequency band, with carrier frequency of 45 GHz. Our model is obtained according to~(\ref{e:extension}).}\label{f:5g}
\end{figure}
We remark that, starting from expression~(\ref{e:pds3gpp}), it is difficult to derive a closed form time-domain model for the PN, while our approach lead to a very simple model for the time-domain PN, that is the sum of two independent AR processes. This is at the expense of an approximation of the PSD. On the other hand, BER simulations of coded systems affected by the PN generated according to our extended model have shown no difference with respect to the case in which the PN follows the 3GPP model.

\section{Conclusions}\label{s:conclusions}
In this paper, we considered the phase noise typical of local oscillators in wireless communications. In particular, we considered both free-running and phase-locked oscillators. Our main contributions are the derivation of analytical models for the power spectral density of the phase noise and of the phasor, described by means of parameters that are directly connected to measurements.
We then considered the discrete-time phase noise channel, and proposed analytical models for the discrete-time phase noise signal, also described by means of parameters linked to measurements. Finally, we evaluated the intersymbol interference and the power loss due to the presence of phase noise, affecting the samples obtained by matched filtering and sampling.

The usefulness of the proposed models and analysis has been demonstrated by using computer simulations, by considering the phase noise typical of satellite communications, i.e., DVB-S2 systems, and the phase noise typical of millimeter wave frequency band, employed in next-generation cellular systems. An important implication of our work is that the proposed models can be used in the design of receiver algorithms and for performance prediction, starting from the knowledge of some fundamental parameters that can be found from measurements of practical oscillators.

\appendices
\section{Experimental measurements}\label{a:mesures}
In this section, we consider experimental measurements of the PN of some oscillators from~\cite{KuLaBa12,ChHeKu17,Pe15}. In Figure~\ref{f:measFreeRunnging}, we show the single-side-band (SSB) PN spectrum of a free-running voltage-controlled oscillator (VCO) designed in~\cite{KuLaBa12}. From the figure, it is possible to see that the spectrum follows a  slope of $-20$ dB/decade on a logarithmic scale. This behaviour is due to integration of white noise inside the oscillator circuitry. Usually, the spectrum presents a flat noise floor at high frequencies, modelling the thermal noise at the oscillator output, but in this experiment it is not visible.
\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/psd2}
\caption{Measured phase noise versus offset frequency for the VCO designed in~\cite{KuLaBa12}.}\label{f:measFreeRunnging}
\end{figure}

We now consider PLL-locked oscillators. They consist of three major functional units, i.e., a phase detector, a loop filter and a VCO. A simplified block diagram of a PLL is shown in Figure~\ref{f:pll}. The PLL generates an output signal whose phase is matched to the phase of a reference signal. This is achieved by comparing the phase of the two signals, by passing the signal proportional to the phase difference to a low-pass filter, and then using the filtered signal to drive the VCO, which tries to  keep the phases matched. PLL systems are in fact composed by two oscillators, the VCO and the oscillator generating the reference signal. Similarly to free-running oscillators, the SSB of PLL-locked oscillators is characterized by a $-20$ dB/decade slope and by the flat noise floor at high frequencies, that mainly come from the PN of the VCO. Moreover, the spectrum of PLL-locked oscillators has a flat part at low frequencies, coming from the low-pass filtering, and it is hence characterized by a loop bandwidth that we have indicated with $f_{\mathrm{3dB}}$. The PN of PLL systems is highly dependent on the VCO, while the oscillator generating the reference signal contributes with a PN at very low frequencies and is often ignored. Figures~\ref{f:measPLL1} and \ref{f:measPLL2} show the measured PN spectra of different PLL-locked oscillators. In both figures, we notice the flat part at low frequencies and the $-20$ db/decade slope. The spectrum in Figure~\ref{f:measPLL1} shows clearly the flat PN at high frequencies, while in Figure~\ref{f:measPLL2}, taken from~\cite{Pe15}, we can see the contribution of the reference oscillator at very low frequencies.

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/PLL}
\caption{Simplified block diagram of a PLL-locked oscillator.}\label{f:pll}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/psd1}
\caption{Measured phase noise spectrum from the oscillator considered in~\cite{ChHeKu17}.}\label{f:measPLL1}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/image001.jpg}
\caption{Measured phase noise for the PLL-locked oscillator designed in~\cite{Pe15}. Figure taken with author's permission. \copyright Persson.}\label{f:measPLL2}
\end{figure}

\section{Proof of theorem~\ref{thm:phasorautocor}}\label{a:thmphasorautocor}
As an help process, we define the process $d_\tau(t)$ of the phase increments in the following way
\begin{equation}
d_\tau(t)=\theta(t)-\theta(t-\tau)\,.
\end{equation}
The process $d_\tau(t)$ is a Gaussian random process~\cite{FoVa88} with zero mean and variance
\begin{align}
\sigma^2_{d_\tau}&=\text{E}\{ (\theta(t)-\theta(t-\tau))^2  \} =\\
&=\text{E}\{ \theta^2(t) + \theta^2(t-\tau) -2 \theta(t)\theta(t-\tau) \} =\\
&=2 R_\theta(0) -2 R_\theta(\tau)=\\
&= \frac{2\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}  (1- e^{-2\pi f_{\mathrm{3dB}}  |\tau|} )  \, ,\label{e:varDelta}
\end{align}
where we have used the expression of the PN autocorrelation given in~(\ref{e:corr}).
This variance increases with $|\tau|$ up to a ceiling at $\frac{2\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}$.

We now consider the autocorrelation of the phasor, which we denote by $R_h(\tau)$. We have
\begin{align}
R_h(\tau)&=\text{E}\{   e^{j\theta(t)} e^{-j\theta(t-\tau)} \} =\\
&= \text{E}\{ e^{jd_\tau(t) }  \} =\\ \label{e:deltaGaus}
& = \int_{-\infty}^{\infty}  e^{jx}  \frac{1}{\sqrt{2\pi } \sigma_{d_\tau} }    e^{-\frac{x^2}{2\sigma^2_{d_\tau}}}  = \\
&= \int_{-\infty}^{\infty}   \frac{1}{\sqrt{2\pi } \sigma_{d_\tau} }    e^{-\frac{x^2-  j 2 x \sigma^2_{d_\tau} }{2\sigma^2_{d_\tau}}}  =\\
&= e^{-\frac{\sigma^2_{d_\tau}}{2}}\int_{-\infty}^{\infty}   \frac{1}{\sqrt{2\pi } \sigma_{d_\tau} }    e^{-\frac{(x-  j x \sigma^2_{d_\tau})^2 }{2\sigma^2_{d_\tau}}}  =\\
&=e^{-\frac{\sigma^2_{d_\tau}}{2}}\label{e:exp_gaus}\, ,
\end{align}
where in~(\ref{e:deltaGaus}) we have used the fact that $d_\tau(t) $ is a zero-mean Gaussian random variable.
Replacing~(\ref{e:varDelta}) in~(\ref{e:exp_gaus}), we obtain~(\ref{e:corr_phasor}).

\section{Proof of theorem~\ref{thm:AR model}}\label{a:thmARmodel}
The continuous-time process $\theta(t)$ with PSD~(\ref{e:psdphi}), with no band limitation, has autocorrelation
\begin{equation}\label{e:corr_psi}
R_\theta(\tau)=\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}  e^{-2\pi f_{\mathrm{3dB}}  |\tau|}\, .
\end{equation}
Therefore, the discrete-time random process $\theta_k$, obtained by sampling the process with sampling period $T_\mathrm{s}$, has autocorrelation
\begin{equation}\label{e:corrthetak}
R_\theta(nT_\mathrm{s})=\text{E} \{ \theta_k \theta_{k-n} \} =  \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}} \, e^{-2\pi f_{\mathrm{3dB}} T_\mathrm{s} |n|}   \, .
\end{equation}
At the same time, being an AR process, the autocorrelation of the process $\theta_k$ can be written in the form
\begin{equation}\label{e:corrpsi}
R_\theta(n)=\text{E} \{ \theta_k \theta_{k-n} \} =\frac{\sigma^2_u}{1-a^2}a^{|n|}\, ,
\end{equation}
from where we see that the parameter $a$ regulates the correlation between the samples $\theta_k$. If we replace the expressions of $a$ and $\sigma^2_u$ in~(\ref{e:corrpsi}), we obtain the autocorrelation given in~(\ref{e:corrthetak}).
Since the PSD of the continuous-time process $\theta(t)$ is not band limited, the PSD of $\theta_k$, that is the repetition of the PSD of $\theta(t)$ with period $1/T_s$, is affected by aliasing. Therefore, the PSD of $\theta_k$ equals~(\ref{e:psdphi}) in the band $ [-1/2T_\mathrm{s}, 1/2T_\mathrm{s}]$, except for the error due to the aliasing, that is
\begin{equation}
    S_\theta(f)=\frac{10^{10}\ell^2_{100}}{f^2_{\mathrm{3dB}}+f^2} +\alpha \, ,
\end{equation}
where $\alpha$ is the aliasing.
 The variance of the error $\alpha$ is given in Lemma~\ref{l:aliasing}. This variance tends to zero when the system bandwidth $1/T_\mathrm{s}$ is much larger than $f_{\mathrm{3dB}}$, that is when $f_{\mathrm{3dB}} T_s\ll 1$.

\section{Proof of theorem~\ref{thm:eta}}\label{a:thmeta}
The proof starts with the computation of $\eta_\text{D}$ and $\eta_\text{ISI}$. Then Lemma~\ref{l:eta} is used to compute $\eta$ as the sum of the previous two terms.

The mean square error between the phasor $e^{j \theta_k}$ and $g_{0,k}$ is given in~(\ref{e:etad}). Expanding~(\ref{e:etad}), we obtain
\begin{equation}\label{e:etaD}
\eta_\text{D}=1 + \gamma_0 - 2 \Re\Big( E_{\theta}\Big \{  e^{-j \theta_k} g_{0,k} \Big\} \Big )\,,
\end{equation}
where $\gamma_0$ is given in~(\ref{e:gamma_l}) with $\ell=0$.
We first consider the computation of $\gamma_0$. We substitute the expression of the pulse $p(t)$ in~(\ref{e:gamma_l}) with $\ell=0$,
\begin{equation}
\gamma_0= \frac{1}{{T_s}^2 }\!\! \!  \int_{-\infty}^{+\infty}   \! \!\! \int_{-\infty}^{+\infty}   \!\!\! \text{sinc}^2\Big(\frac{z}{T_s}  \Big) \text{sinc}^2\Big(\frac{z_1}{T_s}  \Big)  \!R_h(z_1\!\!-\!\!z) d z d z_1, \label{e:gamma0f}
\end{equation}
then we replace $R_h(z_1\! -\!z) $ with
\begin{equation}\nonumber
R_h(z_1\!-\!z) \!\!= \!\!  \mathcal{F}^{-1} \{ S_h (f) e^{-j 2\pi z f} \}\!\! = \!\!\!\! \int_{-\infty}^{+\infty}   \! \!\!  \!S_h (f) e^{j 2\pi z_1 f}\! e^{-j 2\pi z f} \!df
\end{equation}
and rearrange the integrals in~(\ref{e:gamma0f}) in the following way
\begin{align}\nonumber
\gamma_0= &\frac{1}{{T_s}^2 } \!\! \!  \int_{-\infty}^{+\infty}  \Big(  \int_{-\infty}^{+\infty}  \! \!\! \text{sinc}^2\Big(\frac{z}{T_s}  \Big)  e^{-j 2\pi z f}    \!d z  \Big) \cdot \\ \label{e:integral}
&\qquad \cdot  \Big(  \int_{-\infty}^{+\infty}   \!\!\!  \text{sinc}^2\Big(\frac{z_1}{T_s}  \Big) e^{j 2\pi z_1 f} d z_1   \Big)  S_h (f)   d f\,.
\end{align}
The two internal integrals in $z$ and $z_1$ equal the Fourier transform of $\text{sinc}^2$, that is
\begin{equation}
\int_{-\infty}^{+\infty}  \frac{1}{{T_s}^2 } \text{sinc}^2\Big(\frac{t}{T_s}  \Big)  e^{-j 2\pi t f}  d t= \frac{1}{T_s} \Lambda(f T_s),
\end{equation}
where $\Lambda(x)$ is the triangular function, i.e., $\Lambda(x)=1-| x|$, for $-1<x<1$ and zero otherwise. Replacing the expression of $S_h (f)$ given in~(\ref{e:free_phasor_spectrum}) in~(\ref{e:integral}), we finally obtain
\begin{equation}\label{e:gamma0}
\gamma_0=\frac{2}{\pi}  \Big[   \text{tan}^{-1}   \Big( \frac{1}{\rho } \Big) \Big[ 1 -\rho ^2 \Big] -\rho  \text{log}\Big(1+\frac{1}{\rho ^2}\Big) +\rho  \Big]\, ,
\end{equation}
where $\rho =\pi 10^{10} \ell^2_{100} T_s$. We consider the last term in~(\ref{e:etaD}), and compute
\begin{align}\nonumber
E_{\theta}\Big \{ \! e^{-j \theta_k} g_{0,k} \! \Big\}\!\! &= \int_{-\infty}^{+\infty}  \!\!\!\! p^2(z)  E_{\theta}\Big \{ \!  e^{-j \theta(k T_s)} e^{j \theta(z+k T_s)}  \! \Big \} dz=\\\nonumber
&= \frac{1}{T_s} \!\!\int_{-\infty}^{+\infty}  \!\!\!\! \text{sinc}^2\Big(\frac{z}{T_s}\Big)  R_h(z) dz =\\ \nonumber
&=\frac{1}{T_s} \!\!\int_{-\infty}^{+\infty}  \!\!\!\! \text{sinc}^2\Big(\frac{z}{T_s}\Big) e^{-2\pi^2 10^{10}\ell^2_{100} |z| }dz =\\ \label{e:corr_g}
&= \frac{2}{\pi}  \Big[   \text{tan}^{-1}   \Big( \frac{1}{\rho } \Big)-\frac{\rho }{2} \text{log}\Big(1+\frac{1}{\rho ^2}\Big)  \Big]\, .
\end{align}
Using~(\ref{e:gamma0}) and~(\ref{e:corr_g}) in~(\ref{e:etaD}), we obtain equation~(\ref{e:etaDfinal}) and conclude the proof for $\eta_\text{D}$.

Now we move on the computation of $\eta_\text{ISI}$. Using the definition of $\gamma_\ell$ in~(\ref{e:gammal}), we can write
\begin{equation}\label{e:etaisi1}
\eta_\text{ISI}=\sum_{\ell\neq 0} \gamma_\ell\, .
\end{equation}
We show how to compute the sum $\sum_{\ell} \gamma_\ell$ and then we use the expression of $\gamma_0$ reported in~(\ref{e:gamma0}) to compute $\eta_\text{ISI}$. We define $r'(t)$ the received signal without AWGN. Similarly,  we define $y'(t)$ and $y'_k$ the signal after matched filtering and the samples after sampling, respectively, without AWGN. The process $r'(t)$ is cycle-stationary, and can be made wide-sense stationary with the classical approach of introducing a random delay, uniformly distributed in $[-T_s,T_s]$. For this reason, we assume that $r'(t)$ is wide-sense stationary. All the following results can in fact be easily extended to the case where we consider the random delay. Similarly, we assume that the process $y'(t)$, which originates from linear filtering of $r'(t)$, is wide-sense stationary. The power of the process $y'(t)$ can be computed by integrating its PSD, that is
\begin{align}
&P_{y'(t)}= \int_{-\infty}^{+\infty} S_{r'(t)} (f) |P(f)|^2 d f\\
&=T_s\int_{-\frac{1}{2T_s}}^{+\frac{1}{2T_s}} S_{r'(t)} (f) d f \\ \label{e:Sr}
&=\! \frac{T_s}{\pi} \!\!\int_{-\frac{1}{2T_s}}^{+\frac{1}{2T_s}} \!\!\text{tan}^{-1}\!\Big(\frac{f+1/2T_s}{\pi 10^{10} \ell^2_{100}}\Big) \!\!-\!\!\text{tan}^{-1}\!\Big(\frac{f-1/2T_s}{\pi 10^{10} \ell^2_{100} }\Big) d f\\
&=\!\!\frac{2}{\pi} \text{tan}^{-1}\!\Big (\!\frac{1}{\pi 10^{10} \ell^2_{100} T_s} \Big)\!\!-\!\! 10^{10} \ell^2_{100} T_s \text{log}\Big(\!1\!\!+\!\frac{1}{\pi^2 10^{20} \ell^4_{100} T_s^2} \Big)\\
& =\frac{2}{\pi} \Big[ \text{tan}^{-1}\Big (\!\frac{1}{\rho }    \Big) - \frac{\rho }{2} \text{log}\Big(  1 + \frac{1}{\rho ^2}    \Big)  \Big]\label{e:Py}\, .
\end{align}
In~(\ref{e:Sr}), we have used the fact that
the PSD of the received signal is given by the convolution of the PSD of the information carrying  signal and that of the phasor $e^{j\theta(t)}$, since they are independent processes. Since we assume that the process $y'(t)$ is a wide-sense stationary process, and it is bandlimited with bandwidth $1/T_s$, we have the following relationship between the power of the continuous-time process and that of the sequence $y'_k$
\begin{equation}\label{e:Pcd}
P_{y'(t)}= E_{x,\theta}  \Big\{ |y'_k|^2    \Big \}\, .
\end{equation}
The power of the discrete-time signal can be computed starting from~(\ref{e:rec_k}) as
\begin{equation}\label{e:sumgamma1}
E_{x}  \Big\{ |y'_k|^2    \Big \}=   \sum_\ell  \gamma_\ell\, .
\end{equation}
By using (\ref{e:Py}), (\ref{e:Pcd}) and (\ref{e:sumgamma1}), we obtain
\begin{equation}\label{e:sumgamma}
\sum_\ell  \gamma_\ell = \frac{2}{\pi} \Big[ \text{tan}^{-1}\Big (\!\frac{1}{\rho }    \Big) - \frac{\rho }{2} \text{log}\Big(  1 + \frac{1}{\rho ^2}    \Big)  \Big] \, ,
\end{equation}
and, finally, from~(\ref{e:etaisi1}) and~(\ref{e:gamma0}) we get~(\ref{e:etaisifinal}).

Finally, equation~(\ref{e:eta}) is obtained by using~Lemma~\ref{l:eta}.

\section{Phasor PSD: extension to the case $\ell^2_\infty\neq0$}\label{a:linfinity}
In this proof and in the proof in Appendix~\ref{a:linfinity2}, we introduce the bandwidth $B_\theta$ of the white PN as an help parameter. This is necessary in the derivations to avoid dealing with a white process with infinite bandwidth, but the parameter does not appear in the end results. Ideally, $B_\theta$ can be chosen arbitrarily large. On the other hand, we choose $1/T_\mathrm{s}\leq B_\theta<10^{10}$ to limit the variance of the error of the used approximations. This limits the analysis to communication systems with symbol rate $\frac{1}{T_s}<10^{10}$. The extension of the analysis to higher symbol rates must be done without resorting to the approximation of the white PN process with the first order Taylor expansion.

In the following, we mark with a circumflex accent the statistical quantities that refer to the case of $\ell^2_\infty\neq0$.
When we remove the assumption $\ell^2_\infty=0$, the autocorrelation of the PN in~(\ref{e:corr}) becomes
\begin{equation}\label{e_a:corr}
\hat{R}_\theta(\tau)=\frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}  e^{-2\pi f_{\mathrm{3dB}}  |\tau|}+\ell^2_\infty B_\theta\text{sinc}(\tau B_\theta) \, .
\end{equation}
The variance of the help process $d_\tau$ introduced in Appendix~\ref{a:thmphasorautocor} is then
\begin{equation}
\hat{\sigma}^2_{d_\tau}= \frac{2\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}  (1- e^{-2\pi f_{\mathrm{3dB}}  |\tau|}) +  2\ell^2_\infty B_\theta (1-\text{sinc}(\tau B_\theta) ) \label{ea:varDelta}
\end{equation}
and hence the autocorrelation of the phasor given in~(\ref{e:exp_gaus}) becomes
\begin{align}\label{ea:exp_gaus}
\hat{R}_h(\tau)&=e^{- \frac{\pi 10^{10} \ell^2_{100}}{f_{\mathrm{3dB}}}   (1- e^{-2\pi f_{\mathrm{3dB}}  |\tau|} )  - \ell^2_\infty B_\theta (1-\text{sinc}(\tau B_\theta) } \\ \label{ea:exp_gauss2}
&= R_{h}(\tau)  e^{-\ell^2_\infty B_\theta (1-\text{sinc}(\tau B_\theta))}\, ,
\end{align}
where in the last equation we have made explicit the relationship with $R_{h}(\tau)$, the autocorrelation of the phasor when $\ell^2_\infty=0$, whose expression is given in~(\ref{e:corr_phasor}). In practical oscillators, the value $\ell^2_\infty$ is in the order of -120 dB or less, therefore, even for large $B_\theta$, we can use the first order Taylor expansion of the exponential $e^x\simeq 1+x$ in~(\ref{ea:exp_gauss2}), that is
\begin{equation}\label{ea:app}
   e^{-\ell^2_\infty B_\theta (1-\text{sinc}(\tau B_\theta)}\simeq 1-  \ell^2_\infty B_\theta (1-\text{sinc}(\tau B_\theta))\, .
\end{equation}
The variance of the error in the above approximation, denoted by $\sigma^2_\text{e}$, is upper bounded by the second term of the Taylor expansion,
that is
\begin{equation}
   \sigma^2_\text{e}\leq  \ell^4_\infty B_\theta^2 (1-\text{sinc}(\tau B_\theta))^2 /2\leq \ell^4_\infty B_\theta^2\, .
\end{equation}
As stated, this error is low for a large range of values of $B_\theta$.
Using~(\ref{ea:app}) in~(\ref{ea:exp_gauss2}), we obtain
\begin{equation}
\hat{R}_h(\tau)=R_{h}(\tau)  (1 -\ell^2_\infty B_\theta + \ell^2_\infty B_\theta\text{sinc}(\tau B_\theta))\, ,
\end{equation}
and the PSD of the phasor is obtained by taking the Fourier transform of the above autocorrelation, that is
\begin{equation}
    \hat{S}_h(f)=(1 -\ell^2_\infty B_\theta) S_{h}(f)    + S_{h}(f) \otimes \ell^2_\infty \Pi\Big(\frac{f} {B_\theta}\Big)\, .
\end{equation}
In the above equation, the term $1 -\ell^2_\infty B_\theta$ that multiplies $S_{h}(f)$ makes the power of the phasor unitary, while the effect of the flat PN, dominating at high frequencies, is represented by the second term. In the case of free-running oscillator, we have
\begin{align}\nonumber
    S_{h}(f) \otimes  \ell^2_\infty \Pi\Big(\frac{f} {B_\theta}\Big)&= \frac{\ell^2_\infty }{\pi} \Big[ \text{tan}^{-1} \Big(\frac{1}{\pi 10^{10}\ell^2_{100}} \Big(f+\frac{B_\theta}{2}  \Big)\Big)\\ &\!-\!\text{tan}^{-1} \Big(\frac{1}{\pi 10^{10}\ell^2_{100}} \Big(f\!-\!\frac{B_\theta}{2}  \Big)\Big)   \Big]\, .\label{e:flat}
\end{align}
This term is shown in Figure~\ref{f:flat} for $B_\theta=10^6$, $\ell^2_\infty=-120$~dB and $\ell^2_{100}=-90$~dB. We observe that it has a  behavior similar to the flat part of the PN PSD at high frequencies, being it flat in the bandwidth $B_\theta$.
\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{Figures/PSDflat}
\caption{Flat part of the phasor power spectral density given in~(\ref{e:flat}).}\label{f:flat}
\end{figure}

When $B_\theta$ is large, the term~(\ref{e:flat}) is approximately equal to $\ell^2_\infty$. At the same time, being $\ell^2_\infty$ very small, the multiplicative term $1 -\ell^2_\infty B_\theta$ is approximately unitary even for large $B_\theta$. Therefore, we obtain that the PSD of the phasor in the communication band can be written with good approximation as
\begin{equation}
        \hat{S}_h(f)=S_{h}(f)+ \ell^2_\infty\, .
\end{equation}

\section{Discrete-time model for the white PN ($\ell^2_\infty\neq0$)}\label{a:linfinity2}
We here consider the extension of the discrete-time PN model described in Section~\ref{ss:DT model} to account for the case $\ell^2_\infty\neq0$. We assume that the PN is given by the following sum
\begin{equation}
    \theta(t)=\theta^c(t)+\theta^f(t)\, ,
\end{equation}
where $\theta^c(t)$ and $\theta^f(t)$ are two independent processes describing the PN PSD with $-20$~dB/decade slope and the flat part of the PSD, i.e., the first and second term in equation~(\ref{e:psd_l100}), respectively. The discrete-time description of the close-to-carrier PN is given in Section~\ref{ss:DT model}, while here we focus on the term $\theta^f(t)$.
In the following, we use the white PN bandwidth $B_\theta$ introduced in Appendix~\ref{a:linfinity} and we show that, if $B_\theta\geq\frac{1}{T_s}$, we can model the realizations of the process  $\theta^f_k\triangleq\theta^f(k T_s)$ as independent random variables, each with variance $\ell^2_\infty/T_s$.
We consider the following phasor
\begin{align}\label{one}
   e^{j \theta (t)}& = e^{j \theta^c (t)}e^{j \theta^f (t)}\\\label{two}
   &\simeq e^{j \theta^c (t)} +  j e^{j \theta^c (t)} \theta^f (t)\\\label{three}
   &=e^{j \theta^c (t)} +  j \theta^r (t)
\end{align}
where in~(\ref{two}) we have approximated the phasor $e^{j \theta^f (t)}$ with its first order Taylor expansion, while in~(\ref{three}) the process $\theta^r (t)$ is white like $\theta^f (t)$ since it is a random phase rotation of a white process. The process $\theta^f (t)$ has autocorrelation
\begin{equation}
R_{\theta^f}(\tau)=\ell^2_\infty B_\theta\text{sinc}(\tau B_\theta) \, ,
\end{equation}
and the variance of the error of the approximation in~(\ref{two}) is upper bounded by
\begin{equation}
    \sigma^2_\text{exp}\leq\text{E}\{ [\theta^f(t)]^2 \}/2=\ell^2_\infty B_\theta/2\,.
\end{equation}

We now consider the effect of the matched filtering. For simplicity, we assume that $p(t)$ is the normalized sinc function given in~(\ref{e:sinc}), but the analysis can be extended to other kinds of pulse. We define $y'(t)$ the signal after matched filtering
\begin{align}
 y'(t)= &  e^{j \theta^c (t)}\otimes p(t)+ j \theta^r (t)  \otimes p(t)\\
 =& n_c(t)+n_f(t),
\end{align}
where $n_c(t)=e^{j \theta^c (t)}\otimes p(t)$ and $n_f(t)=j \theta^r (t)  \otimes p(t)$.
Since the PN processes are independent, the PSD of $y'(t)$ is the sum of the PSD of $n_c(t)$ and the PSD of $n_f(t)$, which, for $B_\theta\geq\frac{1}{T_s}$, is a rectangle of bandwidth $1/T_s$. Therefore, the signal $y'(t)$ is independent of $B_\theta$, as long as it is larger or equal to the system bandwidth. For this reason, to find a discrete-time model for $\theta^f(t)$, we can assume that $B_\theta=1/T_\mathrm{s}$, and the autocorrelation of the process $\theta^f_k=\theta^f(kT_\mathrm{s})$ becomes
\begin{equation}
    R_\theta^f(m)=\frac{\ell^2_\infty}{ T_s} \delta( m )\, ,
\end{equation}
where $\delta(m)$ is the Kronecker delta. Accordingly, samples $\theta^f_k$ can be generated as  i.i.d. Gaussian random variables with zero mean and variance $\frac{\ell^2_\infty}{ T_s} $.

\section*{Acknowledgement}
We would like to thank Dr. Kuylenstierna for providing the measurement data of Figures~\ref{f:measFreeRunnging} and \ref{f:measPLL1}, and Bj{\"o}rn Persson for giving permission to use Figure~\ref{f:measPLL2} from his Master's thesis.