==== Front
Entropy (Basel)
Entropy (Basel)
entropy
Entropy
1099-4300 MDPI 

33287089
10.3390/e22111324
entropy-22-01324
Article
Extending Fibre Nonlinear Interference Power Modelling to Account for General Dual-Polarisation 4D Modulation Formats
https://orcid.org/0000-0002-8805-7815Liga Gabriele 1* https://orcid.org/0000-0002-2846-7679Barreiro Astrid 1 Rabbani Hami 12 https://orcid.org/0000-0002-2172-3051Alvarado Alex 1 1 Information and Communication Theory Lab, Signal Processing Systems Group, Department of Electrical Engineering, Eindhoven University of Technology, 5600 MB Eindhoven, The Netherlands; a.barreiro.berrio@tue.nl (A.B.); hami.rabbani@email.kntu.ac.ir (H.R.); a.alvarado@tue.nl (A.A.)
2 Department of Electrical Engineering, K. N. Toosi University of Technology, Tehran 1355-16315, Iran
* Correspondence: g.liga@tue.nl
20 11 2020 
11 2020 
22 11 132407 9 2020 27 10 2020 © 2020 by the authors.2020Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).In optical communications, four-dimensional (4D) modulation formats encode information onto the quadrature components of two arbitrary orthogonal states of polarisation of the optical field. Many analytical models available in the optical communication literature allow, within a first-order perturbation framework, the computation of the average power of the nonlinear interference (NLI) accumulated in coherent fibre-optic transmission systems. However, all such models only operate under the assumption of transmitted polarisation-multiplexed two-dimensional (PM-2D) modulation formats, which only represent a limited subset of the possible dual-polarisation 4D (DP-4D) formats. Namely, only those where data transmitted on each polarisation channel are mutually independent and identically distributed. This paper presents a step-by-step mathematical derivation of the extension of existing NLI models to the class of arbitrary DP-4D modulation formats. In particular, the methodology adopted follows the one of the popular enhanced Gaussian noise model, albeit dropping most assumptions on the geometry and statistic of the transmitted 4D modulation format. The resulting expressions show that, whilst in the PM-2D case the NLI power depends only on different statistical high-order moments of each polarisation component, for a general DP-4D constellation, several other cross-polarisation correlations also need to be taken into account.

4D modulation formatsoptical communicationschannel modelling
==== Body
1. Introduction
With the resurgence of polarisation-diverse, optical coherent detection, transmission of information over an optical fibre is typically performed exploiting four degrees of freedom of the optical field: two quadrature components over two orthogonal states of polarisation. The standard approach consists in encoding data independently over the two polarisation channels using the same two-dimensional (2D) modulation format. The resulting four-dimensional (4D) constellation is often referred to as a polarisation-multiplexed 2D (PM-2D) modulation format. The strong point of PM-2D formats is their simplicity of generation and performance analysis: as the two polarisation channels are independent and under the assumption of data-independent cross-polarisation interference in the fibre channel, transmission performance can be evaluated using the 2D component format.

Despite the popularity of PM-2D formats, a substantial amount of research work in the literature has been devoted to more general 4D formats, i.e., 4D constellations which are not necessarily generated as Cartesian products of a component 2D constellation [1,2]. These formats have recently regained attention due to their potential power efficiency, nonlinearity tolerance, and ultimately their still unexplored shaping gains. The reason for that relies on the fact that, by exploiting the full 4D space, sensitivity and other relevant performance metrics such as mutual information or generalised mutual information can be improved compared to traditional PM-2D formats [3,4,5,6,7].

Previous works on optimised 4D modulation formats have either operated under an additive white Gaussian noise channel hypothesis [1,2,3] or exploited some heuristic approaches to derive nonlinearly tolerant formats in the fibre-optic channel [5,6,7]. However, accurately predicting the amount of nonlinear interference generated by transmission of a given constellation in an optical fibre is key to optimising its shape in N dimensions.

Modelling of nonlinear interference (NLI) in optical fibre transmission is quite a mature field of research where an impressive amount of progress was made in the first half of the 2010s, e.g., in [8,9,10,11]. In particular, [10,11] introduced for the first time the possibility of predicting the dependency of nonlinear interference power as a function of the modulation format features, i.e., geometrical shape and statistical properties. Among other assumptions, one underlying key point of all previous models is the transmission of PM-2D modulation formats, where data on the two polarisation channels are assumed to be independent and identically distributed. Under this constraint, one can predict the NLI power using the statistical properties of the 2D component modulation format. It is clear, however, that this approach ceases to be applicable to general DP-4D formats, where a single 2D component format might not even exist.

In this work, we extend the existing analytical expressions for NLI power to account for DP-4D constellations where the two 2D polarisation components are not identically distributed or when there is statistical dependency between them. The undertaken approach is the same as in [9], i.e., a frequency-domain, first-order perturbation study. Unlike [9], no assumptions are made on either the marginal or joint statistics of the two polarisation components of the transmitted 4D constellation (besides being zero-mean). The final expressions reveal the impact of several different (nontrivial) cross-polarisation statistics on the NLI power.

The formulas presented in this work enable an accurate computation of the NLI power for all possible dual-polarisation formats in optical fibre transmission. As a result, a reliable optimisation of both geometry and symbol probability of occurrence of such 4D formats is also enabled for the optical fibre channel.

2. Organisation of the Manuscript and Notation
The manuscript is organised as follows: (i) in Section 3, the investigated system model is described and the model assumptions are presented; (ii) Section 4, Section 5, Section 6, Section 7 and Section 8 are devoted to a step-by-step analytical derivation of the model; and (iii) ultimately, the main model expression is presented in Section 8. In particular: in Section 4, the regular perturbation (RP) solution to the frequency domain Manakov equation is derived for a multi-span fibre system and its power spectral density (PSD) is evaluated in the case of a transmitted periodic signal; in Section 5, the contributions of the different high-order moments and cross-polarisation correlations of the transmitted 4D modulation format are highlighted; finally, Section 8 derives, via Theorem 2, an expression for the PSD as the signal period is extended into infinity. A flowchart of the main derivation steps performed in this work, with their corresponding references in the manuscript, is shown in Figure 1.

Throughout this manuscript, we denote 2D (column) vectors with boldface letters (e.g., a), whereas 2D column vector functions are indicated with boldface capital letters (e.g., E(f,z),E˜(t,z), etc.). For indicating the optical field, the first variable of represents either the time or frequency variable whereas the second one represents the fibre propagation section. An exception is made for the multi-span system case, where second and third variables are assigned to the number of spans and span length, respectively. This highlights the joint dependence of the output optical field on these two variables, as shown later in the paper. F{·}, E{·}, and Re{·} indicate the Fourier transform, the statistical expectation, and the real part operators, respectively. The delta distribution is indicated by δ(·), whereas δk denotes the Kronecker delta defined as
 δk≜1fork=0,0elsewhere. 
Finally, Z, R, and C denote the integer, real, and complex fields, respectively, and j is the imaginary unit.

Figure 1 Flowchart of the analytical derivation in this work.

3. Model Assumptions
3.1. System Model
The baseband equivalent model of the optical fibre system under investigation in this work is shown in Figure 2. The fibre channel is a multi-span fibre system using Erbium-doped fibre amplification (EDFA). In this manuscript, it is assumed that a single-channel signal is transmitted. The transmitter is assumed to generate for each symbol period n the 4D symbol an=[ax,n,ay,n]T, where ax,n,ay,n∈C are complex symbols modulated on two arbitrary orthogonal polarisation states x and y, respectively. The sequence of symbols an for n∈Z is assumed to be a cyclostationary process of period W. The set of random variables (RVs) within each period of such process are also assumed to be statistically independent. Linear modulation with a single, real pulse p(t) on x and y polarisation is adopted. The pulse p(t) with spectrum P(f) is assumed to be strictly band-limited within the range of frequencies [−Rs/2,Rs/2]. As discussed in Section 3.3, the transmitted signal E˜(t,0) is assumed to be periodic with period T, such that
 (1) E˜(t,0)=∑n=−(W−1)/2(W−1)/2anp(t−nTs),for0≤t≤T, 
Ts=1/Rs=T/W represents the symbol period, and Rs is the symbol rate. A schematic representation of the transmitted signal is shown in Figure 3.

Figure 2 System model under investigation in this work which consists of an optical fibre system model and a nonlinear interference (NLI) variance estimation block: the two branches in the NLI variance estimator block indicate alternative ways of estimating the NLI variance ΣNLI.

The signal E˜(t,0) is transmitted over Ns (homogeneous) fibre spans, each of length Ls and each followed by an ideal lumped optical amplifier for which the gain exactly recovers from the span losses. Since in this work we are only concerned about the prediction of NLI arising from the signal–signal nonlinear interactions along the fibre propagation, the optical noise added by the amplifier plays no role in the model and will be entirely neglected. The signal at the channel output E˜(t,Ns,Ls) is ideally compensated for accumulated chromatic dispersion in the link (see Section 4). In the frequency-domain output of the chromatic dispersion compensation (CDC) block E˚(f,Ns,Ls) (Figure 2), we ideally isolate the first-order RP term E1(f,Ns,Ls) (see Section 4) and we compute its PSD S(f,Ns,Ls). The vector of the NLI powers ΣNLI≜[σNLI,x2,σNLI,y2]T for both x and y polarisations is obtained by integrating over the frequency interval [−Rs/2,Rs/2] the NLI PSD weighted by the function |P(f)|2, where P∗(f) is the frequency response of a matched filter (MF) for the system under consideration. As shown in Figure 2, this quantity is equivalent to the variance of the output of the MF followed by symbol-rate sampling, which more naturally arises when assessing the transmission performance of systems employing an MF at the receiver. The model in this manuscript provides an analytical relationship between the statistical features of the transmitted symbols an and ΣNLI.

3.2. DP-4D vs. PM-2D Formats
The model presented in this paper allows for prediction of the NLI for generic 4D real modulation formats. A 4D format is defined as a set
 (2) A≜{a(i)=[ax(i),ay(i)]T∈C2,i=1,2,…,M}, 
where ax and ay are the symbols modulated on two orthogonal polarisation states x and y, respectively, and M is the modulation cardinality. It can be seen that the elements in A are 2D vectors in C as opposed to 4D. This is only due the to baseband-equivalent representation of signals used throughout this paper, while it is common to refer to a modulation format dimensionality based on the real signal dimensions, which justifies the 4D format label.

Two important particular cases of the formats in (2) are (i) the so-called polarisation-multiplexed 2D (PM-2D) modulation formats, which are characterised by A=X2, X∈C, where X represents the 2D component constellation, and (ii) polarisation-hybrid 2D modulation formats characterized by A=X×Y, with X,Y∈C, X≠Y, where X and Y are two distinct component 2D formats in x and y polarisation, respectively. PM-2D formats are the most common ones in optical communications due to their generation’s simplicity. Both PM-2D and polarisation-hybrid 2D formats are often analysed in terms of their 2D polarisation components. This is because A can be factorised in two component formats which are independently encoded. Hence, if the generic transmitted constellation point is regarded as a random variable, in a conventional PM-2D format, the two polarisation components are statistically independent. In the remainder of this paper, no specific assumption on either the geometry or the statistic of the transmitted 4D symbols will be made, except the zero-mean feature E{a(i)}=0.

Figure 3 Schematic representation of the periodic signal assumption, where W symbols are transmitted every T [s], each symbol with a duration of Ts [s]: the periodicity assumption will be lifted in Section 8 by letting Δf→0.

3.3. Transmitted Signal Form
Let E˜(t,z)=E˜x(t,z)ix+E˜y(t,z)iy be the complex envelope of the optical field vector at time t and fibre section z, and let ix, iy denote 2 orthonormal polarisations of the transversal plane of propagation. Let also E(f,z)=Ex(f,z)ix+Ey(f,z)iy be the (vector) Fourier transform of E˜(t,z) defined as
 E(f,z)=F{E˜(t,z)}≜∫−∞∞E˜(t,z)e−j2πftdt. 

Because of the periodicity assumption made in (1) (see Figure 3), we can write E˜(t,0) as
 (3) E˜(t,0)=∑k=−∞∞Ckej2πkΔft, 
where Ck=[Cx,k,Cy,k]T, Cx/y,k are the Fourier series coefficients of E˜(t,0) and Δf=1/T is the frequency spacing of the spectral lines in Ex/y(f,z). Hence, E(f,0) can be then written as
 (4) E(f,0)=∑k=−∞∞Ckδ(f−kΔf). 
Since each component of E˜(t,0) is periodic with period T, we can write
 E˜(t,0)=∑n=−∞∞E^(t−nT,0), 
where, as per assumption in (1), we have
 E^(t,0)≜∑n=−(W−1)/2(W−1)/2anp(t−nTs),for−T2≤t≤T20,otherwise, 
and, W is assumed to be odd without loss of generality. Under the above assumptions, the Fourier coefficients in (3), for k∈Z, are given by
 (5a) Ck=Δf∫−T2T2E^(t,0)e−j2πkΔftdt 
 (5b) =Δf∫−T2T2∑n=−(W−1)/2(W−1)/2anp(t−nTs)e−j2πkΔftdt 
 (5c) =Δf∑n=−(W−1)/2(W−1)/2an∫−T2T2p(t−nTs)e−j2πkΔftdt 
 (5d) ≈Δf∑n=−(W−1)/2(W−1)/2anP(kΔf)e−j2πkΔfnTs 
 (5e) =ΔfP(kΔf)∑n=−(W−1)/2(W−1)/2ane−j2πknW 
 (5f) =ΔfP(kΔf)νk, 
where P(f)≜F{p(t)} and
 (6) νk=[νx,k,νy,k]T=Δf∑n=−(W−1)/2(W−1)/2ane−j2πknW,∀k∈Z, 
are the discrete Fourier transforms of the sequence an, n=0,1,…,W−1. Note that the approximation in (5c)–(5d) is justified only for large enough values of T as
 limT→∞∫−T2T2p(t−nTs)e−j2πkΔftdt=F{p(t−nTs)}|f=kΔf,forn,k∈Z, 
and letting T→∞ will be the approach taken at a later stage in this derivation.

Finally, combining (4) and (5f), we obtain
 (7) E(f,0)=Δf∑k=−∞∞P(kΔf)νkδ(f−kΔf)≈∑k=−(W−1)/2(W−1)/2P(kΔf)νkδ(f−kΔf), 
where the approximate equality on the right-hand side of (7) stems from the fact that p(t) is assumed to be strictly or quasi-strictly band-limited (see Section 3.1). Hence, P(kΔf) is effectively equal to zero for k=−W/2,−W/2+1,…,W/2.

4. PSD of the First-Order NLI for Periodic Transmitted Signals
To find an analytical expression for NLI power, first, a solution as explicit as possible to the Manakov equation [12]
 (8) ∂E˜(t,z)∂z=−α2E˜(t,z)−jβ22∂2E˜(t,z)∂t2+j89γ|E˜(t,z)|2E˜(t,z), 
must be found. Equation (8) describes the propagation of the optical field E˜(t,z) in a single strand of fibre (e.g., a fibre span with no amplifier in the system in Figure 2). In this case, α, β2, and γ representing the attenuation, group velocity dispersion, and nonlinearity coefficients, respectively, can be assumed to be spatially constant. As it is well-known, general closed-form solutions are not available for (8). Like most of the existing NLI power models in the literature, the model derived here operates within a first-order perturbative framework. In particular, a frequency-domain first-order regular perturbation (RP) approach in the γ coefficient is performed [13,14], i.e., the Fourier transform of the solution in (8) is expressed as
 (9) E(f,z)=∑n=0∞γnAn(f,z)≈A0(f,z)+γA1(f,z), 
where
 (10) En(f,z)=γnAn(f,z)forn=0,1,…, 
represents the so-called nth order term of the expansion.

In the following theorem, we present the expressions for E0(f,z), and E1(f,z), when a multiple fibre span system like the one in Figure 2 is considered. These expressions are well-known in the literature (see, e.g., [13]). Nevertheless, we present the proof in Appendix A for completeness.

Theorem 1 (First-order frequency-domain RP solution for a multi-span fibre system).
Let E(f,z) be the solution in the frequency domain of the Manakov equation for the system in Figure 2 with initial condition at distance z=0 given by the transmitted signal E(f,0). Then, the first-order RP solution after Ns spans E(f,Ns,Ls) is given by
 E(f,Ns,Ls)≈E0(f,Ns,Ls)+E1(f,Ns,Ls), 
where the zeroth-order term is given by
 E0(f,Ns,Ls)=E0(f,0)ej2π2f2β2NsLs, 
and the first-order term is
 (11) E1(f,Ns,Ls)=−j89γej2π2f2β2NsLs∫−∞∞∫−∞∞ET(f1,0)E∗(f2,0)E(f−f1+f2,0)η(f1,f2,f,Ns,Ls)df1df2, 
with
 (12) η(f1,f2,f,Ns,Ls)≜1−e−αLsej4π2β2(f−f1)(f2−f1)Lsα−j4π2β2(f−f1)(f2−f1)∑l=1Nse−j4π2β2(l−1)(f−f1)(f2−f1)Ls, 
where Ns and Ls are the number of spans and the span length of each span, respectively.

Proof.  See Appendix A. □

