Generalized Integral Transform Modeling of Complex Coupled-Field Transport in Two-Dimensional Keller–Segel Chemotaxis: Verification, Spectral Convergence, and Self-Similar Localization

Complex engineering systems often involve interacting transport fields in which diffusion, field-driven fluxes, production or decay, and nonlinear feedback jointly control stability, transport, and localization.Although originating in mathematical biology, the Keller-Segel system is a canonical coupled-field transport model with this structure and provides a stringent benchmark for computational methodologies aimed at complex and multiphysics transport problems.A multidimensional hybrid numerical-analytical formulation based on the Generalized Integral Transform Technique (GITT) is developed for a parabolic-parabolic Keller-Segel system with homogeneous Neumann boundary conditions.Spatial averages are treated analytically, zero-mean fields are expanded in normalized Neumann eigenfunctions, and exact cosine-product selection rules generate a sparse nonlinear modal system.The transformed equations are advanced using CVODE, a variable-step, variable-order ordinary differential equation (ODE) solver from the SUNDIALS (Suite of Nonlinear and Differential/Algebraic Equation Solvers) package, with variable-order backward differentiation formulas (BDF), together with a matrix-free scaled preconditioned generalized minimal residual (SPGMR) solver in two dimensions.Systematic one-and two-dimensional convergence studies are performed for chemotactic sensitivities χ = 1-5 and are complemented by independent conservative finite-volume/method-of-lines solutions in physical space.The finite-volume discretizations exhibit second-order self-convergence and closely reproduce the converged GITT fields.The parametric study quantifies the transition from weak relaxation to strong aggregation and the growth of genuinely two-dimensional mixed modes.For χ = 5, successive truncations up to 7,920 spatial modes (15,840 coupled ODEs) resolve a strongly localized regime.Physical and spectral localization scales contract consistently, while normalized profiles collapse onto an approximately invariant shape.Independent late-time extrapolations yield characteristic times clustered near t ≈ 2.184.The results demonstrate a high-resolution hybrid framework for nonlinear coupled-field partial differential equation (PDE) systems with cross-field transport, evolving localization, and demanding multiscale resolution, while providing numerical evidence compatible with possible finite-time singularization without constituting a mathematical proof of blow-up.