# An Alternative Derivation of the Gaussian Noise model

Alberto Bononi and Paolo Serena
Università degli studi di Parma, Parma, Italy
14 Aug. 2022

###### Abstract

By extending the results in _Bononi 2012_ , we provide here a complete
alternative derivation of Turin’s Gaussian Noise (GN) model for dual-
polarization dispersion uncompensated coherent optical links. This paper
contains the lecture notes used by the authors first on July 19, 2012, and
then again on june 6, 2016 at the University of Parma.

## References

  * [1] A. Bononi and P. Serena, "An alternative derivation of Johannisson’s regular perturbation model,” arXiv:1207.4729v1 [physics.optics], 2012.
  * [2] P. Johannisson, M. Karlsson, “Perturbation Analysis of Nonlinear Propagation in a Strongly Dispersive Optical Communication System,” J. Lightw. Technol., vol. 31, pp. 1273-1282, 2013.
  * [3] A. Carena, V. Curri, G. Bosco, P. Poggiolini, and F. Forghieri, “Modeling of the Impact of Non-Linear Propagation Effects in Uncompensated Optical Coherent Transmission Links,” J. Lightw. Technol., vol. 30, pp. 1524-1539, 2012.
  * [4] A. Vannucci, P. Serena, and A. Bononi, "The RP method: a new tool for the iterative solution of the nonlinear Schroedinger equation,” J. Lightw. Technol., vol. 20, pp. 1102–1112, Jul. 2002.
  * [5] A. Bononi, P. Serena, N. Rossi, E. Grellier, and F. Vacondio, "Modeling Nonlinearity in Coherent Transmissions with Dominant Intrachannel-Four-Wave-Mixing, Opt. Express, vol. 20, no. 7, pp. 7777-7791, Mar. 2012.
  * [6] A. Papoulis, “Probability, Random Variables, and Stochastic Processes.” III Ed., McGraw-Hill, 1991.
  * [7] J. W. Goodman, “Statistical Optics”. Wiley, NY: 1985, page 44 eq. (2.8-21)
  * [8] I. S. Reed, “On a Moment Theorem for Complex Gaussian Processes,” IRE Trans. Inform. Theory, vol. 8, pp. 194-195, Apr. 1962\.
  * [9] P. Poggiolini _et al.,_ “A Detailed Analytical Derivation of the GN Model of Non-Linear Interference in Coherent Optical Transmission Systems” arXiv:1209.0394v13 [physics.optics].

##  I Introduction

Goal of this paper is to provide an alternative derivation of the results
appearing in P. Johannisson’s et al. [2] on (a generalization of) the well-
known frequency-domain GN model introduced by Turin’s group in [3], i.e., a
closed-form formula for the power spectral density (PSD) of the nonlinear
interference (NLI) at the output of a dispersion uncompensated (DU) coherent
system. The main proof was already reported in [1], and we here make it more
clear by providing more detail. The novelty here is that we also show that the
removal of the extra “phase term” appearing in the output PSD formula in [2]
is due to a change of variable that tracks the average nonlinear phase and
yields exactly Turin’s GN formula [3].

Ref. [2] first proved that the GN model is indeed a first-order regular
perturbation (RP1) model fed by independent Gaussian spectral lines. Given the
equivalence of RP1 and the Volterra series model proved in [4], these results
also explain the coincidence of the Volterra series and the GN model
approaches to the study of DU coherent system.

##  II Frequency Domain NLI RP1 Solution

We start from the dual-polarization (DP) single-channel first-order Regular
Perturbation (RP1) solution of the Manakov-based dispersion-managed nonlinear
Schroedinger equation (DMNLSE) [5, Appendix 2]:

$\tilde{\boldsymbol{U}}(L,f)=\tilde{\boldsymbol{U}}(0,f)+\tilde{\boldsymbol{U}}_{p}(L,f)$
(1)

where $L$ is the total link length, and the NLI _p_ erturbation field is

$\tilde{\boldsymbol{U}}_{p}(L,f)=-jP_{0}\iint_{-\infty}^{\infty}\mathcal{K}(f_{1}f_{2})\tilde{\boldsymbol{U}}(0,f+f_{1})\tilde{\boldsymbol{U}}^{\dagger}(0,f+f_{1}+f_{2})\tilde{\boldsymbol{U}}(0,f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
(2)

where:

i) boldface fields are 2x1 vectors containing the Fourier transforms (denoted
by a tilde) of the X and Y polarizations in the transmitter polarization frame
of reference; a dagger stands for transposition and conjugation; and the DP
field power is normalized to an arbitrary reference power $P_{0}$ (which in
[2] is chosen as the _per-polarization_ average power);

ii) the un-normalized scalar frequency kernel is defined as (Cfr. [2] eq. (45)
and [5] eq. (25)):

$\mathcal{K}(F)\triangleq\int_{0}^{L}\gamma^{\prime}(s)\mathcal{G}(s)\mbox{e}^{-jC(s)(2\pi)^{2}F}ds$
(3)

where $F=f_{1}f_{2}$ is the product of two frequencies,
$\gamma^{\prime}=\frac{8}{9}\gamma$ with $\gamma$ the fiber nonlinear
coefficient, $\mathcal{G}(s)$ the line power gain from $z=0$ to $z=s$, and
$C(s)=-\int_{0}^{s}\beta_{2}(z)dz$ is the cumulated dispersion in the
transmission fibers (with dispersion coefficient $\beta_{2}$) up to coordinate
$s$. We do not normalize here the frequency axis, as done instead in [2]. In
[5] we use the normalized kernel

$\tilde{\eta}(F)=\frac{\mathcal{K}(F)}{\mathcal{K}(0)}$

which then at $F=0$ equals 1. The nonlinear phase referred to power $P_{0}$ is

$\Phi_{NL}=P_{0}\mathcal{K}(0).$

Hence by multiplying and dividing by $\mathcal{K}(0)$ we can recast (2) as

$\tilde{\boldsymbol{U}}_{p}(L,f)=-j\Phi_{NL}\iint_{-\infty}^{\infty}\tilde{\eta}(f_{1}f_{2})\tilde{\boldsymbol{U}}(0,f+f_{1})\tilde{\boldsymbol{U}}^{\dagger}(0,f+f_{1}+f_{2})\tilde{\boldsymbol{U}}(0,f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
(4)

From (4), the X component of the RP1 solution writes explicitly as

$\displaystyle\frac{\tilde{U}_{x,p}(L,f)}{-j\Phi_{NL}}=$
$\displaystyle\iint_{-\infty}^{\infty}\tilde{\eta}(f_{1}f_{2})\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}+$
$\displaystyle\,\iint_{-\infty}^{\infty}\tilde{\eta}(f_{1}f_{2})\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
(5)

where the first line gives the self-phase modulation (SPM) of X on X, while
the second line gives the intra-channel cross-polarization modulation
(I-XPolM) of Y on X. A perfectly dual expression for component Y is obtained
by exchanging the indices $x$ and $y$.

##  III Gaussian Assumption and Johannisson’s Result

In [2],[3] the key assumption is that the input fields are the sum of
independent spectral lines:

$\displaystyle\tilde{U}_{x}(0,f)$
$\displaystyle=\sqrt{f_{0}}\sum_{k=-\infty}^{\infty}\xi_{k}\sqrt{\hat{G}_{x}(kf_{0})}\delta(f-kf_{0})$
$\displaystyle\tilde{U}_{y}(0,f)$
$\displaystyle=\sqrt{f_{0}}\sum_{k=-\infty}^{\infty}\zeta_{k}\sqrt{\hat{G}_{y}(kf_{0})}\delta(f-kf_{0})$

with $\xi_{k}$ and $\zeta_{k}$ independent identically distributed standard
(i.e. zero-mean unit variance) circular complex Gaussian random variables
(RV). Such signals do have a _per-polarization_ power spectral density
$\hat{G}_{x/y}(f)$ (normalized to $P_{0}$) in the limit $f_{0}\to 0$ [3].
After long statistical averaging calculations, the authors get the power
spectral density of the $\tilde{U}_{x,p}(L,f)$ RV as per eq. (89) of reference
[2]. Note that our PSD $\hat{G}(f)$ is normalized such that $G(f)\equiv
P_{0}\hat{G}(f)$, where $G$ is the un-normalized PSD per polarization. Also,
$P_{x}=P_{0}\int_{-\infty}^{\infty}\hat{G}_{x}(f)df$ and
$P_{y}=P_{0}\int_{-\infty}^{\infty}\hat{G}_{y}(f)df$. We report in our
notation the result in [2]:

$\displaystyle\hat{G}_{x,p}(f)=$ $\displaystyle
P_{0}^{2}\\{2\iint_{-\infty}^{\infty}\left|\mathcal{K}((f_{1}-f)(f_{2}-f))\right|^{2}\hat{G}_{x}(f_{1})\hat{G}_{x}(f_{2})\hat{G}_{x}(f_{1}+f_{2}-f)\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+$
$\displaystyle\iint_{-\infty}^{\infty}\left|\mathcal{K}((f_{1}-f)(f_{2}-f))\right|^{2}\hat{G}_{x}(f_{1})\hat{G}_{y}(f_{2})\hat{G}_{y}(f_{1}+f_{2}-f)\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+$
$\displaystyle\mathcal{K}(0)^{2}\hat{G}_{x}(f)\left(4(\int_{-\infty}^{\infty}\hat{G}_{x}(f)df)^{2}+4\int_{-\infty}^{\infty}\hat{G}_{x}(f)df\int_{-\infty}^{\infty}\hat{G}_{y}(f)df+(\int_{-\infty}^{\infty}\hat{G}_{y}(f)df)^{2}\right)\\}$
(6)

and a dual expression for Y is obtained by swapping $x\leftrightarrow y$.
Recall that $\hat{G}_{x,p}(f)$ is the NLI PSD, normalized by $P_{0}$.

An equivalent form of (6) using the normalized kernel is the following