While Theorem 1 gives an approximation for the field at the output of the fibre, we are interested in the field after ideal CDC (see Figure 2). Ideal CDC ideally removes the exponential ej2β2π2f2NsLs from (11), leading to a first-order term in the RP solution for the system in Figure 2 given by
 (13) E˚1(f,Ns,Ls)=[E˚1,x,E˚1,y]T=−j89γ∫−∞∞∫−∞∞ET(f1,0)E∗(f2,0)E(f−f1+f2,0)·η(f1,f2,f,Ns,Ls)df1df2. 

Substituting the spectrum of the transmitted periodic signal (7) in (13), we obtain, for instance, for the x component in (13),
 (14) E˚1,x(f,Ns,Ls)=−j89γΔf3/2∑k=−∞∞∑m=−∞∞∑n=−∞∞P(kΔf)P∗(mΔf)P(nΔf)νx,kνx,m∗νx,n+νy,kνy,m∗νx,n·∫−∞∞∫−∞∞δ(f1−kΔf)δ(f2−mΔf)δ(f−f1+f2−nΔf)η(f1,f2,f,Ns,Ls)df1df2. 
Although the product of Dirac’s deltas in (14) is not well-defined in the standard distribution theory framework, in this case, such product can be dealt with in the same way as products between distributions and smooth functions. This approach was formalised by Colombeau in his theory of product between distributions [15]. Thus, integrating in f1 and f2, we obtain
 (15) E˚1,x(f,Ns,Ls)=−j89γΔf3/2∑k=−∞∞∑m=−∞∞∑n=−∞∞P(kΔf)P∗(mΔf)P(nΔf)·νx,kνx,m∗νx,n+νy,kνy,m∗νx,nη(kΔf,mΔf,(k−m+n)Δf,Ns,Ls)δ(f−(k−m+n)Δf). 

Setting i=k−m+n and defining
 (16) ηk,m,n≜η(kΔf,mΔf,(k−m+n)Δf,Ns,Ls)=1−e−αLsej4π2Δf2β2(n−m)(m−k)Lsα−j4π2Δf2β2(n−m)(m−k)∑l=1Nse−j4π2Δf2β2(l−1)(n−m)(m−k)Ls 
(15) can be rewritten as
 (17) E˚1,x(f,Ns,Ls)=∑i=−∞∞ciδ(f−iΔf), 
where
 (18) ci≜−j89γΔf3/2∑(k,m,n)∈SiP(kΔf)P∗(mΔf)P(nΔf)νx,kνx,m∗νx,n+νy,kνy,m∗νx,nηk,m,n, 
and
 (19) Si≜{(k,m,n)∈Z3:k−m+n=i}. 

The PSD of the received nonlinear interference (to the 1st-order) is defined as
 S(f,Ns,Ls)=[Sx(f,Ns,Ls),Sy(f,Ns,Ls)]T≜E|E˚1,x(f,Ns,Ls)|2,E|E˚1,y(f,Ns,Ls)|2T. 
For periodic signals, which in the frequency domain can be expressed as in (17), the PSD can be expressed as [16] (Section 4.1.2)
 (20) Sx(f,Ns,Ls)=∑i=−∞∞E{|ci|2}δ(f−iΔf). 

Substituting the expression (18) for ci in (20), we obtain
 Sx(f,Ns,Ls)=892γ2Δf3∑i=−∞∞δ(f−iΔf)E{∑(k,m,n)∈SiP(kΔf)P∗(mΔf)P(nΔf)νx,kνx,m∗νx,n 
 (21a) +νy,kνy,m∗νx,nηk,m,n∑(k′,m′,n′)∈SiP∗(k′Δf)P(m′Δf)P∗(n′Δf)νx,k′∗νx,m′νx,n′∗+νy,k′∗νy,m′νx,n′∗ηk′,m′,n′∗}=892γ2Δf3∑i=−∞∞δ(f−iΔf)E{∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′νx,kνx,m∗νx,nνx,k′∗νx,m′νx,n′∗ 
 (21b) +νx,kνx,m∗νx,nνy,k′∗νy,m′νx,n′∗+νy,kνy,m∗νx,nνx,k′∗νx,m′νx,n′∗+νy,kνy,m∗νx,nνy,k′∗νy,m′νx,n′∗ηk,m,nηk′,m′,n′∗}, 
where we have defined
 (22) Pk,m,n,k′,m′,n′≜P(kΔf)P∗(mΔf)P(nΔf)P∗(k′Δf)P(m′Δf)P∗(n′Δf). 

The following proposition can be used to make (21b) more compact, and in particular, it will be used to group the two inner correlation terms in (21b) (νx,kνx,m∗νx,nνy,k′∗νy,m′νx,n′∗ and νy,kνy,m∗νx,nνx,k′∗νx,m′νx,n′∗).

Proposition 1. For Pk,m,n,k′,m′,n′ in (22), we have
 (23) ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′νx,kνx,m∗νx,nνy,k′∗νy,m′νx,n′∗ηk,m,nηk′,m′,n′∗=(∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′νy,kνy,m∗νx,nνx,k′∗νx,m′νx,n′∗ηk,m,nηk′,m′,n′∗)∗. 


Proof.  See Appendix B. □

Using (23) and (21b) can be written as
 (24) Sx(f,Ns,Ls)=892γ2Δf3∑i=−∞∞δ(f−iΔf)E{∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′·νx,kνx,m∗νx,nνx,k′∗νx,m′νx,n′∗+νy,kνy,m∗νx,nνy,k′∗νy,m′νx,n′∗ηk,m,nηk′,m′,n′∗+2Re{Pk,m,n,k′,m′,n′νx,kνx,m∗νx,nνy,k′∗νy,m′νx,n′∗ηk,n,mηk′,n′,m′∗}}=892γ2Δf3∑i=−∞∞δ(f−iΔf)∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′Eνx,kνx,m∗νx,nνx,k′∗νx,m′νx,n′∗+Eνy,kνy,m∗νx,nνy,k′∗νy,m′νx,n′∗ηk,m,nηk′,m′,n′∗+2Re{Pk,m,n,k′,m′,n′Eνx,kνx,m∗νx,nνy,k′∗νy,m′νx,n′∗ηk,n,mηk′,n′,m′∗}, 
where the real part operator arises from the sum of the complex conjugate terms discussed in the Proposition section (Section 1).

According to (24), calculation of the PSD of the NLI reduces to the computation of a four-dimensional summation (per frequency component iΔf) of three sixth-order correlations of the sequence of random variables νx/y,n,n=0,1,…,W−1. The y-component Sy(f,Ns,Ls) of the PSD can be calculated once Sx(f,Ns,Ls) is obtained by simply swapping the polarisation labels x→y and y→x. This is due to the invariance of the Manakov equation in (8) to such a transformation.

5. Classification of the Modulation-Dependent Contributions in the 6th-Order Frequency-Domain Correlation
In this section, we will break down the frequency-domain sixth-order correlation terms in (24) to highlight different contributions in terms of 4D modulation-dependent cross-polarisation correlations.

5.1. Expansion in Terms of the Stochastic Moments of the Transmitted Modulation Format
To relate the PSD in (24) to the statistical properties of the transmitted modulation format, we replace (6) into (24), obtaining
 (25) Sx(f,Ns,Ls)=892γ2Δf3∑i=−∞∞δ(f−iΔf)∑(k,m,n)∈Si(k′,m′,n′)∈Si[Pk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗·∑i∈{0,1,…,W−1}6Si(k,m,n,k′,m′,n′)+2Re{Pk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗·∑i∈{0,1,…,W−1}6Ti(k,m,n,k′,m′,n′)}], 
where i≜(i1,i2,…,i6),
 (26) Si(k,m,n,k′,m′,n′)≜Δf3Eax,i1ax,i2∗ax,i3ax,i4∗ax,i5ax,i6∗+Eay,i1ay,i2∗ax,i3ay,i4∗ay,i5ax,i6∗·e−j2πW(ki1−mi2+ni3−k′i4+m′i5−n′i6), 
and
 (27) Ti(k,m,n,k′,m′,n′)≜Δf3Eax,i1ax,i2∗ax,i3ay,i4∗ay,i5ax,i6∗e−j2πW(ki1−mi2+ni3−k′i4+m′i5−n′i6). 

The terms Si(k,m,n,k′,m′,n′) and Ti(k,m,n,k′,m′,n′) give rise to several correlations among the transmitted symbols ax,i and ay,j at different time-slots i,j, each weighted by a complex exponential. As discussed in Section 3, in this work, we operate under the assumption that the sequence of vector RVs ai for i=0,1,…,W−1 are independent, identically distributed (i.i.d.), and with E{ai}=E{a}=0. As shown in the following example, this assumption allows us to discard the Si(k,m,n,k′,m′,n′) and Ti(k,m,n,k′,m′,n′) terms which are identically zero for some values of i. Moreover, as it will be shown in Example 2 for all other values of i, Si(k,m,n,k′,m′,n′), and Ti(k,m,n,k′,m′,n′) can be expressed as a product of high-order statistical moments of the RVs ax and ay, which enables a more compact expression for (25).

Example 1. Under the i.i.d. assumption for the sequence of vector RVs ai, i=0,1,…W−1 made in this work, in any of the cases where
 (28) iκ1≠iκ2=iκ3=…=iκ6forκ1,κ2,…,κ6=1,2,…,6;κ1≠κ2≠…≠κ6, 
any of the sixth-order correlations in (26) and (27) degenerate into a product between a first-order moment and a fifth-order correlation. Such a product is equal to zero under our assumption E{ax,i}=E{ax}=0. For example, for i1≠i2=i3=…=i6, we have
 E{ax,i1ax,i2∗ax,i3ax,i4∗ax,i5ax,i6∗}=E{ax,i1}E{|ax,i2|4ax,i2∗}=E{ax}E{|ax|4ax∗}=0. 
From this follows that, for all elements in the set defined in (28), Si(k,m,n,k′,m′,n′)=0 and Ti(k,m,n,k′,m′,n′)=0. This example highlights a zero-contribution region in the 6D space
{0,1,…,W−1}6, as illustrated in Figure 4.

The Si(k,m,n,k′,m′,n′) and Ti(k,m,n,k′,m′,n′) contributions for the set in (28) are identically zero regardless of the values taken by k,m,n,k′,m′, and n′. However, as it will be shown in Section 6, for a specific subset of values k,m,n,k′,m′, and n′, such contributions cancel each other in the inner sums in (25) due to the complex exponential weights.

Example 2. Under the i.i.d. assumption for the sequence of vector RVs ai, i=0,1,…W−1 made in this work, we have that, for all elements in the set {i∈{0,1,…,W−1}6,i1=i2,i3=i4=i5=i6,i1≠i3}
 Si(k,m,n,k′,m′,n′)=Δf3[E{|ax,i1|2}E{|ax,i3|4}+E{|ay,i1|2}E{|ax,i3|2|ay,i3|2}]·e−j2πW((k−m)i1+(n−k′+m′−n′)i3)=Δf3[E{|ax|2}E{|ax|4}+E{|ay|2}E{|ax|2|ay|2}]e−j2πW((k−m)i1+(n−k′+m′−n′)i3),Ti(k,m,n,k′,m′,n′)=Δf3E{|ax,i1|2}E{|ax,i3|2|ay,i3|2}e−j2πW(ki1−mi2+ni3−k′i4+m′i5−n′i6)=Δf3E{|ax|2}E{|ax|2|ay|2}e−j2πW(ki1−mi2+ni3−k′i4+m′i5−n′i6). 
It can be noted that (i) the sixth-order correlation degenerates into a product of marginal high-order moments of ax and ay and into the cross-polarisation correlation E{|ax|2|ay|2} and (ii) all elements within the set in this example contribute to the inner summation in (25) with the same set of moments, cross-polarisation correlations, and products thereof (i.e., E{|ax|2},E{|ax|4},E{|ay|2},andE{|ax|2|ay|2}). This example illustrates how to break down each instance of the contributions Si(k,m,n,k′,m′,n′) and Ti(k,m,n,k′,m′,n′), which will be then added up in Section 6.


In the remainder of this section, we first partition the six-dimensional space i∈{0,1,…,W−1}6 and list all sets corresponding to nonzero elements of Si(k,m,n,k′,m′,n′) and Ti(k,m,n,k′,m′,n′). As shown in Example 2, this will help highlight the contribution of a specific set in terms of high-order moments of the transmitted symbols a in (25). Then, we proceed to list all such contributions.

5.2. Set Partitioning
The six-dimensional space i∈{0,1,…,W−1}6 can be partitioned in different subsets, each one uniquely defined by a partition on the set of indices (i1, i2, i3, i4, i5, and i6). Each partition defines its corresponding subset in {0,1,…,W−1}6 as follows: for each index partition, the indices belonging to the same subset all take the same value, whilst the indices belonging to different subsets have distinct values. This is schematically illustrated in Figure 4. For example, the subset of {0,1,…,W−1}6 labelled by the index partition {(i1,i2),(i3,i4),(i5,i6)} is defined as {i∈{0,1,…,W−1}6:i1=i2,i3=i4,i5=i6,i1≠i3≠i5}. This subset is shown in Figure 4 as part of L1.

In Figure 4, the families of subsets of {0,1,…,W−1}6 labelled Li, i=1,2,…,4, are also highlighted. These families are characterised by subsets sharing the same cardinality of elements associated to their corresponding index partition. For example, in L1, all index partitions are characterised by 3 subsets, each one containing 2 indices. As shown in Example 2, this way of partitioning the set {0,1,…,W−1}6 is useful as it separates out the different contributions of (25) based on the high-order moments of a, as it is highlighted in region L3 of Figure 4.

Since we have 6 different indices, the number of subsets in a partition can vary from 1 to 6. Each of these subsets can contain a number of elements also ranging from 1 to 6. However, the subsets of {0,1,…,W−1}6, where the corresponding index partition has one or more index subsets with only one element, bring no contribution to (25) and thus can be discarded. This is illustrated in Example 1. The above class of index partitions then forms a zero contribution region, as shown in Figure 4. Such a region also includes all subsets where the corresponding index partitions contain 4 or more index subsets, as at least one of these subsets will have to contain only one element.

As shown in Figure 4, by removing the zero contribution region from {0,1,…,W−1}6, only 4 different families of subsets are left:(i) L1={i∈{0,1,…,W−1}6:iκ1=iκ2,iκ3=iκ4,iκ5=iκ6;κ1,κ2,…κ6=1,2,…,6;κ1≠κ2≠κ3≠κ4≠κ5≠κ6}. This set contains all sets of elements where the indices i1,i2,…,i6, can be grouped in 3 pairs. The indices take up the same value within each pair but different values across different pairs. It can be found that this set can be partitioned in 15 different subsets C1(i),i=1,2,…,15, representing all possible distinct ways of pairing the ik indices for k=1,2,…6. These sets are listed in Table A1 in Appendix C, where each column shows a subgroup of indices taking the same value;

(ii) L2:{i∈{0,1,…,W−1}6,iκ1=iκ2=iκ3,iκ4=iκ5=iκ6;κ1,κ2,…,κ6=1,2,…,6;κ1≠κ2≠κ3≠κ4≠κ5≠κ6} which can be broken down in 10 subsets C2(i),i=1,2,…,10, listed in Table A2 in Appendix C. Each index subgroup identifies a triplet of indices assuming the same value;

(iii) L3={i∈{0,1,…,W−1}6:iκ1=iκ2,iκ3=iκ4=iκ5=iκ6;κ1,κ2,…,κ6=1,2,…,6,κ1≠κ2≠κ3≠κ4≠κ5≠κ6} which can be partitioned in 15 subsets C3(i),i=1,2,…,15, listed in Table A3 in Appendix C. Each of the two index subgroups identifies the pair and the quadruple of indices assuming the same value;

(iv) L4:{i∈{0,1,…,W−1}6:i1=i2=i3=i4=i5=i6}.



6. Evaluation of the L-Based Contributions
In this section, we provide three examples for the computation of the contributions of a generic element in L1, L2, and L3. The full list of contributions in these three sets and the ones in L4 are given in Section 6.1, Section 6.2, Section 6.3, Section 6.4.

We label each contribution as Mg(h)(k,m,n,k′,m′,n′) and Ng(h)(k,m,n,k′,m′,n′), where
 (29) Mg(h)(k,m,n,k′,m′,n′)≜∑i∈Cg(h)Si(k,m,n,k′,m′,n′),Ng(h)(k,m,n,k′,m′,n′)≜∑i∈Cg(h)Ti(k,m,n,k′,m′,n′), 
and the subsets Cg(h) are taken from Table A1, Table A2, Table A3 in Appendix C.

Example 3 (Contributions in L1).
M1(1), i.e., one of the 2 contributions for the set C1(1)={i∈{0,1,…,W−1}6:i1=i2,i3=i4,i5=i6,i1≠i3,i1≠i5,i3≠i5} is given by
 (30) M1(1)≜∑i∈C1(1)Si(k,m,n,k′,m′,n′)=Δf3E3|ax|2+E|ay|2|Eaxay∗|2∑i1=0W−1e−j2πW(k−m)i1∑i3≠i1e−j2πW(n−k′)i3·∑i5≠i1,i5≠i3e−j2πW(m′−n′)i5. 

