Coherent propagation of powerful out-of-phase wave beams in linear arrays of weakly coupled cores

Transverse filamentation instability in media with cubic nonlinearity is the most important mechanism restricting the propagation of high-power laser pulses [1]. In fact, this is a fundamental limitation for the generation of high-power laser pulses. Currently, several mechanisms are proposed to bypass this dangerous self-focusing instability. First, the transverse filamentation instability in media with an inertial nonlinear response [2] has no specific transverse scale with the maximum growth rate. Second, the instability is suppressed for partially coherent radiation, if the correlation radius is less than the scale of the maximal growth rate [3–5]. Third, the filamentation instability can be suppressed in the case of few-cycle laser pulses, for which the synchronization between pulses and growing field perturbations can be disturbed [6,7]. This takes place for intense laser pulses of ultrashort durations, when the nonlinear dependence of the group velocity on the light intensity comes into play.

To eliminate filamentation instability, multicore fibers (MCF) can be used to achieve coherent propagation of laser pulses at a power level significantly exceeding the critical self-focusing power in homogeneous media [8–12]. However, there is a specific critical self-focusing power for in-phase wave field distributions in such MCF [13–19]. Nevertheless, a number of interesting results were obtained for such systems: generation of a supercontinuum [20–22], shortening of the laser pulse duration [15,16,21–25], demonstration of the nonlinear switching [26], and formation of light bullets [16,21–24,27–31]. Recently, strongly varying stable nonlinear solutions have been found in MCF from cores located in a ring, allowing to operate with a total power being noticeably larger than the power that can be transmitted in a single core [32]. Moreover, these out-of-phase nonlinear solutions still survive even in a strongly deformed MCF [25].

However, the technological problems in the manufacturing of such MCFs are quite large. At the same time, significant success was achieved in the creation of MCF with quasi-one-dimensional coupled cores [16]. These MCFs consist of a small number of lines with tens or hundreds of cores along. Therefore, the generalization of the methods of coherent radiation transport to such quasi-one-dimensional MCFs seems to be interesting and important.

In this paper, we generalize the out-of-phase distribution of the wave field [32] to the case of MCF that consist of one-dimensional and two-dimensional lattices of cores. Let the fundamental guided modes of optical waveguides oriented parallel to the z-axis be weakly coupled to each other. Then, the propagation of laser radiation in the MCF can be approximately described as a superposition of modes localized in each core. Considering one-dimensional arrays of cores and assuming all cores to be the same, we obtain the system of equations for the field amplitudes in the cores [33,34]

$$\begin{equation} \mathrm{i} \frac{\partial u_n}{\partial z}+u_{n+1}+u_{n-1}+|u_n|^2u_n-\mathrm{i} \gamma u_n=0. \end{equation} \tag{ 1 }$$

The envelope amplitude u_n in the n-th core changes during propagation along the z-axis under the influence of the following factors: the Kerr nonlinearity of a single core, amplification in an active media, and interaction with the nearest neighboring cores. Here, the index n varies from 1 to N, the evolutionary coordinate z is normalized to the coupling coefficient χ between neighboring cores, the envelope complex amplitude of the wave packet is measured in $\sqrt{\chi/\beta}$, β is the nonlinearity coefficient, and $\gamma = \Gamma / \chi$ is the normalized gain ratio. Equation (1) preserves the total power $P = \sum |u_n|^2$ in the absence of amplification $(\gamma = 0)$. The applicability of eqs. (1) is limited by the approximation of the single-mode wave field propagation in each core. It is violated, when the radiation power in any core becomes close to the critical self-focusing power in homogeneous media $\mathcal{P}_{\textrm{cr}}$.

Let us consider the case of no gain $(\gamma=0)$ in this section. The dispersion relation for equation (1) is $k_z = 2\cos \varkappa$ for $\mathrm{e}^{\mathrm{i} k_z z + \mathrm{i} \varkappa n}$ processes. Since we are interested in stationary distributions, we will search for them near the points of zero group velocity: $\varkappa=0$ and $\varkappa=\pi$. They correspond to in-phase $u_n(z) = s(n,z) \mathrm{e}^{2\mathrm{i} z}$ and out-of-phase $u_n(z) = (-1)^n a(n,z) \mathrm{e}^{-2\mathrm{i} z}$ field distributions. Here s >  0 and a >  0 are the envelopes of the desired solutions, $n=0 \ldots N+1$ is the coordinate along the lattice. We specially added two "virtual" points $n = 0$ and $n = N + 1$ at the boundaries, so that the field at these points is zero:

$$\begin{equation} s(0,z) = s(N+1,z) = a(0,z) = a(N+1,z) = 0.\quad \end{equation} \tag{ 2 }$$

Let us analyze a one-dimensional lattice consisting of a large but finite number of cores $N \gg 1$. Assuming that the envelopes of the field are a smooth function

$$\begin{equation} \left|\frac{\partial a}{\partial n}\right| \ll \pi |a|, \qquad\left|\frac{\partial s}{\partial n}\right| \ll \pi |s|, \end{equation} \tag{ 3 }$$

we arrive at the nonlinear Schrödinger equations for in-phase and out-of-phase envelopes

$$\begin{eqnarray} &&\mathrm{i} \frac{\partial s}{\partial z} + \frac{\partial^2 s}{\partial n^2} + |s|^2 s=0, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&\mathrm{i} \frac{\partial a}{\partial z} - \frac{\partial^2 a}{\partial n^2} + |a|^2 a=0. \end{eqnarray} \tag{ 5 }$$

Here and below, we analyze in-phase and out-of-phase solutions in parallel due to their similarity. A distinctive feature of eq. (5) for the out-of-phase envelope is the presence of the negative sign of the second derivative. As a result, the focusing nonlinearity of the medium behaves as the defocusing one for the out-of-phase distribution of the wave field. This suggests the existence of a stable out-of-phase distribution of the wave field, which is not subject to filamentation instability.

We will seek a solution in the form $a = a_0 A(\varkappa n) \mathrm{e}^{-\mathrm{i} \lambda z}$ and $s = a_0 S(\varkappa n) \mathrm{e}^{-\mathrm{i} \lambda z}$, where a₀ is the wave beam amplitude and λ is the longitudinal wave number. Substituting this expression into eqs. (4), (5) we obtain the equation for the envelopes A, S

$$\begin{eqnarray} \frac{\partial^2 S}{\partial n^2} +\lambda S + S^3 &= 0, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} -\frac{\partial^2 A}{\partial n^2} +\lambda A + A^3 &= 0. \end{eqnarray} \tag{ 7 }$$

These equations have a solution in the form $A(\xi)=\mathop{\mathrm{sn}}(\xi,m)$, $S(\xi)=\mathop{\mathrm{cn}}(\xi,m)$ with parameters obtained from the boundary conditions (2):

$$\begin{equation} 2m \varkappa^2=a_0^2, \quad (N+1) \varkappa = 2 K(m), \quad \lambda = (1+m) \varkappa^2.\quad \end{equation} \tag{ 8 }$$

Here, $K(m) =\int_0^{\pi/2} \frac{\textrm{d}\varphi}{\sqrt{1-m \sin^2 \varphi}}$ is the complete elliptic integral of the first kind, the parameter m varies in the range $[0,1]$ and is found from the transcendental equation

$$\begin{equation} m K(m)^2 = \frac{1}{8} (N+1)^2 a_0^2. \end{equation} \tag{ 9 }$$

The graphical solution for the parameter m is shown in fig. 1(a), which has the form of a step with an inflection point $(N+1)a_0 \approx 4$. The asymptotics of this solution at small and large amplitudes have the form