$\displaystyle\hat{G}_{x,p}(f)=$
$\displaystyle\Phi_{NL}^{2}\\{2\iint_{-\infty}^{\infty}\left|\tilde{\eta}(f_{1}f_{2})\right|^{2}\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f+f_{1}+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+$
$\displaystyle\iint_{-\infty}^{\infty}\left|\tilde{\eta}(f_{1}f_{2})\right|^{2}\hat{G}_{x}(f+f_{1})\hat{G}_{y}(f+f_{2})\hat{G}_{y}(f+f_{1}+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+$
$\displaystyle\hat{G}_{x}(f)\left(4(\int_{-\infty}^{\infty}\hat{G}_{x}(f)df)^{2}+4\int_{-\infty}^{\infty}\hat{G}_{x}(f)df\int_{-\infty}^{\infty}\hat{G}_{y}(f)df+(\int_{-\infty}^{\infty}\hat{G}_{y}(f)df)^{2}\right)\\}$
(7)

which better shows the formal parallel with the field equation (5): the field
double integral in $f_{1},f_{2}$ of the product kernel-field-field∗-field
becomes a PSD double integral in $f_{1},f_{2}$ of the product squared kernel
magnitude-PSD-PSD-PSD.

In the next section we provide an alternative proof of (7).

##  IV The New Proof

We here present and comment the results in [1]. We start from (5) and make the
following two assumptions regarding the input X,Y fields $U_{x}(0,t)$,
$U_{y}(0,t)$:

1) they are wide-sense stationary (WSS);

2) they are jointly Gaussian processes.

Regarding assumption 1), we plan to exploit the following extension of result
([6], p. 418, eq. (12-76) ):

_Theorem_ 1

Consider the jointly WSS stochastic processes $x(t)$ and $y(t)$, and let

$\tilde{X}(f)\equiv\mathcal{F}[x(t)]=\int_{-\infty}^{\infty}x(t)\mbox{e}^{-j2\pi
ft}\mbox{d}t$

the Fourier transform of $x$ ( in the mean-square (MS) sense), and
$\tilde{Y}(f)$ is similarly defined. Let their cross power spectral density
(PSD) be
$G_{xy}(f)=\mathcal{F}[R_{xy}(\tau)]=\mathcal{F}[E[x(t+\tau)y^{*}(t)]]$. Then

$E[\tilde{X}(f)\tilde{Y}^{*}(u)]=G_{xy}(f)\delta(u-f)\qquad\square$ (8)

As a byproduct, we also have

$E[\tilde{X}(f)\tilde{X}^{*}(u)]=G_{x}(f)\delta(u-f).$

This theorem thus shows that the Fourier transform of _any MS-integrable WSS
process_ is nonstationary white noise, and thus _the spectral lines of its
Fourier transform are_ _uncorrelated_. This was one of the key assumptions
about the input field in [2],[3], which therefore was tantamount to assuming a
WSS input field in the time domain.

Regarding assumption 2), we plan to exploit the following result, known as the
_complex Gaussian moment theorem_ (CGMT), a generalization to complex
variables of Isserlis theorem [8, 7]:

_Theorem 2_

Let $U_{1},U_{2},\ldots,U_{2k}$ be zero-mean jointly circular complex Gaussian
random variables. Then

$E[U_{1}^{*}U_{2}^{*}...U_{k}^{*}U_{k+1}U_{k+2}...U_{2k}]=\sum_{\pi}E[U_{1}^{*}U_{p}]E[U_{2}^{*}U_{q}].....E[U_{k}^{*}U_{r}]$
(9)

where $\sum_{\pi}$ denotes a summation over the $k!$ possible permutations
$(p,q,...,r)$ of indices $(k+1,k+2,...,2k)\qquad\square$

For instance,

$\displaystyle E[U_{1}^{*}U_{2}^{*}U_{3}^{*}U_{4}U_{5}U_{6}]$
$\displaystyle=E[U_{1}^{*}U_{4}]E[U_{2}^{*}U_{5}]E[U_{3}^{*}U_{6}]$
$\displaystyle+E[U_{1}^{*}U_{4}]E[U_{2}^{*}U_{6}]E[U_{3}^{*}U_{5}]$
$\displaystyle+E[U_{1}^{*}U_{5}]E[U_{2}^{*}U_{4}]E[U_{3}^{*}U_{6}]$
$\displaystyle+E[U_{1}^{*}U_{5}]E[U_{2}^{*}U_{6}]E[U_{3}^{*}U_{4}]$
$\displaystyle+E[U_{1}^{*}U_{6}]E[U_{2}^{*}U_{4}]E[U_{3}^{*}U_{5}]$
$\displaystyle+E[U_{1}^{*}U_{6}]E[U_{2}^{*}U_{5}]E[U_{3}^{*}U_{4}].$ (10)

(scheme: 1,2,3 containing the conjugate terms stay at their place. Terms 4,5,6
get al possible permuted positions).

Let’s now start the new proof. We are interested in the PSD $\hat{G}_{x,p}(f)$
of the NLI field $U_{x,p}(L,t)=\mathcal{F}^{-1}[\tilde{U}_{x,p}(L,f)]$. By
theorem 1 we have:

$E[\tilde{U}_{x,p}(L,f)\tilde{U}_{x,p}^{*}(L,u)]=\hat{G}_{x,p}(f)\delta(u-f).$
(11)

and we wish to get an expression for $\hat{G}_{x,p}(f)$ from the left hand
side, which can be explicitly calculated using (5):

$\displaystyle\frac{E[\tilde{U}_{x,p}(L,f)\tilde{U}_{x,p}^{*}(L,u)]}{\Phi_{NL}^{2}}=$
$\displaystyle
E[\iint_{-\infty}^{\infty}\tilde{\eta}(f_{1}f_{2})[\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})+$
$\displaystyle\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})]\mbox{d}f_{1}\mbox{d}f_{2}\cdot$
$\displaystyle\iint_{-\infty}^{\infty}\tilde{\eta}(f_{3}f_{4})^{*}[\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{x}(0,u+f_{3}+f_{4})\tilde{U}_{x}^{*}(0,u+f_{4})+$
$\displaystyle\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{y}(0,u+f_{3}+f_{4})\tilde{U}_{y}^{*}(0,u+f_{4})]\mbox{d}f_{3}\mbox{d}f_{4}]=$
$\displaystyle\iiiint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}\mbox{d}f_{3}\mbox{d}f_{4}\tilde{\eta}(f_{1}f_{2})\tilde{\eta}(f_{3}f_{4})^{*}\cdot$
$\displaystyle\\{E[\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{x}(0,u+f_{3}+f_{4})\tilde{U}_{x}^{*}(0,u+f_{4})]+$
$\displaystyle
E[\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{y}(0,u+f_{3}+f_{4})\tilde{U}_{y}^{*}(0,u+f_{4})]+$
$\displaystyle
E[\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{y}(0,u+f_{3}+f_{4})\tilde{U}_{y}^{*}(0,u+f_{4})]+$
(12) $\displaystyle
E[\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{x}(0,u+f_{3}+f_{4})\tilde{U}_{x}^{*}(0,u+f_{4})]\\}.$
(13)

The above expectations of circular Gaussian RVs can now be obtained by the
following Theorem 3, which exploits Theorems 1 and 2 and is proved in Appendix
1:

Theorem 3

For jointly stationary circular complex Gaussian zero-mean processes
$A(t),B(t),C(t),D(t),E(t),F(t)$ we have the general formula

$\displaystyle
E\left[\tilde{A}(f+f_{1})\tilde{B}^{*}(f+f_{1}+f_{2})\tilde{C}(f+f_{2})\tilde{D}^{*}(u+f_{3})\tilde{E}(u+f_{3}+f_{4})\tilde{F}^{*}(u+f_{4})\right]$
$\displaystyle=$
$\displaystyle\Bigl{[}G_{ab}(f+f_{1})G_{cd}(f)G_{ef}(f+f_{4})\delta(f_{2})\delta(f_{3})+$
$\displaystyle
G_{ab}(f+f_{1})G_{ed}(f+f_{3})G_{cf}(f)\delta(f_{2})\delta(f_{4})+$
$\displaystyle
G_{cb}(f+f_{2})G_{ad}(f)G_{ef}(f+f_{4})\delta(f_{1})\delta(f_{3})+$
$\displaystyle
G_{cb}(f+f_{2})G_{ed}(f+f_{3})G_{af}(f)\delta(f_{1})\delta(f_{4})+$
$\displaystyle
G_{eb}(f+f_{1}+f_{2})G_{ad}(f+f_{1})G_{cf}(f+f_{2})\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})+$
$\displaystyle
G_{eb}(f+f_{1}+f_{2})G_{cd}(f+f_{2})G_{af}(f+f_{1})\delta(f_{4}-f_{1})\delta(f_{3}-f_{2})\Bigr{]}$
$\displaystyle\cdot\delta(u-f)$ $\displaystyle\qquad\square$ (14)

(scheme: b,d,f containing the conjugate terms stay at their place. Terms a,c,e
get al possible permuted positions. Arguments of terms in products must be the
same, hence the deltas)

We next apply the general formula (14) to the expectations in (13) to get:

First expectation:

$\displaystyle
E\left[\underbrace{\tilde{U}_{x}(0,f+f_{1})}_{A}\underbrace{\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})}_{\mathbf{B}}\underbrace{\tilde{U}_{x}(0,f+f_{2})}_{C}\underbrace{\tilde{U}_{x}^{*}(0,u+f_{3})}_{\mathbf{D}}\underbrace{\tilde{U}_{x}(0,u+f_{3}+f_{4})}_{E}\underbrace{\tilde{U}_{x}^{*}(0,u+f_{4})}_{\mathbf{F}}\right]=$
$\displaystyle\delta(u-f)\Bigl{[}\underbrace{G_{xx}(f+f_{1})}_{G_{ab}(f+f_{1})}\underbrace{G_{xx}(f)}_{G_{cd}(f)}\underbrace{G_{xx}(f+f_{4})}_{G_{ef}(f+f_{4})}\delta(f_{2})\delta(f_{3})+$
$\displaystyle\underbrace{G_{xx}(f+f_{1})}_{G_{ab}(f+f_{1})}\underbrace{G_{xx}(f+f_{3})}_{G_{ed}(f+f_{3})}\underbrace{G_{xx}(f)}_{G_{cf}(f)}\delta(f_{2})\delta(f_{4})+$
$\displaystyle\underbrace{G_{xx}(f+f_{2})}_{G_{cb}(f+f_{2})}\underbrace{G_{xx}(f)}_{G_{ad}(f)}\underbrace{G_{xx}(f+f_{4})}_{G_{ef}(f+f_{4})}\delta(f_{1})\delta(f_{3})+$
$\displaystyle\underbrace{G_{xx}(f+f_{2})}_{G_{cb}(f+f_{2})}\underbrace{G_{xx}(f+f_{3})}_{G_{ed}(f+f_{3})}\underbrace{G_{xx}(f)}_{G_{af}(f)}\delta(f_{1})\delta(f_{4})+$
$\displaystyle\underbrace{G_{xx}(f+f_{1}+f_{2})}_{G_{eb}(f+f_{1}+f_{2})}\underbrace{G_{xx}(f+f_{1})}_{G_{ad}(f+f_{1})}\underbrace{G_{xx}(f+f_{2})}_{G_{cf}(f+f_{2})}\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})+$
$\displaystyle\underbrace{G_{xx}(f+f_{1}+f_{2})}_{G_{eb}(f+f_{1}+f_{2})}\underbrace{G_{xx}(f+f_{2})}_{G_{cd}(f+f_{2})}\underbrace{G_{xx}(f+f_{1})}_{G_{af}(f+f_{1})}\delta(f_{4}-f_{1})\delta(f_{3}-f_{2})\Bigr{]}$
(15)