Since
 (31) ∑k=0W−1ejnk2πW=W,forn=pW,p∈Z0,elsewhere, 
we can compute (30) using the following approach:
we add up the terms for all i1,i3,i5 values including all cases when i1, i3, and i5 are equal among each other. Because of (31), these terms sum up to W3 only when k=m+pW,n=k′+pW,andm′=n′+pW, p∈Z; otherwise, they sum to 0;

we subtract the terms corresponding to the cases: i1=i3,i1≠i5; i1=i5,i1≠i3; and i3=i5,i1≠i3. As an example, the number of terms defined by i1=i3,i1≠i5 is given by the difference between the number of all pairs i1,i5∈{0,1,2,…,W−1} and the number of terms for i1=i5. According to (31), the former terms sum to W2 only for k−m+n−k′=pW,m′−n′=pW, whereas the latter sum to W only for k−m+n−k′+m′−n′=pW, with p∈Z. In all other cases, they all bring zero contribution. Similar results are obtained for i1=i5,i1≠i3 and i3=i5,i1≠i3;

we finally subtract the terms i1=i3=i5, which sum to W only for k−m+n−k′+m′−n′=pW, p∈Z; otherwise, they sum to 0 (see (31)).


Hence, we obtain
 M1(1)=Δf3[E3{|ax|2}+2E2{|ax|2}|E{axay∗}|2+E{|ay|2}|E{axay∗}|2][W3δk−m−pWδn−k′−pWδm′−n′−pW−[W2(δk−m+n−k′−pWδm′−n′−pW+δk−m+m′−n′−pWδn−k′−pW+δm′−n′+n−k′−pWδk−m−pW)−3Wδk−m+m′−n′+n−k′−pW]−Wδk−m+m′−n′+n−k′−pW]=[E3{|ax|2}+2E2{|ax|2}|E{axay∗}|2+E{|ay|2}|E{axay∗}|2][Rs3δk−m+pWδn−k′−pWδm′−n′−pW−Rs2Δf(δk−m+n−k′−pWδm′−n′−pW+δk−m+m′−n′−pWδn−k′−pW+δm′−n′+n−k′−pWδk−m−pW)+2RsΔf2δk−m+m′−n′+n−k′−pW], 
where we have used Rs=WΔf.

The same approach can be followed to compute N1(1), which is, thus, given by
 N1(1)≜∑i∈C1(1)Ti(k,m,n,k′,m′,n′)=E2{|ax|2}|E{axay∗}|2[Rs3δk−m−pWδn−k′−pWδm′−n′−pW−Rs2Δf(δk−m+n−k′−pWδm′−n′−pW+δk−m+m′−n′−pWδn−k′−pW+δm′−n′+n−k′−pWδk−m−pW)+2RsΔf2δk−m+m′−n′+n−k′−pW]. 
All other contributions in L1 can be computed using the approach used in this example.

Example 4 (Contributions in L2).
M2(1), i.e., the contribution for the set C2(1)={i∈{0,1,…,W−1}6:i1=i2=i3,i4=i5=i6,i1≠i4} is given by
 (32) M2(1)=∑i∈C2(1)Si(k,m,n,k′,m′,n′)=Δf3[E{ax,i1|ax,i1|2}E∗{ax,i4|ax,i4|2}+E{ax,i1|ay,i1|2}E∗{ax,i4|ay,i4|2}]·∑i1=0W−1e−j2πW(k−m+n)i1∑i4≠i1e−j2πW(−k′+m′−n′)i4=Δf3[|E{ax|ax|2}|2+|E{ax|ay|2}|2]∑i1=0W−1e−j2πW(k−m+n)i1∑i4≠i1e−j2πW(−k′+m′−n′)i4. 

Following a similar approach as in Example 3, we compute (32) by
adding up the terms for all i1 and i4 values including all cases when i1 and i4 are equal to each other. These terms sum up to W2 only when k−m+n=pW and −k′+m′−n′=pW, with p∈Z; otherwise, they sum to 0;

subtracting the terms corresponding to the cases i1=i4. These terms sum to W only for k−m+n−k′+m′−n′=pW, p∈Z; otherwise, they sum to zero.


We, thus, obtain
 M2(1)=[|E{ax|ax|2}|2+|E{ax|ay|2}|2][Rs2Δfδk−m+n−pWδk′−m′+n′−pW−RsΔf2δk−m+n−k′+m′−n′−pW], 


Following the same approach for N2(1), we have
 N2(1)≜∑i∈C1(1)Ti(k,m,n,k′,m′,n′)=E{ax|ax|2}E{ax|ay|2}[Rs2Δfδk−m+nδk′−m′+n′−pW−RsΔf2δk−m+n−k′+m′−n′−pW]. 
All other contributions in L2 can be computed using the approach used in this example.

Example 5 (Contributions in L3).
M3(3), i.e., the contribution for the values in the set C3(3)={i∈{0,1,…,W−1}6:i1=i4,i2=i3=i5=i6,i1≠i4,i2≠i3,i5≠i6} is given by
 (33) M3(3)=∑i∈C3(3)Si(k,m,n,k′,m′,n′)=Δf3[E{|ax,i1|2}E{|ax,i2|4}+E{|ay,i1|2}E{|ax,i2|2|ay,i2|2}]∑i1=1W−1e−j2πW(k−k′)i1∑i2≠i1e−j2πW(−m+n+m′−n′)i2=Δf3[E{|ax|2}E{|ax|4}+E{|ay|2}E{|ax|2|ay|2}]∑i1=1W−1e−j2πW(k−k′)i1∑i2≠i1e−j2πW(−m+n+m′−n′)i2. 
As in the L2 case described in Example 4, in L3, each subset is characterized by 2 subgroups of indices. Hence, the approach followed to compute (33) is identical to (32) and gives
 M3(3)=[E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}][Rs2Δfδk−k′−pWδm−n−m′+n′−pW−RsΔf2δk−m+n−k′+m′−n′−pW]. 
Similarly,
 N3(3)=E{axay∗}E{ax∗ay|ax|2}[Rs2Δfδk−k′−pWδm−n−m′+n′−pW−RsΔf2δk−m+n−k′+m′−n′−pW]. 
All other contributions in L3 can be computed using the approach used in this example.

As shown in the above examples, each contribution Mg(h) and Ng(h) is nonzero only for a specific set of (k,m,n,k′,m′,n′) values which is spanned by p∈Z. However, the terms (k,m,n,k′,m′,n′) arising for all p≠0 bring a total contribution to (25) that can be considered negligible. This is due to our assumption on P(f) being strictly band-limited (see Section 3.1) and to the magnitude of the functions product ηk,m,nηk′,m′,n′∗ (see definitions (16) and (22)). Thus, in the computations performed in the following subsections, we will restrict ourselves to the case p=0.

6.1. Contributions in L1
In this section, the contributions M1(i), N1(i) for i=1,2,…,15 are computed following Example 3. These contributions are listed in Table 1.

6.2. Contributions in L2
Following Example 4, the contributions M2(h),N2(h),h=1,2,…,10, are computed and listed in Table 2.

6.3. Contributions in L3
Following Example 5, the contributions M3(h),N3(h),h=1,2,…,15, are computed and listed in Table 3.

6.4. Contributions in L4
Since L4 comprises a single subset characterised by the single subgroup of all 6 indices (see Section 5.2), only one pair of contributions M4(1), N4(1) exists, and it is given by
 M4(1)≜∑i∈C4(1)Si(k,m,n,k′,m′,n′)=∑i1=0W−1[E{|ax|6}+E{|ax|2|ay|4}]e−j2πW(k−m)i1=[E{|ax|6}+E{|ax|2|ay|4}]RsΔf2δk−m+n−k′+m′−n′,N4(1)≜∑i∈C4(1)Si(k,m,n,k′,m′,n′)=E{|ax|4|ay|2}RsΔf2δk−m+n−k′+m′−n′. 

7. Sum of All Contributions
In Section 5, we evaluated all contributions Mg(h) and Ng(h) to the PSD in (25). In particular, from (25)–(27), and (29), we have
 (34) Sx(f,Ns,Ls)=892γ2Δf3∑i=−∞∞δ(f−iΔf)∑(k,m,n)∈Si(k′,m′,n′)∈SiP∑g=14∑h=1H(g)Mg(h)+2ReP∑g=14∑h=1H(g)Ng(h), 
where H(g) is the number of subsets in the partitions of Lg, g=1,2,3,4, described in Section 5.2 (H(g)=15,10,15,1 for g=1,2,3,4, resp.) and
 (35) P≜Pk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗. 
In this section, we evaluate ∑g=14∑h=1H(g)Mg(h) and ∑g=14∑h=1H(g)Ng(h) as well as compact the resulting expression as much as possible.

Before we proceed with computing the abovementioned summation, we remove the Kronecker deltas in Mg(h) and Ng(h) corresponding to contributions in the following subspaces: (i) k=m; (ii) n=m; (iii) k′=m′; and (iv) n′=m′. These contributions correspond to so-called bias terms, i.e., they arise from a component of the field E1(f,z) which is fully correlated with the transmitted field E0(t,0). This component, after CDC and MF, only results in a deterministic and static complex scaling of the received constellation, which is typically compensated at the receiver even in the presence of other noise sources in the system. Thus, it does not contribute to the power of the additive zero-mean interference component that we observe at the output of the MF + sampling stage once the received constellation is synchronised (in phase and amplitude) with the transmitted one. A more detailed discussion on these bias terms can be found in [8] (Appendix A), [14] (Appendix C). Moreover, the component δk−m+nδk′−m′+n′ in Mg(h) and Ng(h) is also removed as it only gives nonzero contribution to the PSD for frequency f=iΔf=0; hence, its effect on the total NLI variance vanishes as we let Δf→0 (see Section 8). A total of 23 terms from the last columns of Table 1, Table 2 and Table 3 are thus removed. The remaining contributions are given in Table A4 in Appendix D.

We now compact the contributions in Table A4 by grouping the Kronecker delta products based on each correlation term they multiply. We use three pairs of curly brackets {·} to denote the terms multiplying Rs3, Rs2, and Rs. The list of all Kronecker delta products multiplying each correlation term is shown in Table 4. The correlation terms are divided into intra-polarisation (expectations containing only ax) and cross-polarisation terms (expectations containing both ax and ay). Moreover, the correlations are categorised based on the specific contribution (either M or N) in (34) to which they belong.

As it can be observed in Table 4, each correlation term is associated with different delta functions. To compact these terms, we exploit a property introduced in the following proposition.

Proposition 2. Let D1(k,m,n,k′,m′,n′) and D2(k,m,n,k′,m′,n′) be two Kronecker delta products of the kind shown in Table 4. If
 (36) D1(k,m,n,k′,m′,n′)=D2(n,m,k,n′,m′,k′), 
then
 (37) ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D1(k,m,n,k′,m′,n′)=∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D2(k,m,n,k′,m′,n′). 
This property also holds when applying the transformations k=n, n=k, k′=n′, and n′=k′, individually.

Proof.  See Appendix E. □

The property in (37) allows us to group many of the Kronecker function products in Table 4 under a single term. Namely, the Kronecker delta products in Table 4 can be grouped in subsets that are closed to property (36), since they all result in the same value of summations in (37). In particular, 14 distinct subsets can be identified for the list of Kronecker delta products in Table 4. We label these subsets, which are shown in Table 5, as Dl for l=1,2,…,14.

Summing all the contributions in Table 4, using Proposition 2 for the elements in the subsets listed in Table 4, and finally ordering by Kronecker delta product, we obtain from (34)
 (38) Sx(f,Ns,Ls)=892γ2Δf∑i=−∞∞δ(f−iΔf)·∑(k,m,n)∈Si(k′,m′,n′)∈Si[Rs3Δf2[(a1P+2Re{a1′P})δk−k′δm−m′δn−n′+(a2P+2Re{a2′P})δk−k′δm+n′δn+m′+(a3P+2Re{a3′P})δk+nδm−m′δk′+n′]+Rs2Δf3[(b1P+2Re{b1′P})δk−m−k′δn+m′−n′+(b2P+2Re{b2′P})δk−m+m′δn−k′−n′+(b2∗P+2Re{b3′P})δk+n−k′δm−m′+n′+(b4P+2Re{b4′P})δk+n+m′δm+k′+n′+(c1P+2Re{c1′P})δk+nδm+k′−m′+n′+(c2P+2Re{c2′P})δk−k′δm−n−m′+n′+(c3P+2Re{c3′P})δk+m′δm−n+k′+n′+(c4P+2Re{c4′P})δm−m′δk+n−k′−n′+(c3∗P+2Re{c5′P})δm+k′δk+n+m′−n′+(c1∗P+2Re{c6′P})δk′+n′δk−m+n+m′]+RsΔf4(d1P+2Re{d1′P})δk−m+n−k′+m′−n′], 
where the coefficients multiplying P are listed in Table A5 in Appendix F and where coset leaders in Table 5 have been used.

entropy-22-01324-t004_Table 4Table 4 List of Kronecker delta contributions ordered by the corresponding high-order moment or correlation of the transmitted modulation format.

Correlation Terms	Kronecker Delta Products	

Intra-Polarisation Terms
	

In Mg(h)
	
	

E3{|ax|2}
	{δk−k′δm−m′δn−n′,δk−n′δm−m′δn−k′},{−2δn−k′δk−m+m′−n′,−2δn−n′δk−m−k′+m′,	

	−2δk−k′δm−n−m′+n′,−2δm−m′δk+n−k′−n′,−2δk−n′δm−n+k′−m′},{12δk−m+n−k′+m′−n′}	

E{|ax|2}|E{ax2}|2
	{δk+nδm−m′δk′+n′,δk−k′δm+n′δn+m′,δk+m′δm+k′δn−n′,δk+m′δm+n′δn−k′,	

	δk−n′δm+k′δn+m′},{−δn−k′δk−m+m′−n′,−3δm+k′δk+n+m′−n′,−3δk+nδm+k′−m′+n′	

	
−3δk′+n′δk+n−m+m′,−δm−m′δk+n−k′−n′,−3δm+n′δk+n−k′+m′,−3δn+m′δk−m−k′−n′,
	

	−δk−k′δm−n−m′+n′,−3δk+m′δm−n+k′+n′,−δn−n′δk−m−k′+m′,−δk−n′δm−n+k′−m′,}	

	{18δk−m+n−k′+m′−n′}	

|E{ax|ax|2}|2
	{}, {δk−m−k′δn+m′−n′,δk−m+m′δn−k′−n′,δk−m−n′δn−k′+m′,δk+n−k′δm−m′+n′,	

	δk+n−n′δm+k′−m′,δk−k′+m′δm−n+n′,δk−k′−n′δm−n+m′,δk+m′−n′δm−n+k′},	

	{−9δk−m+n−k′+m′−n′}	

|E{ax3}|2
	{}, {δk+n+m′δm+k′+n′},{−δk−m+n−k′+m′−n′}	

E{|ax|4}E{|ax|2}
	{}, {δk−k′δm−n−m′+n′,δk−n′δm−n+k′−m′,δm−m′δk+n−k′−n′,δn−k′δk−m+m′−n′,	

	
δn−n′δk−m−k′+m′},{−9δk−m+n−k′+m′−n′}
	

E∗{ax2|ax|2}E{ax2}
	{}, {δk+nδm+k′−m′+n′,δk+m′δm−n+k′+n′,δn+m′δk−m−k′−n′},	

	
{−3δk−m+n−k′+m′−n′}
	

E{ax2|ax|2}E∗{ax2}
	{}, {δm+k′δk−n+m′−n′,δm+n′δk+n−k′+m′,δk′+n′δk−m+n+m′},{−3δk−m+n−k′+m′−n′}	

E{|ax|6}
	{}, {}, {δk−m+n−k′+m′−n′}	

Cross-polarisation terms
	

In Mg(h)
	
	

E{|ax|2}E2{|ay|2}
	
{δk−k′δm−m′δn−n′},{−2δn−n′δk−m−k′+m′,−δm−m′δk+n−k′−n′,−δk−k′δm−n−m′+n′},
	

	
{4δk−m+n−k′+m′−n′}
	

E{|ax|2}|E{ay2}|2
	
{δk+m′δm+k′δn−n′},{−δn−n′δk−m−k′+m′,−δm+k′δk+n+m′−n′,−δk+m′δm−n+k′+n′},
	

	
{2δk−m+n−k′+m′−n′}
	

E{|ax|2}E{|ay|4}
	{}, {δn−n′δk−m−k′+m′},{−δk−m+n−k′+m′−n′}	

E{|ax|2ay∗}E{ay|ay|2}
	{}, {δk−m+m′δn−k′−n′,δk−k′+m′δm−n+n′},{−2δk−m+n−k′+m′−n′}	

E{|ax|2|ay|2}E{|ay|2}
	{}, {δk−k′δm−n−m′+n′,δm−m′δk+n−k′−n′},{−4δk−m+n−k′+m′−n′}	

E{axay}E∗{axay|ax|2}
	{}, {δk+m′δm−n+k′+n′,δm+k′δk−n+m′−n′,δn+m′δk−m−k′−n′,δk′+n′δk−m+n+m′},	

	
{−4δk−m+n−k′+m′−n′}
	

E∗{|ax|2ay2}E{ay2}
	{}, {δk+m′δm−n+k′+n′},{−δk−m+n−k′+m′−n′}	

