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

Complex and multiphysics engineering systems are frequently governed by interacting transport fields in which one variable modifies the flux, source, stability, or spatial organization of another. Such structures occur in reactive and chemically influenced transport, thermo-fluid and electrochemical systems, micro-scale transport, and other coupled-field models. Despite the diversity of applications, these systems share recurring numerical difficulties: nonlinear feedback, conservation constraints, stability transitions, strong localization, and the emergence of progressively finer spatial scales. Reliable computation therefore requires not only stable time integration but also transparent control of spatial resolution and independent verification of the interacting fields.

Although introduced in mathematical biology, Keller–Segel chemotaxis provides a particularly transparent representative of this broader class of nonlinear coupled-field transport problems. Cell density evolves by diffusion and by a directed flux driven by gradients of a second field, while the chemoattractant itself diffuses, is produced by the cell population, and decays. The resulting two-way feedback loop couples transport, reaction, and field-induced drift and can trigger instability, pattern formation, and sharp localization. In its biological setting, chemotaxis is the directed migration of cells or microorganisms in response to chemical cues and is involved in processes such as bacterial nutrient search, immune-cell recruitment, tissue organization, wound healing, and developmental cell migration. From a computational perspective, however, the same mathematical structure makes the Keller–Segel system a demanding benchmark for methodologies intended for complex coupled transport and multiphysics models.

The mathematical foundations of continuum chemotaxis modeling can be traced to the work of Patlak [1] and to the seminal models introduced by Keller and Segel [2], [3]. Since then, Keller–Segel-type systems have become a central framework for studying the competition between random cell diffusion and directed chemotactic drift. Reviews by Horstmann [4], [5] and Hillen and Painter [6] document the extensive mathematical development of these models and their many variants. Even in their minimal forms, Keller–Segel systems can exhibit a remarkably rich spectrum of behavior, including aggregation, spatial patterning, spike formation, traveling structures, stabilization, and, for suitable parameter ranges and model classes, finite-time concentration or blow-up phenomena [7], [8], [9].

These same nonlinear mechanisms make Keller–Segel systems demanding computational problems. The numerical representation of aggregation must simultaneously preserve global quantities, maintain physically \sloppy admissible solutions, and resolve progressively steeper spatial gradients. A broad range of discretization strategies has therefore been developed, including finite-volume schemes [10], positivity-preserving finite-volume/finite-difference formulations [11], [12], adaptive moving-mesh methods [13], and multidimensional difference-potential algorithms [14]. Accurate computation becomes particularly delicate when a localized peak narrows rapidly, since a numerically stable time integration does not by itself guarantee adequate spatial resolution. For blow-up-like solutions, the convergence of the computed concentration time is itself a nontrivial issue [15].

A central numerical question in strongly aggregative regimes is therefore not only whether a discrete solution can be advanced in time, but also whether the spatial representation remains resolved as increasingly high wavenumbers become dynamically relevant. This distinction is especially important in two dimensions, where localized concentration can transfer progressively more information toward the upper end of a spectral representation. Under such conditions, low-order observables may appear converged while peak height and width remain sensitive to spatial truncation. Reliable interpretation consequently requires complementary diagnostics that monitor physical-space localization, spectral-space resolution, and agreement with an independently converged numerical methodology.

Despite these advances, the numerical literature considered here remains dominated by physical-space discretizations in which spatial resolution is controlled through mesh refinement. In strongly localized Keller–Segel dynamics, this makes it difficult to monitor directly how nonlinear interactions redistribute content across spatial scales, while independent verification is commonly performed within closely related discretization frameworks. The methodological gap addressed here is therefore a multidimensional transform formulation that (i) treats the spatial means analytically, (ii) evaluates the two-dimensional nonlinear modal couplings exactly rather than by numerical quadrature, (iii) exposes a sparse interaction structure whose truncation can be monitored directly in spectral space, and (iv) is independently verified against a converged conservative physical-space method. To the authors’ knowledge, this combination of exact sparse Generalized Integral Transform Technique (GITT) coupling, systematic one- and two-dimensional modal convergence, independent finite-volume method (FVM)/Method of lines (MOL) verification, and high-resolution physical–spectral localization diagnostics has not previously been developed for the parabolic–parabolic Keller–Segel benchmark considered here.

Within this context, the GITT provides an alternative hybrid numerical–analytical framework. GITT combines eigenfunction expansions with numerical integration of the transformed ordinary differential equation (ODE) system, thereby retaining an analytical representation of the spatial dependence while controlling approximation through modal truncation. Recent applications by independent groups [16], [17] complement the broader computational–analytical GITT framework developed for linear and nonlinear transport problems [18], [19], [20]. This structure is attractive for Keller–Segel systems because spatial convergence can be examined directly in modal space, while nonlinear field coupling can be treated without discretizing the partial differential equation (PDE) itself on a conventional physical-space mesh. More generally, these features are valuable in coupled engineering transport models for which interaction mechanisms and resolution transfer between spatial scales must be monitored explicitly.