where $G_{xx}\equiv\hat{G}_{x}$.

Second expectation:

$\displaystyle
E\left[\underbrace{\tilde{U}_{x}(0,f+f_{1})}_{A}\underbrace{\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})}_{\mathbf{B}}\underbrace{\tilde{U}_{y}(0,f+f_{2})}_{C}\underbrace{\tilde{U}_{x}^{*}(0,u+f_{3})}_{\mathbf{D}}\underbrace{\tilde{U}_{y}(0,u+f_{3}+f_{4})}_{E}\underbrace{\tilde{U}_{y}^{*}(0,u+f_{4})}_{\mathbf{F}}\right]=$
$\displaystyle\delta(u-f)\Bigl{[}\underbrace{G_{xy}(f+f_{1})}_{G_{ab}(f+f_{1})}\underbrace{G_{yx}(f)}_{G_{cd}(f)}\underbrace{G_{yy}(f+f_{4})}_{G_{ef}(f+f_{4})}\delta(f_{2})\delta(f_{3})+$
$\displaystyle\underbrace{G_{xy}(f+f_{1})}_{G_{ab}(f+f_{1})}\underbrace{G_{yx}(f+f_{3})}_{G_{ed}(f+f_{3})}\underbrace{G_{yy}(f)}_{G_{cf}(f)}\delta(f_{2})\delta(f_{4})+$
$\displaystyle\underbrace{G_{yy}(f+f_{2})}_{G_{cb}(f+f_{2})}\underbrace{G_{xx}(f)}_{G_{ad}(f)}\underbrace{G_{yy}(f+f_{4})}_{G_{ef}(f+f_{4})}\delta(f_{1})\delta(f_{3})+$
$\displaystyle\underbrace{G_{yy}(f+f_{2})}_{G_{cb}(f+f_{2})}\underbrace{G_{yx}(f+f_{3})}_{G_{ed}(f+f_{3})}\underbrace{G_{xy}(f)}_{G_{af}(f)}\delta(f_{1})\delta(f_{4})+$
$\displaystyle\underbrace{G_{yy}(f+f_{1}+f_{2})}_{G_{eb}(f+f_{1}+f_{2})}\underbrace{G_{xx}(f+f_{1})}_{G_{ad}(f+f_{1})}\underbrace{G_{yy}(f+f_{2})}_{G_{cf}(f+f_{2})}\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})+$
$\displaystyle\underbrace{G_{yy}(f+f_{1}+f_{2})}_{G_{eb}(f+f_{1}+f_{2})}\underbrace{G_{yx}(f+f_{2})}_{G_{cd}(f+f_{2})}\underbrace{G_{xy}(f+f_{1})}_{G_{af}(f+f_{1})}\delta(f_{4}-f_{1})\delta(f_{3}-f_{2})\Bigr{]}$

where $G_{yy}\equiv\hat{G}_{y}$, and assuming uncorrelated X and Y (i.e.,
$G_{yx}\equiv 0$) we get

$\displaystyle
E\left[\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{y}(0,u+f_{3}+f_{4})\tilde{U}_{y}^{*}(0,u+f_{4})\right]=$
$\displaystyle\delta(u-f)\Bigl{[}G_{yy}(f+f_{2})G_{xx}(f)G_{yy}(f+f_{4})\delta(f_{1})\delta(f_{3})+$
$\displaystyle
G_{yy}(f+f_{1}+f_{2})G_{xx}(f+f_{1})G_{yy}(f+f_{2})\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})\Bigr{]}.$
(16)

Third expectation:

$\displaystyle
E\left[\underbrace{\tilde{U}_{x}(0,f+f_{1})}_{A}\underbrace{\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})}_{\mathbf{B}}\underbrace{\tilde{U}_{x}(0,f+f_{2})}_{C}\underbrace{\tilde{U}_{x}^{*}(0,u+f_{3})}_{\mathbf{D}}\underbrace{\tilde{U}_{y}(0,u+f_{3}+f_{4})}_{E}\underbrace{\tilde{U}_{y}^{*}(0,u+f_{4})}_{\mathbf{F}}\right]=$
$\displaystyle\delta(u-f)\Bigl{[}\underbrace{G_{xx}(f+f_{1})}_{G_{ab}(f+f_{1})}\underbrace{G_{xx}(f)}_{G_{cd}(f)}\underbrace{G_{yy}(f+f_{4})}_{G_{ef}(f+f_{4})}\delta(f_{2})\delta(f_{3})+$
$\displaystyle\underbrace{G_{xx}(f+f_{1})}_{G_{ab}(f+f_{1})}\underbrace{G_{yx}(f+f_{3})}_{G_{ed}(f+f_{3})}\underbrace{G_{xy}(f)}_{G_{cf}(f)}\delta(f_{2})\delta(f_{4})+$
$\displaystyle\underbrace{G_{xx}(f+f_{2})}_{G_{cb}(f+f_{2})}\underbrace{G_{xx}(f)}_{G_{ad}(f)}\underbrace{G_{yy}(f+f_{4})}_{G_{ef}(f+f_{4})}\delta(f_{1})\delta(f_{3})+$
$\displaystyle\underbrace{G_{xx}(f+f_{2})}_{G_{cb}(f+f_{2})}\underbrace{G_{yx}(f+f_{3})}_{G_{ed}(f+f_{3})}\underbrace{G_{xy}(f)}_{G_{af}(f)}\delta(f_{1})\delta(f_{4})+$
$\displaystyle\underbrace{G_{yx}(f+f_{1}+f_{2})}_{G_{eb}(f+f_{1}+f_{2})}\underbrace{G_{xx}(f+f_{1})}_{G_{ad}(f+f_{1})}\underbrace{G_{xy}(f+f_{2})}_{G_{cf}(f+f_{2})}\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})+$
$\displaystyle\underbrace{G_{yx}(f+f_{1}+f_{2})}_{G_{eb}(f+f_{1}+f_{2})}\underbrace{G_{xx}(f+f_{2})}_{G_{cd}(f+f_{2})}\underbrace{G_{xy}(f+f_{1})}_{G_{af}(f+f_{1})}\delta(f_{4}-f_{1})\delta(f_{3}-f_{2})\Bigr{]}$

and assuming uncorrelated X and Y (i.e., $G_{yx}\equiv 0$) we get

$\displaystyle
E\left[\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{y}(0,u+f_{3}+f_{4})\tilde{U}_{y}^{*}(0,u+f_{4})\right]=$
$\displaystyle\delta(u-f)\Bigl{[}G_{xx}(f+f_{1})G_{xx}(f)G_{yy}(f+f_{4})\delta(f_{2})\delta(f_{3})+$
$\displaystyle
G_{xx}(f+f_{2})G_{xx}(f)G_{yy}(f+f_{4})\delta(f_{1})\delta(f_{3})\Bigr{]}.$
(17)

Fourth expectation:

$\displaystyle
E\left[\underbrace{\tilde{U}_{x}(0,f+f_{1})}_{A}\underbrace{\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})}_{\mathbf{B}}\underbrace{\tilde{U}_{y}(0,f+f_{2})}_{C}\underbrace{\tilde{U}_{x}^{*}(0,u+f_{3})}_{\mathbf{D}}\underbrace{\tilde{U}_{x}(0,u+f_{3}+f_{4})}_{E}\underbrace{\tilde{U}_{x}^{*}(0,u+f_{4})}_{\mathbf{F}}\right]=$
$\displaystyle\delta(u-f)\Bigl{[}\underbrace{G_{xy}(f+f_{1})}_{G_{ab}(f+f_{1})}\underbrace{G_{yx}(f)}_{G_{cd}(f)}\underbrace{G_{xx}(f+f_{4})}_{G_{ef}(f+f_{4})}\delta(f_{2})\delta(f_{3})+$
$\displaystyle\underbrace{G_{xy}(f+f_{1})}_{G_{ab}(f+f_{1})}\underbrace{G_{xx}(f+f_{3})}_{G_{ed}(f+f_{3})}\underbrace{G_{yx}(f)}_{G_{cf}(f)}\delta(f_{2})\delta(f_{4})+$
$\displaystyle\underbrace{G_{yy}(f+f_{2})}_{G_{cb}(f+f_{2})}\underbrace{G_{xx}(f)}_{G_{ad}(f)}\underbrace{G_{xx}(f+f_{4})}_{G_{ef}(f+f_{4})}\delta(f_{1})\delta(f_{3})+$
$\displaystyle\underbrace{G_{yy}(f+f_{2})}_{G_{cb}(f+f_{2})}\underbrace{G_{xx}(f+f_{3})}_{G_{ed}(f+f_{3})}\underbrace{G_{xx}(f)}_{G_{af}(f)}\delta(f_{1})\delta(f_{4})+$
$\displaystyle\underbrace{G_{xy}(f+f_{1}+f_{2})}_{G_{eb}(f+f_{1}+f_{2})}\underbrace{G_{xx}(f+f_{1})}_{G_{ad}(f+f_{1})}\underbrace{G_{yx}(f+f_{2})}_{G_{cf}(f+f_{2})}\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})+$
$\displaystyle\underbrace{G_{xy}(f+f_{1}+f_{2})}_{G_{eb}(f+f_{1}+f_{2})}\underbrace{G_{yx}(f+f_{2})}_{G_{cd}(f+f_{2})}\underbrace{G_{xx}(f+f_{1})}_{G_{af}(f+f_{1})}\delta(f_{4}-f_{1})\delta(f_{3}-f_{2})\Bigr{]}$