E{|ax|2ay2}E∗{ay2}
	{}, {δm+k′δk−n+m′−n′},{−δk−m+n−k′+m′−n′}	

E{|ax|2|ay|4}
	{}, {}, {δk−m+n−k′+m′−n′}	

|E{axay∗}|2E{|ay|2}
	
{δk−n′δm−m′δn−k′},{−2δn−k′δk−m+m′−n′,−2δk−n′δm−n+k′−m′,−δk−k′δm−n−m′+n′,
	

	
−δm−m′δk+n−k′−n′},{8δk−m+n−k′+m′−n′}
	

E{axay}E{ax∗ay}E∗{ay2}
	
{δk−n′δm+k′δn+m′},{−2δm+k′δk+n+m′−n′,−δk+nδm+k′−m′+n′,−δn+m′δk−m−k′−n′,
	

	
−δk−n′δm−n+k′−m′},{4δk−m+n−k′+m′−n′}
	

E∗{axay}E{axay∗}E{ay2}
	
{δk+m′δm+n′δn−k′},{−2δk+m′δm−n+k′+n′,−δk′+n′δk−m+n+m′,−δn−k′δk−m+m′−n′,
	

	
−δm+n′δk+n−k′+m′},{4δk−m+n−k′+m′−n′}
	

|E{axay}|2E{|ay|2}
	{δk+nδm−m′δk′+n′,δk−k′δm+n′δn+m′},{−2δk′+n′δk−m+n+m′,−2δk+nδm−m′+k′+n′,	

	−2δn+m′δk−m−k′−n′,−2δm+n′δk+n−k′+m′,−δm−m′δk+n−k′−n′,−δk−k′δm−n−m′+n′},	

	
{8δk−m+n−k′+m′−n′}
	

|E{ax|ay|2}|2
	{}, {δk−m−n′δn−k′+m′,δk+n−k′δm−m′+n′,δk−k′−n′δm−n+m′},{−4δk−m+n−k′+m′−n′}	

|E{axay2}|2
	{}, {δk+n+m′δm+k′+n′},{−δk−m+n−k′+m′−n′}	

|E{ax∗ay2}|2
	{}, {δk+m′−n′δm−n+k′},{−δk−m+n−k′+m′−n′}	

In Ng(h)
	
	

E{|ax|2}|E{axay∗}|2
	{δk−k′δm−m′δn−n′,δk−n′δm−m′δn−k′},{−2δn−k′δk−m+m′−n′,−2δk−k′δm−n−m′+n′,	

	
−2δm−m′δk+n−k′−n′,−δn−n′δk−m−k′+m′,−δk−n′δm−n+k′−m′},{8δk−m+n−k′+m′−n′}
	

E{|ax|2}|E{axay}|2
	
{δk+m′δm+k′δn−n′,δk−n′δm+k′δn+m′},{−2δk′+n′δk−m+n+m′,−2δn+m′δk−m−k′−n′,
	

	−2δk+m′δm−n+k′+n′,−2δm+k′δk+n+m′−n′,−δn−n′δk−m−k′+m′,−δk−n′δm−n+k′−m′},	

	
{8δk−m+n−k′+m′−n′}.
	

E2{|ax|2}E{|ay|2}
	{}, {−δn−n′δk−m−k′+m′,−δk−n′δm−n+k′−m′},{4δk−m+n−k′+m′−n′}	

E{|ax|2}E{|ax|2|ay|2}
	{}, {δk−n′δm−n+k′−m′,δn−n′δk−m−k′+m′},{−4δk−m+n−k′+m′−n′}	

E{|ax|4}E{|ay|2}
	{}, {}, {−δk−m+n−k′+m′−n′}	

E{ax|ax|2}E{ax|ay|2}
	{}, {}, {−δk−m+n−k′+m′−n′}	

E{ay∗|ay|2}E{|ax|2ay}
	{}, {δk−m−k′δn+m′−n′,δk+n−n′δm+k′−m′},{−2δk−m+n−k′+m′−n′}	

E{ax∗|ax|2}E{ax|ay|2}
	{}, {δk−m−n′δn−k′+m′,δk+n−n′δm+k′−m′},{−2δk−m+n−k′+m′−n′}	

|E{|ax|2ay}|2
	{}, {δk−m−k′δn+m′−n′,δk−m+m′δn−k′−n′,δk−k′−n′δm−n+m′,δk+m′−n′δm−n+k′},	

	
{−4δk−m+n−k′+m′−n′}
	

E{axay∗}E{ax∗ay|ax|2}
	{}, {δk−k′δm−n−m′+n′,δm−m′δk+n−k′−n′,δn−k′δk−m+m′−n′},{−4δk−m+n−k′+m′−n′}	

E{|ax|4|ay|2}
	{}, {}, {δk−m+n−k′+m′−n′}	

E{ax2}E∗{axay}E{ax∗ay}
	
{δk+nδm−m′δk′+n′},{−2δk+nδm+k′−m′+n′,−δm+k′δk+n+m′−n′,−δk′+n′δk+n−m+m′,
	

	
−δm−m′δk+n−k′−n′},{4δk−m+n−k′+m′−n′}
	

E∗{ax2}E{axay}E{axay∗}
	
{δk−k′δm+n′δn+m′,δk+m′δm+n′δn−k′},{−2δm+n′δk+n−k′+m′,−δn+m′δk−m−k′−n′,
	

	
−δk−k′δm−n−m′+n′,−δn−k′δk−m+m′−n′,−δk+m′δm−n+k′+n′},{4δk−m+n−k′+m′−n′}
	

|E{ax2}|2E{|ay|2}
	{}, {−δm+n′δk+n−k′+m′,−δk+nδm+k′−m′+n′},{2δk−m+n−k′+m′−n′}	

|E{ax2ay∗}|2
	{}, {δk+n−k′δm−m′+n′},{−δk−m+n−k′+m′−n′}	

|E{ax2ay}|2
	{}, {δk+n+m′δm+k′+n′},{−δk−m+n−k′+m′−n′}	

E{ax2}E∗{ax2|ay|2}
	{}, {δk+nδm+k′−m′+n′,δm+n′δk+n−k′+m′},{−2δk−m+n−k′+m′−n′}	

E{axay}E∗{axay|ay|2}
	{}, {δk+nδm+k′−m′+n′,δn+m′δk−m−k′−n′},{−2δk−m+n−k′+m′−n′}	

E∗{axay}E{axay|ay|2}
	{}, {δm+n′δk+n−k′+m′,δk′+n′δk−m+n+m′},{−2δk−m+n−k′+m′−n′}	

E{ax∗ay}E{axay∗|ay|2}
	{}, {δk−n′δm−n+k′−m′},{−δk−m+n−k′+m′−n′}	

E{axay∗}E{ax∗ay|ay|2}
	{}, {δn−k′δk−m+m′−n′},{−δk−m+n−k′+m′−n′}	
entropy-22-01324-t005_Table 5Table 5 Subsets of Kronecker delta products which are closed to property (36). The terms in boldface are the ones used to group all other elements within each set.

Set Name	Set Elements	

D1
	
δk−k′δm−m′δn−n′,δk−n′δm−m′δn−k′
	

D2
	
δk−k′δm+n′δn+m′,δk+m′δm+k′δn−n′,δk+m′δm+n′δn−k′,δk−n′δm+k′δn+m′
	

D3
	
δk+nδm−m′δk′+n′
	

D4
	
δk−m−k′δn+m′−n′,δk−m−n′δn−k′+m′,δk+m′−n′δm−n+k′,δk−k′+m′δm−n+n′
	

D5
	
δk−m+m′δn−k′−n′,δk−k′−n′δm−n−m′
	

D6
	
δk+n−k′δm−m′+n′,δk+n−n′δm+k′−m′
	

D7
	
δk+n+m′δm+k′+n′
	

D8
	
δk+nδm+k′−m′+n′
	

D9
	
δk−k′δm−n−m′+n′,δk−n′δm−n+k′−m′,δn−k′δk−m+m′−n′,δn−n′δk−m−k′+m′
	

D10
	
δk+m′δm−n+k′+n′,δn+m′δk−m−k′−n′
	

D11
	
δm−m′δk+n−k′−n′
	

D12
	
δm+k′δk+n+m′−n′,δm+n′δk+n−k′+m′
	

D13
	
δk′+n′δk−m+n+m′
	

D14
	
δk−m+n−k′+m′−n′
	
Equation (38) can be further manipulated using the following proposition.

Proposition 3. Let D1(k,m,n,k′,m′,n′) and D2(k,m,n,k′,m′,n′) be two Kronecker delta products of the kind shown in the second column of Table 4. If D1(k,m,n,k′,m′,n′)=D2(k′,m′,n′,k,m,n), then
 (39) ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D1(k,m,n,k′,m′,n′)=∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D2(k,m,n,k′,m′,n′)∗. 


Proof.  See Appendix G. □

Corollary 1. Let D(k,m,n,k′,m′,n′) be a Kronecker delta product for which the following property holds
 D(k,m,n,k′,m′,n′)=D(k′,m′,n′,k,m,n), 
then ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D(k,m,n,k′,m′,n′)∈R.

Proof.  This corollary directly follows from Proposition 3 when D1=D2=D. □

Based on Corollary 1, we obtain from (38)
 Sx(f,Ns,Ls)=892γ2Δf∑i=−∞∞δ(f−iΔf)·∑(k,m,n)∈Si(k′,m′,n′)∈Si[Rs3Δf2[(a1+2Re{a1′})Pδk−k′δm−m′δn−n′+(a2+2Re{a2′})Pδk−k′δm+n′δn+m′+(a3+2Re{a3′})Pδk+nδm−m′δk′+n′]+Rs2Δf3[(b1+2Re{b1′})Pδk−m−k′δn+m′−n′+(b2P+2Re{b2′P})δk−m+m′δn−k′−n′+(b2∗P+2Re{b3′P})δk+n−k′δm−m′+n′+(b4+2Re{b4′})Pδk+n+m′δm+k′+n′+(c1P+2Re{c1′P})δk+nδm+k′−m′+n′+(c2+2Re{c2′})Pδk−k′δm−n−m′+n′+(c3P+2Re{c3′P})δk+m′δm−n+k′+n′+(c4+2Re{c4′})Pδm−m′δk+n−k′−n′+(c5P+2Re{c5′P})δm+k′δk+n+m′−n′+(c6P+2Re{c6′P})δk′+n′δk−m+n+m′]+RsΔf4(d1+2Re{d1′})Pδk−m+n−k′+m′−n′], 
where we have used the fact that the sets Di for i=1,2,3,4,7,9,11,14, are closed to the transformation in Proposition 3. Furthermore, we notice that the set pairs (D5,D6), (D8,D13), and (D10,D12) represent pairs of complementary sets under the transformation in Corollary 1; hence, their elements can be grouped. Consequently,
 (40) Sx(f,Ns,Ls)=892γ2Δf∑i=−∞∞δ(f−iΔf)·Rs3Δf2Φ1Q1+Φ2Q2+Φ3Q3+Rs2Δf3Ψ1Q4+2Re{Ψ2Q5+Ψ3Q5∗}+Ψ4Q6+2Re{Λ1Q7+Λ2Q7∗}+Λ3Q8+2Re{Λ4Q9+Λ5Q9∗}+Λ6Q10+RsΔf4Ξ1Q11, 
where
 (41) Ql≜∑(k,m,n)∈Si(k′,m′,n′)∈SiPD(l)=∑Tl,iPl=1,2,…,11, 
the coefficients Φi, i=1,2,3, Ψi, i=1,…,4, Λi, i=1,…,6, and Ξ1 in (40) are given in Table 6; the sets Si are defined in (19); and D(l) is the coset leaders highlighted in boldface in Table 5 and listed in Table 7 with their corresponding set D. Finally, the sets Tl,i are defined as
 Tl,i≜{(k,m,n,k′,m′,n′)∈{0,1,…,W−1}6:(k,m,n)∈Si,(k′,m′,n′)∈Si,D(l)=1}.  Note how, in the second equality of (41), we have accounted for the multiplication by D(l) by restricting the summation set to Tl,i.

8. Final Result
Equation (40) expresses the NLI PSD for a periodic signal of period T=1/Δf as a function of the statistical moments and cross-polarisation correlations of a generic 4D modulation format. To generalise this result to aperiodic signals, we take the same approach in [8,17], i.e., we let the period T go to infinity or, equivalently, Δf→0 (see Figure 3).

The limit of (40) for Δf→0 is a limit of a distribution (a Dirac’s delta comb) which is parametric in Δf. To rigorously evaluate such a limit, we use Lemma 1 and Theorem 2 presented in the following. In particular, Theorem 2 presents the key result of this work.

Lemma 1 (Dimensionality of the sets Tl,i).
The sets Tl,i, for l=1,2,3, for l=4,…,10, and for l=11 have dimensionalities 2–4, respectively, ∀i∈Z.

Proof.  See Appendix H. □

Theorem 2 (Limit of the distribution Sx(f,Ns,Ls)).
For a generic aperiodic transmitted signal and a fibre transmission system like the one in Figure 2 and under the following assumptions:
i.i.d. sequence of zero-mean input DP-4D symbols an for n∈Z (see Section 3.1)

rectangular (or quasi-rectangular) spectrum of the transmitted pulse p(t) (see Section 3.1)

first-order RP framework for the solution of the Manakov equation in (8)


the NLI PSD S¯x(f,Ns,Ls)≜limΔf→0Sx(f,Ns,Ls), where Sx(f,Ns,Ls) is given in (40), is
 (42) S¯x(f,Ns,Ls)=892γ2Rs3Φ1χ1(f)+Φ2χ2(f)+Φ3χ3(f)+Rs2Ψ1χ4(f)+2Re{Ψ2χ5(f)+Ψ3χ5∗(f)}+Ψ4χ6(f)+2Re{Λ1χ7(f)+Λ2χ7∗(f)}+Λ3χ8(f)+2Re{Λ4χ9(f)+Λ5χ9∗(f)}+Λ6χ10(f)+RsΞ1χ11(f), 
where the coefficients Φi, i=1,2,3, Ψi, i=1,2,…,4, Λi, i=1,2,…,6, and Ξ1 as well as the integrals χi(f), i=1,2,…,11, are given in Table 8. As discussed at the end of Section 4, S¯y(f) can be obtained applying the transformation x→y, y→x to (42).

The NLI power vector ΣNLI can be obtained from the PSDs in x and y as
 (43) ΣNLI≜[σNLI,x2,σNLI,y2]T=∫−∞∞S¯x(f,Ns,Ls)|P(f)|2df,∫−∞∞S¯y(f,Ns,Ls)|P(f)|2dfT, 
where P(f) is the transmitted pulse spectrum.

Proof.  See Appendix I. □

The results in (42) and (43) generalise the formulas presented in [9] for PM-2D modulation formats. In particular, assuming that
ax and ay are statistically independent, which leads to, e.g., E{axay}=E{ax}E{ay}=0,

ax and ay are identically distributed, which, for example, leads to E{|a|2}≜E{|ax|2}=E{|ay|2},

E{ax2}=E{ax3}=E{ay2}=E{ay3}=0, which applies, for instance, to distributions with a certain degree of symmetry,

 it can be seen that (42) reduces to
 S¯x(f,Ns,Ls)=892γ2[Rs3Φ1χ1(f)+Rs2(Λ3χ8(f)+Λ6χ10(f))+RsΞ1χ11(f)], 
with
 Φ1=3E3{|a|2},Λ3=5E{|a|4}E{|a|2}−10E3{|a|2},Λ6=E{|a|4}E{|a|2}−2E3{|a|2},Ξ1=E{|a|6}−9E{|a|4}E{|a|2}+12E3{|a|2}, 
which matches the formulation given in [9] (Equation (41)).

entropy-22-01324-t008_Table 8Table 8 Table of high-order moments, correlation coefficients, and integrals appearing in (42). The function η(f1,f2,f) is defined in (12).

Name	Value	

Correlation coefficients
	

Φ1
	
2E3{|ax|2}+4E{|ax|2}|E{axay∗}|2+E{|ax|2}E2{|ay|2}+|E{axay∗}|2E{|ay|2}
	

Φ2
	
4E{|ax|2}|E{ax2}|2+E{|ax|2}|E{ay2}|2+4E{|ax|2}|E{axay}|2+|E{axay}|2E{|ay|2}


+2Re{E{axay}E{ax∗ay}E∗{ay2}+2E∗{ax2}E{axay}E{axay∗}}
	

Φ3
	
E{|ax|2}|E{ax2}|2+|E{axay}|2E{|ay|2}+2Re{E{ax2}E∗{axay}E{ax∗ay}}
	

Ψ1
	
4|E{ax|ax|2}|2+4|E{|ax|2ay}|2+E{|ax|2ay}E{ay∗|ay|2}+E{|ax|2ay∗}E{ay|ay|2}+|E{ax|ay|2}|2


+|E{ax∗ay2}|2+2Re{E{ax∗|ax|2}E{ax|ay|2}}
	

Ψ2
	
2|E{ax|ax|2}|2+2|E{|ax|2ay}|2+E{|ax|2ay∗}E{ay|ay|2}+|E{ax|ay|2}|2
	

Ψ3
	