The present work treats the parabolic–parabolic Keller–Segel model both as a chemotaxis system and as a canonical nonlinear coupled-field transport problem within a broader multiphysics modeling context. A multidimensional GITT formulation with homogeneous Neumann boundary conditions is developed in which the spatially averaged components are separated and treated analytically, so that only zero-mean fluctuations are represented through the eigenfunction expansion. In two dimensions, a tensor-product Neumann basis is employed and nonlinear modal interactions are evaluated by exact cosine-product selection rules. This produces a highly sparse transformed system despite the large number of formally possible mode interactions. The nonlinear ODEs are integrated using variable-order backward differentiation formulas (BDF) through CVODE (a variable-step, variable-order ODE solver from SUNDIALS), with a matrix-free scaled preconditioned generalized minimal residual (SPGMR) linear solver for the multidimensional calculations. The formulation is therefore designed not only to compute the coupled fields but also to expose modal interaction, stability, scale transfer, and localization in a form suitable for systematic verification.

The numerical assessment is deliberately broader than a single high-sensitivity calculation. First, systematic one-and two-dimensional GITT convergence studies are carried out for $\chi=1,2,3,4,5$, spanning the linear threshold $\chi_c=2$ and several increasingly aggregative regimes. Convergence is evaluated through global norms, selected physical-space locations, modal-tail information, and representative profiles. Second, the GITT solutions are confronted with conservative finite-volume/method-of-lines formulations developed independently in physical space. These finite-volume solutions are themselves subjected to grid-refinement studies, allowing verification to be based on two separately converged spatial methodologies rather than on a single reference computation. Third, the resulting converged solutions are used for a quantitative parametric analysis of the effect of $\chi$ in both one and two dimensions, including the growth of genuinely multidimensional mixed harmonics. The independent FVM/MOL comparisons are restricted to the uniformly resolved verification windows ($t \leq 3$ in one dimension and $t \leq 2$ in two dimensions); the subsequent late-time concentration regime is assessed through systematic GITT truncation refinement and dedicated spectral-resolution diagnostics.

The strongly aggregative $\chi=5$ regime is then followed beyond the uniform parametric window through successive spectral truncations up to 7,920 spatial modes, corresponding to 15,840 coupled ODEs. An outer-shell spectral criterion defines a truncation-dependent resolution frontier. Physical and spectral localization are quantified independently through the full width at half maximum and a root-mean-square wavenumber scale, while normalized profiles test for self-similar collapse. These diagnostics reveal an approximately self-similar concentration regime and provide mutually consistent characteristic-time estimates near $t \approx 2.184$. The computations therefore furnish strong numerical evidence compatible with possible finite-time singularization, while they are not interpreted as a mathematical proof of blow-up of the continuous PDE system.

The remainder of the paper is organized as follows. Section 2 presents the Keller–Segel formulation, spatial-average equations, and zero-mean decomposition. Section 3 develops the one- and two-dimensional GITT formulations and sparse nonlinear modal coupling. Section 4 presents the temporal checks, systematic GITT convergence studies, and independent FVM/MOL verification. Section 5 develops the one- and two-dimensional parametric analysis and quantifies nonlinear multidimensional transfer. Section 6 examines the $\chi=5$ high-resolution concentration regime, including the resolution frontier, localization scales, profile collapse, and characteristic-time extrapolation. Section 7 summarizes the principal findings.

The formulation considered in this work is the parabolic–parabolic minimal Keller–Segel system [2], [3], [6]. Let $\Omega \subset \mathbb{R}^d$ be a bounded spatial domain, with d = 1 or 2 in the numerical studies presented below. The cell density $u(\mathbf{x}, t)$ and the chemoattractant concentration $v(\mathbf{x}, t)$ satisfy

where $d_u>0$ and $d_v>0$ are the cell and chemoattractant diffusivities, respectively, $\chi>0$ is the chemotactic sensitivity, $\alpha>0$ is the chemoattractant production coefficient, and $\beta>0$ is its decay coefficient. The coefficients in Eqs. (1)–(2) are retained explicitly so that the transformed formulation remains general. The numerical examples introduced below are normalized benchmark problems; no specific physical dimensionalization is required for the methodological and convergence analyses pursued here.

The initial conditions are prescribed as:

and homogeneous Neumann boundary conditions are imposed:

where, $\mathbf{n}$ denotes the outward unit normal to $\partial \Omega$. Because Eq. (6) eliminates the normal chemotactic drift and Eq. (5) eliminates the normal diffusive flux, these conditions define a closed system for the cell population. They also imply zero normal diffusive flux of the chemoattractant.

Eq. (1) describes the competition between random cell diffusion and directed chemotactic transport. Eq. (2) accounts for chemoattractant diffusion, production by the cell population, and linear degradation. Although the differential system is compact, the nonlinear transport term in Eq. (1) can generate sharply localized cell-density structures, which is the principal numerical difficulty examined in this work.

