A Self-Consistent Wave-Optics Workflow for Graphene Nonlinearity in a D-Shaped Fiber

Anoma D. McCoy1, Morten Ibsen1
1Optoelectronics Research Centre, University of Southampton, Southampton, UK
Published in 2026

Graphene-coated optical fibers offer a compact route to fiber-based nonlinear switching by combining guided-wave interaction with graphene’s intensity-dependent optical absorption. However, early-phase design of these devices is poorly suited to full 3D optical-frequency simulation: subwavelength meshing makes repeated design sweeps computationally expensive when the goal is identifying trends rather than full 3D device-level fidelity. In this work, we present a computationally efficient 2D COMSOL Multiphysics workflow using LiveLink for MATLAB to estimate the nonlinear modal response of a graphene-coated D-shaped fiber. The modeled structure is a side-polished fiber with a graphene-Ta2O5-graphene stack on the polished flat. Each atomically thin graphene layer is modeled using COMSOL's Transition Boundary Condition (TBC), which captures the material response through a surface conductivity, avoiding explicit volume meshing while retaining its electromagnetic effect. The baseline linear graphene surface conductivity is calculated using the Kubo formulation, while the nonlinear response includes Kerr-based refractive-index changes and saturable absorption driven by the local optical intensity at the graphene surface. To avoid a full 3D nonlinear wave-propagation model, we use a 2D Electromagnetic Waves, Frequency Domain eigenmode study in COMSOL's Wave Optics Module, embedded in a MATLAB self-consistent iteration loop. The workflow starts from a Kubo-only eigenmode solution that defines the baseline guided mode and graphene-field interaction. For each input optical power, MATLAB uses the extracted graphene-surface intensity and absorbance to compute an updated graphene surface conductivity by adding Kerr and saturable-absorption corrections to the Kubo baseline. This updated material state is then pushed back into COMSOL, which re-solves the eigenmode problem with the updated conductivity, and the loop is repeated until the complex modal effective index, neff, converges. For this study the modeled segment length is 250 um, over which the graphene material state is treated as locally uniform; longer devices are represented by cascaded segments, with optical power updated from one segment to the next. To cross-check model consistency, we implemented the nonlinear material update in two ways: an effective nonlinear-index pathway, nNL, and a nonlinear surface-conductivity pathway, sigmaNL. The nNL pathway adds Kerr and saturation effects to the Kubo-derived graphene index before converting the result to the surface conductivity used by COMSOL. The sigmaNL pathway applies the nonlinear correction directly at the conductivity level. Both pathways ultimately update the same boundary-condition parameter in COMSOL: the nonlinear graphene surface conductivity used by the TBC. For the present study, the self-consistent solution converged within three iterations for a tolerance of 5 x 10^-6, requiring ~120-130 seconds per optical-power point versus the multi-day runtimes expected for comparable full 3D optical models. The nNL and sigmaNL pathways produced overlapping trends in both the nonlinear index shift and the effective modal attenuation, as shown in Figure 4, confirming that the two formulations are equivalent routes to the same physical result. As expected for graphene saturable absorption, increasing normalized graphene-surface intensity produces a saturating index shift and decreasing modal attenuation. This workflow provides a fast design-screening method for identifying graphene-induced nonlinear modal trends before committing to more computationally intensive 3D simulations.