E{ax∗|ax|2}E{ax|ay|2}+|E{ax2ay∗}|2
	

Ψ4
	
|E{ax3}|2+2|E{ax2ay}|2+|E{axay2}|2
	

Λ1
	
−3E{|ax|2}|E{ax2}|2+E∗{ax2|ax|2}E{ax2}−|E{ax2}|2E{|ay|2}−2|E{axay}|2E{|ay|2}


+E{ax2}E∗{ax2|ay|2}−2E{ax2}E∗{axay}E{ax∗ay}+E{axay}E∗{axay|ay|2}−E{axay}E{ax∗ay}E∗{ay2}
	

Λ2
	
−2E{|ax|2}|E{axay}|2+E{axay}E∗{axay|ax|2}−E{ax2}E∗{axay}E{ax∗ay}
	

Λ3
	
4E{|ax|4}E{|ax|2}−4E{|ax|2}|E{ax2}|2−8E3{|ax|2}+4E{|ax|2}E{|ax|2|ay|2}


−12E{|ax|2}|E{axay∗}|2−4E{|ax|2}|E{axay}|2−4E2{|ax|2}E{|ay|2}−3E{|ax|2}E2{|ay|2}


−E{|ax|2}|E{ay2}|2+E{|ax|2|ay|2}E{|ay|2}+E{|ax|2}E{|ay|4}−5|E{axay∗}|2E{|ay|2}


−|E{axay}|2E{|ay|2}+2Re{2E{axay∗}E{ax∗ay|ax|2}−E{axay}E{ax∗ay}E∗{ay2}


+E{ax∗ay}E{axay∗|ay|2}−2E∗{ax2}E{axay}E{axay∗}}
	

Λ4
	
−6E{|ax|2}|E{ax2}|2+2E∗{ax2|ax|2}E{ax2}+4E{|ax|2}|E{axay}|2−E{|ax|2}|E{ay2}|2


+E∗{|ax|2ay2}E{ay2}+2E{axay}E∗{ax|ax|2ay}−2|E{axay}|2E{|ay|2}−2E∗{ax2}E{axay}E{axay∗}


+E{axay}E∗{axay|ay|2}−E∗{axay}E{axay∗}E{ay2}−2Re{E∗{axay}E{axay∗}E{ay2}}
	

Λ5
	
−2E{|ax|2}|E{axay}|2+E{axay}E∗{axay|ax|2}−|E{ax2}|2E{|ay|2}−E∗{ax2}E{axay}E{axay∗}


−2Re{E{ax2}E∗{axay}E{ax∗ay}}
	

Λ6
	
−2E3{|ax|2}+E{|ax|4}E{|ax|2}−E{|ax|2}|E{ax2}|2−4E{|ax|2}|E{axay∗}|2−E{|ax|2}E2{|ay|2}


+E{|ax|2|ay|2}E{|ay|2}−|E{axay∗}|2E{|ay|2}−|E{axay}|2E{|ay|2}


+2Re{E{axay∗}E{ax∗ay|ax|2}−E{ax2}E∗{axay}E{ax∗ay}}
	

Ξ1
	
E{|ax|6}−9E{|ax|4}E{|ax|2}+12E3{|ax|2}−2E{|ax|4}E{|ay|2}+E{|ax|2|ay|4}


−8E{|ax|2}E{|ax|2|ay|2}−4E{|ax|2|ay|2}E{|ay|2}+2E{|ax|4|ay|2}−E{|ax|2}E{|ay|4}


+4E{|ax|2}E2{|ay|2}+8E2{|ax|2}E{|ay|2}+18E{|ax|2}|E{ax2}|2−|E{ax3}|2−9|E{ax|ax|2}|2


+2E{|ax|2}|E{ay2}|2−4|E{ax|ay|2}|2−8|E{|ax|2ay}|2+8|E{axay∗}|2E{|ay|2}+8|E{axay}|2E{|ay|2}


−|E{axay2}|2−|E{ax∗ay2}|2+16E{|ax|2}|E{axay∗}|2−2|E{ax2ay∗}|2+16E{|ax|2}|E{axay}|2


+4|E{ax2}|2E{|ay|2}−2|E{ax2ay}|2+2Re{4E{axay}E{ax∗ay}E∗{ay2}−3E{ax2|ax|2}E∗{ax2}


−2E{|ax|2ay}E{ay∗|ay|2}−E{|ax|2ay2}E∗{ay2}−2E{axay}E∗{axay|ay|2}−E{axay∗}E{ax∗ay|ay|2}


−2E{ax∗|ax|2}E{ax|ay|2}−2E{ax2}E∗{ax2|ay|2}−E{ax|ax|2}E{ax|ay|2}−4E{axay∗}E{ax∗ay|ax|2}


−4E{axay}E∗{axay|ax|2}+8E{ax2}E∗{axay}E{ax∗ay}}
	

Integrals
	