For any sufficiently regular scalar field $f(\mathbf{x}, t)$, define its spatial average by:

where, $|\Omega|$ is the measure of the spatial domain. Integrating Eq. (1) over $\Omega$ and applying the divergence theorem together with Eq. (5) and Eq. (6) gives:

Hence,

Likewise, integration of Eq. (2) yields:

whose explicit solution is:

with:

Therefore, the zero spatial mode of the cell density is exactly conserved, whereas the zero mode of the chemoattractant follows the scalar linear problem (Eq. (10)). These components need not be included in the numerically integrated transformed system.

Accordingly, the original fields are decomposed as:

where the fluctuation fields satisfy:

Substitution of Eq. (13) into Eq. (1)–Eq. (2), followed by use of Eq. (8) and Eq. (10), gives the governing equations for the zero-mean fields:

The corresponding initial conditions are:

and the boundary conditions retain the homogeneous Neumann form:

The decomposition (Eq. (13)) is useful for two reasons. First, it treats the spatial averages analytically and therefore preserves the exact mean balances independently of modal truncation. Second, because U and V have zero spatial mean, the constant Neumann eigenfunction can be removed from the numerical eigenfunction expansions without loss of information. This point is exploited explicitly in Section 3.

For the benchmark problems considered below, the model coefficients are normalized as:

while the chemotactic sensitivity $\chi$ is varied. All quantities reported for these numerical benchmarks are expressed in normalized form; no physical units are assigned to the benchmark variables or parameters. For these parameters and for the initial conditions adopted below:

so that Eq. (11) reduces to:

Consequently, Eqs. (15)–(16) become:

The normalized intervals are chosen as $[ 0, \pi]$ in each active coordinate direction. This normalization is especially convenient for the Neumann eigensystem because it yields cosine eigenfunctions with integer wavenumbers and eigenvalues, avoiding an otherwise repetitive factor of $\pi / L$ in the transformed equations. The choice is made solely for analytical and numerical convenience and does not restrict the general integral-transform formulation developed in Section 3.

For verification of the transformed formulation, the one-dimensional problem is posed on:

with

Hence,

with

This case is used only as a compact benchmark for modal convergence, conservation, and verification of the transformed equations.

The principal calculations are performed on the square domain:

with the symmetric initial condition:

The corresponding zero-mean fields initially satisfy:

Homogeneous Neumann conditions are imposed on all four boundaries:

The two-dimensional benchmark is symmetric with respect to interchange of $x$ and $y$. The parametric calculations reported in Section 5 consider $\chi=1,2,3,4,5$, including the linear threshold $\chi_c=2$, with particular emphasis on $\chi=5$, for which the solution develops the strongly localized concentration regime examined through the high-resolution spectral diagnostics of Section 6.

The zero-mean formulation developed in Section 2.2 is particularly convenient for an integral-transform treatment because the constant Neumann eigenfunction is represented exactly by the analytically known spatial averages and can therefore be excluded from the numerical modal system. The GITT solution is presented first in one dimension and then extended to the two-dimensional tensor-product basis used in the principal calculations.

For the normalized interval $\Omega_{1 D}=[ 0, \pi]$, consider the Neumann eigenvalue problem:

An orthonormal set of eigenfunctions is:

where, $\delta_{i j}$ denotes the Kronecker delta. The orthonormality relation is:

Because the fluctuation fields have zero spatial mean, the constant mode has identically zero transformed amplitudes. Since $X_0=1 / \sqrt{\pi}$, the transform/inverse pairs may therefore be written directly for $n \geq 1$:

Multiplying Eqs. (15) and (16) by $X_i(x)$, integrating over $[ 0, \pi]$, using integration by parts and the homogeneous Neumann boundary conditions, and then inserting Eqs. (38) and (40), yields:

where,

For $i, j, k \geq 1$, direct evaluation gives the exact selection rule:

Since $j \geq 1$, the first Kronecker delta vanishes automatically when $i=k$. Eq. (44) shows that only a small subset of the formally possible nonlinear interactions is nonzero.

For the normalized benchmark defined by Eq. (20) and $u_{av, 0}=1$, Eqs. (41)-(42) reduce to:

The initial condition Eq. (27) gives:

For diagnostics expressed in the unnormalized cosine representation, it is convenient to define the physical cosine amplitudes:

For the square domain $\Omega_{2 D}=[ 0, \pi]^2$, the one-dimensional normalized Neumann eigenfunctions are used in tensor-product form:

They satisfy:

Only the global constant mode ($m, n$)=(0, 0) is excluded. Modes with one zero coordinate index, such as ($m$, 0) or (0, $n$), must be retained because they represent legitimate spatial variation along one coordinate direction. Introducing the compound mode index:

the transform/inverse pairs become:

Here and below, the notation $I \neq 0$ is shorthand for exclusion of the single pair $(0,0)$.