and assuming uncorrelated X and Y (i.e., $G_{yx}\equiv 0$) we get

$\displaystyle
E\left[\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})\tilde{U}_{x}^{*}(0,u+f_{3})\tilde{U}_{x}(0,u+f_{3}+f_{4})\tilde{U}_{x}^{*}(0,u+f_{4})\right]=$
$\displaystyle\delta(u-f)\Bigl{[}G_{yy}(f+f_{2})G_{xx}(f)G_{xx}(f+f_{4})\delta(f_{1})\delta(f_{3})+$
$\displaystyle
G_{yy}(f+f_{2})G_{xx}(f+f_{3})G_{xx}(f)\delta(f_{1})\delta(f_{4})\Bigr{]}.$
(18)

Substitution of (15),(16),(17),(18) into (13) finally gives

$\displaystyle\frac{E[\tilde{U}_{x,p}(L,f)\tilde{U}_{x,p}^{*}(L,u)]}{\Phi_{NL}^{2}}=\delta(u-f)\\{\iiiint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}\mbox{d}f_{3}\mbox{d}f_{4}\tilde{\eta}(f_{1}f_{2})\tilde{\eta}(f_{3}f_{4})^{*}\cdot$
$\displaystyle\cdot\\{\Bigl{[}\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f)\hat{G}_{x}(f+f_{4})\delta(f_{2})\delta(f_{3})+\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{3})\hat{G}_{x}(f)\delta(f_{2})\delta(f_{4})+$
$\displaystyle\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f)\hat{G}_{x}(f+f_{4})\delta(f_{1})\delta(f_{3})+\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f+f_{3})\hat{G}_{x}(f)\delta(f_{1})\delta(f_{4})+$
$\displaystyle\\!\\!\\!\\!\\!\\!\hat{\\!\\!\\!G}_{x}(f+f_{1}+f_{2})\\!\hat{G}_{x}(f+f_{1})\negthinspace\hat{G}_{x}(f+f_{2})\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})\\!+\\!\hat{G}_{x}(f+f_{1}+f_{2})\\!\hat{G}_{x}(f+f_{2})\\!\hat{G}_{x}(f+f_{1})\delta(f_{4}-f_{1})\delta(f_{3}-f_{2})\Bigr{]}\\!+$
$\displaystyle\Bigl{[}\hat{G}_{y}(f+f_{2})\hat{G}_{x}(f)\hat{G}_{y}(f+f_{4})\delta(f_{1})\delta(f_{3})+\hat{G}_{y}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{y}(f+f_{2})\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})\Bigr{]}+$
$\displaystyle\Bigl{[}\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f)\hat{G}_{y}(f+f_{4})\delta(f_{2})\delta(f_{3})+\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f)\hat{G}_{y}(f+f_{4})\delta(f_{1})\delta(f_{3})\Bigr{]}+$
$\displaystyle\Bigl{[}\hat{G}_{y}(f+f_{2})\hat{G}_{x}(f)\hat{G}_{x}(f+f_{4})\delta(f_{1})\delta(f_{3})+\hat{G}_{y}(f+f_{2})\hat{G}_{x}(f+f_{3})\hat{G}_{x}(f)\delta(f_{1})\delta(f_{4})\Bigr{]}\\}\\}.$

From (11), the term multiplying $\delta(u-f)$ must be the desired PSD.

Each of the 6 lines in the above quadruple integral provides 2 integrands
whose integral we must calculate:
1.1)

$\displaystyle\iiiint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}\mbox{d}f_{3}\mbox{d}f_{4}\tilde{\eta}(f_{1}f_{2})\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f)\hat{G}_{x}(f+f_{4})\delta(f_{2})\delta(f_{3})$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)\left[\iint\tilde{\eta}(f_{1}f_{2})\hat{G}_{x}(f+f_{1})\delta(f_{2})\mbox{d}f_{1}\mbox{d}f_{2}\right]\left[\iint\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{4})\delta(f_{3})\mbox{d}f_{3}\mbox{d}f_{4}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)\left[\int\hat{G}_{x}(f+f_{1})\underbrace{\left(\int\tilde{\eta}(f_{1}f_{2})\delta(f_{2})\mbox{d}f_{2}\right)}_{\tilde{\eta}(0)\cdot\int\delta(f_{2})\mbox{d}f_{2}}\mbox{d}f_{1}\right]\left[\int\hat{G}_{x}(f+f_{4})\underbrace{\left(\int\tilde{\eta}(f_{3}f_{4})^{*}\delta(f_{3})\mbox{d}f_{3}\right)}_{\tilde{\eta}^{*}(0)\cdot\int\delta(f_{3})\mbox{d}f_{3}}\mbox{d}f_{4}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)|\tilde{\eta}(0)|^{2}\left[\int\hat{G}_{x}(f+f_{1})\mbox{d}f_{1}\right]\left[\int\hat{G}_{x}(f+f_{4})\mbox{d}f_{4}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)|\tilde{\eta}(0)|^{2}\frac{P_{x}}{P_{0}}$

1.2)

$\displaystyle\iiiint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}\mbox{d}f_{3}\mbox{d}f_{4}\tilde{\eta}(f_{1}f_{2})\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{3})\hat{G}_{x}(f)\delta(f_{2})\delta(f_{4})$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)\left[\iint\tilde{\eta}(f_{1}f_{2})\hat{G}_{x}(f+f_{1})\delta(f_{2})\mbox{d}f_{1}\mbox{d}f_{2}\right]\left[\iint\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{3})\delta(f_{4})\mbox{d}f_{3}\mbox{d}f_{4}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)\left[\int\hat{G}_{x}(f+f_{1})\underbrace{\left(\int\tilde{\eta}(f_{1}f_{2})\delta(f_{2})\mbox{d}f_{2}\right)}_{\tilde{\eta}(0)\cdot\int\delta(f_{2})\mbox{d}f_{2}}\mbox{d}f_{1}\right]\left[\int\hat{G}_{x}(f+f_{3})\underbrace{\left(\int\tilde{\eta}(f_{3}f_{4})^{*}\delta(f_{4})\mbox{d}f_{4}\right)}_{\tilde{\eta}^{*}(0)\cdot\int\delta(f_{4})\mbox{d}f_{4}}\mbox{d}f_{3}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)|\tilde{\eta}(0)|^{2}\left[\int\hat{G}_{x}(f+f_{1})\mbox{d}f_{1}\right]\left[\int\hat{G}_{x}(f+f_{3})\mbox{d}f_{3}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)|\tilde{\eta}(0)|^{2}\frac{P_{x}}{P_{0}}$

2.1)

$\displaystyle\iiiint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}\mbox{d}f_{3}\mbox{d}f_{4}\tilde{\eta}(f_{1}f_{2})\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f)\hat{G}_{x}(f+f_{4})\delta(f_{1})\delta(f_{3})$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)\left[\iint\tilde{\eta}(f_{1}f_{2})\hat{G}_{x}(f+f_{2})\delta(f_{1})\mbox{d}f_{1}\mbox{d}f_{2}\right]\left[\iint\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{4})\delta(f_{3})\mbox{d}f_{3}\mbox{d}f_{4}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)\left[\int\hat{G}_{x}(f+f_{2})\underbrace{\left(\int\tilde{\eta}(f_{1}f_{2})\delta(f_{1})\mbox{d}f_{1}\right)}_{\tilde{\eta}(0)}\mbox{d}f_{2}\right]\left[\int\hat{G}_{x}(f+f_{4})\underbrace{\left(\int\tilde{\eta}(f_{3}f_{4})^{*}\delta(f_{3})\mbox{d}f_{3}\right)}_{\tilde{\eta}^{*}(0)}\mbox{d}f_{4}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)|\tilde{\eta}(0)|^{2}\left[\int\hat{G}_{x}(f+f_{2})\mbox{d}f_{2}\right]\left[\int\hat{G}_{x}(f+f_{4})\mbox{d}f_{4}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)|\tilde{\eta}(0)|^{2}\frac{P_{x}}{P_{0}}$

2.2)

$\displaystyle\iiiint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}\mbox{d}f_{3}\mbox{d}f_{4}\tilde{\eta}(f_{1}f_{2})\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f+f_{3})\hat{G}_{x}(f)\delta(f_{1})\delta(f_{4})$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)\left[\iint\tilde{\eta}(f_{1}f_{2})\hat{G}_{x}(f+f_{2})\delta(f_{1})\mbox{d}f_{1}\mbox{d}f_{2}\right]\left[\iint\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{3})\delta(f_{4})\mbox{d}f_{3}\mbox{d}f_{4}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)\left[\int\hat{G}_{x}(f+f_{2})\underbrace{\left(\int\tilde{\eta}(f_{1}f_{2})\delta(f_{1})\mbox{d}f_{1}\right)}_{\tilde{\eta}(0)}\mbox{d}f_{2}\right]\left[\int\hat{G}_{x}(f+f_{3})\underbrace{\left(\int\tilde{\eta}(f_{3}f_{4})^{*}\delta(f_{4})\mbox{d}f_{4}\right)}_{\tilde{\eta}^{*}(0)}\mbox{d}f_{3}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)|\tilde{\eta}(0)|^{2}\left[\int\hat{G}_{x}(f+f_{2})\mbox{d}f_{2}\right]\left[\int\hat{G}_{x}(f+f_{3})\mbox{d}f_{3}\right]$
$\displaystyle=$
$\displaystyle\hat{G}_{x}(f)|\tilde{\eta}(0)|^{2}\frac{P_{x}}{P_{0}}$

3.1)

