The new phenomena arising from Anderson localization in nonlinear optical systems attract significant research interest. Previous experiments have achieved the localization of light waves within disordered media, and extensive theoretical studies show the presence of optical solitons in nonlinear Schr?dinger equations show disordered potentials. However, the presence, stability, and dynamics of certain disordered solitons in cubic-quintic nonlinear dielectric waveguides with disordered potentials remain unexplored. Additionally, the phenomenon and mechanism of competitive cubic-quintic nonlinearity interactions in disordered waveguide arrays are not fully understood. This study provides a new approach for exploring the physical properties of disordered solitons in high-order nonlinear media.

We primarily use the Newton iteration method to solve the nonlinear Schr?dinger equation and obtain soliton solutions. Once the steady-state solutions are determined, verifying their stability is crucial. We employ two methods to assess soliton stability: direct dynamic simulations and linear stability analysis. In linear stability analysis, perturbations are introduced to the steady-state solution, and eigenvalues are calculated by diagonalizing the Hamiltonian matrix. If any eigenvalue has a positive real part, the perturbation grows exponentially, indicating instability. Conversely, if all eigenvalues have non-positive real parts, the solution is considered stable. Soliton propagation is simulated using the stepwise Fourier method, which separates the Hamiltonian into linear and nonlinear components, each treated with different approaches. Finally, we compare transmission simulation results with those from the linear stability analysis.

We first investigate the system’s band structure and the Anderson modes in the band under linear conditions, as illustrated in Fig. 1. Anderson modes can randomly appear at any lattice position, and as the eigenvalue b increases, these modes become less localized. Under nonlinear conditions, we examine disordered solitons originating from Anderson modes and their P-b curves in a semi-infinite energy gap, as shown in Fig. 2. Figs. 2(a) and 2(c) depict the P-b curves corresponding to the first and fourth Anderson modes. Subgraphs in Fig. 2(c) demonstrate that resonant interactions between Anderson modes and disordered solitons cause deviations in the power curve. We then analyze disordered solitons originating from localized modes and their power curves in Fig. 3 and Fig. 5. Under cubic-quintic nonlinearity, the soliton family forms different branches as the power P increases, and the profile of the disordered solitons does not always widen with increasing power. Fig. 4 and Fig. 6 further examine the stability and dynamics of these soliton families through linear stability analysis and numerical transmission simulations, identifying the intervals where stable and unstable disordered solitons exist. Finally, Fig. 8(a) presents the disordered soliton family originating from the Anderson mode in the first bandgap, while Fig. 8(b) shows soliton profile variations with increasing power P.

We begin with Anderson modes and numerically investigate disordered solitons and multistability effects in cubic-quintic nonlinear media. Based on the band structure of disordered waveguide arrays, we find that localized Anderson modes bifurcate into soliton families that can be excited without a threshold and resonate with other Anderson modes when spatially overlapping, leading to soliton envelope broadening. In the semi-infinite gap, we analyze a high-order soliton family originating from localized modes, verifying its stability through linear stability analysis and numerical transmission simulations. Finally, we explore disordered solitons in the first bandgap under self-defocusing conditions, finding that as the beam power P increases, solitons are more likely to occur in both the first band and the first gap.