Applying the integral transformation to Eqs. (15)-(16) gives, for every nonzero mode $I$,

The nonlinear coupling coefficient is:

For the normalized benchmark:

For the symmetric initial condition (Eq. (31)), all transformed cell-density fluctuations initially vanish:

whereas the chemoattractant field contains only the two axis modes:

The subsequent appearance of modes such as (1, 1), (2, 0), and higher harmonics is therefore generated entirely by the nonlinear chemotactic coupling.

For later modal diagnostics, the corresponding physical cosine amplitudes are:

so that $U=\sum a_{m n} \cos (m x) \cos (n y)$ and analogously for $V$, with the (0, 0) term excluded.

The coefficient in Eq. (61) can be reduced to a scalar triple-product integral. Using:

together with the homogeneous Neumann boundary conditions and:

one obtains:

where,

For

the tensor-product structure gives:

with the one-dimensional normalized cosine triple product:

Defining:

the product-to-sum identity yields the exact expression :

Because all indices are nonnegative integers, Eq. (75) is nonzero only when one index equals the sum of the other two, including the appropriate zero-index degeneracies, or when all three indices are zero. The same rule applies independently in the y direction.

Eqs. (69)–(75) eliminate the need for numerical quadrature of the nonlinear modal coefficients. More importantly, they reveal that the three-index interaction structure is extremely sparse: a coupling is retained only when the cosine selection rules are satisfied simultaneously in both coordinate directions. The implementation therefore constructs only admissible nonzero interactions rather than storing or traversing the full formal tensor.

For a square truncation in which $0 \leq m, n \leq N$ and only (0,0) is excluded, the number of retained spatial modes is:

and the coupled first-order system contains:

Accordingly, this quantity defines the dimension of the transformed ODE system. The same modal truncation is adopted for U and V in the present work. Convergence is assessed a posteriori by systematic increases in N and, in the strongly localized regime, by explicit monitoring of the high-wavenumber spectral shell as described in Section 6.

The transformed systems derived in Section 3 are nonlinear initial-value problems whose dimension is controlled by the modal truncation. The numerical strategy adopted here separates temporal accuracy, GITT spatial convergence, and independent physical-space verification. The one-dimensional benchmark is first used for analytical and temporal checks and for a systematic modal-convergence study. The same convergence architecture is then extended to the two-dimensional problem. Finally, both formulations are confronted with conservative finite-volume/method-of-lines (FVM/MOL) solutions obtained directly in physical space.

The modal ordinary differential equations (ODEs) are integrated with CVODE from the SUNDIALS library using variable-order BDF and adaptive temporal error control. The production tolerances are:

where $\Delta t_{\text {max }}$ is the maximum internal CVODE time step.

For the one-dimensional GITT benchmark, the transformed system is sufficiently small for a dense Newton/direct linear solve. For the multidimensional calculations, direct factorization becomes impractical as the number of retained modes increases; the production implementation therefore uses matrix-free SPGMR with maximum Krylov dimension 80 and at most 10 restarts. No preconditioner was used over the truncation ranges reported here. All nonlinear modal coefficients are generated once from the exact selection rules of Section 3.4, and only nonzero interactions are stored.

The independent FVM/MOL calculations use the same CVODE tolerances but a completely different spatial representation. The one-dimensional FVM system is solved with a banded direct linear solver, whereas the two-dimensional FVM system uses matrix-free SPGMR because the finest 161 × 161 grid generates 51,842 coupled ODEs. Table 1 summarizes the principal settings.

Linearization of the one-dimensional benchmark about the homogeneous state provides an analytical verification. A mode with squared wavenumber $q=n^2$ satisfies:

where, $a_n$ and $b_n$ are physical cosine amplitudes of $U$ and $V$. The growth rates are:

and the neutral condition $\sigma_{+}=0$ gives:

Since the lowest nonconstant Neumann mode has $q=1$, the first linear-instability threshold is:

This threshold concerns linear growth of perturbations and is not a blow-up criterion.

For the prescribed initial perturbation, $a_1(0)=0$ and $b_1(0)=\varepsilon$. Defining:

the corresponding linear cell-density amplitude is:

The transformed implementation reproduces this early-time solution before nonlinear modal transfer becomes significant.

Temporal accuracy was checked independently at $N=40$ by integrating the same nonlinear one-dimensional modal system with CVODE/BDF and with classical fourth-order Runge-Kutta using $\Delta t=10^{-4}$. For $\chi=5$ at $t=$ 3:

confirming that temporal errors are negligible on the spatial-convergence scales discussed below.

The systematic one-dimensional convergence study used $\chi$ = 1, 2, 3, 4, 5, the truncations:

and the comparison times $t$ = 0.5, 1.5, 3. The $N$ = 64 solution was used as reference. Relative errors are normalized by the fluctuation fields:

with analogous definitions for $V$.