$\displaystyle\iiiint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}\mbox{d}f_{3}\mbox{d}f_{4}\tilde{\eta}(f_{1}f_{2})\tilde{\eta}(f_{3}f_{4})^{*}\hat{G}_{x}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})$
$\displaystyle=$
$\displaystyle\iint\tilde{\eta}(f_{1}f_{2})\hat{G}_{x}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\left[\iint\tilde{\eta}(f_{3}f_{4})^{*}\delta(f_{3}-f_{1})\delta(f_{4}-f_{2})\mbox{d}f_{3}\mbox{d}f_{4}\right]\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle=$
$\displaystyle\iint\tilde{\eta}(f_{1}f_{2})\hat{G}_{x}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\underbrace{\left[\int\delta(f_{3}-f_{1})\underbrace{\int\tilde{\eta}(f_{3}f_{4})^{*}\delta(f_{4}-f_{2})\mbox{d}f_{4}}_{\tilde{\eta}^{*}(f_{3}f_{2})}\mbox{d}f_{3}\right]}_{\tilde{\eta}^{*}(f_{1}f_{2})}\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle=$
$\displaystyle\iint|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{x}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$

Each pair of delta removes two integrals, so that the PSD turns out to be
(first two lines above produce first 4 lines, third line above produces 5th
line, 4th line above produces 6th and 7th lines, and final line above produces
the last two lines):

$\displaystyle\frac{\hat{G}_{x,p}(f)}{\Phi_{NL}^{2}}$
$\displaystyle=|\tilde{\eta}(0)|^{2}{\color[rgb]{0,1,0}\iintop_{-\infty}^{\infty}\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f)\hat{G}_{x}(f+f_{4})\mbox{d}f_{1}\mbox{d}f_{4}}$
$\displaystyle+|\tilde{\eta}(0)|^{2}{\color[rgb]{0,1,0}\iintop_{-\infty}^{\infty}\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{3})\hat{G}_{x}(f)\mbox{d}f_{1}\mbox{d}f_{3}}$
$\displaystyle+|\tilde{\eta}(0)|^{2}{\color[rgb]{0,1,0}\iintop_{-\infty}^{\infty}\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f)\hat{G}_{x}(f+f_{4})\mbox{d}f_{2}\mbox{d}f_{4}}$
$\displaystyle+|\tilde{\eta}(0)|^{2}{\color[rgb]{0,1,0}\iintop_{-\infty}^{\infty}\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f+f_{3})\hat{G}_{x}(f)\mbox{d}f_{2}\mbox{d}f_{3}}$
$\displaystyle+2\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{x}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+|\tilde{\eta}(0)|^{2}{\color[rgb]{1,0,0}\iintop_{-\infty}^{\infty}\hat{G}_{y}(f+f_{2})\hat{G}_{x}(f)\hat{G}_{y}(f+f_{4})\mbox{d}f_{2}\mbox{d}f_{4}}$
$\displaystyle+\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{y}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{y}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+2\left({\color[rgb]{1,0,1}|\tilde{\eta}(0)|^{2}\iintop_{-\infty}^{\infty}\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f)\hat{G}_{y}(f+f_{4})\mbox{d}f_{1}\mbox{d}f_{4}}\right.+$
$\displaystyle\,\,\,\,\,\,\left.{\color[rgb]{1,0,1}|\tilde{\eta}(0)|^{2}\iintop_{-\infty}^{\infty}\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f)\hat{G}_{y}(f+f_{4})\mbox{d}f_{2}\mbox{d}f_{4}}\right).$

In summary, considering that by construction $\tilde{\eta}(0)=1$, we have:

$\displaystyle\frac{\hat{G}_{x,p}(f)}{\Phi_{NL}^{2}}$
$\displaystyle=2\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{x}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{y}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{y}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+\hat{G}_{x}(f)\left({\color[rgb]{0,1,0}4}(\frac{P_{x}}{P_{0}})^{2}+{\color[rgb]{1,0,1}4}\frac{P_{x}}{P_{0}}\frac{P_{y}}{P_{0}}+{\color[rgb]{1,0,0}1}(\frac{P_{y}}{P_{0}})^{2}\right)$
(19)

which confirms Johannisson’s equation (7) and completes the desired
alternative proof.

##  V Intermediate Summary

We have presented an alternative derivation of Johannisson’s et al. [2] result
based on the new method in [1]. We first remark that both the result in [2]
and our new approach are able to deal with correlated X and Y, although this
feature was not exploited here. Next we note that we did not have to assume
independent input spectral lines: this comes naturally from the stationarity
of the input process. Finally, the truly critical assumption in the model in
[2, 3] is therefore the assumption of Gaussianity at any $z$ during
propagation, which is implicit in the assumption of a Gaussian input process,
and the fact that the “forcing terms” in the RP equation are the linearly
distorted signals at any $z$, which thus remain Gaussian.

Therefore the true limit of the GN model in [2, 3] is that indeed starting
from a non-Gaussian spectrum such as the one of a digitally modulated
signal111Although the authors in [3] present in their Appendix B an appealing
heuristic justification of their Gaussian signal assumption, still their
invoking the central limit theorem at their equation (37) is not rigorous.
They would conclude that any digitally modulated signal with any number of
levels has a Gaussian Fourier transform (which in turn implies the time-domain
signal itself is Gaussian), which is clearly _not_ the case. , it takes some
finite propagation in a non-infinite dispersion line to approximately get both
a Gaussian spectrum and a Gaussian-like time-domain signal.

##  VI Role of the last line in (19)

We now focus the attention of the last “phase term” in the main result of [2],
i.e. our equation (19). We note that
$\int_{-\infty}^{\infty}\hat{G}_{x}(f)df=P_{x}/P_{0}$,
$\int_{-\infty}^{\infty}\hat{G}_{y}(f)df=P_{y}/P_{0}$, hence the last term is

$\displaystyle\left({\color[rgb]{0,1,0}4}(\int_{-\infty}^{\infty}\hat{G}_{x}(f)df)^{2}+{\color[rgb]{1,0,1}4}\int_{-\infty}^{\infty}\hat{G}_{x}(f)df\int_{-\infty}^{\infty}\hat{G}_{y}(f)df+{\color[rgb]{1,0,0}1}(\int_{-\infty}^{\infty}\hat{G}_{y}(f)df)^{2}\right)$
$\displaystyle=$
$\displaystyle(2\frac{P_{x}}{P_{0}})^{2}+2(2\frac{P_{x}}{P_{0}})\frac{P_{y}}{P_{0}}+(\frac{P_{y}}{P_{0}})^{2}$
$\displaystyle=$ $\displaystyle(\frac{2P_{x}+P_{y}}{P_{0}})^{2}$

Going back to un-normalized kernels and PSDs, we can finally write the sought
nonlinear interference PSD as

$\displaystyle G_{x,p}(f)$
$\displaystyle=\mathcal{K}(0)^{2}\\{2\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}G_{x}(f+f_{1}+f_{2})G_{x}(f+f_{1})G_{x}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}G_{y}(f+f_{1}+f_{2})G_{x}(f+f_{1})G_{y}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}\\}$
$\displaystyle+\mathcal{K}(0)^{2}\left(2P_{x}+P_{y}\right)^{2}G_{x}(f)\\}.$
(20)

The last line can be interpreted as the PSD of the first-order regular
perturbation (RP1) solution of a system whose “exact” output is:

$\displaystyle
U_{x}(L,t)=U_{x}(0,t)\mbox{e}^{-j\gamma\int_{0}^{L}\mathcal{G}(s)\mbox{d}s\cdot(2P_{x}+P_{y})}\cong$
$\displaystyle
U_{x}(0,t)(1-j\overbrace{\gamma\int_{0}^{L}\mathcal{G}(s)\mbox{d}s}^{\mathcal{K}(0)}\cdot(2P_{x}+P_{y}))$
$\displaystyle=$ $\displaystyle
U_{x}(0,t)-\underbrace{j\mathcal{K}(0)(2P_{x}+P_{y})U_{x}(0,t)}_{\triangleq
U_{x,p}(L,t)}$

so that the PSD of $U_{x,p}(L,t)$ is
$\mathcal{K}(0)^{2}(2P_{x}+P_{y})^{2}G_{x}(f)$. Note the surprising average
nonlinear phase rotation by an effective power $2P_{x}+P_{y}$ which is larger
than the actual average power in the fiber $P_{x}+P_{y}$! This is one of the
new unexpected features of dispersion uncompensated systems.

In such DU systems we are thus naturally lead to postulate an ansatz of the DP
dispersion managed nonlinear Schroedinger equation (DMNLSE) of the kind

$\tilde{\boldsymbol{A}}(z,f)=\sqrt{P_{0}\mathcal{G}(z)}\mbox{e}^{j\frac{C(z)\omega^{2}}{2}}\left[\begin{array}[]{c}\tilde{U}_{x}(z,f){\color[rgb]{1,0,0}\mbox{e}}^{{\color[rgb]{1,0,0}-j\gamma\int_{0}^{L}\mathcal{G}(s)\mbox{d}s\cdot(2P_{x}+P_{y})}}\\\
\tilde{U}_{y}(z,f){\color[rgb]{1,0,0}\mbox{e}}^{{\color[rgb]{1,0,0}-j\gamma\int_{0}^{L}\mathcal{G}(s)\mbox{d}s\cdot(2P_{y}+P_{x})}}\end{array}\right]$

which we call the dual-polarization enhanced regular perturbation (DP-ERP)
ansatz (from a similar ERP definition in [4]). We prove below as an exercise
that this change of variables removes the unwanted spectral components in the
last line of equation (20). That’s what tacitly also Turin’s group does in
their detailed GN model derivation [9].

##  VII Analysis

The procedure reported below re-derives from scratch the DM-NLSE and the RP1
solution in dual polarization and in the correct phase-reference.

The Manakov Nonlinear Schroedinger equation (M-NLSE) in engineering notation
writes in the frequency domain as [5]

$\displaystyle\frac{\partial\tilde{\boldsymbol{A}}(z,f)}{\partial z}=$
$\displaystyle\frac{g(z)-j\omega^{2}\beta_{2}(z)}{2}\tilde{\boldsymbol{A}}(z,f)+$
(21)
$\displaystyle-j\gamma(z)\iint_{-\infty}^{\infty}\underbrace{\tilde{\boldsymbol{A}}(z,f+f_{1})\tilde{\boldsymbol{A}}^{\dagger}(z,f+f_{1}+f_{2})\tilde{\boldsymbol{A}}(z,f+f_{2})}_{\left[\begin{array}[]{c}(\tilde{A}_{x}^{*}(z,f+f_{1}+f_{2})\tilde{A}_{x}(z,f+f_{2})+\tilde{A}_{y}^{*}\tilde{A}_{y})\tilde{A}_{x}(z,f+f_{1})\\\
(\tilde{A}_{x}^{*}(z,f+f_{1}+f_{2})\tilde{A}_{x}(z,f+f_{2})+\tilde{A}_{y}^{*}\tilde{A}_{y})\tilde{A}_{y}(z,f+f_{1})\end{array}\right]}\mbox{d}f_{1}\mbox{d}f_{2}$
(24)