χ1(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2|P(f1)|2|P(f2)|2|P(f−f1+f2)|2|η(f1,f2,f)|2df1df2
	

χ2(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2|P(f1)|2|P(f2)|2|P(f−f1+f2)|2η(f1,f2,f)η∗(f1,f1−f2−f,f)df1df2
	

χ3(f)
	
|P(f)|2∫−Rs/2Rs/2∫−Rs/2Rs/2|P(f1)|2|P(f2)|2η(f1,−f,f)η∗(f2,−f,f)df1df2
	

χ4(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2P(f1)P∗(f2)P(f−f1+f2)P∗(f1−f2)P(f3)P∗(f−f1+f2+f3)η(f1,f2,f)


·η∗(f1−f2,f3,f)df1df2df3
	

χ5(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2P(f1)P∗(f2)P(f−f1+f2)P∗(f3)P(f2−f1)P∗(f−f1+f2−f3)η(f1,f2,f)


·η∗(f3,f2−f1,f)df1df2df3
	

χ6(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2P(f1)P∗(f2)P(f−f1+f2)P∗(f3)P∗(f+f2)P(f2+f3)η(f1,f2,f)


·η∗(f3,−f−f2,f)df1df2df3
	

χ7(f)
	
P(f)∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2|P(f1)|2P∗(f2)P(f3)P∗(f−f2+f3)η(f1,−f,f)η∗(f2,f3,f)df1df2df3
	

χ8(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2|P(f1)|2P∗(f2)P(f−f1+f2)P(f3)P∗(f−f1+f3)η(f1,f2,f)η∗(f1,f3,f)df1df2df3
	

χ9(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2|P(f1)|2P∗(f2)P(f−f1+f2)P∗(f3)P(f−f1−f3)η(f1,f2,f)η∗(f3,−f1,f)df1df2df3
	

χ10(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2P(f1)|P(f2)|2P(f−f1+f2)P∗(f3)P∗(f+f2−f3)η(f1,f2,f)η∗(f3,f2,f)df1df2df3
	

χ11(f)
	
∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2∫−Rs/2Rs/2P(f1)P∗(f2)P(f−f1+f2)P∗(f3)P(f4)P∗(f−f3+f4)η(f1,f2,f)


·η∗(f3,f4,f)df1df2df3df4
	
9. Discussion and Conclusions
In this paper, we have derived a comprehensive analytical expression for NLI power when a general DP-4D modulation format is transmitted. The transmitted format is only assumed to be zero-mean. The reported result extends the model in [9] by accounting for any constellation geometry and statistic in four dimensions. This extension is performed by lifting two underlying assumptions in [9] (and other existing models): (i) the transmitted formats are PM versions of a 2D format and (ii) some high-order moments of the 2D components of the transmitted modulation format, such as E{ax2} and E{ax3}, are implicitly assumed to be equal to zero.

The presented results are derived in a single-channel transmission scenario. However, as it can be inferred from previous works, extending the expressions to the wavelength-division multiplexing (WDM) case does not lead to a different set of modulation-related statistical quantities in the NLI power expression. An extension of this work to the WDM transmission scenario will be addressed in a future publication.

Future work will also focus on comparing the presented model with possible heuristic extensions of existing PM-2D models to the general DP-4D case, for instance by using the 4D constellation standardised fourth-order moment (or so-called kurtosis). For such a study, a numerical validation of the model via the split-step Fourier method will certainly play a key role. Lastly, 4D constellation shaping in the optical fibre channel arguably represents the most attractive application and future research direction for the model derived in this manuscript. 

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author Contributions
Conceptualization, G.L.; investigation, G.L., A.B., and A.A.; methodology, G.L., A.B., and A.A.; writing—original draft preparation, G.L.; writing—review and editing, G.L., A.B., H.R., and A.A.; visualisation, G.L. and A.B. All authors have read and agreed to the published version of the manuscript.

Funding
The work of G.L. is funded by the EuroTechPostdoc programme under the European Union’s Horizon 2020 research and innovation programme (Marie Skłodowska-Curie grant agreement No. 754462). This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 757791).

Conflicts of Interest
The authors declare no conflict of interest.

Appendix A. Proof of Theorem 1
The Manakov Equation (8) can be written in the frequency domain as
 (A1) ∂E(f,z)∂z=−α2E(f,z)+j4π2f2β22E(f,z)+jγ89F{|E˜(t,z)|2E˜(t,z)}=−α2+j4π2f2β22E(f,z)+jγ89F{|E˜(t,z)|2}∗E(f,z), 
where ∗ denotes a modified convolution operator between a scalar function and a vector function (For a scalar function α and a vector function B=[Bx,By]T, the operator α∗B is defined here as α∗B≜α∗Bx,α∗ByT). Expanding the nonlinear term in (A1), we have
 F{|E˜(t,z)|2}∗E(f,z)=F{E˜x(t,z)E˜x∗(t,z)}+F{E˜y(t,z)E˜y∗(t,z)}∗[Ex(f,z),Ey(f,z)]T, 
which, for instance, for the x component, becomes
 (A2) Ex(f,z)∗Ex∗(−f,z)∗Ex(f,z)+Ey(f,z)∗Ey∗(−f,z)∗Ex(f,z).  Expanding the first term in (A2), we obtain
 (A3) Ex(f,z)∗Ex∗(−f,z)∗Ex(f,z)=∫−∞∞∫−∞∞Ex(f1,z)Ex∗(f1−f2,z)Ex(f−f2,z)df1df2, 
which by substitution f1−f2=f˜2 becomes (for notation’s simplicity, the integration variable f˜2 is relabelled as f2)
 (A4) Ex(f,z)∗Ex∗(−f,z)∗Ex(f,z)=−∫−∞∞∫−∞∞Ex(f1,z)Ex∗(f2,z)Ex(f−f1+f2,z)df1df2. 
Similar to the steps in (A3) and (A4), the second term in (A2) can be found as
 Ey(f,z)∗Ey∗(−f,z)∗Ex(f,z)=−∫−∞∞∫−∞∞Ey(f1,z)Ey∗(f2,z)Ex(f−f1+f2,z)df1df2. 
The x component in (A1) can be then rewritten as
 (A5) ∂Ex(f,z)∂z=−α2+j2π2f2β2Ex(f,z)−j89γ∫−∞∞∫−∞∞Ex(f1,z)Ex∗(f2,z)Ex(f−f1+f2,z)+Ey(f1,z)Ey∗(f2,z)Ex(f−f1+f2,z)df1df2. 

Following the first-order RP approach to finding the solution to the Manakov equation [13], we replace the x component of the first-order expansion in (9) into (A5) and equate terms with the same power of γ. After some algebra and after substituting the An terms with the corresponding En using (10), we find the following set of differential equations
 (A6) ∂E0,x(f,z)∂z=−α2+j2π2f2β2E0,x(f,z), 
 (A7) ∂E1,x(f,z)∂z=−j89γ∫−∞∞∫−∞∞E0,x(f1,z)E0,x∗(f2,z)E0,x(f−f1+f2,z)+E0,y(f1,z)E0,y∗(f2,z)E0,x(f−f1+f2,z)df1df2. 

The zeroth-order term for a single fibre span of length z is given by
 (A8) E0,x(f,z)=E(f,0)e(−α/2+j2π2β2f2)z. 
On the other hand, the first-order term (for the x component) E1,x(f,z), with initial conditions given by the transmitted signal E(f,0), can be found solving the following differential equation
 (A9) ∂E1,x(f,z)∂z=−α2+j2π2β2f2E1,x(f,z)−j89γ∫−∞∞∫−∞∞E0,x(f1,z)E0,x∗(f2,z)E0,x(f−f1+f2,z)+E0,y(f1,z)E0,y∗(f2,z)E0,y(f−f1+f2,z)df1df2. 
The solution to (A9) with initial condition E1,x(f,0)=0 is given by
 (A10a) E1,x(f,z)=−j89γe−αz+j2β2π2f2z∫0zeα2−j2β2π2f2z′∫−∞∞∫−∞∞E0,x(f1,z′)E0,x∗(f2,z′)E0,x(f−f1+f2,z′)+E0,y(f1,z′)E0,y∗(f2,z′)E0,x(f−f1+f2,z′)df1df2dz′ 
 (A10b) =−j89γe−αz+j2β2π2f2z∫0zeα2−j2β2π2f2z′∫−∞∞∫−∞∞Ex(f1,0)Ex∗(f2,0)Ex(f−f1+f2,0)+Ey(f1,0)Ey∗(f2,0)Ex(f−f1+f2,0)e−32αz′ej4π2β22(f12−f22+(f−f1+f2)2)z′df1df2dz′, 
where (A8) was used in the step from (A10a) to (A10b).

The power profile assumed in Section 3.1 for the multi-span optical link exponentially decays with a lumped amplification at the end of each span, which brings the power back to the transmitted level. This behavior leads to a discontinuity in the function α(z) across the interface where an amplifier is located. For such a power profile, we can solve the differential Equations (A6) and (A7) by exploiting the continuity of their coefficients within each span and by imposing the initial conditions at the input of each new fibre span E0,x(f,lLs+)=eαLs/2E0,x(f,lLs−) and E1,x(f,lLs+)=eαLs/2E1,x(f,lLs−), for l=1,2,…,Ns. Here, z=lLs− and z=lLs+ indicate the sections at the input and at the output of the lth amplifier, respectively. Thus, we obtain that the zeroth and first-order term after Ns fibre spans are given by
 (A11) E0,x(f,Ns,Ls)=E(f,0)ej2π2β2f2NsLs, 
 (A12) E1,x(f,Ns,Ls)=−j89γej2β2π2f2NsLs∑l=1Ns∫(l−1)LslLseα2−j2β2π2f2z′·∫−∞∞∫−∞∞E0,x(f1,z′)E0,x∗(f2,z′)E0,x(f−f1+f2,z′)+E0,y(f1,z′)E0,y∗(f2,z′)E0,x(f−f1+f2,z′)df1df2dz′. 
Using (A11) in (A12) and swapping the integral in z′ with the double integral in df1df2, we obtain
 E1,x(f,Ns,Ls)=−j89γej2β2π2f2NsLs∫−∞∞∫−∞∞Ex(f1,0)Ex∗(f2,0)Ex(f−f1+f2,0)+Ey(f1,0)Ey∗(f2,0)Ex(f−f1+f2,0)·∑l=1Ns∫(l−1)LslLse−α+jβ2(f−f1)(f2−f1)z′dz′df1df2=−j89γej2β2π2f2NsLs∫−∞∞∫−∞∞Ex(f1,0)Ex∗(f2,0)Ex(f−f1+f2,0)+Ey(f1,0)Ey∗(f2,0)Ex(f−f1+f2,0)·1−e−αzejβ2(f−f1)(f2−f1)zα−jβ2(f−f1)(f2−f1)∑l=1Nse−j4π2β2(l−1)(f−f1)(f2−f1)Lsdf1df2. 

The y-component of the zeroth-order term E0,y(f,z) and first-order term E1,y(f,z) can be found using the transformation x→y, y→x in (A11) and (A12), respectively. Finally, bringing together the x and y components, we have
 E0(f,Ns,Ls)=E0(f,0)ej2π2f2β2NsLs,E1(f,Ns,Ls)=−j89γej2β2π2f2NsLs∫−∞∞∫−∞∞ET(f1,0)E∗(f2,0)E(f−f1+f2,0)η(f1,f2,f,z)df1df2, 
where η(f,f1,f2,z) is defined in (12), which proves the theorem.

Appendix B. Proof of Proposition 1
Applying the variable transformation k˜=k′,m˜=m′,n˜=n′,k˜′=k,m˜′=m,andn˜′=n to the left-hand side of (23), we obtain
 (A13a) ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′νx,kνx,m∗νx,nνy,k′∗νy,m′νx,n′∗ηk,n,mηk′,n′,m′∗ 
 (A13b) =∑(k˜′,m˜′,n˜′)∈Si(k˜,m˜,n˜)∈SiPk˜′,m˜′,n˜′,k˜,m˜,n˜νx,k˜′νx,m˜′∗νx,n˜′νy,k˜∗νy,m˜νx,n˜∗ηk˜′,m˜′,n˜′ηk˜,m˜,n˜∗ 
 (A13c) =∑(k˜,m˜,n˜)∈Si(k˜′,m˜′,n˜′)∈SiPk˜,m˜,n˜,k˜′,m˜′,n˜′∗(νy,k˜∗νy,m˜νx,n˜∗νx,k˜′νx,m˜′∗νx,n˜′ηk˜,m˜,n˜ηk˜′,m˜′,n˜′∗)∗ 
 (A13d) =(∑(k˜,m˜,n˜)∈Si(k˜′,m˜′,n˜′)∈SiPk˜,m˜,n˜,k˜′,m˜′,n˜′νy,k˜νy,m˜∗νx,n˜νx,k˜′∗νx,m˜′νx,n˜′∗ηk˜,m˜,n˜∗ηk˜′,m˜′,n˜′)∗, 
where in the step between (A13b) and (A13c), we have used the property
 (A14) Pk˜′,m˜′,n˜′,k˜,m˜,n˜=Pk˜,m˜,n˜,k˜′,m˜′,n˜′∗, 
which can be easily verified based on definition (22). Using the relabelling k˜→k,m˜→m,n˜→n,k˜′→k′,m˜′→m′,n˜′→n′ for (A13d), the proposition is proven.

Appendix C. Partition Tables for Sets L1, L2, and L3
entropy-22-01324-t0A1_Table A1Table A1 List of all subsets in L1: for each subset, the index subgroups identify the corresponding pairs of indices assuming the same value.


	Index Subgroup	1	2	3	
Subset Label	
	

C1(1)
	
i1,i2
	
i3,i4
	
i5,i6
	

C1(2)
	
i1,i2
	
i3,i5
	
i4,i6
	

C1(3)
	
i1,i2
	
i3,i6
	
i4,i5
	

C1(4)
	
i1,i3
	
i2,i4
	
i5,i6
	

C1(5)
	
i1,i3
	
i2,i5
	
i4,i6
	

C1(6)
	
i1,i3
	
i2,i6
	
i4,i5
	

C1(7)
	
i1,i4
	
i2,i3
	
i5,i6
	

C1(8)
	
i1,i4
	
i2,i5
	
i3,i6
	

C1(9)
	
i1,i4
	
i2,i6
	
i3,i5
	

C1(10)
	
i1,i5
	
i2,i3
	
i4,i6
	

C1(11)
	
i1,i5
	
i2,i4
	
i3,i6
	

C1(12)
	
i1,i5
	
i2,i6
	
i3,i4
	

C1(13)
	
i1,i6
	
i2,i3
	
i4,i5
	

C1(14)
	
i1,i6
	
i2,i4
	
i3,i5
	

C1(15)
	
i1,i6
	
i2,i5
	
i3,i4
	
entropy-22-01324-t0A2_Table A2Table A2 List of all subsets in L2: for each subset, the index subgroups identify the corresponding triplets of indices assuming the same value.


	Subgroup Index	1	2	
Subset Label	
	

C2(1)
	
i1,i2,i3
	
i4,i5,i6
	

C2(2)
	
i1,i2,i4
	
i3,i5,i6
	

C2(3)
	
i1,i2,i5
	
i3,i4,i6
	

C2(4)
	
i1,i2,i6
	
i3,i4,i5
	

C2(5)
	
i1,i3,i4
	
i2,i5,i6
	

C2(6)
	
i1,i3,i5
	
i2,i4,i6
	

C2(7)
	
i1,i3,i6
	
i2,i4,i5
	

C2(8)
	
i1,i4,i5
	
i2,i3,i6
	

C2(9)
	
i1,i4,i6
	
i2,i3,i5
	

C2(10)
	
i1,i5,i6
	
i2,i3,i4
	
entropy-22-01324-t0A3_Table A3Table A3 List of all subsets in L3: for each subset, the index subgroups identify the corresponding pair and quadruple of indices assuming the same value. The set C3(1) corresponds to the case discussed in Example 2.


	Subgroup Index	1	2	
Subset Label	
	

C3(1)
	
i1,i2
	
i3,i4,i5,i6
	

C3(2)
	
i1,i3
	
i2,i4,i5,i6
	

C3(3)
	
i1,i4
	
i2,i3,i5,i6
	

C3(4)
	
i1,i5
	
i2,i3,i4,i6
	

C3(5)
	
i1,i6
	
i2,i3,i4,i5
	

C3(6)
	
i2,i3
	
i1,i4,i5,i6
	

C3(7)
	
i2,i4
	
i1,i3,i5,i6
	

C3(8)
	
i2,i5
	
i1,i3,i4,i6
	

C3(9)
	
i2,i6
	
i1,i3,i4,i5
	

C3(10)
	
i3,i4
	
i1,i2,i5,i6
	

C3(11)
	
i3,i5
	
i1,i2,i4,i6
	

C3(12)
	
i3,i6
	
i1,i2,i4,i5
	

C3(13)
	
i4,i5
	
i1,i2,i3,i6
	

C3(14)
	
i4,i6
	
i1,i2,i3,i5
	

C3(15)
	
i5,i6
	
i1,i2,i3,i4
	
Appendix D. Mg(h) and Ng(h) Contributions
entropy-22-01324-t0A4_Table A4Table A4 List of Mg(h) and Ng(h) contributions computed in Section 5 without bias terms.


g
	
h
	Corr. Terms in Mg(h)	Corr. Terms in Ng(h)	Delta Products	
1	1	
E3{|ax|2}+|E{axay∗}|2E{|ay|2}
	
E{|ax|2}|E{axay∗}|2
	
−Rs2Δfδn−k′δk−m+m′−n′


+2RsΔf2δk−m+n−k′+m′−n′
	
2	
E{|ax|2}|E{ax2}|2+|E{axay}|2E{|ay|2}
	
E{|ax|2}|E{axay}|2
	
−Rs2Δf(δk′+n′δk−m+n+m′


+δn+m′δk−m−k′−n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
3	
E3{|ax|2}+E{|ax|2}E2{|ay|2}
	
E2{|ax|2}E{|ay|2}
	
−Rs2Δfδn−n′δk−m−k′+m′


+2RsΔf2δk−m+n−k′+m′−n′
	
4	
E{|ax|2}|E{ax2}|2+E{axay}E{ax∗ay}E∗{ay2}
	
E{ax2}E∗{axay}E{ax∗ay}
	
−Rs2Δf(δm+k′δk+n+m′−n′


+δk+nδm+k′−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
5	
E{|ax|2}|E{ax2}|2+|E{axay}|2E{|ay|2}
	
E{ax2}E∗{axay}E{ax∗ay}
	
Rs3δk+nδm−m′δk′+n′


−Rs2Δf(δk′+n′δk−m+n+m′


+δm−m′δk+n−k′−n′


+δk+nδm−m′+k′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
6	
E{|ax|2}|E{ax2}|2+|E{axay}|2E{|ay|2}
	
|E{ax2}|2E{|ay|2}
	
−Rs2Δf(δm+n′δk+n−k′+m′


+δk+nδm+k′−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
7	
E3{|ax|2}+|E{axay∗}|2E{|ay|2}
	
E{|ax|2}|E{axay∗}|2
	
−Rs2Δfδk−k′δm−n−m′+n′


+2RsΔf2δk−m+n−k′+m′−n′
	
8	
E3{|ax|2}+E{|ax|2}E2{|ay|2}
	
E{|ax|2}|E{axay∗}|2
	
Rs3δk−k′δm−m′δn−n′


−Rs2Δf(δn−n′δk−m−k′+m′


+δm−m′δk+n−k′−n′


+δk−k′δm−n−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
9	
E{|ax|2}|E{ax2}|2+|E{axay}|2E{|ay|2}
	
E∗{ax2}E{axay}E{axay∗}
	
Rs3δk−k′δm+n′δn+m′


−Rs2Δf(δn+m′δk−m−k′−n′


+δm+n′δk+n−k′+m′


+δk−k′δm−n−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
10	
E{|ax|2}|E{ax2}|2+E∗{axay}E{axay∗}E{ay2}
	
E{|ax|2}|E{axay}|2
	
−Rs2Δf(δk′+n′δk−m+n+m′


+δk+m′δm−n+k′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
11	
E{|ax|2}|E{ax2}|2+E{|ax|2}|E{ay2}|2
	
E{|ax|2}|E{axay}|2
	
Rs3δk+m′δm+k′δn−n′


−Rs2Δf(δn−n′δk−m−k′+m′


+δm+k′δk+n+m′−n′


+δk+m′δm−n+k′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
12	
E{|ax|2}|E{ax2}|2+E∗{axay}E{axay∗}E{ay2}
	
E∗{ax2}E{axay}E{axay∗}
	
Rs3δk+m′δm+n′δn−k′


−Rs2Δf(δn−k′δk−m+m′−n′


+δm+n′δk+n−k′+m′


+δk+m′δm−n+k′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
13	
E3{|ax|2}+|E{axay∗}|2E{|ay|2}
	
E2{|ax|2}E{|ay|2}
	
−Rs2Δfδk−n′δm−n+k′−m′


+2RsΔf2δk−m+n−k′+m′−n′
	
14	
E{|ax|2}|E{ax2}|2+E{axay}E{ax∗ay}E∗{ay2}
	
E{|ax|2}|E{axay}|2
	
Rs3δk−n′δm+k′δn+m′


−Rs2Δf(δn+m′δk−m−k′−n′


+δm+k′δk+n+m′−n′


+δk−n′δm−n+k′−m′)


+2RsΔf2δk−m+n−k′+m′−n′
	
15	
E3{|ax|2}+|E{axay∗}|2E{|ay|2}
	
E{|ax|2}|E{axay∗}|2
	
Rs3δk−n′δm−m′δn−k′


−Rs2Δf(δn−k′δk−m+m′−n′


+δm−m′δk+n−k′−n′


+δk−n′δm−n+k′−m′)


+2RsΔf2δk−m+n−k′+m′−n′
	
2	1	
|E{ax|ax|2}|2+|E{ax|ay|2}|2
	
E{ax|ax|2}E{ax|ay|2}
	
−RsΔf2δk−m+n−k′+m′−n′
	
2	
|E{ax|ax|2}|2+E{ay∗|ay|2}E{|ax|2ay}
	
|E{|ax|2ay}|2
	
Rs2Δfδk−m−k′δn+m′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
3	
|E{ax|ax|2}|2+E{|ax|2ay∗}E{ay|ay|2}
	
|E{|ax|2ay}|2
	
Rs2Δfδk−m+m′δn−k′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
4	
|E{ax|ax|2}|2+|E{ax|ay|2}|2
	
E{ax∗|ax|2}E{ax|ay|2}
	
Rs2Δfδk−m−n′δn−k′+m′


−RsΔf2δk−m+n−k′+m′−n′
	
5	
|E{ax|ax|2}|2+|E{ax|ay|2}|2
	
|E{ax2ay∗}|2
	
Rs2Δfδk+n−k′δm−m′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
6	
|E{ax3}|2+|E{axay2}|2
	
|E{ax2ay}|2
	
Rs2Δfδk+n+m′δm+k′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
7	
|E{ax|ax|2}|2+E{|ax|2ay}E{ay∗|ay|2}
	
E{ax|ax|2}E{ax∗|ay|2}
	
Rs2Δfδk+n−n′δm+k′−m′


−RsΔf2δk−m+n−k′+m′−n′
	
8	
|E{ax|ax|2}|2+E{|ax|2ay∗}E{ay|ay|2}
	
E{ax|ay|2}E{ax∗|ax|2}
	
Rs2Δfδk−k′+m′δm−n+n′


−RsΔf2δk−m+n−k′+m′−n′
	
9	
|E{ax|ax|2}|2+|E{ax|ay|2}|2
	
|E{|ax|2ay}|2
	
Rs2Δfδk−k′−n′δm−n−m′


−RsΔf2δk−m+n−k′+m′−n′
	
10	
|E{ax|ax|2}|2+|E{ax∗ay2}|2
	
|E{|ax|2ay}|2
	
Rs2Δfδk+m′−n′δm−n+k′


−RsΔf2δk−m+n−k′+m′−n′
	
3	1	
E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}
	
E{|ax|2}E{|ax|2|ay|2}
	
−RsΔf2δk−m+n−k′+m′−n′
	
2	
E∗{ax2|ax|2}E{ax2}+E{axay}E∗{axay|ay|2}
	
E{ax2}E∗{ax2|ay|2}
	
Rs2Δfδk+nδm+k′−m′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
3	
E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}
	
E{axay∗}E{ax∗ay|ax|2}
	
Rs2Δfδk−k′δm−n−m′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
4	
E∗{ax2|ax|2}E{ax2}+E∗{|ax|2ay2}E{ay2}
	
E{axay}E∗{axay|ax|2}
	
Rs2Δfδk+m′δm−n+k′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
5	
E{|ax|4}E{|ax|2}+E{ax∗ay}E{axay∗|ay|2}
	
E{|ax|2}E{|ax|2|ay|2}
	
Rs2Δfδk−n′δm−n+k′−m′


−RsΔf2δk−m+n−k′+m′−n′
	
6	
E{|ax|4}E{|ax|2}+E{axay∗}E{ax∗ay|ay|2}
	
E{|ax|2}E{|ax|2|ay|2}
	
−RsΔf2δk−m+n−k′+m′−n′
	
7	
E{ax2|ax|2}E∗{ax2}+E{|ax|2ay2}E∗{ay2}
	
E∗{axay}E{axay|ax|2}
	
Rs2Δfδm+k′δk−n+m′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
8	
E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}
	
E{ax∗ay}E{axay∗|ax|2}
	
Rs2Δfδm−m′δk+n−k′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
9	
E{ax2|ax|2}E∗{ax2}+E∗{axay}E{axay|ay|2}
	
E∗{ax2}E{ax2|ay|2}
	
Rs2Δfδm+n′δk+n−k′+m′


−RsΔf2δk−m+n−k′+m′−n′
	
10	
E{|ax|4}E{|ax|2}+E{axay∗}E{ax∗ay|ay|2}
	
E{axay∗}E{ax∗ay|ax|2}
	
Rs2Δfδn−k′δk−m+m′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
11	
E∗{ax2|ax|2}E{ax2}+E{axay}E∗{axay|ay|2}
	
E{axay}E∗{axay|ax|2}
	
Rs2Δfδn+m′δk−m−k′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
12	
E{|ax|4}E{|ax|2}+E{|ax|2}E{|ay|4}
	
E{|ax|2}E{|ax|2|ay|2}
	
Rs2Δfδn−n′δk−m−k′+m′


−RsΔf2δk−m+n−k′+m′−n′
	
13	
E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}
	
E{|ax|4}E{|ay|2}
	
−RsΔf2δk−m+n−k′+m′−n′
	
14	
E{ax2|ax|2}E∗{ax2}+E∗{axay}E{axay|ay|2}
	
E∗{axay}E{axay|ax|2}
	
Rs2Δfδk′+n′δk−m+n+m′


−RsΔf2δk−m+n−k′+m′−n′
	
15	
E{|ax|4}E{|ax|2}+E{ax∗ay}E{axay∗|ay|2}
	
E{ax∗ay}E{axay∗|ax|2
	
−RsΔf2δk−m+n−k′+m′−n′
	
4	1	
E{|ax|6}+E{|ax|2|ay|4}
	
E{|ax|4|ay|2}
	
RsΔf2δk−m+n−k′+m′−n′
	
Appendix E. Proof of Proposition 2
Applying the variable transformation k˜=n, m˜=m, n˜=k, k˜′=n′, m˜′=m′, and n˜′=k′ to the right-hand side of (37), we have
 (A15) ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D2(k,m,n,k′,m′,n′)=∑(n˜,m˜,k˜)∈Si(n˜′,m˜′,k˜′)∈SiPn˜,m˜,k˜,n˜′,m˜′,k˜′ηn˜,m˜,k˜ηn˜′,m˜′,k˜′∗D2(n˜,m˜,k˜,n˜′,m˜′,k˜′). 

From definitions (16) and (22), it can be easily verified that Pn,m,k,n′,m′,k′=Pk,m,n,k′,m′,n′ and ηn,m,k=ηk,m,n. Moreover, based on the definition of the set Si in (19), it can be observed that the condition (k,m,n)∈Si is equivalent to (n,m,k)∈Si, i.e., generates the same set of triplets (k,m,n). We can thus write (A15) as
 ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D2(k,m,n,k′,m′,n′)=∑(k˜,m˜,n˜)∈Si(k˜′,m˜′,n˜′)∈SiPk˜,m˜,n˜,k˜′,m˜′,n˜′ηk˜,m˜,n˜ηk˜′,m˜′,n˜′∗D2(n˜,m˜,k˜,n˜′,m˜′,k˜′)=∑(k˜,m˜,n˜)∈Si(k˜′,m˜′,n˜′)∈SiPk˜,m˜,n˜,k˜′,m˜′,n˜′ηk˜,m˜,n˜ηk˜′,m˜′,n˜′∗D1(k˜,m˜,n˜,k˜′,m˜′,n˜′), 
which proves the proposition.

Appendix F. Coefficients in (38)
entropy-22-01324-t0A6_Table A5Table A5 Expressions for the coefficients in (38).

Name	Value	Name	Value	

a1
	
2E3{|ax|2}+E{|ax|2}E2{|ay|2}+|E{axay∗}|2E{|ay|2}
	
a1′
	
2E{|ax|2}|E{axay∗}|2
	

a2
	
4E{|ax|2}|E{ax2}|2+E{|ax|2}|E{ay2}|2


+|E{axay}|2E{|ay|2}+2Re{E{axay}E{ax∗ay}E∗{ay2}}
	
a2′
	
2E{|ax|2}|E{axay}|2


+2E∗{ax2}E{axay}E{axay∗}
	

a3
	
E{|ax|2}|E{ax2}|2+|E{axay}|2E{|ay|2}
	
a3′
	
E{ax2}E∗{axay}E{ax∗ay}
	

b1
	
4|E{ax|ax|2}|2+E{|ax|2ay}E{ay∗|ay|2}


+E{|ax|2ay∗}E{ay|ay|2}+|E{ax|ay|2}|2+|E{ax∗ay2}|2
	
b1′
	
E{ax∗|ax|2}E{ax|ay|2}


+2|E{|ax|2ay}|2
	

b2
	
2|E{ax|ax|2}|2+E{|ax|2ay∗}E{ay|ay|2}+|E{ax|ay|2}|2
	
b2′
	
2|E{|ax|2ay}|2
	

	
b3′
	
E{ax∗|ax|2}E{ax|ay|2}+|E{ax2ay∗}|2
	

b4
	
|E{ax3}|2+|E{axay2}|2
	
b4′
	
|E{ax2ay}|2
	

c1
	
−3E{|ax|2}|E{ax2}|2+E∗{ax2|ax|2}E{ax2}


−E{axay}E{ax∗ay}E∗{ay2}−2|E{axay}|2E{|ay|2}


+E{axay}E∗{axay|ay|2}
	
c1′
	
−2E{ax2}E∗{axay}E{ax∗ay}


−|E{ax2}|2E{|ay|2}


+E{ax2}E∗{ax2|ay|2}
	

c2
	
−8E3{|ax|2}−4E{|ax|2}|E{ax2}|2+4E{|ax|4}E{|ax|2}


−3E{|ax|2}E2{|ay|2}−E{|ax|2}|E{ay2}|2


+E{|ax|2}E{|ay|4}+E{|ax|2|ay|2}E{|ay|2}


−5|E{axay∗}|2E{|ay|2}−|E{axay}|2E{|ay|2}


+2Re{E{ax∗ay}E{axay∗|ay|2}−E{axay}E{ax∗ay}E∗{ay2}}


	
c2′
	
−6E{|ax|2}|E{axay∗}|2


−2E{|ax|2}|E{axay}|2


−2E2{|ax|2}E{|ay|2}


+2E{|ax|2}E{|ax|2|ay|2}


+2E{axay∗}E{ax∗ay|ax|2}


−2E∗{ax2}E{axay}E{axay∗}
	

c3
	
−6E{|ax|2}|E{ax2}|2+2E∗{ax2|ax|2}E{ax2}


−E{|ax|2}|E{ay2}|2+E∗{|ax|2ay2}E{ay2}


−E∗{axay}E{axay∗}E{ay2}−2|E{axay}|2E{|ay|2}


+E{axay}E∗{axay|ay|2}−2Re{E∗{axay}E{axay∗}E{ay2}}
	
c3′
	
4E{|ax|2}|E{axay}|2


+2E{axay}E∗{ax|ax|2ay}


−2E∗{ax2}E{axay}E{axay∗}
	

c4
	
−2E3{|ax|2}−E{|ax|2}|E{ax2}|2−E{|ax|2}E2{|ay|2}


+E{|ax|2|ay|2}E{|ay|2}−|E{axay∗}|2E{|ay|2}


−|E{axay}|2E{|ay|2}+E{|ax|4}E{|ax|2}
	
c4′
	
−2E{|ax|2}|E{axay∗}|2


+E{axay∗}E{ax∗ay|ax|2}


−E{ax2}E∗{axay}E{ax∗ay}
	

	
c5′
	
−2E{|ax|2}|E{axay}|2


+E{axay}E∗{axay|ax|2}


−E∗{ax2}E{axay}E{axay∗}


−|E{ax2}|2E{|ay|2}


−2Re{E{ax2}E∗{axay}E{ax∗ay}}
	

	
c6′
	
−2E{|ax|2}|E{axay}|2


+E{axay}E∗{axay|ax|2}


−E{ax2}E∗{axay}E{ax∗ay}
	

d1
	12E3{|ax|2}+18E{|ax|2}|E{ax2}|2−|E{ax3}|2

−9|E{ax|ax|2}|2−9E{|ax|4}E{|ax|2}

+E{|ax|6}+4E{|ax|2}E2{|ay|2}+2E{|ax|2}|E{ay2}|2

−E{|ax|2}E{|ay|4}−4E{|ax|2|ay|2}E{|ay|2}

−4|E{ax|ay|2}|2+8|E{axay∗}|2E{|ay|2}

+8|E{axay}|2E{|ay|2}−|E{axay2}|2−|E{ax∗ay2}|2

+E{|ax|2|ay|4}+2Re{4E{axay}E{ax∗ay}E∗{ay2}

−E{axay∗}E{ax∗ay|ay|2}−3E{ax2|ax|2}E∗{ax2}

−2E{axay}E∗{axay|ay|2}−2E{|ax|2ay}E{ay∗|ay|2}

−E{|ax|2ay2}E∗{ay2}}	
d1′
	
8E{|ax|2}|E{axay∗}|2−|E{ax2ay∗}|2


+8E{|ax|2}|E{axay}|2+4E2{|ax|2}E{|ay|2}


−4E{|ax|2}E{|ax|2|ay|2}−E{|ax|4}E{|ay|2}


−2E{ax∗|ax|2}E{ax|ay|2}


−2E{ax2}E∗{ax2|ay|2}−E{ax|ax|2}E{ax|ay|2}


−4|E{|ax|2ay}|2+E{|ax|4|ay|2}


−4E{axay∗}E{ax∗ay|ax|2}


−4E{axay}E∗{axay|ax|2}


+2|E{ax2}|2E{|ay|2}−|E{ax2ay}|2


+8Re{E{ax2}E∗{axay}E{ax∗ay}}
	
Appendix G. Proof of Proposition 3
Since D1(k,m,n,k′,m′,n′)=D2(k′,m′,n′,k,m,n), the left-hand side of (39) can be written as
 (A16) ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D1(k,m,n,k′,m′,n′)=∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D2(k′,m′,n′,k,m,n). 

Using the change of variables k˜=k′, m˜=m′, n˜=n′, k˜′=k, m˜′=m, and n˜′=n, the right-hand side of (A16) can be equivalently expressed as
 (A17a) ∑(k,m,n)∈Si(k′,m′,n′)∈SiPk,m,n,k′,m′,n′ηk,m,nηk′,m′,n′∗D2(k′,m′,n′,k,m,n) 
 (A17b) =∑(k˜′,m˜′,n˜′)∈Si(k˜,m˜,n˜)∈SiPk˜′,m˜′,n˜′,k˜,m˜,n˜ηk˜,m˜,n˜∗ηk˜′,m˜′,n˜′D2(k˜,m˜,n˜,k˜′,m˜′,n˜′) 
 (A17c) =∑(k˜,m˜,n˜)∈Si(k˜′,m˜′,n˜′)∈SiPk˜,m˜,n˜,k˜′,m˜′,n˜′∗ηk˜,m˜,n˜∗ηk˜′,m˜′,n˜′D2(k˜,m˜,n˜,k˜′,m˜′,n˜′) 
 (A17d) =(∑(k˜,m˜,n˜)∈Si(k˜′,m˜′,n˜′)∈SiPk˜,m˜,n˜,k˜′,m˜′,n˜′ηk˜,m˜,n˜ηk˜′,m˜′,n˜′∗D2(k˜,m˜,n˜,k˜′,m˜′,n˜′))∗, 
where in the step from (A17b) to (A17c), we have used (A14). Equation (A17d) is identical to the right-hand side of (39) up to the variable relabelling k˜→k,m˜→m,n˜→n,k˜′→k′,m˜′→m′,n˜′→n′, which proves the proposition.

Appendix H. Proof of Lemma 1
To prove the statement about dimensionality of the sets Tl,i, we take as an example the cases for l=1,2,3. In these instances, the sets Tl,i∀i∈Z are identified by 5 linear constraints on the set of variables (k,m,n,k′,m′,n′)∈{0,1,…,W−1}6 given by (i) the 2 linearly independent constraints, (k,m,n)∈Si and (k′,m′,n′)∈Si and (ii) the 3 linearly independent constraints induced by the condition D(l)=1 for l=1,2,3 (see Table 7). Then, let Al≜[a1T;a2T;…;a5T] be a 5×6 matrix in which the rows ak, k=1,2,…,5, describe each of these 5 linear combinations, x≜[k,m,n,k′,m′,n′], and yi=[i,i,0,0,0]. Thus, the set Tl,i can be equivalently defined as
 (A18) Tl,i={x∈{0,1,…,W−1}6:Alx=yi}. 

From (A18), it can be seen that Tl,i is a vector space for which the number of dimensions is given by
 (A19) dim{Tl,i}=6−rank(Al). 

Due to the construction of the delta products D(l),l∈{1,2,3}, it can be shown that the rows of Al are linearly dependent under the relationship a1−a2=±a3±a4±a5. Hence, ∀l∈{1,2,3} and i∈Z, we have rank(Al)=4. As a result, from (A19), dim{Tl,i}=2,∀l∈{1,2,3}, and i∈Z.

For l∈{4,5,…,10}, we have that Tl,i is identified by 4 linear constraints, 2 of them related to the Si set and 2 related to the condition D(l)=1. Furthermore, it can be seen that a1−a2=±a3±a4, hence leading to rank (Al)=3,∀l∈{4,5,…,10} and dim{Tl,i}=3. Finally, based on similar arguments, one can show that rank (Al)=2 for l=11 and dim{Tl,i}=4, which proves the lemma.

Appendix I. Proof of Theorem 2
The limit of a sequence of distributions f(Δf) is defined as the distribution f˜ such that [18] (Section 2.2)
 (A20) 〈f˜,ψ〉=limΔf→0〈f(Δf),ψ〉,∀ψ, 
where
 (A21) 〈f,ψ〉≜∫−∞∞f(f)ψ(f)df 
denotes the functional corresponding to the distribution f applied to a generic test function ψ [18] (Section 1.1). In particular, the delta distribution centered in f0 is defined as
 (A22) 〈δf0,ψ〉≜∫−∞∞δ(f−f0)ψ(f)df=ψ(f0). 
Based on (A21), we have for the distribution Sx(f,Ns,Ls) in (40),
 (A23a) 〈Sx(f,Ns,Ls),ψ〉=892γ2Δf[Rs3Δf2(Φ1∫−∞∞∑i=−∞∞∑T1,iPδ(f−iΔf)ψ(f)df+…+Φ3∫−∞∞∑i=−∞∞∑T3,iPδ(f−iΔf)ψ(f)df)+Rs2Δf3(Ψ1∫−∞∞∑i=−∞∞∑T4,iPδ(f−iΔf)ψ(f)df+…+Λ6∫−∞∞∑i=−∞∞∑T10,iPδ(f−iΔf)ψ(f)df)+RsΔf4Ξ1∫−∞∞∑i=−∞∞∑T11,iPδ(f−iΔf)ψ(f)df] 
 (A23b) =892γ2[Rs3Δf3(Φ1∑i=−∞∞∑T1,iPψ(iΔf)+…+Φ3∑i=−∞∞∑T3,iPψ(iΔf))+Rs2Δf4(Ψ1∑i=−∞∞∑T4,iPψ(iΔf)+…+Λ6∑i=−∞∞∑T10,iPψ(iΔf))+RsΔf5Ξ1∑i=−∞∞∑T11,iPψ(iΔf)], 
where we have used (A22) in the step between (A23a) and (A23b).

Now, we want to show that all the terms in (A23b) are multidimensional Riemann sums, which then will converge to multidimensional integrals in the limit for Δf→0. From (16), (22), and (35), it can be seen that the terms Pψ(iΔf) are samples on a multidimensional grid of step Δf of the multivariate function
 (A24) P˜(f1,f2,f3,f1′,f2′,f3′)≜P(f1)P∗(f2)P(f3)P∗(f1′)P(f2′)P∗(f3′)η(f1,f2,f3,Ns,Ls)η∗(f1′,f2′,f3′,Ns,Ls),f1,f2,f3,f1′,f2′,f3′∈R. 
Moreover, Δft(l), which represents the power of Δf multiplying the lth element in (A23b), where
 (A25) t(l)≜3forl=1,2,3;4forl=4,…,10;5forl=11,  is a measure of the t(l)th dimensional hypercube in Rt(l) for which the side measures Δf. Hence, to prove that each term in (A23b) converges to a sum of multiple integrals of the multivariate functions Pl˜ψ(f), we simply need to show that the dimensionality of the summation sets, i.e., Z×Tl,i, is equal to t(l), i.e., dim{Tl,i}=t(l)−1, for l=1,…,11, and ∀i∈Z. This can be easily verified comparing Lemma 1 to (A25). Defining the subspaces of R6
 (A26) Ql(f)≜{(f1,f2,f3,f1′,f2′,f3′)∈R6∩Gl:f1−f2+f3=f,f1′−f2′+f3′=f}, 
where Gl is the set defined by the condition D(l)=1, and by replacing the variables (k,m,n,k′,m′,n′)→(f1,f2,f3,f1′,f2′,f3′), we have
 (A27a) limΔf→0〈Sx(f,Ns,Ls),ψ〉=892γ2[Rs3(Φ1∫−∞∞∫...∫Q1(f)P˜(f1,…,f3′)ψ(f)df1…df3′df+…+Φ3∫−∞∞∫...∫Q3(f)P˜(f1,…,f3′)ψ(f)df1…df3′df)+Rs2(Ψ1∫−∞∞∫...∫Q4(f)P˜(f1,…,f3′)ψ(f)df1…df3′df+…+Λ6∫−∞∞∫...∫Q10(f)P˜(f1,…,f3′)ψ(f)df1…df3′df)+RsΞ1∫−∞∞∫...∫Q11(f)P˜(f1,…,f3′)ψ(f)df1…df3′df] 
 (A27b) =892γ2[Rs3(Φ1∫−∞∞∫−∞∞∫−∞∞P˜1(f1,f2,f)ψ(f)df1df2df+…+Φ3∫−∞∞∫−∞∞∫−∞∞P˜3(f1,f2,f)ψ(f)df1df2df)+Rs2(Ψ1∫...∫R4P˜4(f1,…,f3,f)ψ(f)df1…df3df+…+Λ6∫...∫R4P˜10(f1,…,f3,f)ψ(f)df1…df3df)+RsΞ1∫...∫R5P˜11(f1,…,f4,f)ψ(f)df1…df4df]. 
In the steps from (A27a) to (A27b), we have replaced in each integral the function P˜ with its constrained instance over Ql(f)
 P˜l≜P˜(f1,f2,f3,f1′,f2′,f3′)|(f1,f2,f3,f1′,f2′,f3′)∈Ql(f), 
and explicitly expressed the dimensionality of the integrals based on the dimension of their corresponding integration domains Ql(f). By construction (see (A26)), dim{Ql(f)}=t(l)−1, for l=1,…,11.

Finally, using (A20) and comparing definition (A21) with (A27b), we obtain
 S¯x(f,Ns,Ls)=limΔf→0Sx(f,Ns,Ls)=892γ2[Rs3(Φ1∫−∞∞∫−∞∞P˜1(f1,f2,f)df1df2+…+Φ3∫−∞∞∫−∞∞P˜3(f1,f2,f)df1df2)+Rs2(Ψ1∫−∞∞∫−∞∞∫−∞∞P˜4(f1,f2,f3,f)df1df2df3+…+Λ6∫−∞∞∫−∞∞∫−∞∞P˜10(f1,f2,f3,f)df1df2df3)+RsΞ1∫...∫R4P˜11(f1,…,f4,f)df1…df4], 
which, defining
 (A28) χl(f)≜∫R2P˜l(f1,f2,f)df1df2=∫−Rs2Rs2∫−Rs2Rs2P˜l(f1,f2,f)df1df2,l=1,2,3;∫R3P˜l(f1,f2,f3,f)df1df2df3=∫−Rs2Rs2∫−Rs2Rs2∫−Rs2Rs2P˜l(f1,f2,f3,f)df1df2df3,l=4,…,10;∫R4P˜11(f1,…,f4,f)df1…df4=∫...∫R4P˜11(f1,…,f4,f)df1…df4,l=11, 
with R4≜[−Rs/2,Rs/2]4, proves the theorem. The second equalities in (A28) are justified by the form of the functions P˜l (see Table 8) which, due to the assumption of strictly bandlimited pulses, have limited support within the hybercube [−Rs/2,Rs/2]t(l)−1. To derive the explicit expressions for χl(f) in Table 8, we used (A28) and the property P(−f)=P∗(f), which stems from the fact that p(t) is assumed to be a real-valued function (see Section 3.1).

Figure 4 Venn diagram of the partition on the 6D space i=(i1,i2,i3,i4,i5,i6)∈{0,1,…,W−1}6 discussed in Section 5.2

entropy-22-01324-t001_Table 1Table 1 List of contributions M1(h) and N1(h) for i=1,2,…,15.


h
	Corr. Terms in M1(h)	Corr. Terms in N1(h)	Delta Products	
1	
E3{|ax|2}


+|E{axay∗}|2E{|ay|2}
	
E{|ax|2}|E{axay∗}|2
	
Rs3δk−mδn−k′δm′−n′−Rs2Δf(δm′−n′δk−m+n−k′


+δn−k′δk−m+m′−n′+δk−mδm′−n′+n−k′)


+2RsΔf2δk−m+n−k′+m′−n′
	
2	
E{|ax|2}|E{ax2}|2


+|E{axay}|2E{|ay|2}
	
E{|ax|2}|E{axay}|2
	
Rs3δk−mδn+m′δk′+n′−Rs2Δf(δk′+n′δk−m+n+m′


+δn+m′δk−m−k′−n′+δk−mδn−k′+m′−n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
3	
E3{|ax|2}


+E{|ax|2}E2{|ay|2}
	
E2{|ax|2}E{|ay|2}
	
Rs3δk−mδn−n′δk′−m′−Rs2Δf(δk′−m′δk−m+n−n′


+δn−n′δk−m−k′+m′+δk−mδn−n′−k′+m′)


+2RsΔf2δk−m+n−k′+m′−n′
	
4	
E{|ax|2}|E{ax2}|2


+E{axay}E{ax∗ay}E∗{ay2}
	
E{ax2}E∗{axay}E{ax∗ay}
	
Rs3δk+nδm+k′δm′−n′−Rs2Δf(δm′−n′δk+n−m−k′


+δm+k′δk+n+m′−n′+δk+nδm+k′−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
5	
E{|ax|2}|E{ax2}|2


·|E{axay}|2E{|ay|2}
	
E{ax2}E∗{axay}E{ax∗ay}
	
Rs3δk+nδm−m′δk′+n′−Rs2Δf(δk′+n′δk+n−m+m′


+δm−m′δk+n−k′−n′+δk+nδm−m′+k′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
6	
E{|ax|2}|E{ax2}|2


+|E{axay}|2E{|ay|2}
	
|E{ax2}|2E{|ay|2}
	
Rs3δk+nδm+n′δk′−m′−Rs2Δf(δk′−m′δk+n−m−n′


+δm+n′δk+n−k′+m′+δk+nδm+k′−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
7	
E3{|ax|2}


+|E{axay∗}|2E{|ay|2}
	
E{|ax|2}|E{axay∗}|2
	
Rs3δk−k′δm−nδm′−n′−Rs2Δf(δm′−n′δk−k′−m+n


+δm−nδk−k′+m′−n′+δk−k′δm−n−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
8	
E3{|ax|2}


+E{|ax|2}E2{|ay|2}
	
E{|ax|2}|E{axay∗}|2
	
Rs3δk−k′δm−m′δn−n′−Rs2Δf(δn−n′δk−m−k′+m′


+δm−m′δk+n−k′−n′+δk−k′δm−n−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
9	
E{|ax|2}|E{ax2}|2


+|E{axay}|2E{|ay|2}
	
E∗{ax2}E{axay}E{axay∗}
	
Rs3δk−k′δm+n′δn+m′−Rs2Δf(δn+m′δk−m−k′−n′


+δm+n′δk+n−k′+m′+δk−k′δm−n−m′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
10	
E{|ax|2}|E{ax2}|2


+E∗{axay}E{axay∗}E{ay2}
	
E{|ax|2}|E{axay}|2
	
Rs3δk+m′δm−nδk′+n′−Rs2Δf(δk′+n′δk−m+n+m′


+δm−nδk−k′+m′−n′+δk+m′δm−n+k′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
11	
E{|ax|2}|E{ax2}|2


+E{|ax|2}|E{ay2}|2
	
E{|ax|2}|E{axay}|2
	
Rs3δk+m′δm+k′δn−n′−Rs2Δf(δn−n′δk−m−k′+m′


+δm+k′δk+n+m′−n′+δk+m′δm−n+k′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
12	
E{|ax|2}|E{ax2}|2


+E∗{axay}E{axay∗}E{ay2}
	
E∗{ax2}E{axay}E{axay∗}
	
Rs3δk+m′δm+n′δn−k′−Rs2Δf(δn−k′δk−m+m′−n′


+δm+n′δk+n−k′+m′+δk+m′δm−n+k′+n′)


+2RsΔf2δk−m+n−k′+m′−n′
	
13	
E3{|ax|2}


+|E{axay∗}|2E{|ay|2}
	
E2{|ax|2}E{|ay|2}
	
Rs3δk−n′δm−nδk′−m′−Rs2Δf(δk′−m′δk−m+n−n′


+δm−nδk−k′+m′−n′+δk−n′δm−n+k′−m′)


+2RsΔf2δk−m+n−k′+m′−n′
	
14	
E{|ax|2}|E{ax2}|2


+E{axay}E{xa∗ay}E∗{ay2}
	
E{|ax|2}|E{axay}|2
	
Rs3δk−n′δm+k′δn+m′−Rs2Δf(δn+m′δk−m−k′−n′


+δm+k′δk+n+m′−n′+δk−n′δm−n+k′−m′)


+2RsΔf2δk−m+n−k′+m′−n′
	
15	
E3{|ax|2}


+|E{axay∗}|2E{|ay|2}
	
E{|ax|2}|E{axay∗}|2
	
Rs3δk−n′δm−m′δn−k′−Rs2Δf(δn−k′δk−m+m′−n′


+δm−m′δk+n−k′−n′+δk−n′δm−n+k′−m′)


+2RsΔf2δk−m+n−k′+m′−n′
	
entropy-22-01324-t002_Table 2Table 2 List of contributions M2(h) and N2(h) for i=1,2,…,10.


h
	Corr. Terms in M2(h)	Corr. Terms in N2(h)	Delta Products	
1	
|E{ax|ax|2}|2+|E{ax|ay|2}|2
	
E{ax|ax|2}E{ax|ay|2}
	
Rs2Δfδk−m+nδk′−m′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
2	
|E{ax|ax|2}|2+E{ay∗|ay|2}E{|ax|2ay}
	
|E{|ax|2ay}|2
	
Rs2Δfδk−m−k′δn+m′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
3	
|E{ax|ax|2}|2+E{|ax|2ay∗}E{ay|ay|2}
	
|E{|ax|2ay}|2
	
Rs2Δfδk−m+m′δn−k′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
4	
|E{ax|ax|2}|2+|E{ax|ay|2}|2
	
E{ax∗|ax|2}E{ax|ay|2}
	
Rs2Δfδk−m−n′δn−k′+m′


−RsΔf2δk−m+n−k′+m′−n′
	
5	
|E{ax|ax|2}|2+|E{ax|ay|2}|2
	
|E{ax2ay∗}|2
	
Rs2Δfδk+n−k′δm−m′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
6	
|E{ax3}|2+|E{axay2}|2
	
|E{ax2ay}|2
	
Rs2Δfδk+n+m′δm+k′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
7	
|E{ax|ax|2}|2+E{|ax|2ay}E{ay∗|ay|2}
	
E{ax|ax|2}E{ax∗|ay|2}
	
Rs2Δfδk+n−n′δm+k′−m′


−RsΔf2δk−m+n−k′+m′−n′
	
8	
|E{ax|ax|2}|2+E{|ax|2ay∗}E{ay|ay|2}
	
E{ax|ay|2}E{ax∗|ax|2}
	
Rs2Δfδk−k′+m′δm−n+n′


−RsΔf2δk−m+n−k′+m′−n′
	
9	
|E{ax|ax|2}|2+|E{ax|ay|2}|2
	
|E{|ax|2ay}|2
	
Rs2Δfδk−k′−n′δm−n−m′


−RsΔf2δk−m+n−k′+m′−n′
	
10	
|E{ax|ax|2}|2+|E{ax∗ay2}|2
	
|E{|ax|2ay}|2
	
Rs2Δfδk+m′−n′δm−n+k′


−RsΔf2δk−m+n−k′+m′−n′
	
entropy-22-01324-t003_Table 3Table 3 List of contributions M3(h) and N3(h) for i=1,2,…,15.


h
	Corr. Terms in M3(h)	Corr. Terms in N3(h)	Delta Products	
1	
E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}
	
E{|ax|2}E{|ax|2|ay|2}
	
Rs2Δfδk−mδn−k′+m′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
2	
E∗{ax2|ax|2}E{ax2}+E{axay}E∗{axay|ay|2}
	
E{ax2}E∗{ax2|ay|2}
	
Rs2Δfδk+nδm+k′−m′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
3	
E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}
	
E{axay∗}E{ax∗ay|ax|2}
	
Rs2Δfδk−k′δm−n−m′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
4	
E∗{ax2|ax|2}E{ax2}+E∗{|ax|2ay2}E{ay2}
	
E{axay}E∗{axay|ax|2}
	
Rs2Δfδk+m′δm−n+k′+n′


−RsΔf2δk−m+n−k′+m′−n′
	
5	
E{|ax|4}E{|ax|2}+E{ax∗ay}E{axay∗|ay|2}
	
E{|ax|2}E{|ax|2|ay|2}
	
Rs2Δfδk−n′δm−n+k′−m′


−RsΔf2δk−m+n−k′+m′−n′
	
6	
E{|ax|4}E{|ax|2}+E{axay∗}E{ax∗ay|ay|2}
	
E{|ax|2}E{|ax|2|ay|2}
	
Rs2Δfδm−nδk−k′+m′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
7	
E{ax2|ax|2}E∗{ax2}+E{|ax|2ay2}E∗{ay2}
	
E∗{axay}E{axay|ax|2}
	
Rs2Δfδm+k′δk−n+m′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
8	
E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}
	
E{ax∗ay}E{axay∗|ax|2}
	
Rs2Δfδm−m′δk+n−k′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
9	
E{ax2|ax|2}E∗{ax2}+E∗{axay}E{axay|ay|2}
	
E∗{ax2}E{ax2|ay|2}
	
Rs2Δfδm+n′δk+n−k′+m′


−RsΔf2δk−m+n−k′+m′−n′
	
10	
E{|ax|4}E{|ax|2}+E{axay∗}E{ax∗ay|ay|2}
	
E{axay∗}E{ax∗ay|ax|2}
	
Rs2Δfδn−k′δk−m+m′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
11	
E∗{ax2|ax|2}E{ax2}+E{axay}E∗{axay|ay|2}
	
E{axay}E∗{axay|ax|2}
	
Rs2Δfδn+m′δk−m−k′−n′


−RsΔf2δk−m+n−k′+m′−n′
	
12	
E{|ax|4}E{|ax|2}+E{|ax|2}E{|ay|4}
	
E{|ax|2}E{|ax|2|ay|2}
	
Rs2Δfδn−n′δk−m−k′+m′


−RsΔf2δk−m+n−k′+m′−n′
	
13	
E{|ax|4}E{|ax|2}+E{|ax|2|ay|2}E{|ay|2}
	
E{|ax|4}E{|ay|2}
	
Rs2Δfδk′−m′δk−m+n−n′


−RsΔf2δk−m+n−k′+m′−n′
	
14	
E{ax2|ax|2}E∗{ax2}+E∗{axay}E{axay|ay|2}
	
E∗{axay}E{axay|ax|2}
	
Rs2Δfδk′+n′δk−m+n+m′


−RsΔf2δk−m+n−k′+m′−n′
	
15	
E{|ax|4}E{|ax|2}+E{ax∗ay}E{axay∗|ay|2}
	
E{ax∗ay}E{axay∗|ax|2
	
Rs2Δfδm′−n′δk−m+n−k′


−RsΔf2δk−m+n−k′+m′−n′
	
entropy-22-01324-t006_Table 6Table 6 Correlation coefficients in (40): the values of a1,a1′,b1,b1′,…d1′ are given in Table A5.

Name	Value	Name	Value	

Φ1
	
a1+2Re{a1′}
	
Λ1
	
c1+c1′
	

Φ2
	
a2+2Re{a2′}
	
Λ2
	
c6′
	

Φ3
	
a3+2Re{a3′}
	
Λ3
	
c2+2Re{c2′}
	

Ψ1
	
b1+2Re{b1′}
	
Λ4
	
c3+c3′
	

Ψ2
	
b2+b2′
	
Λ5
	
c5′
	

Ψ3
	
b3′
	
Λ6
	
c4+2Re{c4′}
	

Ψ4
	
b4+2Re{b4′}
	
Ξ1
	
d1+2Re{d1′}
	
entropy-22-01324-t007_Table 7Table 7 Delta products D(l) in the Ql terms, l=1,2,…,11, in (41) with their corresponding D set in Table 5.


l
	
D(l)
	Set D	
1	
δk−k′δm−m′δn−n′
	
D1
	
2	
δk−k′δm+n′δn+m′
	
D2
	
3	
δk+nδm−m′δk′+n′
	
D3
	
4	
δk−m−k′δn+m′−n′
	
D4
	
5	
δk−m+m′δn−k′−n′
	
D5
	
6	
δk+n+m′δm+k′+n′
	
D7
	
7	
δk+nδm+k′−m′+n′
	
D8
	
8	
δk−k′δm−n−m′+n′
	
D9
	
9	
δk+m′δm−n+k′+n′
	
D10
	
10	
δm−m′δk+n−k′−n′
	
D11
	
11	
δk−m+n−k′+m′−n′
	
D14
==== Refs
References
1. Agrell E.  Karlsson M.   Power-efficient modulation formats in coherent transmission systems J. Lightwave Technol. 2009 27 5115 5126 10.1109/JLT.2009.2029064 
2. Karlsson M.  Agrell E.   Which is the most power-efficient modulation format in optical links? Opt. Express 2009 17 10814 10819 10.1364/OE.17.010814 19550481 
3. Alvarado A.  Agrell E.   Four-dimensional coded modulation with bit-wise decoders for future optical communications J. Lightwave Technol. 2015 33 1993 2003 10.1109/JLT.2015.2396118 
4. Eriksson T.A.  Fehenberger T.  Andrekson P.A.  Karlsson M.  Hanik N.  Agrell E.   Impact of 4D channel distribution on the achievable rates in coherent optical communication experiments J. Lightwave Technol. 2016 34 2256 2266 10.1109/JLT.2016.2528550 
5. Kojima K.  Yoshida T.  Koike-Akino T.  Millar D.S.  Parsons K.  Pajovic M.  Arlunno V.   Nonlinearity-tolerant four-dimensional 2A8PSK family for 5–7 bits/symbol spectral efficiency J. Lightwave Technol. 2017 35 1383 1391 10.1109/JLT.2017.2662942 
6. Chen B.  Chigo O.  Hafermann H.  Alvarado A.   Polarization-ring-switching for nonlinearity-tolerant geometrically-shaped four-dimensional formats maximizing generalized mutual information J. Lightwave Technol. 2019 37 3579 3591 10.1109/JLT.2019.2918072 
7. Chen B.  Alvarado A.  van der Heide S.  van den Hout M.  Hafermann H.  Okonkwo C.   Analysis and experimental demonstration of orthant-symmetric four-dimensional 7 bit/4D-sym modulation for optical fiber communication arXiv 2020 2003.12712v2 
8. Poggiolini P.  Bosco G.  Carena A.  Curri V.  Jiang Y.  Forghieri F.   A detailed analytical derivation of the GN model of non-linear interference in coherent optical transmission systems arXiv 2012 1209.0394v13 
9. Carena A.  Bosco G.  Curri V.  Jiang Y.  Poggiolini P.  Forghieri F.   On the accuracy of the GN-model and on analytical correction terms to improve it arXiv 2014 1401.6946v7 
10. Mecozzi A.  Essiambre R.   Nonlinear Shannon limit in pseudolinear coherent systems J. Lightwave Technol. 2012 30 2011 2024 10.1109/JLT.2012.2190582 
11. Dar R.  Feder M.  Mecozzi A.  Shtaif M.   Properties of nonlinear noise in long, dispersion-uncompensated fiber links Opt. Express 2013 21 25685 25699 10.1364/OE.21.025685 24216794 
12. Marcuse D.  Menyuk C.R.  Wai P.K.A.   Application of the Manakov-PMD equation to studies of signal propagation in optical fibers with randomly varying birefringence J. Lightwave Technol. 1997 15 1735 1745 10.1109/50.622902 
13. Vannucci A.  Serena P.  Member S.  Bononi A.   The RP method: A new tool for the iterative solution of the nonlinear Schrödinger equation J. Lightwave Technol. 2002 20 1102 1112 10.1109/JLT.2002.800376 
14. Johannisson P.  Karlsson M.   Perturbation analysis of nonlinear propagation in a strongly dispersive optical communication system J. Lightwave Technol. 2013 31 1273 1282 10.1109/JLT.2013.2246543 
15. Colombeau J.F.   New Generalized Functions and Multiplication of Distributions North-Holland, Elsevier Science Publishers B.V. Amsterdam, The Netherlands 1984 
16. Proakis J.G.  Manolakis D.G.   Digital Signal Processing: Principles, Algorithms, and Applications 4th ed. Prentice-Hall, Inc. Upper Saddle River, NJ, USA 2006 
17. Carena A.  Curri V.  Bosco G.  Poggiolini P.  Forghieri F.   Modeling of the impact of nonlinear propagation effects in uncompensated optical coherent transmission links J. Lightwave Technol. 2012 30 1524 1539 10.1109/JLT.2012.2189198 
18. Strichartz R.   A Guide to Distribution Theory and Fourier Transforms CRC-Press Boca Raton, FL, USA 1994