Table 2 shows the most demanding case, $\chi$ = 5 at $t$ = 3. The cell-density error decreases by more than seven orders of magnitude between $N$ = 8 and $N$ = 48, while the chemoattractant field reaches the temporal-error floor earlier. The last-four-mode amplitude at $N$ = 64 is approximately $8.0 \times 10^{-16}$ for $U$, making this truncation a conservative reference over $0 \leq t \leq 3$. The complete one-dimensional modal-convergence trends for $\chi$ = 1 – 5 are shown in Figure 1.

The corresponding direct GITT-FVM/MOL profile comparison for $\chi$ = 5 at $t$ = 3 is shown in Figure 2.

The analytical zero-mode separation also provides a stringent conservation check. For every one-dimensional case and sampled time, the reconstructed means satisfy $u_{a v}=v_{a v}=1$ to approximately machine precision.

An independent physical-space verification was performed with a conservative FVM/MOL discretization of:

On a uniform grid, the cell-density flux at an interior face is approximated by:

with the diffusive chemoattractant flux treated analogously and zero boundary fluxes imposed directly. The nested sequence $N_x$ = 101, 201, 401, 801, 1601 was used, with $N_x$ = 1601 as the refined FVM solution.

Reference-free successive-grid comparisons gave clean second-order convergence; averaged over all $\chi$ and times, the observed orders were 2.00008 and 2.00007 for the $L_{\infty}$ and $L_2$ errors of $U$, respectively. For the hardest case, $\chi=5$ and $t=3$, the direct GITT-FVM comparison gives:

Table 3 shows that the independent-method agreement is not restricted to the most nonlinear case: both $U$ and $V$ remain closely matched over the complete parametric set. The $\chi=5$ case is the most demanding because the cell-density peak is both largest and narrowest.

Because the FVM discretization is demonstrably second order, Richardson extrapolation is applied as follows:

where the subscripts denote the number of finite-volume intervals $M$. This extrapolation reduces these relative $L_{\infty}$ discrepancies to $8.25 \times 10^{-9}$ for $U$ and $4.02 \times 10^{-9}$ for $V$. This reduction shows that the residual direct-method difference is dominated by the expected finite-volume spatial truncation error. Additional one-dimensional convergence times, pointwise comparisons, FVM self-convergence data, and detailed error distributions are provided in the Supplementary Material.

Two structural tests were applied before the multidimensional convergence study. First, a y-independent initial condition:

produced no modes with nonzero $y$ index, demonstrating exact dimensional reduction to the one-dimensional transformed system. Second, the symmetric two-dimensional initial condition preserved:

to roundoff while initially absent mixed modes such as (1,1) became nonzero. This simultaneously verifies coordinate-exchange symmetry and active multidimensional nonlinear transfer.

The systematic two-dimensional convergence study considered:

and the modal truncations:

with $N=40$ as reference. Global errors were evaluated on a common 101 × 101 reconstruction grid. Table 4 reports the discriminating $\chi=5, t=2$ case. The relative $L_{\infty}$ error of $U$ falls from $5.08 \times 10^{-2}$ at $N=8$ to $4.77 \times 10^{-7}$ at $N=32$, while the $N=40$ outer-shell fractions are only $6.02 \times 10^{-9}$ for $U$ and $1.69 \times 10^{-11}$ for $V$. The complete two-dimensional modal-convergence trends for $\chi=1-5$ are summarized in Figure 3.

The corresponding GITT–FVM/MOL comparison along the peak-containing cut $y$ = 0 is shown in Figure 4.

The independent two-dimensional FVM discretizes the conservative fluxes:

on nodal control volumes. Interior nodes use full square volumes, edge nodes use half volumes, and corner nodes use quarter volumes; normal boundary fluxes are set identically to zero. At each interior face, normal derivatives are evaluated by centered differences and the cell density multiplying the chemotactic gradient is the arithmetic average of the two adjacent nodal values, i.e., the direct two-dimensional analogue of Eq. (90). The nested grid sequence contains $21^2, 41^2, 81^2$, and $161^2$ nodes, corresponding at the finest level to 51,842 coupled ODEs.

Reference-free self-convergence gives mean observed orders $p$ = 2.0031 and 2.0013 for the $L_{\infty}$ and $L_2$ errors of $U$, and $p$ = 1.9999 and 1.9992 for $V$. For $\chi$ = 5, $t$ = 2, the finest triplet $M$ = 40,80,160 gives $p_{\infty, U}$ = 1.9987 and $p_{2, U}$ = 1.9995. All 20 FVM integrations terminate normally, with no linear or nonlinear convergence failures. For fieldwise comparison, the refined $M=160\,\, \mathrm{FVM}$ solution is bilinearly evaluated on the same $101 \times 101$ physical grid used to reconstruct the GITT solution; these two grids are not nested. On this common comparison grid, the direct $N$ = 40 GITT versus $M=160\,\, \mathrm{FVM}$ differences for $\chi$ = 5, $t$ = 2 are:

The complete parametric comparison at $t$ = 2 is summarized in Table 5. The discrepancies increase systematically with $\chi$, as expected from the sharpening cell-density structure, but remain small for both dependent variables even in the most demanding case.

The remaining discrepancy is consistent with the second-order FVM truncation. Because Richardson extrapolation on the non-nested 101 × 101 comparison grid would also contain interpolation error, the extrapolation is assessed instead at grid-aligned extrema, for which no interpolation is required. In the same case:

while Richardson extrapolation from $M$ = 80 and 160 yields:

reducing the absolute discrepancy from $3.48 \times 10^{-3}$ to $2.13 \times 10^{-6}$. Similarly, the Richardson errors are $3.73 \times 10^{-10}$ for $u_{\text {min}}$ and $1.67 \times 10^{-7}$ for $v_{\text {max}}$. The close agreement between the independently converged spectral and finite-volume solutions establishes a strong numerical verification of the two-dimensional GITT formulation.

The complete multi- $\chi$ convergence sequences at additional times, pointwise comparisons, diagonal cuts, FVM self-convergence results, and field-difference maps are collected in the Supplementary Material. These verification results provide the numerical basis for the parametric study of Section 5 and the high-resolution analysis of Section 6.

The convergence studies of Section 4 establish $N$ = 64 as a conservative one-dimensional production truncation over $0 \leq t \leq 3$ and $N$ = 40 as a conservative two-dimensional reference over $0 \leq t \leq 2$. These uniformly resolved windows are used here to examine the dependence of the solution on:

including the linear threshold $\chi_c$ = 2 identified in Eq. (82). The high-resolution continuation of the most aggregative case, $\chi$ = 5, is treated separately in Section 6.

The one-dimensional profiles exhibit a systematic transition from relaxation toward the homogeneous state to increasingly strong aggregation. Table 6 summarizes the field extrema at $t$ = 3. For $\chi$ = 1, the solution remains close to the unit mean. At the threshold $\chi$ = 2 , the disturbance remains modest but no longer exhibits the same decay tendency. For $\chi$ = 3 and 4, the contrast between peak and depleted region grows progressively, and for $\chi$ = 5 the cell-density maximum reaches 4.25956 while the minimum falls to 0.382827.

Figure 5 shows both dependent variables at $t=3$. Increasing $\chi$ amplifies the cell-density peak near $x=0$ and deepens the depleted region toward $x=\pi$ while preserving the spatial mean. The chemoattractant responds over a broader spatial scale and remains substantially smoother than $u$, although its amplitude also increases systematically with $\chi$. Additional times are reported in the Supplementary Material.

The same progression is amplified in two dimensions because the nonlinear chemotactic flux acts simultaneously along both coordinates. The symmetric initial condition implies:

and the maximum remains at the Neumann corner (0,0) throughout the resolved parametric interval. Table 7 summarizes the two-dimensional extrema at $t$ = 2.

Figure 6 shows both $u(x, 0, t)$ and $v(x, 0, t)$ at $t=2$ along the cut that passes directly through the corner peak. The increase from $\chi=4$ to 5 is particularly pronounced in $u$, indicating the onset of the rapidly concentrating regime that is later followed with much higher spectral resolution, whereas $v$ retains a broader and smoother spatial distribution. Diagonal cuts, additional contour maps, and the complete set of data-driven three-dimensional surfaces are included in the Supplementary Material.

The two-dimensional character of the solution is quantified directly by modes that are absent initially. For $\chi$ = 5, Table 8 lists the leading axis mode $a_{10}$, the second axis harmonic $a_{20}$, and the mixed mode $a_{11}$. The mixed contribution grows much faster than a perturbative correction: at $t$ = 1.5, $a_{11} / a_{10} \simeq 0.69$, and by $t$ = 2:

so the mixed (1, 1) amplitude has exceeded the leading axis amplitude. This provides a direct spectral measure of genuinely multidimensional nonlinear transfer.

A representative high-resolution snapshot beyond the uniform parametric window is shown in Figure 7. The snapshot corresponds to $\chi$ = 5, $N$ = 72, and $t$ = 2.164, still below the operational resolution frontier $t_{\star}(72) \simeq 2.1651$. Both fields attain their maxima at $(x, y)$ = (0, 0). The cell density is already sharply concentrated:

whereas the chemoattractant remains substantially smoother:

The progression from the uniformly resolved parametric interval to the $t=2.164$ snapshot makes clear that peak amplitude alone is insufficient to establish spatial resolution. The peak height increases while its width contracts and spectral content migrates toward higher wavenumbers. Section 6 therefore treats the $\chi=5$ late-time regime separately, using successive truncations, outer-shell diagnostics, physical localization scales, and normalized profile collapse.