where $\omega=2\pi f$, and $f$ is the frequency normalized to the baud rate.
The input modulated field $\tilde{\boldsymbol{A}}_{0}(f)$ may be pre-chirped
to give
$\tilde{\boldsymbol{A}}(0,f)\equiv\tilde{\boldsymbol{A}}_{0}(f)\mbox{e}^{j\frac{\omega^{2}}{2}\xi_{pre}}$,
where $\xi_{pre}$ is the pre-compensation. Now make the change of variable

$\left[\begin{array}[]{c}\tilde{A}_{x}\\\
\tilde{A}_{y}\end{array}\right](z,f)=\sqrt{P_{0}}\left[\begin{array}[]{c}\tilde{U}_{x}(z,f)\mbox{e}^{\frac{\ln\mathcal{G}(z)+jC(z)\omega^{2}}{2}-j\int_{0}^{z}\gamma(s)\mathcal{G}(s)\mbox{d}s\cdot(2P_{x}+P_{y})}\\\
\tilde{U}_{y}(z,f)\mbox{e}^{\frac{\ln\mathcal{G}(z)+jC(z)\omega^{2}}{2}-j\int_{0}^{z}\gamma(s)\mathcal{G}(s)\mbox{d}s\cdot(2P_{y}+P_{x})}\end{array}\right]$
(25)

where $P_{0}$ is a reference normalizing power,
$\mathcal{G}(z)=\mbox{e}^{\int_{0}^{z}g(s)\mbox{d}s}$ is the power gain at
$z$, and $C(z)\triangleq\xi_{pre}-\int_{0}^{z}\beta_{2}(s)\mbox{d}s$ is the
_cumulated dispersion_ up to $z$. The change can also be written more
explicitly as

$\displaystyle\tilde{A}_{x}(z,f)$
$\displaystyle=\sqrt{P_{0}\mathcal{G}(z)}\tilde{U}_{x}(z,f)\mbox{e}^{\psi_{x}(f^{2})}$
$\displaystyle\tilde{A}_{y}(z,f)$
$\displaystyle=\sqrt{P_{0}\mathcal{G}(z)}\tilde{U}_{y}(z,f)\mbox{e}^{\psi_{y}(f^{2})}$
(26)

where we defined

$\displaystyle\psi_{x}(f^{2})$ $\displaystyle=j\frac{C(z)(2\pi
f)^{2}}{2}-j\phi_{x}$ $\displaystyle\psi_{y}(f^{2})$
$\displaystyle=j\frac{C(z)(2\pi f)^{2}}{2}-j\phi_{y}$

$\displaystyle\phi_{x}(z)$
$\displaystyle=\int_{0}^{z}\gamma(s)\mathcal{G}(s)\mbox{d}s\cdot(2P_{x}+P_{y})$
$\displaystyle\phi_{y}(z)$
$\displaystyle=\int_{0}^{z}\gamma(s)\mathcal{G}(s)\mbox{d}s\cdot(2P_{y}+P_{x})$

Differentiating Eq. (25) gives

$\frac{\partial\tilde{A}_{x}}{\partial
z}=\sqrt{P_{0}\mathcal{G}(z)}\left(\frac{\partial\tilde{U}_{x}}{\partial
z}+[\frac{g(z)-j\omega^{2}\beta_{2}(z)}{2}-j\underbrace{\gamma(z)\mathcal{G}(z)(2P_{x}+P_{y})}_{\phi_{X}^{\prime}(z)}]\tilde{U}_{x}\right)\mbox{e}^{j\psi_{x}(\omega^{2})}$
(27)

and a dual equation holds for $\tilde{A}_{y}$. Substituting into Eq. (21) we
get for the X component:

$\displaystyle\sqrt{P_{0}\mathcal{G}(z)}\left(\frac{\partial\tilde{U}_{x}}{\partial
z}+[\frac{g(z)-j\omega^{2}\beta_{2}(z)}{2}+j\gamma(z)\mathcal{G}(z)]\tilde{U}_{x}\right)\mbox{e}^{\psi_{x}(f^{2})}$
$\displaystyle=\frac{g(z)-j\omega^{2}\beta_{2}(z)}{2}\tilde{U}_{x}(z,f)\sqrt{P_{0}\mathcal{G}(z)}\mbox{e}^{\psi_{x}(f^{2})}$
$\displaystyle-j\gamma(z)\iint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}P_{0}\mathcal{G}(z)$
$\displaystyle(\mbox{e}^{\psi_{x}((f+f_{1}+f_{2})^{2})*}\tilde{U}_{x}^{*}(z,f+f_{1}+f_{2})\mbox{e}^{\psi_{x}((f+f_{2})^{2}}\tilde{U}_{x}(z,f+f_{2})$
$\displaystyle+\mbox{e}^{\psi_{y}((f+f_{1}+f_{2})^{2})*}\tilde{U}_{y}^{*}\mbox{e}^{\psi_{y}((f+f_{2})^{2})}\tilde{U}_{y})$
$\displaystyle\cdot$
$\displaystyle\sqrt{P_{0}\mathcal{G}(z)}\mbox{e}^{\psi_{x}((f+f_{1})^{2})}\tilde{U}_{x}(z,f+f_{1}).$

Dividing both sides by $\sqrt{P_{0}\mathcal{G}(z)}e^{\psi_{x}(f^{2})}$ and
simplifying, one gets:

$\displaystyle\frac{\partial\tilde{U}_{x}}{\partial
z}+j\gamma(z)\mathcal{G}(z)(2P_{x}+P_{y})\tilde{U}_{x}$
$\displaystyle=-j\gamma(z)\iint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}P_{0}\mathcal{G}(z)\mbox{e}^{-\psi_{x}(f^{2})}$
$\displaystyle(\mbox{e}^{\psi_{x}((f+f_{1}+f_{2})^{2})*}\tilde{U}_{x}^{*}(z,f+f_{1}+f_{2})\mbox{e}^{\psi_{x}((f+f_{2})^{2})}\tilde{U}_{x}(z,f+f_{2})$
$\displaystyle+\mbox{e}^{\psi_{y}((f+f_{1}+f_{2})^{2})*}\tilde{U}_{y}^{*}\mbox{e}^{\psi_{y}((f+f_{2})^{2})}\tilde{U}_{y})$
$\displaystyle\cdot$
$\displaystyle\mbox{e}^{\psi_{x}((f+f_{1})^{2}]}\tilde{U}_{x}(z,f+f_{1}).$

Now let’s collect all the exponents of $\mbox{e}^{(.)}$ together in the
r.h.s.. For the first term we have the exponent

$\displaystyle-\psi_{x}(f^{2})+\psi_{x}((f+f_{1}+f_{2})^{2})*+\psi_{x}((f+f_{2})^{2})+\psi_{x}((f+f_{1})^{2})$
$\displaystyle=$
$\displaystyle-[j\frac{C(z)\omega^{2}}{2}-j\phi_{x}]-[j\frac{C(z)(\omega+\omega_{1}+\omega_{2})^{2}}{2}-j\phi_{x}]+[j\frac{C(z)(\omega+\omega_{2})^{2}}{2}-j\phi_{x}]+[j\frac{C(z)(\omega+\omega_{1})^{2}}{2}-j\phi_{x}]$
$\displaystyle=$
$\displaystyle\frac{C(z)}{2}(-\omega^{2}-(\omega+\omega_{1}+\omega_{2})^{2}+(\omega+\omega_{2})^{2}+(\omega+\omega_{1})^{2})$
$\displaystyle=$ $\displaystyle-C(z)\omega_{1}\omega_{2}$

where we used the relation

$[\omega^{2}-(\omega+\omega_{1})^{2}-(\omega+\omega_{1})^{2}+(\omega+\omega_{1}+\omega_{2})^{2}]=2\omega_{1}\omega_{2}.$
(28)

For the second term we have the exponent

$\displaystyle-\psi_{x}(f^{2})+\psi_{y}((f+f_{1}+f_{2})^{2})*+\psi_{y}((f+f_{2})^{2})+\psi_{x}((f+f_{1})^{2})$
$\displaystyle=$
$\displaystyle-[j\frac{C(z)\omega^{2}}{2}-j\phi_{x}]-[j\frac{C(z)(\omega+\omega_{1}+\omega_{2})^{2}}{2}-j\phi_{y}]+[j\frac{C(z)(\omega+\omega_{2})^{2}}{2}-j\phi_{y}]+[j\frac{C(z)(\omega+\omega_{1})^{2}}{2}-j\phi_{x}]$
$\displaystyle=$
$\displaystyle\frac{C(z)}{2}(-\omega^{2}-(\omega+\omega_{1}+\omega_{2})^{2}+(\omega+\omega_{2})^{2}+(\omega+\omega_{1})^{2})$
$\displaystyle=$ $\displaystyle-C(z)\omega_{1}\omega_{2}$

Hence finally we get the Manakov DMNLSE (M-DMNLSE) for X:

$\displaystyle\frac{\partial\tilde{U}_{x}(z,f)}{\partial z}$
$\displaystyle{\color[rgb]{1,0,0}=+j\underbrace{\gamma(z)\mathcal{G}(z)(2P_{x}+P_{y})}_{\phi_{x}^{\prime}(z)}\tilde{U}_{x}(z,f)}-j\gamma(z)\mathcal{G}(z)P_{0}\iint_{-\infty}^{\infty}\mbox{e}^{-jC(z)\omega_{1}\omega_{2}}$
$\displaystyle(\tilde{U}_{x}^{*}(z,f+f_{1}+f_{2})\tilde{U}_{x}(z,f+f_{2})$
$\displaystyle+\tilde{U}_{y}^{*}(z,f+f_{1}+f_{2})\tilde{U}_{y}(z,f+f_{2}))\tilde{U}_{x}(z,f+f_{1})\mbox{d}f_{1}\mbox{d}f_{2}$
(29)

and swapping indices x and y we get the dual equation for Y. In red I have
highlighted the main difference from the traditional result.