$$\begin{equation} m \approx \left\{\begin{array}{l@{\quad}l} (N+1)^2 a_0^2/2 \pi^2, & (N+1) a_0 \ll 4,\\ 1 - 16\, \mathrm{e}^{-(N+1) a_0/\sqrt{2}}, & (N+1) a_0 \gg 4. \end{array}\right. \end{equation} \tag{ 10 }$$

It is seen that the asymptotics for the solution of the parameter m are in good agreement with the graphic solution. Accordingly, for the transverse wave number ϰ , we obtain approximations

$$\begin{equation} \varkappa = \frac{2 K(m)}{N+1} \approx \left\{\begin{array}{l@{\quad}l} \displaystyle\frac{\pi}{N+1}, & (N+1) a_0 \ll 4, \\ a_0/\sqrt{2}, & (N+1) a_0 \gg 4. \end{array}\right. \end{equation} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

Asymptotics of the parameters $\varkappa \approx \pi/(N+1)$, $m \approx 0$ for small amplitudes give the same solutions similar to the field in a plane waveguide [35]

$$\begin{equation} s \approx a \approx a_0 \sin \frac {\pi n} {N + 1}, \qquad a_0 \ll \frac {4} {N + 1}. \end{equation} \tag{ 12 }$$

Asymptotics of the parameters $\varkappa \approx a_0/\sqrt{2}$, $m \approx 1-16 \mathrm{e}^{-(N+1)\varkappa}$ for large amplitudes give fundamentally different solutions

$$\begin{eqnarray} \hspace{-25pt}s &\approx a_0 \delta_{n,\frac{N+1}{2}},\hspace{20pt} a_0 \gg \frac{4}{N+1}, \end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray} \hspace{-25pt}a &\approx a_0, \hspace{49pt} a_0 \gg \frac{4}{N+1}. \end{eqnarray} \tag{ 13b }$$

$\delta_{ij}$

The in-phase distribution of wave beams tends to localize in one core with an increase in amplitude (13). Amplitude distributions are shown in fig. 2(a). This is a discrete analogue of the NSE soliton [33]. We note that the applicability conditions (3) of the shortened solution for s are violated at large amplitudes. Nevertheless, the tendency toward localization and self-isolation of the field in individual core is maintained [6,33].

Zoom In Zoom Out Reset image size

On the contrary, the out-of-phase distribution of wave beams evenly fills all the optical cores with an increase in amplitude (13a). This effect can be used to transport strongly powerful laser beams, whose power can significantly exceed the critical self-focusing power in single core. Amplitude distributions are shown in fig. 2(b). Thus, the most useful stationary solution for the complex envelopes of the wave field has the form

$$\begin{equation} u_n \approx (-1)^n a_0 \mathrm{e}^{\mathrm{i} \lambda z} \cdot \left\{\begin{array}{l@{\quad}l} \sin \displaystyle\frac{\pi n}{N+1}, & (N+1) a_0 \ll 4, \\ 1, & (N+1) a_0 \gg 4. \end{array}\right. \end{equation} \tag{ 14 }$$

Solution (14) was obtained within the framework of the approximate equation (4). The discreteness presented in eq. (1) can violate the stability of this approximate solution and lead to stochastic dynamics, if the difference in the nonlinear phase change $\Delta \varphi$ between adjacent cores at the coupling length between them (equal to 1 for eq. (1)) will exceed π [17,36]. In our case, the phase change difference is equal to the difference in intensities $\Delta \varphi = \max_n | |u_{n+1}|^2-|u_n|^2 | = |u_2|^2-|u_1|^2$ near the boundary, where the gradients are maximal. The value of $\Delta \varphi$ can become somehow noticeable only for the field amplitudes $u_n \sim a_0 > 1$. So, it suffices to estimate the phase change difference in the region $(N + 1) a_0 \gg 4$. Then the series expansion of the solution (14) at the boundary yields

$$\begin{eqnarray} \Delta \varphi &\le |u_2|^2-|u_1|^2 = a_0^2 [\mathop{\mathrm{sn}}(2\varkappa, m)^2-\mathop{\mathrm{sn}}(\varkappa, m)^2]\approx \nonumber\\ & a_0^2 [\tanh^2(\sqrt{2} a_0)-\tanh^2(a_0/\sqrt{2})] < 0.8. \end{eqnarray} \tag{ 15 }$$

Here, we used approximations $\mathop{\mathrm{sn}}(x, m\to 1) \approx \tanh(x)$ and $\varkappa\approx a_0/\sqrt{2}$ for $(N+1)a_0 \gg 4$ (see fig. 1). The maximum phase change difference $\Delta \varphi \simeq 0.8 \ll \pi$ is achieved at the field amplitude $a_0 = u_* \approx 1.852$. So, reaching the threshold of stochastic dynamics on solution (14) is impossible (fig. 1(c)). The reason is as follows. Despite the large values of ϰ for $a_0 \gg 4/(N+1)$, the entire difference in the wave field distribution a(n) from 0 to a₀ falls into the "virtual" region $0 \le n<1$, and there are no sharp gradients of the wave field in the physically correct region $1 \le n \le N$.

The knowledge of out-of-phase solution (14) allows us to predict the nonlinear dynamics of the wave field in the active MCF when the power and, correspondingly, the maximum amplitude $u_{\textrm{max}}$ increase on the scale of the gain length being much longer than the coupling one $(\gamma \ll 1)$. As the wave beam power increases exponentially, we should expect adiabatically smooth equalization of the amplitudes in all cores for $u_{\textrm{max}}> u_{\textrm{th}} \equiv 4/(N + 1)$. Adiabaticity will not be violated if the gain ratio is less than the difference in the propagation constants between the out-of-phase mode and other MCF modes. The propagation constants in the MCF are $h_n = 2 \cos \frac {\pi n}{N+1}$. The smallest of them, $h_N = 2\cos \frac{\pi N}{N+1}$, corresponds to the found out-of-phase distribution, eq. (14). Then, the condition of the amplification adiabaticity takes the form

$$\begin{equation} \gamma \ll \sin \frac{3\pi}{N+1} \sin \frac{\pi}{N+1} \mathop{\approx}\limits_{N \gg 1} \frac{3\pi^2}{(N+1)^2}. \end{equation} \tag{ 16 }$$

It can be seen from the expression that the limiting gain, at which there is no excitation of other modes, decreases quite rapidly as the number of the cores increases.

Figure 3 shows the evolution of the spatial distribution of the wave beam in the active MCF with the gain ratio $\gamma = 0.01$. The calculations were performed for MCF with 9 cores (black curves) and 49 cores (blue curves) in the framework of eq. (1). The field amplitude is normalized to the maximum value $u_{\textrm{max}}$. It is seen that the beam distribution is close to $\sin(\frac {\pi n}{N + 1})$ up to the amplitudes $u_{\textrm{max}}\lesssim u_{\textrm{th}}=4/(N+1)$. The adiabatic equalization of the amplitude distribution is observed with a further increase in the field amplitude during amplification. When the field amplitude is $u_n\gtrsim u_*\approx 1.852 > u_{\textrm{th}}$, the beam distribution over all cores becomes completely uniform. Thus, the results of numerical simulations confirm the stability of the found solution (14) and the adiabatic equalization of the beam distribution from the inhomogeneous to uniform ones with an increase in the total power.

Zoom In Zoom Out Reset image size

Out-of-phase solution (14) can be generalized to the case of the lattice $2\times N$. In this case, the system of equations for the field amplitudes takes the form

$$\begin{eqnarray} & \mathrm{i} \dfrac{\partial u_{1,n}}{\partial z}+u_{1,n+1}+u_{1,n-1}+u_{2,n}\nonumber\\ & +\mu(u_{2,n-1}+u_{2,n+1})+|u_{1,n}|^2u_{1,n}={0}, \end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray} &\mathrm{i} \dfrac{\partial u_{2,n}}{\partial z}+u_{2,n+1}+u_{2,n-1}+u_{1,n}\nonumber\\ &+\mu(u_{1,n-1}+u_{1,n+1})+|u_{2,n}|^2u_{2,n}={0} . \end{eqnarray} \tag{ 17b }$$

Neglecting the diagonal coupling $(\mu = 0)$ and making the change of variables $u_{k,n} = (-1)^k g_n \exp(-\mathrm{i} z)$, we arrive at the equation

$$\begin{equation} \mathrm{i} \frac{\partial g_n}{\partial z} + g_{n-1} + g_{n+1} + |g_{n}|^2 g_{n}={0}, \end{equation} \tag{ 18 }$$

which matches eq. (1) with $\gamma=0$ exactly. Thus, found solution (14) for a one-dimensional lattice can be easily generalized for the $2 \times N$ lattice.

The presence of a weak diagonal coupling $(\mu \neq 0)$ allows the development of filamentation instability. Let us determine the conditions for the development of this instability. It has been shown above that the amplitude distribution in the cores is uniform (eq. (14)) in the case of $u>u_*>u_{\textrm{th}}$. Then, the unperturbed solution with the amplitude a₀ can be written in the form $u_{k,n}^{(s)}=a_0(-1)^{k+n}\mathrm{e}^{\mathrm{i} hz}$, where $h=3-a_0^2-2\mu$. Linearization of eqs. (17) for a wave field of the form $u_{1,n}=a_0[(-1)^n+\delta a\mathrm{e}^{\mathrm{i} \lambda z + \mathrm{i} \varkappa n}]\mathrm{e}^{\mathrm{i} h z}$, $u_{2,n}=a_0[(-1)^{n+1}+\delta b\mathrm{e}^{\mathrm{i} \lambda z + \mathrm{i} \varkappa n}]\mathrm{e}^{\mathrm{i} h z}$ with small perturbations $\delta a, \delta b \ll 1$ gives the dispersion relations

$$\begin{eqnarray*} \lambda_1^2&= 4(1+\cos\varkappa)(1-\mu)[a_0^2+(1-\mu)(1+\cos\varkappa)], \\ \lambda_2^2&= 4[\cos\varkappa(1+\mu)+2-\mu][a_0^2+2-\mu+(1+\mu)\cos\varkappa] \end{eqnarray*}$$

from the condition for the existence of nonzero perturbations. Analysis of the dispersion relations shows the presence of filamentation instability $(\mathrm {Im} \, \lambda_{1,2} <0)$ for the diagonal coefficient $\mu \ge 1/2$. This corresponds to a rather large coupling coefficient, which violates the approximation of weakly coupled cores that we use. Therefore, the found out-of-phase distribution of the wave field in the two-dimensional lattice $2 \times N$ is not subject to filamentation instability.

To confirm the existence of the found out-of-phase mode in the case of a two-dimensional lattice and adiabatic equalization of amplitude distribution with the increasing power, detailed numerical studies were carried out within the framework of a more general model described by the nonlinear unidirectional wave equation [32]

$$\begin{eqnarray} \mathrm{i} \frac{\partial \mathcal{E}}{\partial z} =\,& \sqrt{k_0^2 n_0^2 + \Delta_\perp} \mathcal{E} + k_0 n_2 |\mathcal{E}|^2 \mathcal{E} \nonumber\\ & {}+ k_0 \delta n\, U(x,y) \mathcal{E} + \mathrm{i} \Gamma \mathcal{E}, \end{eqnarray} \tag{ 19 }$$

with the potential $U = \sum_n \exp[-\frac{((x-x_n)^2+(y-y_n)^2)^2}{r_n^4}]$. Here, x_n , y_n , r_n are the position and radius of the cores, n₀ is the refractive index of the cladding, $\delta n$ is the difference between the refractive indices of the cores and the medium, and n₂ is the nonlinear refractive index. The operator $\sqrt{k_0^2 n_0^2 + \Delta_\perp}$ can be easily calculated in Fourier space and allows one to properly describe wave fields with transverse scales of order of the wavelength by taking into account spherical aberrations.

The calculations were performed for the wavelength $\lambda = 1.03\ \mu \text{m}$ for a silica fiber (refractive index $n_0 = 1.45$ and the critical self-focusing power 4 MW). The difference in the refractive indices between the cores and the cladding was $\delta n = 0.005$. The diameters of the cores with a super-Gaussian profile were $6\ \mu \text{m}$, and the distance between them was $12\ \mu \text{m}$. The fields in each core were specified in the form of a Gaussian function closest to the field of the fundamental mode of an individual core.

Figure 4 shows the results of numerical simulations in the active MCF for a two-dimensional lattice $2 \times 9$ with a gain ratio of $\Gamma = 0.02\ \text{cm}^{-1}$. The initial wave field have a low power $\mathcal{P} = 10^{-2} \mathcal{P}_{\textrm{cr}}$ with amplitudes in each of the cores determined by the found solution $u_{k,n} \propto (-1)^{k+n} \sin(\frac{\pi n} {10})$. The initial noise level was about 3%. It is seen that the field amplitudes in all the cores become uniform with the growth of the total power, which confirms the stability of the found solution with respect to the filamentation instability. Numerical results show the possibility of transporting up to $0.5 \mathcal{P}_{\textrm{cr}}$ in each core, even if the wave field narrowing in the cores becomes noticeable (fig. 4). For an MCF of $2 \times 9$ cores, this is about 35 MW of total power. Thus, the results of numerical simulation confirm the stability of the solutions found for the two-dimensional matrix and adiabatic equalization of amplitudes in the found out-of-phase mode with the increasing total power in the active MCF. Along with this, simulation was carried out for a larger number of optical cores. For example, according to (16), one should use $\Gamma = 0.001\ \text{cm}^{-1}$ for an MCF of $2\times 49$ cores, which give a quasi-uniform distribution of wave field intensities on the length of less than 1 m. In this case, the total power will be about 200 MW.

Zoom In Zoom Out Reset image size

The propagation of laser beams in a multi-core fiber consisting of $1 \times N$ and $2 \times N$ lattices of weakly coupled cores is considered. Exact stable nonlinear solutions are found for the out-of-phase mode, which describe the coherent propagation of wave beams in such fibers. Equalization of the wave field amplitudes in all cores with an increase in the total radiation power is shown. This makes it possible to operate with laser radiation with a power 2N times greater than the critical self-focusing power in homogeneous media. The stability of the solutions found is studied both with respect to filamentation instability and with respect to weak amplification in the active medium. The analytical results are confirmed by direct numerical simulations within the framework of the wave equation.

This work was supported by the Center of Excellence "Center of Photonics" funded by the Ministry of Science and Higher Education of the Russian Federation, contract No. 075-15-2020-906.