The parametric calculations of Section 5 show that the $\chi$ = 5 solution enters a rapidly concentrating regime after the initially well-resolved interval. The numerical question then changes qualitatively. It is no longer sufficient to verify convergence of a few low modes or successful completion of the time integrator: the narrowing cell-density structure requires progressively higher spatial wavenumbers. This section therefore combines modal convergence, an outer-shell resolution criterion, physical- and spectral-space localization measures, and scaled-profile analysis. The purpose is to characterize the concentration process without identifying numerical loss of resolution with blowup of the continuous PDE.

For a square truncation $0 \leq m, n \leq N$, excluding only the analytically separated (0,0) mode, the number of retained spatial modes is:

The high-resolution sequence was extended up to $N$ = 88, corresponding to $N_s$ = 7920 spatial modes and 15840 coupled transformed ODEs. At this truncation, the analytically selected sparse nonlinear interaction list contains $1.3795 \times 10^8$ nonzero couplings.

To detect when the truncation boundary begins to carry a non-negligible fraction of the cell-density fluctuation, define the outer shell:

and the modal-norm fraction:

where, $a_{m n}$ is the physical cosine amplitude defined in Eq. (66). The shell fraction is therefore an amplitude-based operational diagnostic; unlike the spectral moment introduced below, it is not a Parseval norm in the orthonormal transform coefficients. The operational resolution criterion adopted here is:

This threshold is not a mathematical error bound. It is a conservative diagnostic indicating that the physical cosine amplitudes on the outer retained shell are still dynamically small relative to the complete retained amplitude set.

Figure 8 illustrates the behavior for $N$ = 80,84, and 88 . At $N$ = 88, the shell fraction increases from $2.71 \times 10^{-4}$ at $t$ = 2.166 to $5.16 \times 10^{-4}$ at $t=$ 2.168 and $1.002 \times 10^{-3}$ at $t$ = 2.170. Thus, even as the concentration accelerates, increasing $N$ moves the threshold-crossing point systematically forward in time.

A direct field-convergence comparison provides a complementary measure. Using $N=88$ as the reference solution, the $N$ = 84 relative $L_{\infty}$ error of the cell-density fluctuation is $4.21 \times 10^{-4}$ at $t=2.166$ and $7.37 \times 10^{-4}$ at $t=$ 2.168, increasing to $1.30 \times 10^{-3}$ at $t=2.170$. In contrast, the leading modal amplitudes remain converged to approximately $10^{-10}-10^{-9}$ even at later times. This separation is important: convergence of the low modes alone does not guarantee convergence of a sharply localized physical-space peak.

The temporal solver itself remains robust throughout this sequence. The $N$ = 88 calculation reaches $t$ = 2.174 with 648 accepted internal steps and 4312 Krylov iterations, with no nonlinear or linear convergence failures. The observed restriction near the concentration regime is therefore spectral rather than a breakdown of the BDF-SPGMR time integrator.

For each truncation, define $t_{\star}(N)$ as the time at which $F_A^{\text {shell }}=10^{-3}$. When the threshold lies between two sampled times $\left(t_1, F_1\right)$ and $\left(t_2, F_2\right)$, logarithmic interpolation is used:

Table 9 consolidates the complete high-resolution sequence.

The monotonic increase of $t_{\star}$ with $N$ shows that the apparent numerical frontier is controlled by spectral resolution. Fitting the empirical form:

to the complete $N=40-88$ sequence gives:

The limiting value $t_{\infty}$ is an extrapolation of the numerical resolution frontier, it is not, by itself, an estimate or proof of the blow-up time of the continuous PDE. The truncation-dependent resolution frontier and its empirical fit are shown in Figure 9.

The physical concentration width is measured relative to the excess above the conserved mean. Define:

The peak occurs at the Neumann corner $(x, y)=(0,0)$. Therefore, the physical width is evaluated by even reflection across the boundary: full width at half maximum $(\mathrm{FWHM})_x$ is twice the distance from the corner to the half-excess crossing along $y=0$, and $\mathrm{FWHM}_y$ is defined analogously. The geometric measure:

is used below. Symmetry gives $\mathrm{FWHM}_x=\mathrm{FWHM}_y$ to numerical precision.

An independent spectral localization scale is obtained directly from the normalized transform coefficients:

As the peak narrows, $k_{\text {rms}}$ increases and $\ell_{\text {spec}}$ decreases.

The $N=88$ values illustrate the contraction clearly. Between $t$ = 2.166 and $t$ = 2.168, both still below the adopted $N$ = 88 shell threshold, $u_{\max }$ increases from 182.98 to $213.78, \mathrm{FWHM}_{\text {geom }}$ decreases from 0.14338 to 0.13158 , and $\ell_{\text {spec }}$ decreases from 0.07741 to 0.07085 . The later points continue the same trend but progressively cross the spectral-resolution criterion. Figures~\ref{fig10} and \ref{fig11} show the systematic behavior with increasing truncation.

A broader, densely sampled $N$ = 80 calculation over the resolved interval $2.140 \leq t \leq 2.166$ provides a particularly useful comparison between the two length scales. Over this interval, the peak excess:

increases from 56.47 to 181.79 , whereas $\mathrm{FWHM}_{\text {geom }}$ decreases from 0.27412 to 0.14337 and $\ell_{\text {spec }}$ from 0.15097 to 0.07741 . At the same time,

a variation of only about 2\%. Thus, the physical-space and spectral-space localization measures contract nearly proportionally.

To examine whether the localization is accompanied by an approximately invariant shape, the boundary-axis profiles are rescaled according to:

By construction, $\widehat{U}(0, t)=1$ and $\widehat{U}(1, t)=1 / 2$.

Figure 12 shows the profiles at fourteen resolved times from $t$ = 2.140 to $t$ = 2.166. Despite the more than threefold increase of the peak excess and nearly twofold contraction of the width, the rescaled profiles almost collapse onto a single curve. Over $0 \leq \xi \leq 4$, the maximum spread among these normalized profiles is only $5.612 \times 10^{-3}$. This constitutes a direct shape-based indication that the concentration is approaching an approximately self-similar regime rather than merely developing an arbitrary narrow spike.

The collapse is not exact. A local effective amplitude–width exponent:

varies over the resolved window from approximately 1.78 to 1.82 for the physical width and from 1.74 to 1.76 for the spectral scale. Moreover, the diagnostic product $A \mathrm{FWHM}_{\text {geom}}^2$ decreases from 4.24 to 3.74. The observed regime is therefore more accurately described as approximately self-similar than as an exact mass-preserving similarity solution with a fixed exponent.

The simultaneous growth of $A(t)$, contraction of two independent localization scales, profile collapse, and systematic displacement of the spectral-resolution frontier strongly suggest a singularizing concentration process. To test whether these independent diagnostics are compatible with a common finite characteristic time, the resolved $N=80$ interval $2.140 \leq t \leq 2.166$ was fitted to the relations below. Here, $p$, $q$, and $r$ denote fitted exponents and should not be confused with the modal indices used in Section 3:

Only points satisfying Eq. (110) were included. The resulting diagnostic estimates are summarized in Table 10, together with the independent extrapolation of the resolution frontier.

All four independent extrapolations fall within the narrow interval:

which motivates the characteristic scale:

The agreement is noteworthy because the four diagnostics probe different aspects of the calculation: peak amplitude, a physical-space width, a Parseval-consistent spectral length, and the truncation-dependent loss-of-resolution frontier.

Nevertheless, Eq. (125) must be interpreted with care. The computations establish a pronounced blow-up-like concentration regime and provide strong numerical evidence compatible with finite-time singularization, but they do not prove that the solution of the continuous parabolic-parabolic Keller-Segel problem becomes unbounded at that time. A mathematical blow-up statement would require analysis beyond numerical extrapolation, and the fitted exponents may still exhibit pre-asymptotic drift as the concentration scale contracts. The present conclusion is therefore restricted to the numerical evidence: the resolved solution approaches an approximately self-similar, rapidly localizing regime whose independent diagnostic time scales consistently cluster near $t \approx 2.184$.

This work developed a multidimensional GITT formulation for a parabolic–parabolic Keller–Segel system, with analytical separation of the spatial means and exact sparse cosine-product coupling in two dimensions. The resulting modal system preserves the mean balances explicitly, exposes nonlinear transfer among retained modes, and provides direct control of spatial truncation within a hybrid numerical–analytical framework.

Systematic one- and two-dimensional convergence studies, together with independent conservative FVM/MOL calculations, verified the transformed solutions over the uniformly resolved comparison windows. The finite-volume discretizations exhibited clean second-order self-convergence, and close agreement with GITT was obtained for both the cell-density and chemoattractant fields. Across $\chi=1-5$, the parametric analysis captured the transition from weak relaxation to strong aggregation and quantified the growth of genuinely two-dimensional mixed modes.

For $\chi=5$, successive refinements up to 7,920 spatial modes allowed the rapidly concentrating regime to be followed beyond the uniform parametric window. Physical width, spectral localization, outer-shell resolution diagnostics, and normalized-profile collapse consistently indicated an approximately self-similar concentrating state, with independent extrapolations clustering near $t \approx 2.184$. These observations are compatible with possible finite-time singularization but do not constitute a mathematical proof of blow-up; the late-time estimates remain numerical and may still contain pre-asymptotic drift.

Taken together, the results show that GITT can serve both as a high-accuracy solver and as a spectral diagnostic framework for nonlinear coupled-field transport. The combination of exact modal coupling, independent physical-space verification, and complementary localization measures provides a transferable computational strategy for other complex transport systems. Beyond chemotaxis, the mathematical structure investigated here, involving cross-gradient transport, coupled source and decay mechanisms, stability transitions, and progressive localization, is shared by a broader class of complex multiphysics engineering systems. Future developments may address adaptive modal enrichment, preconditioning at larger transformed dimensions, and extensions to Keller–Segel variants or additional interacting fields.