Inter-core crosstalk in weakly coupled MCFs with arbitrary core layout and the effect of bending and twisting on the coupling coefficient

We investigate the bending and twisting-induced longitudinal variation of the inter-core coupling coefficient (ICCC) and its effect on inter-core crosstalk (ICXT) in weakly coupled multi-core fibers (MCFs) with an arbitrary core layout. An analytical discrete changes model (DCM) for ICXT field propagation under those conditions is proposed for the first time, providing very fast and rather accurate mean ICXT power estimates. The analytical mean ICXT power estimates are validated through numerical simulation. It is predicted that the mean ICXT power between adjacent cores of the outer ring of the 19-core MCF can be more than 10 dB higher than the one between adjacent cores of the inner ring. It is also predicted that the difference between the mean ICXT power of cores in the inner and outer rings can be much smaller by decreasing the core pitch and increasing the bending radius. This behavior is attributed to the ICXT dependence on the bending and twisting-induced longitudinal variation of the ICCCs. In particular, larger bending and twisting-induced fluctuations of the ICCCs along the longitudinal coordinate are observed in the cores of the outer ring, but the fluctuations become smaller for smaller core pitches and larger bending radii. Furthermore, it is shown that, if the ICCCs’ longitudinal variation is neglected, the mean ICXT power estimates between two adjacent cores are very similar despite the location of those cores. This means that neglecting the longitudinal variation of the ICCCs can lead to misleading estimates of the mean ICXT power, with an error exceeding 15 dB. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Transmission capacity demands are pushing channel capacity closer to Shannon's limit [1]. Mode and spatial division multiplexing techniques have been proposed to address such demands using, for instance, multi-core fibers (MCFs) [2,3]. In particular, single-mode homogeneous MCFs, characterized by cores with similar physical properties, have been used for high transmission capacity [4,5], and networking [6]. However, the inter-core crosstalk (ICXT) can compromise the quality of service, limiting the transmission distance and the use of higher order modulation formats [4].
In order to optimize the design of high transmission capacity in weakly coupled MCF-based systems, network planning and routing, rigorous ICXT models that provide fast ICXT power estimates are required. Models for ICXT estimation and MCF design optimization have been proposed in recent years [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. The discrete changes models (DCMs) proposed in [9][10][11][12][13][14][15][16] for weakly coupled MCF transmission rely on the fact that, when the MCF operates in the phase-matching region, the ICXT can be accurately modeled from the contribution of the so-called phase-matching points (PMPs), i.e. the points along the longitudinal coordinate for which the difference between the effective refractive index of the cores is null. The random fluctuations of the fiber parameters are emulated by random phase shifts (RPSs), introduced in the middle point between consecutive PMPs. Approximating the ICXT as the sum of the contributions associated with all the PMPs enables very fast ICXT estimates, compared to models that rely on solving numerically a system of differential equations for the longitudinal evolution of the electric field. The analytical models proposed in [13,14] rely on such approximation, providing very fast mean ICXT power estimates, yet assuming that one of the cores aligns with the center of the cladding. As a consequence, the models are constrained to a very specific type of MCF structure. Other computationally efficient models, relying on analytical estimates of the ICXT, have been proposed in [17][18][19][20]. However, these models assume that the inter-core coupling coefficient (ICCC) remains constant along the longitudinal coordinate. The models proposed in [22,24,25] take into consideration the longitudinal evolution of the ICCC, and the models proposed in [23][24][25] take into consideration non-linear transmission effects. However, these models rely on the numerical computation of the ICXT field and no analytical solution is provided. As a consequence, the ICXT estimates take a long time to obtain.
In this paper, the ICXT in weakly coupled single-mode MCFs with arbitrary core layout and including the dependence of the coupling coefficient on bending and twisting is studied. Our analysis considers the effect of the deformation of core shape caused by bending and twisting on the effective refractive index but neglects its effect on the birefringence and polarization rotation. In the interest of reducing the computational time needed to estimate the ICXT, taking into consideration the longitudinal bending and twisting perturbations (modeled as in [26]) of the ICCCs, a DCM for the ICXT is proposed. To the authors' knowledge, this is the first analytical model that takes into consideration the longitudinal bending and twisting perturbations of the ICCCs on the ICXT, and, at the same time, assumes an arbitrary core layout. Owing to these considerations, fast mean ICXT power estimates can be obtained for different MCF structures. That allows to easily compare their performance and decide which ones performs better.

ICXT field propagation and inter-core coupling coefficient
In this section, the coupled-mode equation (CME) used to model the ICXT field propagation in weakly coupled MCFs is presented. A step-index single-mode MCF, without trenches enveloping the cores, is considered throughout this paper. The cores of the MCF may have the same effective refractive indexes (perfectly homogeneous MCFs) or slightly different effective refractive indexes (real homogeneous MCFs). A brief overview of the ICCC expressions needed to model the ICXT field propagation is also presented.
The longitudinal evolution of the ICXT field, A n (z), in a weakly coupled MCF with N c cores and lossless linear transmission is given by [10,13,14,17,24] where θ m (z) = k 0 ∫ z 0 n e f f ,m (z )dz and θ n (z) = k 0 ∫ z 0 n e f f ,n (z )dz are the accumulated phases due to propagation from the input of cores m and n, respectively, up to the z coordinate, with n e f f ,m (z ) and n e f f ,n (z ) denoting the longitudinal variation of the effective refractive index of cores m and n, respectively. k 0 = 2π/λ is the wavenumber, with wavelength λ. κ nm (z) denotes the longitudinally-varying ICCC defined as κ nm (z) = (κ nm (z) + κ mn (z))/2 [17], where κ nm (z) and κ mn (z) are the longitudinally-varying ICCCs from core m to n, and from core n to m, respectively. The average of κ nm (z) and κ mn (z) is considered in Eq. (1) so that power conservation is ensured at the MCF output [17].
The formal definition of the ICCC from core m to n is given by [14,24,27] where NA n = n 2 n − n 2 cl denotes the numerical aperture of core n, n n is the refractive index of core n and n cl is the refractive index of the cladding. S n is the cross-sectional area of core n and S ∞ is the infinite cross-sectional area. In the presence of bending and twisting effects, the mode-fields are deformed relative to the case of absence of bending and twisting [28]. However, as R b /a > 10000 [28] (with a bending radius R b > 4 cm and a core radius a ≈ 4 µm) for single-mode MCFs, the field deformation is expected to be low. In the following, we neglect the effect of mode-field deformation described in [28] and the functions F i (r, ξ) (for i ∈ {m, n}), with r and ξ the radial and angular polar coodinates, are defined as in Eqs. (3) and (4) in [14].
An approximated closed-form expression of Eq. (2), generalized to MCFs whose cores have different refractive indexes and core radii, was proposed in [14]: Λnm am a m [a n w m J 0 (u n )I 1 ( a n w m a m )+a m u n J 1 (u n )I 0 ( a n w m a m )] where I l (x) and K l (x) are the l−th order modified Bessel functions of the first and second kind, respectively, and J l (x) is the l−th order Bessel function of the first kind. Υ nm = Γ * n Γ m /|Γ 2 n |, with Γ n and Γ m given by Eq. (7) in [14], and Γ * n denotes the complex conjugate of Γ n . w n and w m are given by Eq. (6) in [14], and u n and u m are given by Eq. (5) in [14]. Λ nm is the core pitch, and a n and a m are the radii of cores n and m, respectively.
As detailed in [14], the formal definition of κ mn can be obtained by replacing the index n by m, and vice-versa, in Eq. (2). By performing the same substitutions in Eq. (3), an approximated closed-form expression of κ mn is obtained. In section 3, the bending and twisting induced longitudinal variation of the ICCCs is described.

Longitudinal variation of the accumulated phases and ICCCs for an arbitrary core layout
In this section, the bending and twisting-induced longitudinal evolution of the accumulated phases and ICCCs, for an arbitrary core location within the cladding, is described.
As an example, Fig. 1 illustrates a 19-core double-ring MCF cross-section [5], where the origin of the system of coordinates aligns with the center of the cladding. The core pitch between any two adjacent cores m and n is denoted by Λ nm . The core location is defined in Cartesian coordinates, relative to the reference (center of the cladding). As an example, in Fig. 1, core 3 is identified as core m and core 2 is identified as core n. The Cartesian coordinates of the center of core m are denoted by x m and y m , and the ones of core n by x n and y n .
In the following, the bending and twisting-induced longitudinal evolution of the effective refractive index of the cores, which in turn affect the longitudinal evolution of the accumulated phases and ICCCs, is detailed. The longitudinal evolution of the effective refractive index of the cores is given by [8] n e f f ,m (z) = n (int) e f f ,m 1 + n e f f ,n (z) = n (int) e f f ,n 1 + x n cos(γz) − y n sin(γz) where n (int) e f f ,m and n (int) e f f ,n are the effective intrinsic refractive index of cores m and n, respectively, i. e. the effective refractive index of the cores in the absence of bending and twisting perturbations, and γ is the twisting rate. The difference between the effective refractive index of the cores is obtained from Eqs. (4) and (5), i.e. δn eq,nm (z) = n e f f ,n (z) − n e f f ,m (z) [13]. It is noted that the models proposed earlier by the authors in [13,14] assumed that the center of one of the cores aligns with the center of the cladding. As a consequence, these models correspond to the particular case of x n = y n = 0 and Λ nm = x 2 m + y 2 m , or x m = y m = 0 and Λ nm = x 2 n + y 2 n . Thus, the framework considered in this paper allows for a much more general study of the impact of the bending and twisting perturbations, modeled as in [26], on the ICXT.
Substituting Eqs. where ρ m = k 0 n (int) e f f ,m /(γR b ) and ρ n = k 0 n (int) e f f ,n /(γR b ). Following the coupled local-mode theory [25], the ICCCs with longitudinal dependence induced by bending and twisting, κ nm (z) and κ mn (z), can be obtained by substituting n e f f ,m and n e f f ,n by n e f f ,m (z) and n e f f ,n (z), respectively, in the parameters of Eq. (2) or (3). The longitudinal bending and twisting effects on n e f f ,m (z) and n e f f ,n (z) lead to the variation of κ nm (z) and κ mn (z) around the intrinsic ICCCs, κ (int) nm and κ (int) mn , respectively, for which n e f f ,m (z) = n (int) and n e f f ,n (z) = n (int) e f f ,n , i.e. the intrinsic ICCCs are the ones obtained in the absence of bending and twisting effects. The logarithm of the ICCC variation along the longitudinal coordinate can be well-approximated by a linear combination of the effective refractive index variation of the two cores, as discussed in further detail in section 5.1. This consideration generalizes the one presented in [14], in which one of the two cores is located at the center of the cladding, and simplifies the derivation of the mean ICXT power expression, detailed in section 4.2. We can then write ln(κ nm (z)) ≈ ln(κ for ∆n e f f ,n (z) = 0) and in the neighborhood of n (int) e f f ,n (and for ∆n e f f ,m (z) = 0), respectively. We stress that, in [14], core n is aligned with the origin of the coordinate system. Thus, s (m) nm is given by s nm derived in the Appendix B of [14], and, since κ mn can be obtained from κ nm by replacing the indexes n by m and vice-versa, as discussed in section 2, s (n) nm can be obtained from the expression of s mn , derived in Appendix C of [14], by substituting the indexes n by m, and vice-versa. Similarly, we can write ln(κ mn (z)) ≈ ln(κ (int) mn ∆n e f f ,n (z), where s (m) mn and s (n) mn are the slopes of ln(κ mn (z)) in the neighborhoods of n (int) e f f ,m (and for ∆n e f f ,n (z) = 0) and n (int) e f f ,n (and for ∆n e f f ,m (z) = 0), respectively. The analytical expression of s (m) mn is the same as the one of s mn provided in Appendix C of [14], and the one of s (n) mn can be obtained from the one of s nm (provided in Appendix B of [14]), by replacing the indexes n by m, and vice-versa.
Applying exponential to both sides of ln(κ nm (z))≈ ln(κ (int) nm ∆n e f f ,n (z) and ln(κ mn (z)) ≈ ln(κ (int) mn ∆n e f f ,n (z) yields nm ∆n e f f , m (z)+s (n) nm ∆n e f f , n (z) mn ∆n e f f , m (z)+s (n) mn ∆n e f f , n (z) .
In this work, κ (int) nm and κ (int) mn , used in the analytical model, are obtained from the closed-form expression of the ICCCs, Eq. (3), and so are the parameters s (n) nm , s (m) nm , s (n) mn and s (m) mn . In doing so, a good compromise between accuracy of the analytical mean ICXT power estimates, short computational times, and simplicity of the model is achieved. In particular, computing the intrinsic ICCCs from Eq. (2) takes longer than using Eq. (3), owing to the necessity to compute the integrals of Eq. (2) numerically. Computing the parameters s (n) nm , s (m) nm , s (n) mn and s (m) mn from Eq. (2) would be very complex. Thus, we compute such parameters from Eq. (3), using the analytical expressions presented in Appendix B of [14].

DCM for an arbitrary core layout
Many of the previous versions of the DCM [12][13][14] assume that one of the cores, considered in the evaluation of the ICXT, is aligned with the origin of the system of coordinates (center of the cladding). This alignment constrains the applicability of the DCM, since, in general, we may be interested in evaluating the ICXT between two cores that are not aligned with the center of the cladding. On the other hand, although the model proposed in [8] considers an arbitrary core layout, the longitudinal evolution of the ICCCs is neglected. For these reasons, we propose a novel DCM that accounts for both the bending and twisting-induced longitudinal variation of the ICCCs and accumulated phases, and an arbitrary core layout. In section 4.1, the DCM and the mean ICXT power estimation using the discrete changes coefficient are presented. In section 4.2, the derivation of the novel discrete changes coefficient is detailed.

DCM and mean ICXT power estimation from the discrete changes coefficient
Following the discussion in sections II-A and IV-A of [13,14], respectively, we consider that: (i) the interfered core n has no input signal; (ii) the crosstalk at the MCF output is low; (iii) the MCF is comprised of 2 cores for the sake of simplicity. Under these conditions, the ICXT field at the output of core n can be written as where A m (0) is the ICXT field at the input of the interfering core m, N = γL/π is the number of PMPs [13], with L the MCF length, φ nm,k is the RPS introduced at the k-th middle point between consecutive PMPs, modeled as in [13], and K nm,k is the discrete changes coefficient of the k-th fiber segment given by [14] K nm,k =(K nm,k +K mn,k )/2 (10) where z k−1 and z k are the longitudinal coordinates denoting the borders of the k-th segment. The analytical expressions of z k and z k−1 are found in the Appendix of [13]. From Eq. (9), the normalized mean ICXT power at the MCF output is given by [13,14] XT

Discrete changes coefficient for an arbitrary core layout
In this section, the derivation of the discrete changes coefficient, used to estimate the mean ICXT power for an arbitrary core layout, is provided. We start by substituting Eqs. (7) and (8) in Eqs. (11) and (12), respectively. Then, Eqs. (11) and (12) are replaced in the discrete changes coefficient shown in Eq. (10). In a first step, we include κ nm (z), shown in Eq. (7), onK nm,k , shown in Eq. (11). Then, we discuss the analytical derivation ofK mn,k , which is similar to the one ofK nm,k .
Replacing Eqs. (4) and (5) in Eq. (7) yields where where The effect of the difference of phase propagation from the MCF input up to the coordinate z k is separated in the definition of the discrete changes coefficient. By doing so, the effect of the field propagation from the MCF input is made explicit Setting with tan −1 (x) the inverse tangent function, we can write the exponent of the first exponential in the integrand of Eq. (17) as Then, the Jacobi-Anger expansion, shown in Eq. (23) of [14], is applied to the first exponential of the integrand of Eq. (17), yielding where ς nm = ∆β mn /γ. Following the discussion in section II-B of [13], the contribution of the discrete changes coefficient, given by Eq. (19), to the mean ICXT power can be separated in two: the contribution of the min-max segments and of the max-min segments. A min-max segment is one such that δn eq,nm (z) starts at a minimum point (occuring at z = z k−1 ) and ends at the next maximum (at z = z k ) [13]. Conversely, a max-min segment is one such that δn eq,nm (z) starts at a maximum point (occuring at z = z k−1 ) and ends at the next minimum (at z = z k ). Equations (42) and (43) in [13] are used to substitute the terms z k − z k−1 and z k + z k−1 in Eq. (19) in order to compute the discrete changes coefficient of the min-max segments and max-min segments. The twisting rate offset, ϕ, shown in Eq. (42) of [13] is given by The discrete changes coefficient of the min-max segments is then given bỹ and, in the case of the max-min segments, It is stressed that Eqs. (14)- (21) only account for the contribution of κ nm (z), shown in Eq. (11), to the discrete changes coefficient, shown in Eq. (10). However, as it can be understood from Eqs. (10) and (12), the contribution of κ mn (z) is also needed. The derivation and subsequent analytical expressions of the min-max and max-min segments ofK mn ,K mn,min−max andK mn,max−min , respectively, are similar to the ones ofK nm . The difference relies on the parameters of the ICCCs (s (m) nm , s (n) nm and κ (int) nm ), and on the parameters that depend on them, i. e. A nm , B nm and α nm . In particular,K mn,min−max is obtained by substituting κ (int) nm , s (m) nm and s (n) nm by κ (int) mn , s (m) mn and s (n) mn , respectively, in the parameters A nm , B nm and α nm shown in Eq. (20). Thus, forK mn,min−max : x n − ρ m y m + ρ n y n , and (iii) α nm becomes α mn = tan −1 (B mn /A mn ). κ (int) mn can be obtained by substituting all the terms that depend on n e f f ,m and n e f f ,n by n (int) e f f ,m and n (int) e f f ,n , respectively, in the expression of κ mn . The analytical expression ofK mn,max−min can be obtained by performing the substitutions previously described in Eq. (21).
From Eq. (13), considering a large number of segments, and noting that |K nm | = |K nm | (see Eq. (18)), the mean ICXT power can be written as a function of discrete changes coefficient of the min-max and max-min segments where K nm,min−max = (K nm,min−max +K mn,min−max )/2 Substituting Eqs. (20) and (21), and the analogous expressions ofK mn,min−max andK mn,max−min described previously, in Eq. (23), and then substituting Eq. (23) in Eq. (22) yields where and represent the contribution of the min-max and max-min segments, respectively, to the mean ICXT power. Although it was concluded in [13,14] that the contribution of the max-min and min-max segments to the mean ICXT power is the same, which could significantly simplify Eq. (24) since both contributions would yield the same value, this was confirmed mainly when one of the cores is aligned with the reference, i.e. x m = y m = 0 or x n = y n = 0. In section 5.2, the influence of the ICCC longitudinal variation on the max-min and min-max segment contribution to the mean ICXT power is discussed.

ICCCs variation with the effective refractive index for an arbitrary core layout
In this section, the discrepancy between the estimates of the ICCC expressions, and the ICCC variation, as a function of the bending and twisting-induced longitudinal variation of the effective refractive indexes, is studied numerically. The refractive index of cores m and n considered in this study are n m = 1.4453 and n n = n m (1 + ∆n (N ) ), respectively, where ∆n (N ) = 0.02% is the normalized difference of refractive index of the cores, ∆n (N ) = (n n − n m )/n m . These, and the remaining MCF parameters used in this section, are presented in Table 1. To compute the intrinsic effective refractive index of cores m and n, n (int) e f f ,m and n (int) e f f ,n , respectively, i.e. the effective refractive indexes in the absence of bending and twisting effects, the dispersion equation is solved for the HE 11 mode [27]: J 0 (u i )/(u i J 1 (u i )) = K 0 (w i )/(w i K 1 (w i )), i ∈ {m, n}, with u i and w i given by Eqs. (5) and (6) in [14], respectively. For the input parameters shown in Table 1 Similarly, the amplitude of the variation of n e f f ,n (z) along z is maximized by increasing x 2 n + y 2 n and minimizing R b . Choosing a small bending radius, e.g. R b = 0.1 m, and Cartesian coordinates of cores m and n located in the outer ring of the 19-core MCF shown in Fig. 1, e.g., x m = 2Λ nm , y m = 0, x n = 3Λ nm /2, y n = √ 3Λ nm /2, yields −0.09% < ∆n (N ) e f f ,m (z) < 0.09% and −0.078% < ∆n (N ) e f f ,n (z) < 0.078%, for a large core pitch Λ nm = 45 µm. To show the ICCCs variation over a wide range of the effective refractive index variation, we consider −0.1% ≤ ∆n (N ) e f f ,i (z) ≤ 0.1%, i ∈ {m, n}. Figure 2 shows the contours of the relative error, in percentage, between the expressions of κ nm , as a function of ∆n (N ) e f f ,m (z) and ∆n (N ) e f f ,n (z). In particular, Fig. 2(a) shows the relative error, in percentage, between κ nm obtained from the formal definition (see Eq. (2)) and the closed-form expression (see Eq. (3) Table 1. the relative error, in percentage, between the formal definition and the linearized expression, i.e. 100 × |k (Eq.(2)) nm − k (Eq. (7)) nm |/k (Eq.(2)) nm . The small relative error observed in Fig. 2(a) shows that Eq. (3) is very accurate over a wide range of the effective refractive index variation. The larger relative error shown in Fig. 2(b) is attributed to approximating the logarithm of the longitudinal variation of the ICCC as a linear combination of the effective refractive index variation, as explained in section 3. Lastly, the error observed in Fig. 2(c) is a combination of the ones shown in Figs. 2(a) and 2(b). Since the relative error shown in Fig. 2(a) is very small compared to the one of Fig. 2(b), the contour of Fig. 2(c) is very similar to the one of Fig. 2(b), i.e. the relative error observed for large variations of ∆n (N ) e f f ,m (z) and ∆n (N ) e f f ,n (z), is mostly imposed by the linearization error.
Since κ mn can be obtained from κ nm while replacing the indexes n by m, and vice-versa, the contours of κ mn are similar to the ones of κ nm , with the axis switched. Thus, the accuracy of κ nm is more sensitive to ∆n (N ) e f f ,m than to ∆n (N ) e f f ,n , whereas the accuracy of κ mn is more sensitive to ∆n  Table 1, similar conclusions were also drawn for other values of ∆n (N ) , a n and a m , typical of real homogeneous MCFs. The impact of the ICCC variation on the mean ICXT power is studied in section 5.3.

Analysis of max-min and min-max segment contribution to the mean ICXT power
In this section, the max-min and min-max segment contribution to the mean ICXT power is studied. As discussed in section 4.2, the mean ICXT power can be estimated analytically from Eq. (24), where the first and second addends of Eq. (24) represent the contribution of the max-min and min-max segments to the mean ICXT power. In previous studies, it was concluded that |K nm,min−max | 2 = |K nm,max−min | 2 , and thus the contribution of the max-min and min-max segments to the mean ICXT power was the same [13,14]. However, it will be shown that, under some circumstances, the bending and twisting effects lead to |K nm,min−max | 2 |K nm,max−min | 2 . Figure 4 shows the parameters used to obtain Fig. 4 are presented in Table 1, except that a n = a m = 4 µm, and y m = 65 µm and y n = 20 µm are considered. Figure 4 Fig. 4(a) reveals that the integral of (κ nm (z) + κ mn (z))/2, with respect to z, is the same within a max-min and min-max segment. Thus, the contribution of the max-min and min-max segments to the mean ICXT power is expected to be similar. Indeed, it was observed that the max-min and min-max segment contribution to the mean ICXT power is 50% each. In contrast, the integral of (κ nm (z) + κ mn (z))/2 is not the same within a max-min and min-max segment, in the cases shown in Figs. 4(b) and 4(c) (and especially in Fig. 4(c)). Thus, the contribution of the max-min and min-max segments to the mean ICXT power is expected to be different, especially in the case of Fig. 4(c). Indeed, the contribution of the max-min and min-max segments is 94% and 6%, respectively, in the case of Fig. 4(b), and 99.7% and 0.3%, respectively, in the case of Fig. 4(c). It is therefore concluded that, in order to observe similar contribution of the max-min and min-max segments to the mean ICXT power, the longitudinal variation of (κ nm (z) + κ mn (z))/2 and n e f f ,n (z) − n e f f ,m (z) have to be in-phase (Fig. 4(a)), or in phase opposition. Otherwise, the integral of (κ nm (z) + κ mn (z))/2 in a max-min segments differs from the one of a min-max segment (Figs. 4(b) and 4(c)), and thus the contribution of the max-min and min-max segments to the mean ICXT power is not the same. Furthermore, the Fig. 5. Most of the mean ICXT power estimates are obtained considering an MCF structure obeying the 19-core MCF design (cases 1-4). Cases 1-3 consider a single interfering core, and case 4 considers multiple interfering cores. A variant of case 1, which does not obey the 19-core MCF design, is also studied. The Cartesian coordinates of cores m (x m and y m ) and core n (x n and y n ) are detailed for each case. results of Figs. 4(a)-4(c) show that the amplitude of (κ nm (z) + κ mn (z))/2 is larger when the cores are located farther away from the center of the cladding. This behavior is attributed to the larger bending and twisting influence on the ICXT of those cores.

Mean ICXT power estimates in a 19-core double-ring MCF structure
In this section, the mean ICXT power is studied for different core locations within the cladding. The mean ICXT power estimates obtained with the analytical model proposed in section 4.2 are compared with the ones obtained through numerical simulation. The simulated mean ICXT power estimates, which can take many hours to compute, are obtained from the numerical computation of Eq. (1), with a step of 10 −4 m. Random fluctuations of the MCF parameters, emulated by RPSs, are inserted at every middle point between consecutive PMPs [13,14]. The RPSs are modeled similarly to [13,14]. Figure 5 shows five cases for which the mean ICXT power is assessed. Cases 1-4 consider that the location of cores m and n obeys the 19-core MCF structure illustrated in Fig. 1, where the core pitch is the same for all adjacent cores, Λ nm . While cases 1-3 consider a single interfering core, case 4 considers multiple interfering cores. For the sake of comparison, a variant of case 1, whose core locations do not fit in the 19-core MCF design, is also studied. Figure 6 shows the mean ICXT power, as a function of ∆n (N ) = (n n − n m )/n m , for cases 1-3 and for the variant of case 1, with Λ nm = 45 µm. The parameters considered while obtaining the results of Fig. 6 are: R b = 0.1 m, γ = 4π rad/m, a n = a m = 4 µm and L = 200 m. The remaining parameters can be found in Table 1. In Fig. 6(a), the analytical estimates of mean ICXT power, obtained from Eq. (24), with ν integers ranging from −25000 to 25000, are compared with the  Fig. 6(a) allows to conclude that: (i) the mean ICXT power estimates are very different from case-to-case. This behavior follows the discussion provided in section 5.2 (see Fig. 4), where it is noted that the amplitude of the longitudinally-varying ICCCs, affecting the electric field propagation, depends on the core location. When comparing the cases shown in Fig. 6(a) where both cores are at the same distance from the cladding center, i.e. case 2 (with cores m and n at a distance Λ nm = 45 µm from the cladding center), and the variant of case 1, where cores m and n are at a distance 2Λ nm = 90 µm from the cladding center, a similar variation of the mean ICXT power with ∆n (N ) is observed; (ii) the mean ICXT power in cores located in the outer ring of the 19-core MCF (case 1) is much higher than the one of the inner ring (case 2), exceeding 10 dB. Thus, the cores located in the outer ring are more susceptible to ICXT, owing to bending and twisting perturbations; (iii) the analytical mean ICXT power estimates agree satisfactorily with the simulation ones. When the cores are located closer to the cladding center, e.g. case 2, the discrepancy between the simulation and analytical estimates does not exceed 0.6 dB. However, the discrepancy becomes larger in case 1 and variant of case 1, not exceeding 2.2 dB. The larger discrepancies observed in these cases are mostly attributed to the error committed in the linearization of the ICCCs, which is larger for larger longitudinal variations of the effective refractive indexes (as addressed in section 5.1) and, equivalently, when the cores are located farther away from the cladding center. Nevertheless, the accuracy of the analytical mean ICXT power estimates is satisfactory, and the variation of the mean ICXT power with ∆n (N ) is similar when comparing the simulated and analytical estimates. Furthermore, the analytical mean ICXT power estimates are obtained very quickly when compared to the simulated ones: the simulated mean ICXT power estimates can take many hours to obtain, owing to the need to solve Eq. (1) with a very small step (10 −4 m); the analytical mean ICXT power estimates take a few seconds to obtain. In Fig. 6(b), the analytical results with and without considering the z−variation of the ICCCs are compared (the z−variation of the accumulated phases is considered in all cases). The analytical mean ICXT power estimates without z−variation of the ICCCs are obtained similarly to the ones that consider such dependence, with the difference that s (m) nm = s (n) nm = s (m) mn = s (n) mn = 0, i.e. the parameters that model the z-dependence of the ICCCs are null. As it can be seen in Fig. 6(b), the mean ICXT power estimates that do not consider the z-dependence of the ICCCs are very similar, regardless the case considered. Thus, it is concluded that the mean ICXT power estimates can be very sensitive to the ICCCs longitudinal variation. In the cases considered in Fig. 6(b), the error committed in the mean ICXT power estimates by neglecting the ICCC longitudinal variation can exceed 15 dB. Figure 7 shows the mean ICXT power as a function of ∆n (N ) , for cases 1, 3 and 4. The results corresponding to cases 1 and 3 have been shown previously in Fig. 6(a), and are shown again to enable a better understanding of the results of case 4. In particular, case 4 corresponds to the combination of cases 1 and 3, as it can be understood from Fig. 5. Thus, it is expected that the mean ICXT power estimates of case 4 are the ones of cases 1 and 3 added together (in linear units). Since the mean ICXT power estimates of case 1 are much higher than the ones of case 3, the mean ICXT power estimates of case 4 are expected to be similar to the ones of case 1, yet slightly higher. The results of Fig. 7 validate this expectation.

Mean ICXT power dependence on the bending and twisting parameters
In this section, the mean ICXT power dependence on the bending and twisting parameters is analyzed. Figure 8 shows the mean ICXT power as a function of the bending radius. In particular, Fig. 8(a) shows the mean ICXT power estimates obtained through numerical simulation and analytically from Eq. (24). Cases 1 and 2 are the ones shown in Fig. 5. The parameters used to obtain the results of Fig. 8(a) are the same as the ones used in Fig. 6(a), with ∆n (N ) = 0.02%. A range of bending radii between 0.1 m and 0.25 m is shown since R b = 0.1 m is considered a small bending radius and the critical bending radius, R pk (see Eq. (17) in [8]), i.e. the maximum bending radius for which the DCM is valid, is around 0.27 m, with the considered parameters. Figure 8(a) reveals good agreement between the analytical and simulation estimates of the mean ICXT power. The worst case of discrepancy between analytical and simulation results occurs for case 1 and R b = 0.1 m. This behavior is attributed to the following: (i) a smaller bending radius leads to a larger amplitude variation of the effective refractive index of the cores, as it can be concluded from Eqs. (4) and (5). As a consequence, the amplitude variation of the ICCCs is also larger, and the error committed in the linearization of the logarithm of the ICCCs increases. The discussion of the linearization error follows the one provided in section 5.1; (ii) when the cores are located farther away from the center of the cladding, i.e. case 1, larger amplitude variation of the ICCCs is observed. Thus, the linearization error is larger in case 1. In Fig. 8(b), the analytical estimates of mean ICXT power already shown in Fig. 8(a) with z-dependent ICCCs are compared with the ones without z-dependent ICCCs. The results of Fig. 8(b) agree with expectations, since the largest improvements of the mean ICXT power estimates, achieved by considering z-dependent ICCCs, occur when the amplitude variation of the ICCCs is high, i.e. for the smallest bending radii and case 1. Figure 9 shows the mean ICXT power as a function of the twisting rate (see Fig. 9(a)) and core pitch (see Fig. 9(b)). The parameters used in Fig. 9(a) are the same as the ones of Fig. 8, with R b = 0.1 m. The twisting rate considered in our analysis assumes typical values (0.1 ≤ γ/(2π) ≤ 10), within the range considered by other authors [10]. The analytical mean ICXT power estimates are shown to agree satisfactorily with the simulation ones. The largest discrepancy between simulation and analytical mean ICXT power estimates occurs for case 1. This behavior was already observed and explained in Fig. 8(a), for R b = 0.1 m. It was concluded from further testing that this discrepancy becomes smaller for larger bending radii (0.15 m and 0.25 m were tested), which is also observed in Fig. 8(a).
The results of Fig. 9(b) use the same parameters as the ones of Fig. 9(a), with γ = 4π rad/m. The curves (c1) and (c2) shown in Fig. 9(b) correspond to ∆n (N ) = −0.026% and ∆n (N ) = 0% (perfectly homogeneous MCF), respectively. Figure 9(b) reveals that smaller differences between the mean ICXT power of cases 1 and 2 tend to occur for smaller core pitches. These results are expected since, as it can be seen in Fig. 5, smaller core pitches reduce the distance of the cores to the center of the MCF and between the rings of the 19-core MCF. It was concluded from Fig. 9(b) that, for ∆n (N ) = −0.026%, the difference between the analytical mean ICXT power of cases 1 and 2 ranges from approximately 0.4 dB (for Λ nm = 30 µm) to 6.3 dB (for Λ nm = 45 µm). For ∆n (N ) = 0%, the range is from approximately 3.3 dB (for Λ nm = 30 µm) to 9.9 dB (for Λ nm = 45 µm). The accuracy of the model and improvement of the mean ICXT power estimates enabled by considering z-dependent ICCCs are expected to be similar for fiber lengths longer than 200 m. Noting that the dependence of the analytical mean ICXT power estimates on the fiber length results from the number of PMPs, N = γL/π (see Eqs. (24)-(26)), the mean ICXT power is expected to increase linearly with the fiber length, with and without z-dependent ICCCs.

Conclusion
We have investigated the bending and twisting-induced longitudinal variation of ICCC and its effect on the ICXT in weakly coupled MCFs with an arbitrary core layout. An analytical DCM for the ICXT field propagation under those conditions has been proposed for the first time, providing very fast and rather accurate mean ICXT power estimates. The analytical mean ICXT power estimates have been validated through numerical simulation. It has been predicted that the mean ICXT power between adjacent cores of the outer ring of the 19-core MCF can be more than 10 dB higher than the one between adjacent cores of the inner ring. It has also been predicted that the difference between the mean ICXT power of cores in the inner and outer rings can be much smaller by decreasing the core pitch and increasing the bending radius. This behavior has been attributed to the ICXT dependence on the bending and twisting-induced longitudinal variation of the ICCCs. In particular, larger bending and twisting-induced fluctuations of the ICCCs along the longitudinal coordinate are observed in the cores of the outer ring, but the fluctuations become smaller for smaller core pitches and larger bending radii. Furthermore, it has been shown that, if the ICCCs longitudinal variation is neglected, the mean ICXT power estimates between two adjacent cores are very similar despite the location of those cores. This means that neglecting the longitudinal variation of the ICCCs can lead to misleading estimates of the mean ICXT power, with an error exceeding 15 dB.
For future work, the analytical DCM may be improved to consider other MCF profiles, e.g. trench-assisted MCFs. The impact of the birefringence caused by the deformation of the core shape on the ICXT in weakly coupled MCFs with an arbitrary core layout may also be assessed.