###  VII-A Dual Polarization first-order Enhanced Regular Perturbation (DP-
ERP1)

The output field $\tilde{\boldsymbol{U}}(L,f)$ is obtained as the following
formal explicit solution of (29):

$\displaystyle\tilde{U}_{x}(L,f)-\tilde{U}_{x}(0,f)$
$\displaystyle=j\int_{0}^{L}\phi_{x}^{\prime}(z)\tilde{U}_{x}(s,f)\mbox{d}s$
$\displaystyle-
jP_{0}\iint_{-\infty}^{\infty}\int_{0}^{L}\gamma(s)\mathcal{G}(s)\mbox{e}^{-jC(s)\omega_{1}\omega_{2}}\cdot$
$\displaystyle(\tilde{U}_{x}^{*}(s,f+f_{1}+f_{2})\tilde{U}_{x}(s,f+f_{2})$
$\displaystyle+\tilde{U}_{y}^{*}(s,f+f_{1}+f_{2})\tilde{U}_{y}(s,f+f_{2}))\tilde{U}_{x}(s,f+f_{1})\mbox{d}s\mbox{d}f_{1}\mbox{d}f_{2}.$
(30)

For systems where the spectrum of the propagating normalized signal
$\tilde{\boldsymbol{U}}(z,\omega)$ does not significantly vary along $z$ we
see that the exact solution (30) may be approximated by the first-order
regular perturbation (DP-ERP1) solution, which we write for the X component:

$\displaystyle\tilde{U}_{x}(L,f)=$
$\displaystyle\tilde{U}_{x}(0,f)+j\underbrace{\phi_{x}(L)}_{(2P_{x}+P_{y})\mathcal{K}(0)}\tilde{U}_{x}(0,f)$
$\displaystyle-
jP_{0}\iint_{-\infty}^{\infty}\mathcal{K}(f_{1}f_{2})(\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})+$
$\displaystyle\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2}))\tilde{U}_{x}(0,f+f_{1})\mbox{d}f_{1}\mbox{d}f_{2}$
(31)

and a dual expression holds for Y by swapping indices $x\leftrightarrow y$,
where we defined

$\mathcal{K}(F)\triangleq\int_{0}^{L}\gamma(s)\mathcal{G}(s)\mbox{e}^{-jC(s)(2\pi)^{2}F}ds$
(32)

as the un-normalized frequency kernel, a function of only the product
$F=f_{1}f_{2}$. We may also define the (cumulated) nonlinear phase as

$\Phi_{NL}=P_{0}\mathcal{K}(0)\qquad[rad]$ (33)

and the normalized (scalar) frequency kernel, or briefly the frequency kernel
as

$\displaystyle\tilde{\eta}(F)=\frac{\mathcal{K}(F)}{\mathcal{K}(0)}=\frac{\int_{0}^{L}\gamma(s)\mathcal{G}(s)\mbox{e}^{-jC(s)(2\pi)^{2}F}ds}{\int_{0}^{L}\gamma(s)\mathcal{G}(s)ds}$
(34)

so that one finally gets a convenient form of the RP1 perturbation as:

$\frac{\tilde{U}_{x,p}(L,f)}{-j\Phi_{NL}}=-A(f)+B(f)$

with

$A(f)\equiv\hat{P}_{Tx}\tilde{U}_{x}(0,f)$

and

$\displaystyle B(f)\equiv$
$\displaystyle\iint_{-\infty}^{\infty}\mathcal{\tilde{\eta}}(f_{1}f_{2})(\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})+$
$\displaystyle\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2}))\tilde{U}_{x}(0,f+f_{1})\mbox{d}f_{1}\mbox{d}f_{2}$

where $\hat{P}_{x}=P_{x}/P_{0}$, $\hat{P}_{y}=P_{y}/P_{0}$, and
$\hat{P}_{Tx}=(2\hat{P}_{x}+\hat{P}_{y})$.

Now we must repeat the long calculations used to verify Joannisson’s result in
[1], and check if indeed the NLI power spectral density obtained from the DP-
ERP1 solution does not have the unwanted “phase” components, third line in
(20). We need

$\frac{E[\tilde{U}_{x,p}(L,f)\tilde{U}_{x,p}^{*}(L,u)]}{\Phi_{NL}^{2}}=E[A(f)A^{*}(u)]-2\mbox{Re}(E[B(f)A^{*}(u)])+E[B(f)B^{*}(u)]$

and in [1] we already calculated $E[B(f)B^{*}(u)]$, while

$\displaystyle E[A(f)A^{*}(u)]$
$\displaystyle=\hat{P}_{Tx}^{2}E[\tilde{U}_{x}(0,f)\tilde{U}_{x}^{*}(0,u)]$
$\displaystyle=\hat{P}_{Tx}^{2}\hat{G}_{x}(f)\delta(f-u)$

where we used Theorem 1 in [1]. We finally have to evaluate

$\displaystyle E[B(f)A^{*}(u)]$
$\displaystyle=\iint_{-\infty}^{\infty}\mbox{d}f_{1}\mbox{d}f_{2}\mathcal{\tilde{\eta}}(f_{1}f_{2})\cdot$
$\displaystyle\hat{P}_{Tx}\\{E[\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,u)]+$
$\displaystyle
E[\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,u)]\\}.$

For 4 generic complex Gaussian RVs, theorem 2 in [1] specializes to

$E[U_{1}^{*}U_{2}^{*}U_{3}U_{4}]=E[U_{1}^{*}U_{3}]E[U_{2}^{*}U_{4}]+E[U_{1}^{*}U_{4}]E[U_{2}^{*}U_{3}].$

Hence

$\displaystyle
E[\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,u)]$
$\displaystyle=$
$\displaystyle\underbrace{E[\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{2})]}_{\hat{G}_{x}(f+f_{2})\delta(f_{1})}\underbrace{E[\tilde{U}_{x}^{*}(0,u)\tilde{U}_{x}(0,f+f_{1})]}_{\hat{G}_{x}(f+f_{1})\delta(u-f-
f_{1})}$ $\displaystyle+$
$\displaystyle\underbrace{E[\tilde{U}_{x}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{1})]}_{\hat{G}_{x}(f+f_{1})\delta(f_{2})}\underbrace{E[\tilde{U}_{x}^{*}(0,u)\tilde{U}_{x}(0,f+f_{2})]}_{\hat{G}_{x}(f+f_{2})\delta(u-f-
f_{2})}$ $\displaystyle=$
$\displaystyle\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f)\delta(f_{1})\delta(u-f)+\hat{G}_{x}(f+f_{1})\delta(f_{2})\hat{G}_{x}(f)\delta(u-f)$

and

$\displaystyle
E[\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})\tilde{U}_{x}(0,f+f_{1})\tilde{U}_{x}^{*}(0,u)]$
$\displaystyle=$
$\displaystyle\underbrace{E[\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{y}(0,f+f_{2})]}_{\hat{G}_{y}(f+f_{2})\delta(f_{1})}\underbrace{E[\tilde{U}_{x}^{*}(0,u)\tilde{U}_{x}(0,f+f_{1})]}_{\hat{G}_{x}(f+f_{1})\delta(u-f-
f_{1})}$ $\displaystyle+$
$\displaystyle\underbrace{E[\tilde{U}_{y}^{*}(0,f+f_{1}+f_{2})\tilde{U}_{x}(0,f+f_{1})]}_{0}\underbrace{E[\tilde{U}_{x}^{*}(0,u)\tilde{U}_{y}(0,f+f_{2})]}_{0}$
$\displaystyle=$
$\displaystyle\hat{G}_{y}(f+f_{2})\hat{G}_{x}(f)\delta(f_{1})\delta(u-f)$ .

Thus

$\displaystyle E[B(f)A^{*}(u)]=$
$\displaystyle\delta(u-f)\iint_{-\infty}^{\infty}\mathcal{\tilde{\eta}}(f_{1}f_{2})\hat{P}_{Tx}[\hat{G}_{x}(f+f_{2})\hat{G}_{x}(f)\delta(f_{1})+$
$\displaystyle+$
$\displaystyle\hat{G}_{x}(f+f_{1})\delta(f_{2})\hat{G}_{x}(f)+\hat{G}_{y}(f+f_{2})\hat{G}_{x}(f)\delta(f_{1})]\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle=$
$\displaystyle\delta(u-f)\hat{G}_{x}(f)\hat{P}_{Tx}\mathcal{\tilde{\eta}}(0)(2\int_{-\infty}^{\infty}\hat{G}_{x}(f_{2})\mbox{d}f_{2}+\int_{-\infty}^{\infty}\hat{G}_{y}(f_{2})\mbox{d}f_{2})$
$\displaystyle=$
$\displaystyle\delta(u-f)\hat{G}_{x}(f)\hat{P}_{Tx}(2\hat{P}_{x}+\hat{P}_{y}).$

Finally from [1] we already know that

$\displaystyle E[B(f)B^{*}(u)]$
$\displaystyle=\delta(u-f)\\{2\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{x}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{y}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{y}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+\hat{G}_{x}(f)(2\hat{P}_{x}+\hat{P}_{y})^{2}\\}.$

Therefore, getting rid of $\delta(u-f)$ we have

$\displaystyle\frac{\hat{G}_{x,p}(f)}{\Phi_{NL}^{2}}$
$\displaystyle{\color[rgb]{1,0,0}=\\{\hat{P}_{Tx}^{2}\hat{G}_{x}(f)}$
$\displaystyle{\color[rgb]{1,0,0}-2\hat{G}_{x}(f)\hat{P}_{Tx}(2\hat{P}_{x}+\hat{P}_{y})}$
$\displaystyle{\color[rgb]{1,0,0}+\hat{G}_{x}(f)(2\hat{P}_{x}+\hat{P}_{y})^{2}\\}}$
$\displaystyle+2\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{x}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{x}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$
$\displaystyle+\iintop_{-\infty}^{\infty}|\tilde{\eta}(f_{1}f_{2})|^{2}\hat{G}_{y}(f+f_{1}+f_{2})\hat{G}_{x}(f+f_{1})\hat{G}_{y}(f+f_{2})\mbox{d}f_{1}\mbox{d}f_{2}$

and we recognize that the red curly bracket is zero since by definition
$\hat{P}_{Tx}\equiv(2\hat{P}_{x}+\hat{P}_{y}).$ Hence indeed the DP-ERP ansatz
removes the unwanted “phase” terms in the NLI PSD.

##  VIII Conclusions

This report is conceptually very important, as it explains the following
apparent contradiction about the RP1 solution of the DMNLSE:

It is known from [4] that the RP1 solution largely overestimates the true
DMNLSE solution, and holds only at extremely small nonlinear phases. How can
it then accurately reproduce the nonlinear interference (NLI) power spectral
density, as claimed by Turin’s group [3]?

The explanation is the following: if instead of the RP1 solution we use the
ERP1 solution [4], and as extended here to dual polarization and called the
DP-ERP1 solution, then the numerical predictions are very close to the true
field, and thus also the NLI variance predictions are close to real. Moreover,
the DP-ERP1 predicted NLI PSD is exactly that reported by Turin’s group [3],
without the puzzling “phase term” in the RP1 result in [2].

## Appendix 1

In this Appendix we prove Theorem 3 in the text. Assuming jointly stationary
circular complex Gaussian zero-mean processes $A(t),B(t),C(t),D(t),E(t),F(t)$,
we have by using (8) in (10):

$\displaystyle T\triangleq
E\left[\tilde{A}(f+f_{1})\tilde{B}^{*}(f+f_{1}+f_{2})\tilde{C}(f+f_{2})\tilde{D}^{*}(u+f_{3})\tilde{E}(u+f_{3}+f_{4})\tilde{F}^{*}(u+f_{4})\right]$
$\displaystyle=$
$\displaystyle\underbrace{E[\tilde{B}^{*}(f+f_{1}+f_{2})\tilde{A}(f+f_{1})]}_{G_{ab}(f+f_{1})\delta(f_{2})}\underbrace{E[\tilde{D}^{*}(u+f_{3})\tilde{C}(f+f_{2})]}_{G_{cd}(f+f_{2})\delta(u+f_{3}-f-f_{2})}\underbrace{E[\tilde{F}^{*}(u+f_{4})\tilde{E}(u+f_{3}+f_{4})]}_{G_{ef}(u+f_{3}+f_{4})\delta(f_{3})}$
$\displaystyle+$
$\displaystyle\underbrace{E[\tilde{B}^{*}(f+f_{1}+f_{2})\tilde{A}(f+f_{1})]}_{G_{ab}(f+f_{1})\delta(f_{2})}\underbrace{E[\tilde{D}^{*}(u+f_{3})\tilde{E}(u+f_{3}+f_{4})]}_{G_{ed}(u+f_{3}+f_{4})\delta(f_{4})}\underbrace{E[\tilde{F}^{*}(u+f_{4})\tilde{C}(f+f_{2})]}_{G_{cf}(f+f_{2})\delta(u+f_{4}-f-f_{2})}$
$\displaystyle+$
$\displaystyle\underbrace{E[\tilde{B}^{*}(f+f_{1}+f_{2})\tilde{C}(f+f_{2})]}_{G_{cb}(f+f_{2})\delta(f_{1})}\underbrace{E[\tilde{D}^{*}(u+f_{3})\tilde{A}(f+f_{1})]}_{G_{ad}(f+f_{1})\delta(u+f_{3}-f-f_{1})}\underbrace{E[\tilde{F}^{*}(u+f_{4})\tilde{E}(u+f_{3}+f_{4})]}_{G_{ef}(u+f_{3}+f_{4})\delta(f_{3})}$
$\displaystyle+$
$\displaystyle\underbrace{E[\tilde{B}^{*}(f+f_{1}+f_{2})\tilde{C}(f+f_{2})]}_{G_{cb}(f+f_{2})\delta(f_{1})}\underbrace{E[\tilde{D}^{*}(u+f_{3})\tilde{E}(u+f_{3}+f_{4})]}_{G_{ed}(u+f_{3}+f_{4})\delta(f_{4})}\underbrace{E[\tilde{F}^{*}(u+f_{4})\tilde{A}(f+f_{1})]}_{G_{af}(f+f_{1})\delta(u+f_{4}-f-f_{1})}$
$\displaystyle+$
$\displaystyle\underbrace{E[\tilde{B}^{*}(f+f_{1}+f_{2})\tilde{E}(u+f_{3}+f_{4})]}_{G_{eb}(u+f_{3}+f_{4})\delta(f+f_{1}+f_{2}-u-f_{3}-f_{4})}\underbrace{E[\tilde{D}^{*}(u+f_{3})\tilde{A}(f+f_{1})]}_{G_{ad}(f+f_{1})\delta(u+f_{3}-f-f_{1})}\underbrace{E[\tilde{F}^{*}(u+f_{4})\tilde{C}(f+f_{2})]}_{G_{cf}(f+f_{2})\delta(u+f_{4}-f-f_{2})}$
$\displaystyle+$
$\displaystyle\underbrace{E[\tilde{B}^{*}(f+f_{1}+f_{2})\tilde{E}(u+f_{3}+f_{4})]}_{G_{eb}(u+f_{3}+f_{4})\delta(f+f_{1}+f_{2}-u-f_{3}-f_{4})}\underbrace{E[\tilde{D}^{*}(u+f_{3})\tilde{C}(f+f_{2})]}_{G_{cd}(f+f_{2})\delta(u+f_{3}-f-f_{2})}\underbrace{E[\tilde{F}^{*}(u+f_{4})\tilde{A}(f+f_{1})]}_{G_{af}(f+f_{1})\delta(u+f_{4}-f-f_{1})}$
$\displaystyle=$

thus

$\displaystyle T$
$\displaystyle=G_{ab}(f+f_{1})G_{cd}(f+f_{2})G_{ef}(u+f_{3}+f_{4})\delta(f_{2})\delta(f_{3})\delta(u+f_{3}-f-f_{2})$
$\displaystyle+G_{ab}(f+f_{1})G_{ed}(u+f_{3}+f_{4})G_{cf}(f+f_{2})\delta(f_{2})\delta(f_{4})\delta(u+f_{4}-f-f_{2})$
$\displaystyle+G_{cb}(f+f_{2})G_{ad}(f+f_{1})G_{ef}(u+f_{3}+f_{4})\delta(f_{1})\delta(f_{3})\delta(u+f_{3}-f-f_{1})$
$\displaystyle+G_{cb}(f+f_{2})G_{ed}(u+f_{3}+f_{4})G_{af}(f+f_{1})\delta(f_{1})\delta(f_{4})\delta(u+f_{4}-f-f_{1})$
$\displaystyle+G_{eb}(u+f_{3}+f_{4})G_{ad}(f+f_{1})G_{cf}(f+f_{2})\cdot$
$\displaystyle\cdot\delta(u+f_{4}-f-f_{2})\delta(u+f_{3}-f-f_{1})\delta(f+f_{1}+f_{2}-u-f_{3}-f_{4})$
$\displaystyle+G_{eb}(u+f_{3}+f_{4})G_{cd}(f+f_{2})G_{af}(f+f_{1})\cdot$
$\displaystyle\cdot\delta(u+f_{4}-f-f_{1})\delta(u+f_{3}-f-f_{2})\delta(f+f_{1}+f_{2}-u-f_{3}-f_{4}).$

Now we use the sampling property of the delta to write, e.g. for the first
line where $f_{2}=0$ and $f_{3}=0$,

$G_{ab}(f+f_{1})G_{cd}(f)G_{ef}(u+f_{4})\delta(f_{2})\delta(f_{3})\delta(u-f)$

and e.g. for the last line where $u+f_{4}=f+f_{1}$ and $u+f_{3}=f+f_{2}$ which
we add up to get

${\color[rgb]{0,1,0}u+f_{3}+f_{4}=(f-u)+f+f_{1}+f_{2}}$

whence

${\color[rgb]{1,0,0}f+f_{1}+f_{2}-u-f_{3}-f_{4}}={\color[rgb]{1,0,0}u-f}$

so that the last line writes as

$\displaystyle
G_{eb}({\color[rgb]{0,1,0}u+f_{3}+f_{4}})G_{cd}(f+f_{2})G_{af}(f+f_{1})\delta(u+f_{4}-f-f_{1})\delta(u+f_{3}-f-f_{2})\delta({\color[rgb]{1,0,0}f+f_{1}+f_{2}-u-f_{3}-f_{4}})=$
$\displaystyle
G_{eb}({\color[rgb]{0,1,0}(f-u)+f+f_{1}+f_{2}})G_{cd}(f+f_{2})G_{af}(f+f_{1})\cdot\delta(u+f_{4}-f-f_{1})\delta(u+f_{3}-f-f_{2})\delta({\color[rgb]{1,0,0}u-f})\stackrel{{\scriptstyle(\mbox{use\,}u=f)}}{{=}}$
$\displaystyle
G_{eb}({\color[rgb]{0,1,0}f+f_{1}+f_{2}})G_{cd}(f+f_{2})G_{af}(f+f_{1})\delta(f_{4}-f_{1})\delta(f_{3}-f_{2})\delta({\color[rgb]{1,0,0}u-f}).$

We therefore get

$\displaystyle T$
$\displaystyle=G_{ab}(f+f_{1})G_{cd}(f)G_{ef}(f+f_{4})\delta(f_{2})\delta(f_{3})\delta(u-f)$
$\displaystyle+G_{ab}(f+f_{1})G_{ed}(f+f_{3})G_{cf}(f)\delta(f_{2})\delta(f_{4})\delta(u-f)$
$\displaystyle+G_{cb}(f+f_{2})G_{ad}(f)G_{ef}(f+f_{4})\delta(f_{1})\delta(f_{3})\delta(u-f)$
$\displaystyle+G_{cb}(f+f_{2})G_{ed}(f+f_{3})G_{af}(f)\delta(f_{1})\delta(f_{4})\delta(u-f)$
$\displaystyle+G_{eb}(f+f_{1}+f_{2})G_{ad}(f+f_{1})G_{cf}(f+f_{2})\delta(f_{4}-f_{2})\delta(f_{3}-f_{1})\delta(u-f)$
$\displaystyle+G_{eb}(f+f_{1}+f_{2})G_{cd}(f+f_{2})G_{af}(f+f_{1})\delta(f_{4}-f_{1})\delta(f_{3}-f_{2})\delta(u-f).$

whence the final form (14) given in Theorem 3.