Rationally Enriched Chebyshev Trunk Bases for DeepONet Surrogates of High P\'eclet Entrance Transport

# Rationally Enriched Chebyshev Trunk Bases for DeepONet Surrogates of High Péclet Entrance Transport
Source: [https://arxiv.org/html/2608.19658](https://arxiv.org/html/2608.19658)
Mingeun ChoiEmail:[mingeun\.choi@gatech\.edu](mailto:[email protected])Affiliation:George W\. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, 30332, GA, USASatish KumarEmail:[satish\.kumar@me\.gatech\.edu](mailto:[email protected])Corresponding author:Corresponding author\.Affiliation:George W\. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, 30332, GA, USA

###### Abstract

This study demonstrates a rationally enriched Chebyshev \(REC\) trunk for deep operator network \(DeepONet\) surrogate models of singularly perturbed and high\-Péclet transport problems whose solution profiles are characterized by thin localized boundary or wall layers\. The REC trunk combines Chebyshev polynomial dictionary elements with rational dictionary elements constructed using the adaptive Antoulas–Anderson \(AAA\) algorithm\. Over five independent training runs, the resulting REC\-trunk DeepONet is evaluated against a vanilla DeepONet and a Chebyshev\-trunk DeepONet whose prescribed dictionary consists only of Chebyshev polynomials across three problems whose singular perturbation parameters are diffusion\-to\-advection ratios: a singularly perturbed scalar boundary\-value problem \(BVP\), the thermal entrance problem with a prescribed wall temperature, and the concentration entrance problem with an absorbing wall\. Across the held\-out test profiles, the REC\-trunk DeepONet improves over the vanilla DeepONet and remains comparable to the Chebyshev\-trunk DeepONet in predicting the scalar profile, with its clearest advantage over the Chebyshev\-trunk DeepONet appearing when the perturbation parameter lies between1\.00×10−41\.00\\times 10^\{\-4\}and1\.78×10−41\.78\\times 10^\{\-4\}, where it reduces the profile\-error metrics by up to19\.5%19\.5\\,\\%relative to the Chebyshev\-trunk DeepONet\. In predicting the wall\-normal temperature and concentration profiles, the REC\-trunk DeepONet reduces the profile\-error metrics by up to60\.2%60\.2\\,\\%and32\.2%32\.2\\,\\%relative to the vanilla and Chebyshev\-trunk DeepONets, respectively, while suppressing artificial near\-wall oscillations as the Péclet or mass\-transfer Péclet number ranges from10210^\{2\}to10410^\{4\}\.

###### Keywords:

DeepONet , Rationally enriched Chebyshev trunk , High\-Péclet transport

## 1Introduction

Transport problems are formulated as partial differential equations \(PDEs\), in which advection, diffusion, and reaction can interact across disparate spatial and temporal scales[courant2008methods](https://arxiv.org/html/2608.19658#bib.bib1)\. When diffusion is weak relative to advection, the associated boundary\-value problem \(BVP\) can become singularly perturbed because the small perturbation parameter \(or diffusion\-to\-advection ratio\) multiplies the highest spatial derivative\. Such problems often develop boundary or interior layers whose widths are far smaller than the domain length[roos2008robust](https://arxiv.org/html/2608.19658#bib.bib2)\. In resolving these localized layers, standard numerical schemes suffer from spurious oscillations or require a prohibitively dense grid[roos2008robust](https://arxiv.org/html/2608.19658#bib.bib2)\. Thus, accurate approximation of solution profiles requires techniques such as fitted finite differences, stabilized finite element methods, and layer\-adapted meshes[roos2008robust](https://arxiv.org/html/2608.19658#bib.bib2), while repeated layer\-resolving simulations across parameter ranges remain computationally expensive\. To reduce this burden, recent machine\-learning \(ML\) studies have developed physics\-informed neural\-network \(PINN\)\-based models for singularly perturbed problems by incorporating asymptotic decompositions, stretched variables, or parameter continuation[arzani2023theory](https://arxiv.org/html/2608.19658#bib.bib3),[zhang2024multi](https://arxiv.org/html/2608.19658#bib.bib4),[cao2023physics](https://arxiv.org/html/2608.19658#bib.bib5)\. However, these models rely on problem\-specific structures chosen from the differential equation, the small parameter, or the layer behavior of the target problem\.

Deep operator networks \(DeepONets\) have emerged as alternative surrogates for singularly perturbed problems because they are designed to learn operators mapping sampled input functions and problem parameters to solution profiles through a branch–trunk representation[lu2021learning](https://arxiv.org/html/2608.19658#bib.bib6),[wang2021learning](https://arxiv.org/html/2608.19658#bib.bib7)\. A pioneering application showed that a DeepONet can approximate sharp\-gradient solutions by evaluating the training loss at layer\-adapted Shishkin points, which depend on the perturbation parameter, rather than by modifying the network representation itself[du2024approximation](https://arxiv.org/html/2608.19658#bib.bib8)\. Subsequent studies then incorporated such layer\-resolving priors within the operator\-learning representation\. Prandtl–Van Dyke DeepONet \(PVD\-ONet\) decomposes the solution into outer, inner, and matching components and assigns these components to multiple DeepONet modules organized by Prandtl and Van Dyke matching[sun2025pvd](https://arxiv.org/html/2608.19658#bib.bib9)\. The enriched finite element operator network \(eFEONet\) augments the finite\-element Galerkin ansatz with singular\-perturbation corrector functions\. The learned coefficients include both nodal finite\-element coefficients and coefficients of the added layer\-corrector basis[lee2025efeo](https://arxiv.org/html/2608.19658#bib.bib10)\. Physics\-informed adaptive\-scale DeepONet \(PAS\-Net\) augments the trunk input with prescribed or learnable locally rescaled coordinates centered at reference points, thereby changing the coordinate representation seen by the trunk[mou2025pas](https://arxiv.org/html/2608.19658#bib.bib11)\. However, these approaches remain problem\-specific, since the solution decomposition, finite element enrichment, or coordinate rescaling must be chosen for the particular layer being represented\.

Parallel efforts have sought to reduce this dependence by expressing the output through prescribed or data\-derived bases and their coefficients\. Proper orthogonal decomposition DeepONet \(POD\-DeepONet\) computes POD modes from the training outputs and uses the resulting modes as the trunk, leaving the branch network to predict modal coefficients[lu2022fair](https://arxiv.org/html/2608.19658#bib.bib12)\. Spectral coefficient learning via operator network \(SCLON\) predicts coefficients in orthogonal expansions such as Fourier or Legendre bases and applies this coefficient\-space formulation to parametric PDEs ranging from singularly perturbed convection\-diffusion equations to Navier–Stokes flows[choi2023spectral](https://arxiv.org/html/2608.19658#bib.bib13)\. The orthogonal polynomial neural operator \(OPNO\) builds neural operators around orthogonal\-polynomial representations on bounded domains and treats Dirichlet, Neumann, and Robin boundary conditions within that polynomial framework[liu2024render](https://arxiv.org/html/2608.19658#bib.bib14)\. Spectral\-embedded DeepONet \(SEDONet\) converts raw coordinate trunk inputs into values of Chebyshev polynomials before a trainable trunk network, improving bounded\-domain DeepONet approximation for sharp gradients, boundary layers, and nonperiodic structures[abid2025sedonet](https://arxiv.org/html/2608.19658#bib.bib15)\. Nevertheless, these spectral and polynomial bases are global representations over the entire domain, and, thus, a thin localized layer can require many degrees of freedom unless the basis is enriched with dictionary elements adapted to the inner scale\.

This inner\-scale representation issue extends beyond idealized scalar BVPs because engineering transport often contains analogous high\-Péclet convection\-diffusion layer structure\. In biomedical engineering, surface\-based biosensors and related biofluidic capture systems involve convective delivery, diffusion, and surface binding near reactive walls[squires2008making](https://arxiv.org/html/2608.19658#bib.bib16)\. In chemical engineering, reactive microchannels and electrochemical chips generate wall\-normal concentration fields governed by convective\-diffusive delivery to reactive interfaces[gervais2006mass](https://arxiv.org/html/2608.19658#bib.bib17),[chevalier2021semianalytical](https://arxiv.org/html/2608.19658#bib.bib18)\. In aerospace engineering, high\-speed wall\-bounded flows involve wall\-normal heat and mass transfer in cooled and transpiration\-cooled boundary layers[sescu2019transpiration](https://arxiv.org/html/2608.19658#bib.bib19),[hillcoat2025transpiration](https://arxiv.org/html/2608.19658#bib.bib20), hypersonic boundary layers with strong near\-wall thermal gradients[xu2022hypersonic](https://arxiv.org/html/2608.19658#bib.bib21), and reacting boundary layers with finite\-rate wall chemistry[passiatore2021finite](https://arxiv.org/html/2608.19658#bib.bib22),[perakis2021recombination](https://arxiv.org/html/2608.19658#bib.bib23)\. Across these examples, the common mathematical feature is a wall\-normal temperature or concentration profile shaped by strong axial transport, transverse diffusion, and, in reactive cases, surface kinetics or wall absorption\. Thus, entrance heat and mass transfer provide benchmark problems that retain the relevant wall\-normal layer geometry under prescribed\-wall, absorbing\-wall, or Robin\-type surface conditions[shah1978laminar](https://arxiv.org/html/2608.19658#bib.bib24),[haase2015graetz](https://arxiv.org/html/2608.19658#bib.bib25),[popel1978mass](https://arxiv.org/html/2608.19658#bib.bib26),[debarnot2018graetz](https://arxiv.org/html/2608.19658#bib.bib27),[aquino2024equilibrium](https://arxiv.org/html/2608.19658#bib.bib28)while allowing accurate numerical reference profiles and controlled variation of Péclet number, axial location, inlet profile, and wall absorption strength\. However, the use of inner\-scale\-enriched DeepONet trunk bases for such high\-Péclet wall\-normal profile reconstruction remains underexplored\.

This paper introduces the rationally enriched Chebyshev \(REC\) trunk as a prescribed output\-coordinate dictionary for DeepONet that combines Chebyshev polynomial dictionary elements with rational dictionary elements without adding trainable trunk parameters\. This REC\-trunk DeepONet is compared over five independent training runs with a vanilla DeepONet and a Chebyshev\-trunk DeepONet, whose prescribed trunk dictionary consists only of Chebyshev polynomial dictionary elements, across three bounded\-domain benchmark problems whose solution profiles contain one\-sided localized layers\. The first problem is the singularly perturbed scalar BVP, serving as the simplest bounded\-interval formulation with an exponentially thin outflow layer[du2024approximation](https://arxiv.org/html/2608.19658#bib.bib8)\. The second problem is the thermal entrance problem with a prescribed wall temperature, where the temperature profile develops under hydrodynamically developed laminar duct flow[shah1978laminar](https://arxiv.org/html/2608.19658#bib.bib24),[haase2015graetz](https://arxiv.org/html/2608.19658#bib.bib25)\. The third problem is the concentration entrance problem with an absorbing wall, where the concentration profile develops under the same hydrodynamically developed laminar duct flow[popel1978mass](https://arxiv.org/html/2608.19658#bib.bib26),[debarnot2018graetz](https://arxiv.org/html/2608.19658#bib.bib27)\. Results demonstrate that DeepONet with the REC trunk remains globally competitive on the singularly perturbed scalar BVP and delivers its clearest advantage on the thermal and concentration entrance problems, where wall\-attached layers control the profile geometry, by answering the following four questions:

1. 1\.How does the REC trunk perform on the singularly perturbed scalar BVP in the small\-perturbation\-parameter regime?
2. 2\.How does the REC trunk perform on the thermal entrance problem in the high\-Péclet regime, when the learned output is a wall\-normal profile at a queried axial location?
3. 3\.How does the REC trunk perform on the concentration entrance problem in the high\-Péclet regime, under Robin wall kinetics and across different Damköhler numbers?
4. 4\.How reproducible are the gains obtained with the REC trunk across repeated training runs and previously unseen solution profiles?

## 2Methods

### 2\.1Benchmark problems and operator\-learning formulation

Figure 1:Computational domains and boundary conditions for \(a\) Problem 1 and \(b\) Problems 2 and 3\. The full channel shown in \(b\) is the symmetric extension of the half\-channel computational domainy∈\[0​;​1\]y\\in\[0\\mathord\{\\mathchar 59\\relax\}1\]\.Figure[1](https://arxiv.org/html/2608.19658#S2.F1)illustrates the computational domains and boundary conditions for the three benchmark problems considered in this study\.Problem 1examines the singularly perturbed scalar BVP posed on the one\-dimensional domain shown in Fig\.[1](https://arxiv.org/html/2608.19658#S2.F1)\(a\) and defined by the equation[du2024approximation](https://arxiv.org/html/2608.19658#bib.bib8)

−ε​d2​ud​x2​\(x\)\+d​ud​x​\(x\)\+u⁡\(x\)=q⁡\(x\);x∈\(0​;​1\);\-\\varepsilon\\frac\{d^\{2\}u\}\{dx^\{2\}\}\(x\)\+\\frac\{du\}\{dx\}\(x\)\+u\(x\)=q\(x\)\\mathchar 59\\relax\\qquad x\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)\\mathchar 59\\relax\(1\)with homogeneous Dirichlet boundary conditions

u⁡\(0\)=u⁡\(1\)=0\.u\(0\)=u\(1\)=0\.\(2\)
Here,xxdenotes the dimensionless spatial coordinate on the bounded interval\(0​;​1\)\(0\\mathord\{\\mathchar 59\\relax\}1\),u⁡\(x\)u\(x\)denotes the dimensionless scalar solution field,q⁡\(x\)q\(x\)denotes the dimensionless randomized source term, andε\\varepsilondenotes the singular perturbation parameter, given by the diffusion\-to\-advection ratio after the constant positive convective velocity is normalized to unity\. The reaction coefficient is also fixed at unity\. These choices allow the effects ofε\\varepsilonandq⁡\(x\)q\(x\)on the solution profiles to be examined without simultaneous changes in convection or reaction\. Sinceε\\varepsilonmultiplies the highest\-order diffusion term, decreasing its value weakens diffusion relative to convection, producing a thinner right\-endpoint exponential layer[roos2008robust](https://arxiv.org/html/2608.19658#bib.bib2)\. The parameter range is sampled asε∈\[10−4​;​10−2\]\\varepsilon\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\], isolating the strongly perturbed regime targeted throughout this study\.

Problem 2addresses the thermal entrance problem with a prescribed wall temperature, known as the Graetz problem, posed on the two\-dimensional domain shown in Fig\.[1](https://arxiv.org/html/2608.19658#S2.F1)\(b\) and governed by the equation[shah1978laminar](https://arxiv.org/html/2608.19658#bib.bib24),[haase2015graetz](https://arxiv.org/html/2608.19658#bib.bib25)

w\(y\)∂Θ∂x\(x;y\)=1Pe∂2Θ∂y2\(x;y\);\(x;y\)∈\(0;1\]×\(0;1\);w\(y\)\\,\\frac\{\\partial\\Theta\}\{\\partial x\}\(x\\mathord\{\\mathchar 59\\relax\}y\)=\\frac\{1\}\{\\mathrm\{Pe\}\}\\,\\frac\{\\partial^\{2\}\\Theta\}\{\\partial y^\{2\}\}\(x\\mathord\{\\mathchar 59\\relax\}y\)\\mathchar 59\\relax\\qquad\(x\\mathord\{\\mathchar 59\\relax\}y\)\\in\(0\\mathord\{\\mathchar 59\\relax\}1\]\\times\(0\\mathord\{\\mathchar 59\\relax\}1\)\\mathchar 59\\relax\(3\)with a parabolic velocity profile for hydrodynamically fully developed flow

w⁡\(y\)=1−y2;w\(y\)=1\-y^\{2\}\\mathchar 59\\relax\(4\)centerline symmetry

∂Θ∂y​\(x​;​0\)=0;\\frac\{\\partial\\Theta\}\{\\partial y\}\(x\\mathord\{\\mathchar 59\\relax\}0\)=0\\mathchar 59\\relax\(5\)a homogeneous Dirichlet condition at the wall

Θ⁡\(x​;​1\)=0\.\\Theta\(x\\mathord\{\\mathchar 59\\relax\}1\)=0\.\(6\)The inlet condition is

Θ⁡\(0​;​y\)=Θin​\(y\)\.\\Theta\(0\\mathord\{\\mathchar 59\\relax\}y\)=\\Theta\_\{\\mathrm\{in\}\}\(y\)\.\(7\)
Here,xxdenotes the dimensionless axial coordinate along the channel,yydenotes the dimensionless wall\-normal coordinate across the channel half\-width,Θ⁡\(x​;​y\)\\Theta\(x\\mathord\{\\mathchar 59\\relax\}y\)denotes the dimensionless temperature difference relative to the prescribed wall temperature, andw⁡\(y\)w\(y\)denotes the prescribed dimensionless axial velocity profile\. Because the boundary conditions are already defined, the surrogate directly maps the randomized inlet profileΘin​\(y\)\\Theta\_\{\\mathrm\{in\}\}\(y\)to the wall\-normal temperature profileΘ⁡\(x​;​y\)\\Theta\(x\\mathord\{\\mathchar 59\\relax\}y\)evaluated at a queried axial locationx∈\[0\.05​;​1\.0\]x\\in\[0\.05\\mathord\{\\mathchar 59\\relax\}1\.0\]\. The thermal Péclet numberPe=U​L/α\\mathrm\{Pe\}=UL/\\alphameasures axial advection relative to transverse thermal diffusion, whereUUdenotes the characteristic axial velocity,LLdenotes the reference length scale, andα\\alphadenotes the thermal diffusivity\. SamplingPe−1∈\[10−4​;​10−2\]\\mathrm\{Pe\}^\{\-1\}\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\], which is the singular perturbation parameter analogous toε\\varepsilonin Problem 1, isolates the high\-Péclet regime characterized by increasingly thin wall\-attached thermal layers\.

Problem 3investigates the concentration entrance problem with an absorbing wall, posed on the same domain used in Problem 2 and expressed by the equation[popel1978mass](https://arxiv.org/html/2608.19658#bib.bib26)

w\(y\)∂c∂x\(x;y\)=1Pem∂2c∂y2\(x;y\);\(x;y\)∈\(0;1\]×\(0;1\);w\(y\)\\,\\frac\{\\partial c\}\{\\partial x\}\(x\\mathord\{\\mathchar 59\\relax\}y\)=\\frac\{1\}\{\\mathrm\{Pe\}\_\{m\}\}\\,\\frac\{\\partial^\{2\}c\}\{\\partial y^\{2\}\}\(x\\mathord\{\\mathchar 59\\relax\}y\)\\mathchar 59\\relax\\qquad\(x\\mathord\{\\mathchar 59\\relax\}y\)\\in\(0\\mathord\{\\mathchar 59\\relax\}1\]\\times\(0\\mathord\{\\mathchar 59\\relax\}1\)\\mathchar 59\\relax\(8\)with the same parabolic velocity profile as in Problem 2

w⁡\(y\)=1−y2;w\(y\)=1\-y^\{2\}\\mathchar 59\\relax\(9\)the same centerline symmetry condition

∂c∂y​\(x​;​0\)=0;\\frac\{\\partial c\}\{\\partial y\}\(x\\mathord\{\\mathchar 59\\relax\}0\)=0\\mathchar 59\\relax\(10\)the inlet condition

c⁡\(0​;​y\)=cin​\(y\);c\(0\\mathord\{\\mathchar 59\\relax\}y\)=c\_\{\\mathrm\{in\}\}\(y\)\\mathchar 59\\relax\(11\)and a Robin condition at the wall[debarnot2018graetz](https://arxiv.org/html/2608.19658#bib.bib27)

−∂c∂y​\(x​;​1\)=Da​c​\(x​;​1\)\.\-\\frac\{\\partial c\}\{\\partial y\}\(x\\mathord\{\\mathchar 59\\relax\}1\)=\\mathrm\{Da\}\\,c\(x\\mathord\{\\mathchar 59\\relax\}1\)\.\(12\)
Here,c⁡\(x​;​y\)c\(x\\mathord\{\\mathchar 59\\relax\}y\)denotes the dimensionless concentration field,cin​\(y\)c\_\{\\mathrm\{in\}\}\(y\)denotes the randomized inlet concentration profile,Pem\\mathrm\{Pe\}\_\{m\}denotes the mass\-transfer Péclet number, andDa\\mathrm\{Da\}denotes the Damköhler number, representing the dimensionless wall\-absorption strength in a partially absorbing channel[aquino2024equilibrium](https://arxiv.org/html/2608.19658#bib.bib28)\. Under this nondimensionalization, both spatial coordinates are scaled by the channel half\-widthLL\. The dimensionless groups arePem=U​L/D\\mathrm\{Pe\}\_\{m\}=UL/DandDa=k​L/D\\mathrm\{Da\}=kL/D, whereUUdenotes the characteristic axial velocity,DDdenotes the molecular diffusivity, andkkdenotes a first\-order surface uptake coefficient with units of velocity\. SamplingPem−1∈\[10−4​;​10−2\]\\mathrm\{Pe\}\_\{m\}^\{\-1\}\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\], analogous toPe−1\\mathrm\{Pe\}^\{\-1\}in Problem 2, isolates the high\-Péclet mass\-transfer regime characterized by increasingly thin wall\-attached concentration layers\. To examine the effect of wall\-absorption strength, Problem 3 is evaluated at three wall\-absorption strengths atDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0\. Similar to the thermal Graetz model, the surrogate maps the randomized inlet condition to the wall\-normal concentration profilec⁡\(x​;​y\)c\(x\\mathord\{\\mathchar 59\\relax\}y\)evaluated at a queried axial locationx∈\[0\.05​;​1\.0\]x\\in\[0\.05\\mathord\{\\mathchar 59\\relax\}1\.0\]\.

### 2\.2Operator learning formulation and baseline trunks

For each benchmark problem, the learning task is to approximate a profile\-valued operator

𝒢⁡\(𝐬​;​μ\)=u⁡\(⋅,𝐬​;​μ\);\\mathcal\{G\}\(\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)=u\(\\cdot;\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)\\mathchar 59\\relax\(13\)where𝐬\\mathbf\{s\}denotes the sensor values of the problem\-specific input function andμ\\mudenotes the problem parameters\. In Problem 1,𝐬\\mathbf\{s\}contains samples of the source termq⁡\(x\)q\(x\), whereas in Problems 2 and 3 it contains samples of the inlet temperature profileΘ⁡\(0​;​y\)\\Theta\(0\\mathord\{\\mathchar 59\\relax\}y\)and the inlet concentration profilec⁡\(0​;​y\)c\(0\\mathord\{\\mathchar 59\\relax\}y\), respectively\. For each query coordinateξ\\xi, the predicted profile value is written in the standard DeepONet form

𝒢⁡\(𝐬​;​μ\)​\(ξ\)≈u^​\(ξ,𝐬​;​μ\)=𝐛⁡\(𝐬​;​μ\)⋅ϕ⁡\(ξ\)\.\\mathcal\{G\}\(\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)\(\\xi\)\\approx\\widehat\{u\}\(\\xi;\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)=\\mathbf\{b\}\(\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)\\cdot\\boldsymbol\{\\phi\}\(\\xi\)\.\(14\)
Here,𝐛⁡\(𝐬​;​μ\)=\[b1​\(𝐬​;​μ\)​;​…​;​bp​\(𝐬​;​μ\)\]\\mathbf\{b\}\(\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)=\[b\_\{1\}\(\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)\\mathord\{\\mathchar 59\\relax\}\\ldots\\mathord\{\\mathchar 59\\relax\}b\_\{p\}\(\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)\]denotes the branch coefficient vector andϕ⁡\(ξ\)=\[ϕ1​\(ξ\)​;​…​;​ϕp​\(ξ\)\]\\boldsymbol\{\\phi\}\(\\xi\)=\[\\phi\_\{1\}\(\\xi\)\\mathord\{\\mathchar 59\\relax\}\\ldots\\mathord\{\\mathchar 59\\relax\}\\phi\_\{p\}\(\\xi\)\]denotes the trunk feature vector, whosekkth components arebk​\(𝐬​;​μ\)b\_\{k\}\(\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)andϕk​\(ξ\)\\phi\_\{k\}\(\\xi\), respectively\. Specifically,ξ\\xirepresents the interval coordinatexxin Problem 1 and the wall\-normal coordinateyyin Problems 2 and 3\. The parameter vectorμ\\mucollects the sample\-dependent parameters that condition the operator output, corresponding toε\\varepsilonin Problem 1,\(Pe−1;x\)\(\\mathrm\{Pe\}^\{\-1\}\\mathchar 59\\relax x\)in Problem 2, and\(Pem−1;x\)\(\\mathrm\{Pe\}\_\{m\}^\{\-1\}\\mathchar 59\\relax x\)in Problem 3\. For Problem 3, separate operator surrogates are constructed for fixed values ofDa\\mathrm\{Da\}instead of including the Damköhler number inμ\\mu\. For notational brevity, the dependence on\(𝐬​;​μ\)\(\\mathbf\{s\}\\mathord\{\\mathchar 59\\relax\}\\mu\)is omitted below whenever the input profile and parameter value are fixed\.

To isolate the effect of trunk\-basis design, the surrogate models compared within each benchmark use the same branch\-network architecture, namely a multilayer perceptron \(MLP\) with Gaussian error linear unit \(GELU\) activations and Xavier initialization\. This MLP comprises three hidden layers of width 256, and its output dimension is fixed top=129p=129to match the dimension of each compared trunk\. Before entering the branch network, the first component of the parameter vectorμ\\muis encoded through its base\-10 logarithm, with any remaining components appended in raw form\.

Within this shared framework, three distinct trunk designs are compared\. The branch network provides the sample\- and parameter\-dependent coefficients, while the trunk representation defines the corresponding output approximation family\. The first baseline model is the standard coordinate\-trunk DeepONet, referred to as theVanillamodel, with a trunk MLP consisting of three hidden layers of width128128and an output dimension ofp=129p=129[lu2021learning](https://arxiv.org/html/2608.19658#bib.bib6)\. Unlike the Chebyshev and REC trunks introduced in the following, this baseline uses a learned trunk MLP that takes both the profile coordinate and the encoded parameter vector as inputs\. Thus, for this baseline, Eq\. \([14](https://arxiv.org/html/2608.19658#S2.E14)\) is understood with the coordinate\-only trunk vectorϕ⁡\(ξ\)\\boldsymbol\{\\phi\}\(\\xi\)replaced by a learned, parameter\-conditioned trunk vectorϕvan​\(ξ​;​μ\)\\boldsymbol\{\\phi\}^\{\\mathrm\{van\}\}\(\\xi\\mathord\{\\mathchar 59\\relax\}\\mu\)\.

The second baseline model is the Chebyshev\-trunk DeepONet, labeled as theChebyshevmodel\. While prior approaches, such as SEDONet, process a Chebyshev feature vector through an additional trainable network[abid2025sedonet](https://arxiv.org/html/2608.19658#bib.bib15), the present baseline uses the Chebyshev polynomials directly as the trunk dictionary\. This formulation isolates the effect of the prescribed trunk dictionary without introducing an auxiliary trainable transformation\. For this baseline, the trunk vectorϕ⁡\(ξ\)\\boldsymbol\{\\phi\}\(\\xi\)in Eq\. \([14](https://arxiv.org/html/2608.19658#S2.E14)\) is the Chebyshev dictionary

ϕcheb​\(ξ\)=\[ϕ1cheb​\(ξ\)​;​…​;​ϕ129cheb​\(ξ\)\]=\[T0​\(2​ξ−1\)​;​…​;​T128​\(2​ξ−1\)\];\\begin\{split\}\\boldsymbol\{\\phi\}^\{\\mathrm\{cheb\}\}\(\\xi\)=&\\big\[\\phi^\{\\mathrm\{cheb\}\}\_\{1\}\(\\xi\)\\mathord\{\\mathchar 59\\relax\}\\ldots\\mathord\{\\mathchar 59\\relax\}\\phi^\{\\mathrm\{cheb\}\}\_\{129\}\(\\xi\)\\big\]\\\\ =&\\big\[T\_\{0\}\(2\\xi\-1\)\\mathord\{\\mathchar 59\\relax\}\\ldots\\mathord\{\\mathchar 59\\relax\}T\_\{128\}\(2\\xi\-1\)\\big\]\\mathchar 59\\relax\\end\{split\}\(15\)whereTjT\_\{j\}denotes the Chebyshev polynomial of the first kind of degreejj\.

### 2\.3Proposed REC trunk

The proposed third model is the REC\-trunk DeepONet, denoted as theRECmodel, motivated by rational discretizations for singularly perturbed BVPs, where layer\-oriented rational mappings resolve thin layers with fewer global polynomial degrees[wang2010rational](https://arxiv.org/html/2608.19658#bib.bib29)\. Unlike the Chebyshev trunk in Eq\. \([15](https://arxiv.org/html/2608.19658#S2.E15)\), whose trunk dictionary consists entirely of Chebyshev polynomial elements, the REC trunk combines a low\-degree outer Chebyshev subdictionary for the smooth outer field with an inner rational subdictionary for the localized layer\. The rational dictionary elements are constructed by the adaptive Antoulas–Anderson \(AAA\) algorithm and remain unchanged during DeepONet training[nakatsukasa2018aaa](https://arxiv.org/html/2608.19658#bib.bib30)\.

The inner rational subdictionary is derived from a canonical exponentially decaying layer family on the unit interval,

ψδ​\(ζ\)=exp⁡\(−1−ζδ\);ζ∈\[0​;​1\]​;δ∈\[10−4​;​10−2\];\\psi\_\{\\delta\}\(\\zeta\)=\\exp\\\!\\left\(\-\\frac\{1\-\\zeta\}\{\\delta\}\\right\)\\mathchar 59\\relax\\qquad\\zeta\\in\[0\\mathord\{\\mathchar 59\\relax\}1\]\\mathord\{\\mathchar 59\\relax\}\\qquad\\delta\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\]\\mathchar 59\\relax\(16\)whereζ\\zetadenotes the dictionary coordinate andδ\\deltadenotes the prototype decay parameter, whose sampled range matches the numerical range ofε\\varepsilon,Pe−1\\mathrm\{Pe\}^\{\-1\}, andPem−1\\mathrm\{Pe\}\_\{m\}^\{\-1\}used in the three benchmark problems defined in Section[2\.1](https://arxiv.org/html/2608.19658#S2.SS1)\. This family represents a one\-sided localized layer attached toζ=1\\zeta=1, corresponding to the right endpointx=1x=1in Problem 1 and to the wall endpointy=1y=1in Problems 2 and 3\.

The construction proceeds by choosing129−Nout129\-N\_\{\\mathrm\{out\}\}logarithmically spaced valuesδk∈\[10−4​;​10−2\]\\delta\_\{k\}\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\], one for each rational dictionary element\. For each sampled valueδk\\delta\_\{k\}, the profileψδk\\psi\_\{\\delta\_\{k\}\}is evaluated on a10251025\-point Chebyshev–Lobatto grid inζ∈\[0​;​1\]\\zeta\\in\[0\\mathord\{\\mathchar 59\\relax\}1\]\. The AAA algorithm approximates each sampled profile by a barycentric rational function with at most1212retained terms and tolerances of10−810^\{\-8\}for both approximation error and support\-point matching checks\. The resulting rational approximant is used as one rational dictionary element and is written as

rk​\(ζ\)=∑j=1mkwk​;​j​fk​;​jζ−zk​;​j∑j=1mkwk​;​jζ−zk​;​j;r\_\{k\}\(\\zeta\)=\\frac\{\\sum\_\{j=1\}^\{m\_\{k\}\}\\dfrac\{w\_\{k\\mathord\{\\mathchar 59\\relax\}j\}f\_\{k\\mathord\{\\mathchar 59\\relax\}j\}\}\{\\zeta\-z\_\{k\\mathord\{\\mathchar 59\\relax\}j\}\}\}\{\\sum\_\{j=1\}^\{m\_\{k\}\}\\dfrac\{w\_\{k\\mathord\{\\mathchar 59\\relax\}j\}\}\{\\zeta\-z\_\{k\\mathord\{\\mathchar 59\\relax\}j\}\}\}\\mathchar 59\\relax\(17\)wherekkindexes the rational dictionary element associated with the sampled layer profileψδk\\psi\_\{\\delta\_\{k\}\},zk​;​jz\_\{k\\mathord\{\\mathchar 59\\relax\}j\}denotes thejjth support point selected by AAA from the10251025Chebyshev–Lobatto points,fk​;​j=ψδk​\(zk​;​j\)f\_\{k\\mathord\{\\mathchar 59\\relax\}j\}=\\psi\_\{\\delta\_\{k\}\}\(z\_\{k\\mathord\{\\mathchar 59\\relax\}j\}\)denotes the sampled layer value atzk​;​jz\_\{k\\mathord\{\\mathchar 59\\relax\}j\},wk​;​jw\_\{k\\mathord\{\\mathchar 59\\relax\}j\}denotes the corresponding barycentric weight, andmkm\_\{k\}denotes the number of retained support points for thekkth rational dictionary element[nakatsukasa2018aaa](https://arxiv.org/html/2608.19658#bib.bib30)\. The stored support points, sampled values, and barycentric weights define the rational dictionary elements used by the REC trunk\.

When inserted into the DeepONet trunk, these rational dictionary elements are evaluated at the profile coordinateξ\\xi\. For the three benchmarks,ξ=x\\xi=xin Problem 1 andξ=y\\xi=yin Problems 2 and 3, withx​;​y∈\[0​;​1\]x\\mathord\{\\mathchar 59\\relax\}y\\in\[0\\mathord\{\\mathchar 59\\relax\}1\]\. No additional coordinate map is introduced, and the same dictionary elements are evaluated withζ=ξ\\zeta=\\xi\. Thus, the trunk dictionaryϕ⁡\(ξ\)\\boldsymbol\{\\phi\}\(\\xi\)in Eq\. \([14](https://arxiv.org/html/2608.19658#S2.E14)\) for REC trunk can be written as

ϕREC​\(ξ\)=\[ϕ1REC​\(ξ\)​;​…​;​ϕ129REC​\(ξ\)\]=\[T0​\(2​ξ−1\)​;​…​;​TNout−1​\(2​ξ−1\)​;​r1​\(ξ\)​;​…​;​r129−Nout​\(ξ\)\];\\begin\{split\}\\boldsymbol\{\\phi\}^\{\\mathrm\{REC\}\}\(\\xi\)=&\\big\[\\phi^\{\\mathrm\{REC\}\}\_\{1\}\(\\xi\)\\mathord\{\\mathchar 59\\relax\}\\ldots\\mathord\{\\mathchar 59\\relax\}\\phi^\{\\mathrm\{REC\}\}\_\{129\}\(\\xi\)\\big\]\\\\ =&\\big\[T\_\{0\}\(2\\xi\-1\)\\mathord\{\\mathchar 59\\relax\}\\ldots\\mathord\{\\mathchar 59\\relax\}T\_\{N\_\{\\mathrm\{out\}\}\-1\}\(2\\xi\-1\)\\mathord\{\\mathchar 59\\relax\}r\_\{1\}\(\\xi\)\\mathord\{\\mathchar 59\\relax\}\\ldots\\mathord\{\\mathchar 59\\relax\}r\_\{129\-N\_\{\\mathrm\{out\}\}\}\(\\xi\)\\big\]\\mathchar 59\\relax\\end\{split\}\(18\)where the firstNoutN\_\{\\mathrm\{out\}\}Chebyshev dictionary elements form the outer Chebyshev subdictionary, and the remaining\(129−Nout\)\(129\-N\_\{\\mathrm\{out\}\}\)rational dictionary elements form the inner rational subdictionary, yielding a total trunk dimension of129129\. The main discussion in Section[3](https://arxiv.org/html/2608.19658#S3)focuses on the REC trunk withNout=16N\_\{\\mathrm\{out\}\}=16, while Section[3\.5](https://arxiv.org/html/2608.19658#S3.SS5)analyzes additional REC trunks withNout=33N\_\{\\mathrm\{out\}\}=33,6565, and9797to examine how the outer\-inner split affects the approximation behavior\. Appendix[A](https://arxiv.org/html/2608.19658#A1)gives an idealized representation argument for this outer\-inner split\. The Chebyshev expansion ofψδ\\psi\_\{\\delta\}shows that decreasingδ\\deltashifts non\-negligible spectral weight toward higher polynomial degrees, whereas the inner rational subdictionary provides coverage across the sampled rangelog10⁡δ∈\[−4​;−2\]\\log\_\{10\}\\delta\\in\[\-4\\mathord\{\\mathchar 59\\relax\}\-2\]\. Under the idealized decomposition in Eq\. \([51](https://arxiv.org/html/2608.19658#A1.E51)\), the REC trunk assigns the smooth outer field and the thin boundary\-attached layer to separate parts of the trunk dictionary, rather than imposing an explicit asymptotic formula\.

### 2\.4Input sampling and reference profile generation

The training, validation, and test datasets are constructed by pairing sampled problem inputs with their corresponding numerical reference profiles\. In Problem 1, each input consists of a randomized source term and the perturbation parameterε\\varepsilon\. In Problems 2 and 3, each input consists of a randomized inlet profile, the inverse\-Péclet\-type transport parameter, and the queried axial locationxx\. The source term in Problem 1 and the inlet profiles in Problems 2 and 3 are sampled at fixed Chebyshev–Lobatto sensor nodes, which include the interval endpoints and cluster near the boundaries\. The resulting sensor vector is supplied to the branch network together with the encoded parameter values, and the corresponding reference profile is evaluated at fixed Chebyshev–Lobatto output nodes\. In all three benchmarks, the source or inlet function is sampled at129129Chebyshev–Lobatto sensor locations, which include the interval endpoints and cluster toward the boundaries\. The target profile is evaluated at257257Chebyshev–Lobatto output locations, providing denser resolution near the endpoint or wall region where the localized layer develops\.

The randomized input functions are constructed as finite smooth expansions rather than taken from external data\. In Problem 1, the source term is formed from four sine modes and two localized Gaussian components,

q⁡\(x\)=∑m=14am\(q\)​sin⁡\(m​π​x\)\+∑j=12a~j\(q\)​exp⁡\[−12​\(x−xj\(q\)ℓj\(q\)\)2\]\.q\(x\)=\\sum\_\{m=1\}^\{4\}a\_\{m\}^\{\(q\)\}\\sin\(m\\pi x\)\+\\sum\_\{j=1\}^\{2\}\\widetilde\{a\}\_\{j\}^\{\(q\)\}\\exp\\\!\\left\[\-\\frac\{1\}\{2\}\\left\(\\frac\{x\-x\_\{j\}^\{\(q\)\}\}\{\\ell\_\{j\}^\{\(q\)\}\}\\right\)^\{2\}\\right\]\.\(19\)The modal and Gaussian amplitudes are drawn independently according toam\(q\)​;​a~j\(q\)∼𝒩⁡\(0​;​1\)a\_\{m\}^\{\(q\)\}\\mathord\{\\mathchar 59\\relax\}\\widetilde\{a\}\_\{j\}^\{\(q\)\}\\sim\\mathcal\{N\}\(0\\mathord\{\\mathchar 59\\relax\}1\), while the Gaussian centers and widths are sampled asxj\(q\)∼𝒰⁡\(0​;​1\)x\_\{j\}^\{\(q\)\}\\sim\\mathcal\{U\}\(0\\mathord\{\\mathchar 59\\relax\}1\)andℓj\(q\)∼𝒰⁡\(0\.03​;​0\.20\)\\ell\_\{j\}^\{\(q\)\}\\sim\\mathcal\{U\}\(0\.03\\mathord\{\\mathchar 59\\relax\}0\.20\)\. The perturbation parameter is sampled by drawinglog10⁡ε\\log\_\{10\}\\varepsilonuniformly from\[−4​;−2\]\[\-4\\mathord\{\\mathchar 59\\relax\}\-2\]\.

For Problems 2 and 3, the inlet function is generated from a positive baseline, four cosine modes, and two localized Gaussian components\. This common inlet representation is denoted bygin​\(y\)g\_\{\\mathrm\{in\}\}\(y\), wheregin=Θing\_\{\\mathrm\{in\}\}=\\Theta\_\{\\mathrm\{in\}\}for Problem 2 andgin=cing\_\{\\mathrm\{in\}\}=c\_\{\\mathrm\{in\}\}for Problem 3\. Before imposing the lower bound, the raw inlet function is

graw​\(y\)=1\+∑m=14am\(in\)​cos⁡\(\(m−1\)​π​y\)\+∑j=12a~j\(in\)​exp⁡\[−12​\(y−yj\(in\)ℓj\(in\)\)2\]\.g\_\{\\mathrm\{raw\}\}\(y\)=1\+\\sum\_\{m=1\}^\{4\}a\_\{m\}^\{\(\\mathrm\{in\}\)\}\\cos\(\(m\-1\)\\pi y\)\+\\sum\_\{j=1\}^\{2\}\\widetilde\{a\}\_\{j\}^\{\(\\mathrm\{in\}\)\}\\exp\\\!\\left\[\-\\frac\{1\}\{2\}\\left\(\\frac\{y\-y\_\{j\}^\{\(\\mathrm\{in\}\)\}\}\{\\ell\_\{j\}^\{\(\\mathrm\{in\}\)\}\}\\right\)^\{2\}\\right\]\.\(20\)The cosine amplitudes are drawn according toam\(in\)∼𝒩⁡\(0​;​0\.182\)a\_\{m\}^\{\(\\mathrm\{in\}\)\}\\sim\\mathcal\{N\}\(0\\mathord\{\\mathchar 59\\relax\}0\.18^\{2\}\), the Gaussian amplitudes asa~j\(in\)∼𝒩⁡\(0​;​0\.122\)\\widetilde\{a\}\_\{j\}^\{\(\\mathrm\{in\}\)\}\\sim\\mathcal\{N\}\(0\\mathord\{\\mathchar 59\\relax\}0\.12^\{2\}\), the Gaussian centers asyj\(in\)∼𝒰⁡\(0​;​1\)y\_\{j\}^\{\(\\mathrm\{in\}\)\}\\sim\\mathcal\{U\}\(0\\mathord\{\\mathchar 59\\relax\}1\), and the Gaussian widths asℓj\(in\)∼𝒰⁡\(0\.05​;​0\.18\)\\ell\_\{j\}^\{\(\\mathrm\{in\}\)\}\\sim\\mathcal\{U\}\(0\.05\\mathord\{\\mathchar 59\\relax\}0\.18\)\. To exclude sign\-changing inlet temperature differences in Problem 2, for which heating and cooling are equivalent up to a global sign change, and negative inlet concentrations in Problem 3, the inlet function used to generate each reference solution is defined as

gin​\(y\)=max⁡\{graw​\(y\)​;​10−3\}\.g\_\{\\mathrm\{in\}\}\(y\)=\\max\\\{g\_\{\\mathrm\{raw\}\}\(y\)\\mathord\{\\mathchar 59\\relax\}10^\{\-3\}\\\}\.\(21\)For Problems 2 and 3, the transport parameter is sampled by drawinglog10⁡Pe−1\\log\_\{10\}\\mathrm\{Pe\}^\{\-1\}andlog10⁡Pem−1\\log\_\{10\}\\mathrm\{Pe\}\_\{m\}^\{\-1\}, respectively, uniformly from\[−4​;−2\]\[\-4\\mathord\{\\mathchar 59\\relax\}\-2\]\. The queried axial location is sampled independently asx∼𝒰⁡\(0\.05​;​1\.0\)x\\sim\\mathcal\{U\}\(0\.05\\mathord\{\\mathchar 59\\relax\}1\.0\)\.

The target profiles are generated by problem\-specific deterministic numerical solvers\. Problem 1 is solved using an upwind finite\-difference discretization on a piecewise\-uniform Shishkin\-type mesh whose transition point depends onε\\varepsilonand which places a finer subgrid near the outflow boundaryx=1x=1where the endpoint layer forms\. For Problems 2 and 3, the entrance\-transport equations are solved by marching in the axial coordinate from the randomized inlet profile to the queried locationxx\. At each axial step, an implicit finite\-difference discretization of the transverse diffusion operator on a uniform grid yields a tridiagonal linear system for the updated wall\-normal profile\. The centerline Neumann condition is imposed aty=0y=0, while the wall condition aty=1y=1is imposed as a Dirichlet condition for Problem 2 and as a Robin condition for Problem 3\.

For the convergence check, the same sampled source term and perturbation parameter in Problem 1 are solved on the40964096\-cell Shishkin mesh used for dataset generation and an81928192\-cell Shishkin mesh and compared after linear interpolation onto the same257257Chebyshev–Lobatto output locations\. For Problems 2 and 3, the same inlet profile, inverse\-Péclet\-type parameter, and queried axial location are solved using the dataset resolution of257257wall\-normal nodes with400400axial steps per unit length and using513513wall\-normal nodes with800800axial steps per unit length, again compared on the same257257Chebyshev–Lobatto output locations\. Across3232randomly sampled checks for each benchmark problem, the mean relative differences between the two numerical profiles are4\.60×10−44\.60\\times 10^\{\-4\}for Problem 1,1\.97×10−41\.97\\times 10^\{\-4\}for Problem 2, and1\.06×10−41\.06\\times 10^\{\-4\},1\.09×10−41\.09\\times 10^\{\-4\}, and1\.43×10−41\.43\\times 10^\{\-4\}for Problem 3 atDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0, respectively\. The corresponding maximum relative differences are9\.30×10−49\.30\\times 10^\{\-4\},6\.64×10−46\.64\\times 10^\{\-4\},4\.30×10−44\.30\\times 10^\{\-4\},4\.53×10−44\.53\\times 10^\{\-4\}, and5\.26×10−45\.26\\times 10^\{\-4\}\.

### 2\.5Training and evaluation

To quantify model accuracy, four profile\-based error metrics are defined on the common discrete evaluation grid\. Letuudenote the numerical reference profile at the queried condition, letu^\\widehat\{u\}denote the corresponding model prediction, and define the profile error ase=u^−ue=\\widehat\{u\}\-u\. The relative discrete profile error is

E2rel=‖e‖2‖u‖2;E\_\{2\}^\{\\mathrm\{rel\}\}=\\frac\{\\\|e\\\|\_\{2\}\}\{\\\|u\\\|\_\{2\}\}\\mathchar 59\\relax\(22\)and the pointwise maximum error is

E∞=‖e‖∞\.E\_\{\\infty\}=\\\|e\\\|\_\{\\infty\}\.\(23\)Here,∥⋅∥2\\\|\\cdot\\\|\_\{2\}denotes the Euclidean norm of the profile values on the fixed output grid, and∥⋅∥∞\\\|\\cdot\\\|\_\{\\infty\}denotes the maximum absolute value over the same grid\.

The layer\-focused maximum error measures the largest absolute discrepancy inside a problem\-dependent layer stripΩs\\Omega\_\{s\},

Emaxlayer=‖e‖∞​;​Ωs\.E\_\{\\max\}^\{\\mathrm\{layer\}\}=\\\|e\\\|\_\{\\infty\\mathord\{\\mathchar 59\\relax\}\\Omega\_\{s\}\}\.\(24\)For Problem 1, the layer strip is the right\-endpoint boundary\-layer region

Ωℓ=\{xj:xj≥1−C​ε​\|log⁡ε\|\};\\Omega\_\{\\ell\}=\\\{x\_\{j\}:\\ x\_\{j\}\\geq 1\-C\\,\\varepsilon\|\\log\\varepsilon\|\\\}\\mathchar 59\\relax\(25\)whereCCdenotes a layer\-width factor\. This definition selects the grid points in the interval\[1−C​ε​\|log⁡ε\|​;​1\]\\left\[1\-C\\varepsilon\\left\|\\log\\varepsilon\\right\|\\mathord\{\\mathchar 59\\relax\}1\\right\]adjacent to the outflow boundaryx=1x=1, where the outflow Dirichlet condition is satisfied through an endpoint layer under positive convection\. In Problem 1,Ωs=Ωℓ\\Omega\_\{s\}=\\Omega\_\{\\ell\}andC=5C=5\. For Problems 2 and 3, the layer strip is the near\-wall region adjacent to the controlled or absorbing wall,

Ωw=\{yj:yj≥1−C​η​x\};\\Omega\_\{w\}=\\\{y\_\{j\}:\\ y\_\{j\}\\geq 1\-C\\sqrt\{\\eta x\}\\\}\\mathchar 59\\relax\(26\)whereη=Pe−1\\eta=\\mathrm\{Pe\}^\{\-1\}in Problem 2,η=Pem−1\\eta=\\mathrm\{Pe\}\_\{m\}^\{\-1\}in Problem 3, andCCdenotes a layer\-width factor\. In these problems,Ωs=Ωw\\Omega\_\{s\}=\\Omega\_\{w\}andC=4C=4\. The near\-wall region in Problems 2 and 3 is used only to define a common near\-wall error measure for comparing the three surrogates and is not intended as an asymptotic estimate of the physical Graetz\-layer thickness\.

The layer\-aware error is defined problem\-dependently to emphasize the localized layer,

ELA=\{\[‖e‖Q2\+ε​‖d​ed​x‖Q2\+‖e‖Q​;​Ωℓ2‖u‖Q2\+ε​‖d​ud​x‖Q2\+‖u‖Q​;​Ωℓ2\]1/2;for Problem 1;\[‖e‖2​;​Ωw2‖u‖2​;​Ωw2\]1/2;for Problems 2 and 3\.E\_\{\\mathrm\{LA\}\}=\\begin\{cases\}\\left\[\\displaystyle\\frac\{\\\|e\\\|\_\{Q\}^\{2\}\+\\varepsilon\\left\\\|\\dfrac\{de\}\{dx\}\\right\\\|\_\{Q\}^\{2\}\+\\\|e\\\|\_\{Q\\mathord\{\\mathchar 59\\relax\}\\Omega\_\{\\ell\}\}^\{2\}\}\{\\\|u\\\|\_\{Q\}^\{2\}\+\\varepsilon\\left\\\|\\dfrac\{du\}\{dx\}\\right\\\|\_\{Q\}^\{2\}\+\\\|u\\\|\_\{Q\\mathord\{\\mathchar 59\\relax\}\\Omega\_\{\\ell\}\}^\{2\}\}\\right\]^\{1/2\}\\mathchar 59\\relax&\\text\{for Problem~1\}\\mathchar 59\\relax\\\\\[8\.61108pt\] \\left\[\\displaystyle\\frac\{\\\|e\\\|\_\{2\\mathord\{\\mathchar 59\\relax\}\\Omega\_\{w\}\}^\{2\}\}\{\\\|u\\\|\_\{2\\mathord\{\\mathchar 59\\relax\}\\Omega\_\{w\}\}^\{2\}\}\\right\]^\{1/2\}\\mathchar 59\\relax&\\text\{for Problems~2 and~3\}\.\\end\{cases\}\(27\)Here,∥⋅∥Q\\\|\\cdot\\\|\_\{Q\}denotes the trapezoidal\-quadrature norm on the fixed output grid\. The restricted norms∥⋅∥Q​;​Ωℓ\\\|\\cdot\\\|\_\{Q\\mathord\{\\mathchar 59\\relax\}\\Omega\_\{\\ell\}\}and∥⋅∥2​;​Ωw\\\|\\cdot\\\|\_\{2\\mathord\{\\mathchar 59\\relax\}\\Omega\_\{w\}\}are evaluated over the corresponding layer strips\. For Problem 1, the derivatives inELAE\_\{\\mathrm\{LA\}\}are evaluated on the same nonuniform output grid using one\-sided two\-point differences at the endpoints and centered two\-point differences at the interior nodes, with the same rule applied to the numerical reference and all model predictions\.

All three models are trained with the same mean\-squared\-error \(MSE\) objective and the same training procedure\. Each model uses30003000training samples and500500validation samples, and performance is evaluated on500500test profiles not used during training or validation\. The Adam optimizer is used with learning rate5×10−45\\times 10^\{\-4\}, batch size6464, weight decay10−610^\{\-6\}, gradient clipping at5\.05\.0, and an exponential learning\-rate scheduler with decay factor0\.9950\.995, for at most250250epochs\. This is followed by a limited\-memory Broyden–Fletcher–Goldfarb–Shanno \(L\-BFGS\) quasi\-Newton refinement with at most8080iterations and strong\-Wolfe line search\. The checkpoint used for testing is selected by the lowest validation value ofEmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}among the models obtained during Adam optimization and after L\-BFGS refinement\. Each trunk design is trained in five independent runs on the same data split, with model comparisons made seed by seed\. Thus, the same500500test profiles are evaluated for each of the five trained instances, yielding25002500error evaluations per model and benchmark\.

## 3Results and discussion

### 3\.1Singularly perturbed scalar BVP

Figure 2:Ratios ofE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}overε∈\[10−4​;​10−2\]\\varepsilon\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\]for the singularly perturbed scalar BVP atNout=16N\_\{\\mathrm\{out\}\}=16\.Figure[2](https://arxiv.org/html/2608.19658#S3.F2)shows the ratios of the error from theRECmodel to the error from theChebyshevmodel and to the error from theVanillamodel overε∈\[10−4​;​10−2\]\\varepsilon\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\]forE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}\. The interval\[10−4​;​10−2\]\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\]is divided into eight equal\-width bins on a logarithmic scale, with nominal bin edges1\.00×10−41\.00\\times 10^\{\-4\},1\.78×10−41\.78\\times 10^\{\-4\},3\.16×10−43\.16\\times 10^\{\-4\},5\.62×10−45\.62\\times 10^\{\-4\},1\.00×10−31\.00\\times 10^\{\-3\},1\.78×10−31\.78\\times 10^\{\-3\},3\.16×10−33\.16\\times 10^\{\-3\},5\.62×10−35\.62\\times 10^\{\-3\}, and1\.00×10−21\.00\\times 10^\{\-2\}\. For each test profile in a given bin, the error from theRECmodel is divided by the error from theChebyshevmodel or theVanillamodel for the same metric\. The marker is placed at the geometric center of the bin and gives the median of these ratios\. The shaded band gives the interquartile range \(IQR\), from the 25th to the 75th percentile\. Relative to theVanillamodel, the median ratios for theRECmodel remain below one in all eightε\\varepsilonbins and for all four metrics, with the largest median value reaching0\.9770\.977only forELAE\_\{\\mathrm\{LA\}\}in the largest\-ε\\varepsilonbin\. Relative to theChebyshevmodel, the median ratios remain near unity across mostε\\varepsilonbins, whereas the smallestε\\varepsilonbin gives ratios below unity for all four metrics, with median values0\.8800\.880,0\.8850\.885,0\.8910\.891, and0\.8050\.805forE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}, respectively\. These trends indicate that, for the singularly perturbed scalar BVP, theRECmodel clearly separates from theVanillamodel across the tested parameter range, while its advantage over theChebyshevmodel is concentrated where the perturbation parameter is smallest\.

Table 1:Comparison of theRECmodel with theVanillaandChebyshevmodels for the singularly perturbed scalar BVP atNout=16N\_\{\\mathrm\{out\}\}=16\.Table[1](https://arxiv.org/html/2608.19658#S3.T1)tabulates error ratios for the singularly perturbed scalar BVP atNout=16N\_\{\\mathrm\{out\}\}=16for theRECmodel against both theVanillaandChebyshevmodels\. Here, the second column gives the ratio of mean errors over the five independent training runs, where each run error is averaged over the500500test profiles\. The third column gives the ratio of mean errors over the five independent training runs, where each run error is averaged over the first three parameter\-bin means using the same parameter bins as Fig\.[2](https://arxiv.org/html/2608.19658#S3.F2)\. The fourth column counts how many of the five independent training runs give lower error for theRECmodel than for the denominator model\. The fifth column gives the percentage of the25002500profile comparisons, from the five independent training runs and500500test profiles per run, in which theRECmodel has lower error than the denominator model\. Relative to theVanillamodel, the second\-column ratios are below0\.3350\.335for all four metrics, and the third\-column ratios decrease further to at most0\.1460\.146\. The fourth and fifth columns show the same separation across all five independent training runs and at least92\.1%92\.1\\%of the25002500profile comparisons\. Relative to theChebyshevmodel, the second column remains close to unity forE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\}, andEmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, whileELAE\_\{\\mathrm\{LA\}\}stays above unity\. In the third column, however,E2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\}, andEmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}decrease to0\.9050\.905,0\.9700\.970, and0\.9770\.977, respectively\. These results indicate that, for the singularly perturbed scalar BVP, theRECmodel is most useful when the endpoint layer inu⁡\(x\)u\(x\)becomes thinnest, while theChebyshevmodel remains competitive when the full test set also includes larger perturbation parameters\.

The REC reductions appear mainly in the first three parameter bins, while the full\-test ratios relative to theChebyshevmodel remain near unity or above unity, consistent with the representation argument in Appendix[A](https://arxiv.org/html/2608.19658#A1)\. The endpoint layer in Problem 1 belongs to the same exponential family used to construct the rational dictionary elementsrkr\_\{k\}in Eq\. \([18](https://arxiv.org/html/2608.19658#S2.E18)\), while the same family also admits the exact Chebyshev expansion in Eq\. \([32](https://arxiv.org/html/2608.19658#A1.E32)\) with the active polynomial degree scale in Eq\. \([38](https://arxiv.org/html/2608.19658#A1.E38)\)\. Thus, the layer alignment of the rational dictionary elements does not by itself require theRECmodel to reduce every metric relative to theChebyshevmodel on the singularly perturbed scalar BVP\.

### 3\.2Thermal entrance problem

Figure 3:Ratios ofE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}overPe−1∈\[10−4​;​10−2\]\\mathrm\{Pe\}^\{\-1\}\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\]for the thermal entrance problem atNout=16N\_\{\\mathrm\{out\}\}=16\.Figure[3](https://arxiv.org/html/2608.19658#S3.F3)illustrates the ratios of the error from theRECmodel to the error from theChebyshevmodel and to the error from theVanillamodel for the thermal entrance problem, calculated in the same way as Fig\.[2](https://arxiv.org/html/2608.19658#S3.F2), withPe−1\\mathrm\{Pe\}^\{\-1\}replacingε\\varepsilonand with the independently sampledxxvalues retained within eachPe−1\\mathrm\{Pe\}^\{\-1\}bin\. Relative to theVanillamodel, the median ratios remain below unity for all eightPe−1\\mathrm\{Pe\}^\{\-1\}bins and all four metrics, with the largest median value equal to0\.8100\.810forEmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}\. Relative to theChebyshevmodel, the median ratios also remain below unity for all bins and metrics, with the largest median value equal to0\.9180\.918forELAE\_\{\\mathrm\{LA\}\}\. These results indicate that, for the thermal entrance problem, theRECmodel improves the reconstruction of the wall\-normal temperature profileθ⁡\(y\)\\theta\(y\)across the testedPe−1\\mathrm\{Pe\}^\{\-1\}range and sampled entrance locationsxx\.

Table 2:Comparison of theRECmodel with theVanillaandChebyshevmodels for the thermal entrance problem atNout=16N\_\{\\mathrm\{out\}\}=16\.Table[2](https://arxiv.org/html/2608.19658#S3.T2)summarizes error ratios for the thermal entrance problem atNout=16N\_\{\\mathrm\{out\}\}=16for theRECmodel against both theVanillaandChebyshevmodels, similar to Table[1](https://arxiv.org/html/2608.19658#S3.T1), but with the first three parameter\-bin means computed using the same parameter bins as Fig\.[3](https://arxiv.org/html/2608.19658#S3.F3)\. Relative to theVanillamodel, all four second\-column ratios are below0\.5580\.558, and all four third\-column ratios are below0\.4350\.435\. The fourth column shows that theRECmodel gives lower mean error than theVanillamodel in all five independent training runs for every metric\. The fifth column remains high as well, ranging from80\.6%80\.6\\%to98\.5%98\.5\\%of the25002500profile comparisons\. Relative to theChebyshevmodel, all four second\-column ratios and all four third\-column ratios remain below unity, with values between0\.6900\.690and0\.7730\.773\. The fourth column again shows that theRECmodel gives lower mean error than theChebyshevmodel in all five independent training runs for all four metrics, and the fifth column remains between74\.7%74\.7\\%and90\.9%90\.9\\%\. These results indicate that, for the thermal entrance problem, theRECmodel provides a more reliable representation of the wall\-normal temperature profileθ⁡\(y\)\\theta\(y\)across sampledPe−1\\mathrm\{Pe\}^\{\-1\}values and entrance locationsxx, where the profile must combine a smooth outer region with a wall\-attached thermal layer\.

### 3\.3Concentration entrance problem

Figure 4:Ratios ofE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}overPem−1∈\[10−4​;​10−2\]\\mathrm\{Pe\}\_\{m\}^\{\-1\}\\in\[10^\{\-4\}\\mathord\{\\mathchar 59\\relax\}10^\{\-2\}\]for the concentration entrance problem atNout=16N\_\{\\mathrm\{out\}\}=16andDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0\.Figure[4](https://arxiv.org/html/2608.19658#S3.F4)reports the ratios of the error from theRECmodel to the error from theChebyshevmodel and to the error from theVanillamodel for the concentration entrance problem atDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0, calculated in the same way as Fig\.[2](https://arxiv.org/html/2608.19658#S3.F2), withPem−1\\mathrm\{Pe\}\_\{m\}^\{\-1\}replacingε\\varepsilonand with the independently sampledxxvalues retained within eachPem−1\\mathrm\{Pe\}\_\{m\}^\{\-1\}bin\. Relative to theVanillamodel, the median ratios forE2relE\_\{2\}^\{\\mathrm\{rel\}\}andE∞E\_\{\\infty\}remain below unity for all three Damköhler numbers and all eightPem−1\\mathrm\{Pe\}\_\{m\}^\{\-1\}bins, with largest median values0\.6880\.688and0\.7450\.745, respectively\. ForEmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}andELAE\_\{\\mathrm\{LA\}\}, the median ratios remain below unity in all bins except the sixthPem−1\\mathrm\{Pe\}\_\{m\}^\{\-1\}bin atDa=0\.1\\mathrm\{Da\}=0\.1and1\.01\.0\. Relative to theChebyshevmodel, all median ratios remain below unity for all three Damköhler numbers, all eight bins, and all four metrics, with the largest median value equal to0\.9140\.914forELAE\_\{\\mathrm\{LA\}\}\. These results indicate that, for the concentration entrance problem, theRECmodel improves the reconstruction of the wall\-normal concentration profilec⁡\(y\)c\(y\)across sampled entrance locationsxxand across wall\-absorption regimes ranging from weak uptake atDa=0\.1\\mathrm\{Da\}=0\.1to stronger wall\-adjacent concentration gradients atDa=10\.0\\mathrm\{Da\}=10\.0\.

Table 3:Comparison of theRECmodel with theVanillaandChebyshevmodels for the concentration entrance problem atNout=16N\_\{\\mathrm\{out\}\}=16\.Metric ratioFull test setFirst three parameter binsSeeds with lower REC errorProfiles with lower REC errorProblem 3 \(Da=0\.1\\mathrm\{Da\}=0\.1\)E2rel\|REC/E2rel\|Vanilla\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.4530\.4530\.2910\.2915/55/594\.9%94\.9\\%E∞\|REC/E∞\|Vanilla\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.4760\.4760\.3290\.3295/55/590\.9%90\.9\\%Emaxlayer\|REC/Emaxlayer\|Vanilla\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.7990\.7990\.5850\.5855/55/564\.2%64\.2\\%ELA\|REC/ELA\|Vanilla\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.8270\.8270\.5780\.5784/54/563\.8%63\.8\\%E2rel\|REC/E2rel\|Chebyshev\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.7610\.7610\.6830\.6835/55/587\.7%87\.7\\%E∞\|REC/E∞\|Chebyshev\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.7200\.7200\.6610\.6615/55/590\.4%90\.4\\%Emaxlayer\|REC/Emaxlayer\|Chebyshev\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.7880\.7880\.7840\.7845/55/577\.7%77\.7\\%ELA\|REC/ELA\|Chebyshev\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.8950\.8950\.9430\.9435/55/565\.5%65\.5\\%Problem 3 \(Da=1\.0\\mathrm\{Da\}=1\.0\)E2rel\|REC/E2rel\|Vanilla\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.4630\.4630\.3040\.3045/55/593\.5%93\.5\\%E∞\|REC/E∞\|Vanilla\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.4840\.4840\.3390\.3395/55/590\.8%90\.8\\%Emaxlayer\|REC/Emaxlayer\|Vanilla\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.7570\.7570\.5640\.5645/55/569\.3%69\.3\\%ELA\|REC/ELA\|Vanilla\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.8290\.8290\.6200\.6204/54/568%68\\%E2rel\|REC/E2rel\|Chebyshev\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.7700\.7700\.7020\.7025/55/587\.8%87\.8\\%E∞\|REC/E∞\|Chebyshev\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.7200\.7200\.6700\.6705/55/591\.6%91\.6\\%Emaxlayer\|REC/Emaxlayer\|Chebyshev\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.7740\.7740\.7820\.7825/55/579\.5%79\.5\\%ELA\|REC/ELA\|Chebyshev\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.8910\.8910\.9540\.9545/55/568\.8%68\.8\\%Problem 3 \(Da=10\.0\\mathrm\{Da\}=10\.0\)E2rel\|REC/E2rel\|Vanilla\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.4330\.4330\.3070\.3075/55/596\.1%96\.1\\%E∞\|REC/E∞\|Vanilla\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.4670\.4670\.3500\.3505/55/592\.7%92\.7\\%Emaxlayer\|REC/Emaxlayer\|Vanilla\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.5840\.5840\.3880\.3885/55/579\.3%79\.3\\%ELA\|REC/ELA\|Vanilla\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{Vanilla\}\}0\.6520\.6520\.5180\.5185/55/579\.2%79\.2\\%E2rel\|REC/E2rel\|Chebyshev\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{2\}^\{\\mathrm\{rel\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.7600\.7600\.7540\.7545/55/587\.9%87\.9\\%E∞\|REC/E∞\|Chebyshev\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\infty\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.6980\.6980\.7050\.7055/55/590\.3%90\.3\\%Emaxlayer\|REC/Emaxlayer\|Chebyshev\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\max\}^\{\\mathrm\{layer\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.6780\.6780\.6810\.6815/55/582\.9%82\.9\\%ELA\|REC/ELA\|Chebyshev\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{REC\}\}/\\left\.E\_\{\\mathrm\{LA\}\}\\right\|\_\{\\mathrm\{Chebyshev\}\}0\.7100\.7100\.6560\.6565/55/572\.8%72\.8\\%Table[3](https://arxiv.org/html/2608.19658#S3.T3)lists error ratios for the concentration entrance problem atNout=16N\_\{\\mathrm\{out\}\}=16andDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0for theRECmodel against both theVanillaandChebyshevmodels, similar to Table[1](https://arxiv.org/html/2608.19658#S3.T1), but with the first three parameter\-bin means computed using the same parameter bins as Fig\.[4](https://arxiv.org/html/2608.19658#S3.F4)\. Relative to theVanillamodel, the second\-column ratios are below unity for all three Damköhler numbers and all four metrics, and the third\-column ratios are smaller than the corresponding second\-column ratios in every row\. The fourth column shows that theRECmodel gives lower mean error than theVanillamodel in all five independent training runs for every metric exceptELAE\_\{\\mathrm\{LA\}\}atDa=0\.1\\mathrm\{Da\}=0\.1and1\.01\.0, where this occurs in four of the five independent training runs, while the fifth column remains above63\.8%63\.8\\%across all rows\. Relative to theChebyshevmodel, all second\-column and third\-column ratios remain below unity forDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0\. The fourth column shows that theRECmodel gives lower mean error than theChebyshevmodel in all five independent training runs for each of the four metrics at each tested Damköhler number, and the fifth column ranges from65\.5%65\.5\\%to91\.6%91\.6\\%\. These results indicate that, for the concentration entrance problem, theRECmodel remains effective for reconstructing the wall\-normal concentration profilec⁡\(y\)c\(y\)as the absorbing\-wall condition changes the near\-wall layer from weakly developed to sharply localized\.

### 3\.4Representative profile behavior and local smoothness

Figure 5:Full representative solution profiles for the three problems atNout=16N\_\{\\mathrm\{out\}\}=16\. \(a\) Scalar profileu⁡\(x\)u\(x\)onx∈\(0​;​1\)x\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)forε=1\.01×10−4\\varepsilon=1\.01\\times 10^\{\-4\}\. \(b\) Temperature profileθ⁡\(y\)\\theta\(y\)ony∈\(0​;​1\)y\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)forPe−1=1\.24×10−4\\mathrm\{Pe\}^\{\-1\}=1\.24\\times 10^\{\-4\}andx=6\.70×10−2x=6\.70\\times 10^\{\-2\}\. \(c\)–\(e\) Concentration profilesc⁡\(y\)c\(y\)ony∈\(0​;​1\)y\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)forDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0, where theDa=0\.1\\mathrm\{Da\}=0\.1case usesPem−1=1\.24×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.24\\times 10^\{\-4\}andx=6\.70×10−2x=6\.70\\times 10^\{\-2\}, while theDa=1\.0\\mathrm\{Da\}=1\.0and10\.010\.0cases usePem−1=1\.15×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.15\\times 10^\{\-4\}andx=8\.32×10−2x=8\.32\\times 10^\{\-2\}\.
Figure[5](https://arxiv.org/html/2608.19658#S3.F5)shows representative cases selected from the held\-out test profiles to visualize the profile behavior underlying the error statistics reported in Sections[3\.1](https://arxiv.org/html/2608.19658#S3.SS1)–[3\.3](https://arxiv.org/html/2608.19658#S3.SS3)\. The figure compares the numerical reference with solution profiles predicted by theVanilla,Chebyshev, andRECmodels atNout=16N\_\{\\mathrm\{out\}\}=16across the three benchmark problems introduced in Section[2\.1](https://arxiv.org/html/2608.19658#S2.SS1)\. For Problems 1 and 2, Figs\.[5](https://arxiv.org/html/2608.19658#S3.F5)\(a\) and[5](https://arxiv.org/html/2608.19658#S3.F5)\(b\) show the scalar solutionu⁡\(x\)u\(x\)overx∈\(0​;​1\)x\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)and the wall\-normal temperature profileθ⁡\(y\)\\theta\(y\)overy∈\(0​;​1\)y\\in\(0\\mathord\{\\mathchar 59\\relax\}1\), respectively, withε=1\.01×10−4\\varepsilon=1\.01\\times 10^\{\-4\}in Fig\.[5](https://arxiv.org/html/2608.19658#S3.F5)\(a\) andPe−1=1\.24×10−4\\mathrm\{Pe\}^\{\-1\}=1\.24\\times 10^\{\-4\}andx=6\.70×10−2x=6\.70\\times 10^\{\-2\}in Fig\.[5](https://arxiv.org/html/2608.19658#S3.F5)\(b\)\. In Fig\.[5](https://arxiv.org/html/2608.19658#S3.F5)\(a\), using theRECmodel decreasesE2relE\_\{2\}^\{\\mathrm\{rel\}\}from7\.90×10−27\.90\\times 10^\{\-2\}with theVanillamodel and2\.49×10−22\.49\\times 10^\{\-2\}with theChebyshevmodel to1\.87×10−21\.87\\times 10^\{\-2\}, and decreasesE∞E\_\{\\infty\}from1\.63×10−11\.63\\times 10^\{\-1\}and3\.03×10−23\.03\\times 10^\{\-2\}to2\.15×10−22\.15\\times 10^\{\-2\}\. In Fig\.[5](https://arxiv.org/html/2608.19658#S3.F5)\(b\),E2relE\_\{2\}^\{\\mathrm\{rel\}\}decreases from4\.45×10−24\.45\\times 10^\{\-2\}with theVanillamodel and2\.41×10−22\.41\\times 10^\{\-2\}with theChebyshevmodel to1\.63×10−21\.63\\times 10^\{\-2\}, andE∞E\_\{\\infty\}decreases from1\.59×10−11\.59\\times 10^\{\-1\}and8\.01×10−28\.01\\times 10^\{\-2\}to5\.91×10−25\.91\\times 10^\{\-2\}\. These reductions show that theRECmodel is much closer to the numerical reference than theVanillamodel and still improves over theChebyshevmodel in the representative scalar and thermal entrance profiles\.

For Problem 3, Figs\.[5](https://arxiv.org/html/2608.19658#S3.F5)\(c\),[5](https://arxiv.org/html/2608.19658#S3.F5)\(d\), and[5](https://arxiv.org/html/2608.19658#S3.F5)\(e\) show the wall\-normal concentration profilesc⁡\(y\)c\(y\)overy∈\(0​;​1\)y\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)as the wall\-absorption strengthDa\\mathrm\{Da\}increases from0\.10\.1to10\.010\.0, withPem−1=1\.24×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.24\\times 10^\{\-4\}andx=6\.70×10−2x=6\.70\\times 10^\{\-2\}forDa=0\.1\\mathrm\{Da\}=0\.1, andPem−1=1\.15×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.15\\times 10^\{\-4\}andx=8\.32×10−2x=8\.32\\times 10^\{\-2\}forDa=1\.0\\mathrm\{Da\}=1\.0and10\.010\.0\. The wall\-endpoint drop of the numerical reference, measured byc⁡\(y=0\.99\)−c⁡\(y=1\)c\(y=0\.99\)\-c\(y=1\), increases from9\.44×10−49\.44\\times 10^\{\-4\}atDa=0\.1\\mathrm\{Da\}=0\.1to1\.43×10−21\.43\\times 10^\{\-2\}atDa=1\.0\\mathrm\{Da\}=1\.0and1\.19×10−11\.19\\times 10^\{\-1\}atDa=10\.0\\mathrm\{Da\}=10\.0, showing the progressive formation of a sharper wall\-adjacent concentration layer asDa\\mathrm\{Da\}increases\. Relative to theVanillamodel, theRECmodel decreasesE2relE\_\{2\}^\{\\mathrm\{rel\}\}from1\.13×10−21\.13\\times 10^\{\-2\},2\.74×10−22\.74\\times 10^\{\-2\}, and2\.48×10−22\.48\\times 10^\{\-2\}to5\.63×10−35\.63\\times 10^\{\-3\},1\.22×10−21\.22\\times 10^\{\-2\}, and1\.82×10−21\.82\\times 10^\{\-2\}, respectively, and decreasesE∞E\_\{\\infty\}from2\.30×10−22\.30\\times 10^\{\-2\},7\.00×10−27\.00\\times 10^\{\-2\}, and7\.64×10−27\.64\\times 10^\{\-2\}to1\.56×10−21\.56\\times 10^\{\-2\},3\.10×10−23\.10\\times 10^\{\-2\}, and5\.69×10−25\.69\\times 10^\{\-2\}\. Relative to theChebyshevmodel,E2relE\_\{2\}^\{\\mathrm\{rel\}\}decreases from1\.17×10−21\.17\\times 10^\{\-2\},1\.88×10−21\.88\\times 10^\{\-2\}, and2\.16×10−22\.16\\times 10^\{\-2\}to the sameRECvalues\. ForDa=0\.1\\mathrm\{Da\}=0\.1and1\.01\.0,E∞E\_\{\\infty\}decreases from2\.83×10−22\.83\\times 10^\{\-2\}to1\.56×10−21\.56\\times 10^\{\-2\}and from5\.88×10−25\.88\\times 10^\{\-2\}to3\.10×10−23\.10\\times 10^\{\-2\}, while theDa=10\.0\\mathrm\{Da\}=10\.0value remains5\.69×10−25\.69\\times 10^\{\-2\}for both models\. Even atDa=0\.1\\mathrm\{Da\}=0\.1, where the wall\-endpoint drop is weak, theChebyshevmodel still introduces visible near\-wall oscillation, whereas theRECmodel follows the nearly flat reference profile more closely\. These results indicate that theRECmodel predicts the numerical reference substantially more accurately than theVanillamodel and also improves the prediction relative to theChebyshevmodel by reducing near\-wall oscillatory deviations, rather than only by resolving an extremely thin localized layer\.

Figure 6:Near\-boundary and near\-wall solution profiles for the three problems atNout=16N\_\{\\mathrm\{out\}\}=16\. \(a\) Scalar profileu⁡\(x\)u\(x\)near the right endpoint forε=1\.01×10−4\\varepsilon=1\.01\\times 10^\{\-4\}\. \(b\) Temperature profileθ⁡\(y\)\\theta\(y\)near the wall forPe−1=1\.24×10−4\\mathrm\{Pe\}^\{\-1\}=1\.24\\times 10^\{\-4\}andx=6\.70×10−2x=6\.70\\times 10^\{\-2\}\. \(c\)–\(e\) Concentration profilesc⁡\(y\)c\(y\)near the absorbing wall forDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0, where theDa=0\.1\\mathrm\{Da\}=0\.1case usesPem−1=1\.24×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.24\\times 10^\{\-4\}andx=6\.70×10−2x=6\.70\\times 10^\{\-2\}, while theDa=1\.0\\mathrm\{Da\}=1\.0and10\.010\.0cases usePem−1=1\.15×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.15\\times 10^\{\-4\}andx=8\.32×10−2x=8\.32\\times 10^\{\-2\}\. Each comparison includes the numerical reference, Vanilla, Chebyshev, and REC\.
Figure[6](https://arxiv.org/html/2608.19658#S3.F6)enlarges the right\-endpoint layer inu⁡\(x\)u\(x\), the near\-wall part of the temperature profileθ⁡\(y\)\\theta\(y\), and the near\-wall part of the concentration profilesc⁡\(y\)c\(y\)atDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0for the five representative cases shown in Fig\.[5](https://arxiv.org/html/2608.19658#S3.F5)\. In these enlarged regions, the predictions from theChebyshevmodel often retain local oscillatory deviations from the numerical reference, whereas the predictions from theRECmodel reduce these deviations while preserving the nearby smooth profile shape\. To quantify this local oscillatory behavior, the second\-difference roughnessW2W\_\{2\}is evaluated for the same representative profiles\. Profiles with stronger oscillatory variation have larger second differences, whileW2W\_\{2\}values closer to the numerical\-reference value indicate closer agreement in local smoothness\. For a discrete profile vector𝐯=\(v1​;​…​;​vm\)\\mathbf\{v\}=\(v\_\{1\}\\mathord\{\\mathchar 59\\relax\}\\ldots\\mathord\{\\mathchar 59\\relax\}v\_\{m\}\)evaluated at them=257m=257Chebyshev–Lobatto output locations defined in Section[2\.4](https://arxiv.org/html/2608.19658#S2.SS4), withvj=u⁡\(xj\)v\_\{j\}=u\(x\_\{j\}\)for Problem 1,vj=θ⁡\(yj\)v\_\{j\}=\\theta\(y\_\{j\}\)for Problem 2, andvj=c⁡\(yj\)v\_\{j\}=c\(y\_\{j\}\)for Problem 3, wherejjindexes the output locations,

W2​\(𝐯\)=∑j=2m−1\|vj\+1−2​vj\+vj−1\|;W\_\{2\}\(\\mathbf\{v\}\)=\\sum\_\{j=2\}^\{m\-1\}\\left\|v\_\{j\+1\}\-2v\_\{j\}\+v\_\{j\-1\}\\right\|\\mathchar 59\\relax\(28\)where\|⋅\|\|\\cdot\|denotes the scalar absolute value\.

Table 4:W2W\_\{2\}values for the numerical reference and for the predictions from theChebyshevandRECmodels in the five representative cases shown in Fig\.[6](https://arxiv.org/html/2608.19658#S3.F6), computed overx∈\(0​;​1\)x\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)for Problem 1 and overy∈\(0​;​1\)y\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)for Problems 2 and 3\.ProblemParameterW2\|Reference\\left\.W\_\{2\}\\right\|\_\{\\mathrm\{Reference\}\}W2\|Chebyshev\\left\.W\_\{2\}\\right\|\_\{\\mathrm\{Chebyshev\}\}W2\|REC\\left\.W\_\{2\}\\right\|\_\{\\mathrm\{REC\}\}Problem 1ε=1\.01×10−4\\varepsilon=1\.01\\times 10^\{\-4\}1\.87×10−11\.87\\times 10^\{\-1\}9\.72×10−19\.72\\times 10^\{\-1\}2\.07×10−12\.07\\times 10^\{\-1\}Problem 2Pe−1=1\.24×10−4\\mathrm\{Pe\}^\{\-1\}=1\.24\\times 10^\{\-4\}x=6\.70×10−2x=6\.70\\times 10^\{\-2\}1\.29×10−11\.29\\times 10^\{\-1\}8\.23×10−18\.23\\times 10^\{\-1\}1\.35×10−11\.35\\times 10^\{\-1\}Problem 3\(Da=0\.1\\mathrm\{Da\}=0\.1\)Pem−1=1\.24×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.24\\times 10^\{\-4\}x=6\.70×10−2x=6\.70\\times 10^\{\-2\}5\.15×10−25\.15\\times 10^\{\-2\}8\.71×10−18\.71\\times 10^\{\-1\}6\.78×10−26\.78\\times 10^\{\-2\}Problem 3\(Da=1\.0\\mathrm\{Da\}=1\.0\)Pem−1=1\.15×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.15\\times 10^\{\-4\}x=8\.32×10−2x=8\.32\\times 10^\{\-2\}6\.39×10−26\.39\\times 10^\{\-2\}1\.43×1001\.43\\times 10^\{0\}9\.26×10−29\.26\\times 10^\{\-2\}Problem 3\(Da=10\.0\\mathrm\{Da\}=10\.0\)Pem−1=1\.15×10−4\\mathrm\{Pe\}\_\{m\}^\{\-1\}=1\.15\\times 10^\{\-4\}x=8\.32×10−2x=8\.32\\times 10^\{\-2\}9\.22×10−29\.22\\times 10^\{\-2\}1\.63×1001\.63\\times 10^\{0\}9\.63×10−29\.63\\times 10^\{\-2\}Table[4](https://arxiv.org/html/2608.19658#S3.T4)shows that theW2W\_\{2\}values from theRECmodel remain closer to the numerical\-reference values than theW2W\_\{2\}values from theChebyshevmodel in all five representative cases\. For the thermal entrance problem,W2W\_\{2\}decreases by83\.6%83\.6\\%from theChebyshevmodel to theRECmodel\. For the concentration entrance problem, the corresponding decreases are92\.2%92\.2\\%,93\.5%93\.5\\%, and94\.1%94\.1\\%atDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0, respectively, giving larger decreases than in the thermal entrance problem\. More importantly, for the singularly perturbed scalar BVP, although theRECmodel does not substantially improve the scalar\-profile prediction accuracy relative to theChebyshevmodel, it decreasesW2W\_\{2\}by78\.7%78\.7\\%\. In addition, for all held\-out test profiles over the five independent training runs, the medianW2W\_\{2\}decreases from theChebyshevmodel to theRECmodel by64\.0%64\.0\\%for the singularly perturbed scalar BVP and by81\.4%81\.4\\%for the thermal entrance problem\. For the concentration entrance problem, the corresponding median decreases are87\.2%87\.2\\%,87\.6%87\.6\\%, and85\.9%85\.9\\%atDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0, respectively\. These results indicate that theRECmodel suppresses the artificial local reversals introduced by theChebyshevmodel near the boundary or wall, thereby giving a near\-boundary or near\-wall profile shape closer to the numerical\-reference behavior across both the five representative cases and all held\-out test profiles\.

### 3\.5Dependence on outer Chebyshev subdictionary size

Sections[3\.1](https://arxiv.org/html/2608.19658#S3.SS1)–[3\.4](https://arxiv.org/html/2608.19658#S3.SS4)show that, first, compared with theVanillamodel, theRECmodel gives lower errors in predictingu⁡\(x\)u\(x\),θ⁡\(y\)\\theta\(y\), andc⁡\(y\)c\(y\), even though theVanillamodel uses a learned, parameter\-conditioned trunk MLP\. This result indicates that the learned coordinate trunk does not by itself provide the same layer\-aligned output functions as the prescribed REC dictionary for the tested bounded\-domain solution profiles\. Compared with theChebyshevmodel, theRECmodel changes the prescribed trunk dictionary while keeping the branch network, trunk dimension, data splits, loss function, optimizer, training procedure, and evaluation protocol fixed, as described in Section[2](https://arxiv.org/html/2608.19658#S2)\. Therefore, the lower errors in the thermal and concentration entrance problems, together with the smallerW2W\_\{2\}values in the representative profiles, support replacing part of the global Chebyshev dictionary by inner\-scale rational dictionary elements in the output\-coordinate representation, rather than adding trainable trunk expressivity, changing the optimization procedure, or imposing a problem\-specific decomposition, rescaling, or sampling rule\. However, the results in Sections[3\.1](https://arxiv.org/html/2608.19658#S3.SS1)–[3\.4](https://arxiv.org/html/2608.19658#S3.SS4)are obtained with the specific splitNout=16N\_\{\\mathrm\{out\}\}=16and do not yet show whether the lower errors in predictingθ⁡\(y\)\\theta\(y\)andc⁡\(y\)c\(y\)can be attributed to using a larger number of inner rational dictionary elements within thep=129p=129trunk dictionary, or whether similar errors would be obtained with a larger outer Chebyshev subdictionary\. They also do not show whether the lower errors relative to theVanillamodel remain whenNoutN\_\{\\mathrm\{out\}\}is increased to3333,6565, and9797\.

Figure 7:Ratios ofE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}from theRECmodel to the corresponding errors from theChebyshevmodel forNout=16N\_\{\\mathrm\{out\}\}=16,3333,6565, and9797\.Figure 8:Ratios ofE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}from theRECmodel to the corresponding errors from theVanillamodel forNout=16N\_\{\\mathrm\{out\}\}=16,3333,6565, and9797\.Figures[7](https://arxiv.org/html/2608.19658#S3.F7)and[8](https://arxiv.org/html/2608.19658#S3.F8)accordingly compare the error from theRECmodel against the errors from theChebyshevandVanillamodels asNoutN\_\{\\mathrm\{out\}\}is varied over1616,3333,6565, and9797\. For each metric and each problem, the solid line gives the median ratio over the five independent training runs and test profiles, and the shaded band gives the IQR\. In Fig\.[7](https://arxiv.org/html/2608.19658#S3.F7), for the singularly perturbed scalar BVP, the ratios remain close to unity and can exceed unity forEmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}andELAE\_\{\\mathrm\{LA\}\}, consistent with Table[1](https://arxiv.org/html/2608.19658#S3.T1), where theChebyshevmodel gives lower errors than theRECmodel on these two metrics\. For the thermal and concentration entrance problems in Fig\.[7](https://arxiv.org/html/2608.19658#S3.F7), the smallest ratios occur atNout=16N\_\{\\mathrm\{out\}\}=16, and increasingNoutN\_\{\\mathrm\{out\}\}moves the ratios toward unity\. In Fig\.[8](https://arxiv.org/html/2608.19658#S3.F8), the ratios remain below unity for the thermal and concentration entrance problems and remain below unity for most singularly perturbed scalar BVP cases\. Thus, the weaker improvement over theChebyshevmodel at largerNoutN\_\{\\mathrm\{out\}\}supports the interpretation that retaining more inner rational dictionary elements is more favorable for predictingθ⁡\(y\)\\theta\(y\)andc⁡\(y\)c\(y\)than increasing the number of outer Chebyshev dictionary elements\. The ratios in Fig\.[8](https://arxiv.org/html/2608.19658#S3.F8)also show that the lower errors relative to theVanillamodel are retained whenNoutN\_\{\\mathrm\{out\}\}is increased\.

### 3\.6Practical implications

Figure 9:Averaged values ofE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}for theVanilla,Chebyshev, andRECmodels atNout=16N\_\{\\mathrm\{out\}\}=16\.Figure[9](https://arxiv.org/html/2608.19658#S3.F9)compares the averaged values ofE2relE\_\{2\}^\{\\mathrm\{rel\}\},E∞E\_\{\\infty\},EmaxlayerE\_\{\\max\}^\{\\mathrm\{layer\}\}, andELAE\_\{\\mathrm\{LA\}\}for theVanilla,Chebyshev, andRECmodels atNout=16N\_\{\\mathrm\{out\}\}=16\. For each model and each problem, the plotted value is the corresponding error metric averaged over the five independent training runs and the test profiles\. These averaged error values summarize the observations in Sections[3\.1](https://arxiv.org/html/2608.19658#S3.SS1)–[3\.4](https://arxiv.org/html/2608.19658#S3.SS4)\. For the singularly perturbed scalar BVP, theRECmodel predicts the scalar profileu⁡\(x\)u\(x\)more accurately than theVanillamodel and remains comparable to theChebyshevmodel\. For the thermal and concentration entrance problems, theRECmodel gives more accurate predictions of the temperature profileθ⁡\(y\)\\theta\(y\)and the concentration profilec⁡\(y\)c\(y\)than both comparison models\. TheNoutN\_\{\\mathrm\{out\}\}sweep in Section[3\.5](https://arxiv.org/html/2608.19658#S3.SS5)further shows that, in these two entrance problems, the lower errors in predictingθ⁡\(y\)\\theta\(y\)andc⁡\(y\)c\(y\)are strengthened by using more inner rational dictionary elements within the REC dictionary\. Therefore, the practical implication is clearest for REC\-type trunk dictionaries in operator surrogates that predict wall\-normal temperature or concentration profiles with a smooth outer part and a near\-wall layer\.

This advantage is expected to be particularly relevant to high\-Péclet entrance\-region heat and mass transfer problems in which repeated evaluations are required as inlet profiles or wall conditions change\. For example, surface\-based biosensors and biosensors based on surface capture require concentration profiles governed by convection, diffusion, reaction, and binding near reactive surfaces[squires2008making](https://arxiv.org/html/2608.19658#bib.bib16)\. Similarly, microfluidic systems with surface reactions and microfluidic electrochemical chips with electrodes on the side walls involve concentration fields formed by laminar convective\-diffusive transport and electrochemical reaction at an electrode interface[gervais2006mass](https://arxiv.org/html/2608.19658#bib.bib17),[chevalier2021semianalytical](https://arxiv.org/html/2608.19658#bib.bib18)\. In high\-speed flow applications, high\-speed boundary layers with wall transpiration, transpiration cooling through a porous wall, hypersonic transitional and turbulent boundary layers, hypersonic turbulent boundary layers with finite\-rate chemical reactions, and reacting boundary layers with recombination reactions also require accurate near\-wall temperature or species profiles[sescu2019transpiration](https://arxiv.org/html/2608.19658#bib.bib19),[hillcoat2025transpiration](https://arxiv.org/html/2608.19658#bib.bib20),[xu2022hypersonic](https://arxiv.org/html/2608.19658#bib.bib21),[passiatore2021finite](https://arxiv.org/html/2608.19658#bib.bib22),[perakis2021recombination](https://arxiv.org/html/2608.19658#bib.bib23)\.

## 4Conclusion

This study examined trunk\-basis design in DeepONet surrogates for singularly perturbed and high\-Péclet transport operators on bounded domains, with particular emphasis on wall\-normal profile reconstruction in entrance\-region heat and mass transfer\. The proposed REC trunk dictionary combines a low\-degree outer Chebyshev subdictionary with an inner rational subdictionary constructed before DeepONet training from a canonical exponentially decaying layer family\. The key conclusions are as follows:

1. 1\.On the singularly perturbed scalar BVP, theRECmodel gives its clearest reductions in the smallest\-parameter regime, while remaining comparable to theChebyshevmodel over the full test set\.
2. 2\.On the thermal entrance problem, theRECmodel improves the reconstruction of the wall\-normal temperature profileθ⁡\(y\)\\theta\(y\)relative to both theVanillaandChebyshevmodels across the tested high\-Péclet regime and sampled entrance locationsxx\.
3. 3\.On the concentration entrance problem, theRECmodel improves the reconstruction of the wall\-normal concentration profilec⁡\(y\)c\(y\)acrossDa=0\.1\\mathrm\{Da\}=0\.1,1\.01\.0, and10\.010\.0, showing that the improvement persists as the absorbing\-wall condition changes the near\-wall layer\.
4. 4\.The reductions obtained with theRECmodel are consistent across the five independent training runs and the held\-out test profiles, indicating that the observed improvements are not restricted to a single training initialization or a small set of selected profiles\.

## Acknowledgments

This work was supported by the Center for Heterogeneous Integration of Micro Electronic Systems \(CHIMES\), one of the seven centers sponsored by the Semiconductor Research Corporation \(SRC\) and the Defense Advanced Research Projects Agency \(DARPA\) under the Joint University Microelectronics Program 2\.0 \(JUMP 2\.0\)\.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work presented in this paper\.

## Data availability

The code and data used in this work are available upon request\.

## Appendix ATheoretical motivation for the REC trunk dictionary

This appendix develops an idealized representation argument for the REC trunk, rather than a convergence theorem for the trained DeepONet\. The argument uses the canonical exponential layer family introduced in Eq\. \([16](https://arxiv.org/html/2608.19658#S2.E16)\), whose sampled instances define the inner rational subdictionary in the REC dictionary construction before training\. Treating the same family as an idealized target clarifies how polynomial and rational trunks represent a thin boundary\-attached layer, and why the REC approximation space combining a low\-degree outer Chebyshev subdictionary with the inner rational subdictionary is structurally compatible with the decomposition into a smooth outer solution and a one\-sided layer\. Define

s=2​ζ−1∈\[−1​;​1\];β=12​δ;s=2\\zeta\-1\\in\[\-1\\mathord\{\\mathchar 59\\relax\}1\]\\mathchar 59\\relax\\qquad\\beta=\\frac\{1\}\{2\\delta\}\\mathchar 59\\relax\(29\)wheressdenotes the shifted coordinate andβ\\betadenotes the rescaled layer parameter\. Then

ψδ​\(ζ\)=e−β​eβ​s\.\\psi\_\{\\delta\}\(\\zeta\)=e^\{\-\\beta\}e^\{\\beta s\}\.\(30\)Using the modified\-Bessel Jacobi\-Anger identity[dlmf1035](https://arxiv.org/html/2608.19658#bib.bib31),

eβ​s=I0​\(β\)\+2​∑k=1∞Ik​\(β\)​Tk​\(s\);e^\{\\beta s\}=I\_\{0\}\(\\beta\)\+2\\sum\_\{k=1\}^\{\\infty\}I\_\{k\}\(\\beta\)T\_\{k\}\(s\)\\mathchar 59\\relax\(31\)whereIkI\_\{k\}denotes the modified Bessel function of the first kind of orderkkandTkT\_\{k\}denotes the Chebyshev polynomial of the first kind of degreekk, one obtains the exact Chebyshev expansion

ψδ​\(ζ\)=e−β​I0​\(β\)\+2​e−β​∑k=1∞Ik​\(β\)​Tk​\(2​ζ−1\)\.\\psi\_\{\\delta\}\(\\zeta\)=e^\{\-\\beta\}I\_\{0\}\(\\beta\)\+2e^\{\-\\beta\}\\sum\_\{k=1\}^\{\\infty\}I\_\{k\}\(\\beta\)\\,T\_\{k\}\(2\\zeta\-1\)\.\(32\)Hence the Chebyshev coefficients are

c0​\(δ\)=e−β​I0​\(β\);ck​\(δ\)=2​e−β​Ik​\(β\);k≥1;c\_\{0\}\(\\delta\)=e^\{\-\\beta\}I\_\{0\}\(\\beta\)\\mathchar 59\\relax\\qquad c\_\{k\}\(\\delta\)=2e^\{\-\\beta\}I\_\{k\}\(\\beta\)\\mathchar 59\\relax\\quad k\\geq 1\\mathchar 59\\relax\(33\)whereck​\(δ\)c\_\{k\}\(\\delta\)denotes the coefficient ofTk​\(2​ζ−1\)T\_\{k\}\(2\\zeta\-1\)\. For the degree\-nnChebyshev truncation

Sn​ψδ​\(ζ\)=c0​\(δ\)\+∑k=1nck​\(δ\)​Tk​\(2​ζ−1\);S\_\{n\}\\psi\_\{\\delta\}\(\\zeta\)=c\_\{0\}\(\\delta\)\+\\sum\_\{k=1\}^\{n\}c\_\{k\}\(\\delta\)\\,T\_\{k\}\(2\\zeta\-1\)\\mathchar 59\\relax\(34\)wherenndenotes the highest retained Chebyshev degree, the uniform truncation error satisfies

‖ψδ−Sn​ψδ‖∞≤2​e−β​∑k=n\+1∞Ik​\(β\);\\\|\\psi\_\{\\delta\}\-S\_\{n\}\\psi\_\{\\delta\}\\\|\_\{\\infty\}\\leq 2e^\{\-\\beta\}\\sum\_\{k=n\+1\}^\{\\infty\}I\_\{k\}\(\\beta\)\\mathchar 59\\relax\(35\)where∥⋅∥∞\\\|\\cdot\\\|\_\{\\infty\}denotes the supremum norm on\[0​;​1\]\[0\\mathord\{\\mathchar 59\\relax\}1\]\. This shows that, asδ\\deltadecreases, appreciable Chebyshev weight shifts to higher polynomial degrees\. A low\-degree polynomial trunk therefore becomes progressively less well matched to this canonical boundary layer family\.

A more quantitative scale estimate follows from the local central\-limit asymptotic for modified Bessel functions[athreya1987bessel](https://arxiv.org/html/2608.19658#bib.bib32),

e−βIk\(β\)∼12​π​βexp\(−k22​β\);β→∞;kβ→γ∈\[0;∞\);e^\{\-\\beta\}I\_\{k\}\(\\beta\)\\sim\\frac\{1\}\{\\sqrt\{2\\pi\\beta\}\}\\exp\\\!\\left\(\-\\frac\{k^\{2\}\}\{2\\beta\}\\right\)\\mathchar 59\\relax\\qquad\\beta\\to\\infty\\mathchar 59\\relax\\qquad\\frac\{k\}\{\\sqrt\{\\beta\}\}\\to\\gamma\\in\[0\\mathord\{\\mathchar 59\\relax\}\\infty\)\\mathchar 59\\relax\(36\)whereγ\\gammadenotes the limiting scaled polynomial degree\. This identifiesk=O⁡\(β\)k=O\(\\sqrt\{\\beta\}\)as the central coefficient scale\. Sinceβ=1/\(2​δ\)\\beta=1/\(2\\delta\), this gives, fork≥1k\\geq 1,

ck​\(δ\)∼2​δπ​exp⁡\(−k2​δ\);c0​\(δ\)∼δπ\.c\_\{k\}\(\\delta\)\\sim 2\\sqrt\{\\frac\{\\delta\}\{\\pi\}\}\\exp\(\-k^\{2\}\\delta\)\\mathchar 59\\relax\\qquad c\_\{0\}\(\\delta\)\\sim\\sqrt\{\\frac\{\\delta\}\{\\pi\}\}\.\(37\)Thus, the envelope of the Chebyshev coefficients has a Gaussian scale inkk\. Ifρ∈\(0​;​1\)\\rho\\in\(0\\mathord\{\\mathchar 59\\relax\}1\)denotes a fixed relative decay level of this envelope andKρ​\(δ\)K\_\{\\rho\}\(\\delta\)denotes the polynomial degree at whichexp⁡\(−k2​δ\)=ρ\\exp\(\-k^\{2\}\\delta\)=\\rho, then

Kρ​\(δ\)=\|log⁡ρ\|δ\.K\_\{\\rho\}\(\\delta\)=\\sqrt\{\\frac\{\|\\log\\rho\|\}\{\\delta\}\}\.\(38\)Therefore, the active polynomial degree scale is proportional toδ−1/2\\delta^\{\-1/2\}asδ→0\\delta\\to 0\. For example, takingρ=e−1\\rho=e^\{\-1\}givesKρ​\(10−4\)=100K\_\{\\rho\}\(10^\{\-4\}\)=100, which explains why a degree128128Chebyshev trunk can remain a strong baseline for the pure canonical exponential layer\. The inner rational subdictionary is therefore not motivated by a claim that it must dominate a sufficiently high degree Chebyshev expansion on this scalar layer alone\. Rather, it is motivated as a compact layer\-aligned inner rational subdictionary that can complement a low\-degree outer Chebyshev subdictionary in the profile\-valued entrance transport problems studied in the main text\.

Next define the logarithmic parameter

θ=log10⁡δ;\\theta=\\log\_\{10\}\\delta\\mathchar 59\\relax\(39\)whereθ\\thetadenotes the log scaled prototype parameter\. Then, withδ=10θ\\delta=10^\{\\theta\},

∂θψ10θ​\(ζ\)=\(ln⁡10\)​\(1−ζδ\)​exp⁡\(−1−ζδ\)\.\\partial\_\{\\theta\}\\psi\_\{10^\{\\theta\}\}\(\\zeta\)=\(\\ln 10\)\\left\(\\frac\{1\-\\zeta\}\{\\delta\}\\right\)\\exp\\\!\\left\(\-\\frac\{1\-\\zeta\}\{\\delta\}\\right\)\.\(40\)If

t=1−ζδ≥0;t=\\frac\{1\-\\zeta\}\{\\delta\}\\geq 0\\mathchar 59\\relax\(41\)wherettdenotes the stretched layer coordinate associated with the canonical exponential profile, then

∂θψ10θ​\(ζ\)=\(ln⁡10\)​t​e−t\.\\partial\_\{\\theta\}\\psi\_\{10^\{\\theta\}\}\(\\zeta\)=\(\\ln 10\)\\,te^\{\-t\}\.\(42\)Sincet​e−t≤e−1te^\{\-t\}\\leq e^\{\-1\}for allt≥0t\\geq 0,

supθ‖∂θψ10θ‖∞≤ln⁡10e;\\sup\_\{\\theta\}\\left\\\|\\partial\_\{\\theta\}\\psi\_\{10^\{\\theta\}\}\\right\\\|\_\{\\infty\}\\leq\\frac\{\\ln 10\}\{e\}\\mathchar 59\\relax\(43\)where the supremum is taken over the log parameter range considered for the layer family\.

Let\{ri\}i=1M\\\{r\_\{i\}\\\}\_\{i=1\}^\{M\}denote the rational dictionary elements, and letδi=10θi\\delta\_\{i\}=10^\{\\theta\_\{i\}\}denote the prototype parameter value associated withrir\_\{i\}\. Suppose that

‖ψδi−ri‖∞≤τi;\\\|\\psi\_\{\\delta\_\{i\}\}\-r\_\{i\}\\\|\_\{\\infty\}\\leq\\tau\_\{i\}\\mathchar 59\\relax\(44\)whereτi\\tau\_\{i\}denotes a uniform approximation error bound for the rational dictionary elementrir\_\{i\}on\[0​;​1\]\[0\\mathord\{\\mathchar 59\\relax\}1\]\. It should not be interpreted as the raw AAA stopping tolerance unless that tolerance has been separately verified to bound the uniform approximation error\.

Then, for anyδ=10θ\\delta=10^\{\\theta\},

‖ψδ−ri‖∞≤τi\+ln⁡10e​\|θ−θi\|=τi\+ln⁡10e​\|log10⁡δ−log10⁡δi\|\.\\\|\\psi\_\{\\delta\}\-r\_\{i\}\\\|\_\{\\infty\}\\leq\\tau\_\{i\}\+\\frac\{\\ln 10\}\{e\}\\,\|\\theta\-\\theta\_\{i\}\|=\\tau\_\{i\}\+\\frac\{\\ln 10\}\{e\}\\left\|\\log\_\{10\}\\delta\-\\log\_\{10\}\\delta\_\{i\}\\right\|\.\(45\)Therefore,

infv∈ℛM‖ψδ−v‖∞≤min1≤i≤M⁡\(τi\+ln⁡10e​\|log10⁡δ−log10⁡δi\|\);\\inf\_\{v\\in\\mathcal\{R\}\_\{M\}\}\\\|\\psi\_\{\\delta\}\-v\\\|\_\{\\infty\}\\leq\\min\_\{1\\leq i\\leq M\}\\left\(\\tau\_\{i\}\+\\frac\{\\ln 10\}\{e\}\\left\|\\log\_\{10\}\\delta\-\\log\_\{10\}\\delta\_\{i\}\\right\|\\right\)\\mathchar 59\\relax\(46\)where

ℛM=span⁡\{r1​;​…​;​rM\}\.\\mathcal\{R\}\_\{M\}=\\mathrm\{span\}\\\{r\_\{1\}\\mathord\{\\mathchar 59\\relax\}\\dots\\mathord\{\\mathchar 59\\relax\}r\_\{M\}\\\}\.\(47\)Here,ℛM\\mathcal\{R\}\_\{M\}denotes the rational subspace,MMdenotes the number of rational dictionary elements, andvvdenotes a candidate approximating function inℛM\\mathcal\{R\}\_\{M\}\.

Consequently, ifℐθ=\[θmin​;​θmax\]\\mathcal\{I\}\_\{\\theta\}=\[\\theta\_\{\\min\}\\mathord\{\\mathchar 59\\relax\}\\theta\_\{\\max\}\]denotes the sampled log parameter interval, and if

hθ=supθ∈ℐθmin1≤i≤M⁡\|θ−θi\|h\_\{\\theta\}=\\sup\_\{\\theta\\in\\mathcal\{I\}\_\{\\theta\}\}\\min\_\{1\\leq i\\leq M\}\|\\theta\-\\theta\_\{i\}\|\(48\)denotes the fill distance of the sampled log grid, then

infv∈ℛM‖ψ10θ−v‖∞≤τmax\+ln⁡10e​hθ;θ∈ℐθ;\\inf\_\{v\\in\\mathcal\{R\}\_\{M\}\}\\\|\\psi\_\{10^\{\\theta\}\}\-v\\\|\_\{\\infty\}\\leq\\tau\_\{\\max\}\+\\frac\{\\ln 10\}\{e\}h\_\{\\theta\}\\mathchar 59\\relax\\qquad\\theta\\in\\mathcal\{I\}\_\{\\theta\}\\mathchar 59\\relax\(49\)whereτmax=max1≤i≤M⁡τi\\tau\_\{\\max\}=\\max\_\{1\\leq i\\leq M\}\\tau\_\{i\}denotes the largest assumed uniform approximation error bound over the inner rational subdictionary\. This is a log parameter coverage statement for the idealized layer family, not a statement about the trained DeepONet optimization error\.

Finally, define

𝒫m−1=span⁡\{T0​\(2​ζ−1\)​;​…​;​Tm−1​\(2​ζ−1\)\};ℋm​;​M=𝒫m−1\+ℛM;\\mathcal\{P\}\_\{m\-1\}=\\mathrm\{span\}\\\{T\_\{0\}\(2\\zeta\-1\)\\mathord\{\\mathchar 59\\relax\}\\dots\\mathord\{\\mathchar 59\\relax\}T\_\{m\-1\}\(2\\zeta\-1\)\\\}\\mathchar 59\\relax\\qquad\\mathcal\{H\}\_\{m\\mathord\{\\mathchar 59\\relax\}M\}=\\mathcal\{P\}\_\{m\-1\}\+\\mathcal\{R\}\_\{M\}\\mathchar 59\\relax\(50\)wheremmdenotes the number of outer Chebyshev dictionary elements andℋm​;​M\\mathcal\{H\}\_\{m\\mathord\{\\mathchar 59\\relax\}M\}denotes the corresponding REC approximation space\.

For theRECmodel used in the main comparison,m=16m=16andM=113M=113\. The Chebyshev trunk with the same nominal dimension corresponds to𝒫128\\mathcal\{P\}\_\{128\}, which is obtained by takingm=129m=129in the definition of𝒫m−1\\mathcal\{P\}\_\{m\-1\}\. In general,ℋ16​;​113\\mathcal\{H\}\_\{16\\mathord\{\\mathchar 59\\relax\}113\}and𝒫128\\mathcal\{P\}\_\{128\}are different approximation spaces, and neither space comparison alone implies that one trained DeepONet must outperform the other\. The following estimate is therefore an existence bound for the REC approximation space, not a dominance theorem relative to the full Chebyshev space\.

For a target field of the form

f⁡\(ζ\)=q⁡\(ζ\)\+A​ψδ​\(ζ\);f\(\\zeta\)=q\(\\zeta\)\+A\\psi\_\{\\delta\}\(\\zeta\)\\mathchar 59\\relax\(51\)whereffdenotes a target profile,qqdenotes a smooth outer component, andAAdenotes a layer amplitude, let

ηm​\(q\):=infπ∈𝒫m−1‖q−π‖∞;\\eta\_\{m\}\(q\):=\\inf\_\{\\pi\\in\\mathcal\{P\}\_\{m\-1\}\}\\\|q\-\\pi\\\|\_\{\\infty\}\\mathchar 59\\relax\(52\)whereπ\\pidenotes a candidate polynomial in𝒫m−1\\mathcal\{P\}\_\{m\-1\}, andηm​\(q\)\\eta\_\{m\}\(q\)denotes the best uniform approximation error of the outer componentqqby the outer Chebyshev subdictionary\. Then

infv∈ℋm​;​M‖f−v‖∞≤ηm​\(q\)\+\|A\|​min1≤i≤M⁡\(τi\+ln⁡10e​\|log10⁡δ−log10⁡δi\|\)\.\\inf\_\{v\\in\\mathcal\{H\}\_\{m\\mathord\{\\mathchar 59\\relax\}M\}\}\\\|f\-v\\\|\_\{\\infty\}\\leq\\eta\_\{m\}\(q\)\+\|A\|\\min\_\{1\\leq i\\leq M\}\\left\(\\tau\_\{i\}\+\\frac\{\\ln 10\}\{e\}\\left\|\\log\_\{10\}\\delta\-\\log\_\{10\}\\delta\_\{i\}\\right\|\\right\)\.\(53\)Although this bound does not prove that the trainedRECmodel must produce lower errors than the trainedVanillaorChebyshevmodels in the comparisons of Section[3](https://arxiv.org/html/2608.19658#S3), it shows that, under the idealized decompositionf=q\+A​ψδf=q\+A\\psi\_\{\\delta\}, the REC trunk is structurally aligned with a smooth outer component and a thin boundary\-attached layer, while the inner rational subdictionary provides an error estimate controlled by the spacing of the sampled prototype parameter nodes on thelog10⁡δ\\log\_\{10\}\\deltaaxis for the canonical layer family\.

The calculation also clarifies the interpretation of the singularly perturbed scalar BVP comparison in Section[3\.1](https://arxiv.org/html/2608.19658#S3.SS1), where theRECmodel improves mainly in the first three parameter bins but does not reduce every full\-test metric relative to theChebyshevmodel\. For the canonical exponential layer, the active polynomial degree scale is proportional toδ−1/2\\delta^\{\-1/2\}, indicating that a degree128128Chebyshev trunk remains an effective representation forδ≥10−4\\delta\\geq 10^\{\-4\}\. Thus, the main role of the AAA\-constructed inner rational subdictionary is not to replace high\-degree Chebyshev approximation in this ideal scalar case, but to supply an explicitly layer\-aligned component within a trunk dictionary having the same total number of dictionary elements when the target profiles combine smooth outer behavior with wall\-attached transport layers\.

## References

- \[1\]R\. Courant, D\. Hilbert, Methods of Mathematical Physics: Partial Differential Equations, Wiley, New York, 1989\.[https://doi\.org/10\.1002/9783527617234](https://doi.org/10.1002/9783527617234)\.
- \[2\]H\.\-G\. Roos, M\. Stynes, L\. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection\-Diffusion\-Reaction and Flow Problems, second ed\., Springer Series in Computational Mathematics, vol\. 24, Springer, Berlin, 2008\.[https://doi\.org/10\.1007/978\-3\-540\-34467\-4](https://doi.org/10.1007/978-3-540-34467-4)\.
- \[3\]A\. Arzani, K\.W\. Cassel, R\.M\. D’Souza, Theory\-guided physics\-informed neural networks for boundary layer problems with singular perturbation, J\. Comput\. Phys\. 473 \(2023\) 111768\.[https://doi\.org/10\.1016/j\.jcp\.2022\.111768](https://doi.org/10.1016/j.jcp.2022.111768)\.
- \[4\]L\. Zhang, G\. He, Multi\-scale\-matching neural networks for thin plate bending problem, Theor\. Appl\. Mech\. Lett\. 14 \(1\) \(2024\) 100494\.[https://doi\.org/10\.1016/j\.taml\.2024\.100494](https://doi.org/10.1016/j.taml.2024.100494)\.
- \[5\]F\. Cao, F\. Gao, X\. Guo, D\. Yuan, Physics\-informed neural networks with parameter asymptotic strategy for learning singularly perturbed convection\-dominated problem, Comput\. Math\. Appl\. 150 \(2023\) 229–242\.[https://doi\.org/10\.1016/j\.camwa\.2023\.09\.030](https://doi.org/10.1016/j.camwa.2023.09.030)\.
- \[6\]L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, G\.E\. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nat\. Mach\. Intell\. 3 \(3\) \(2021\) 218–229\.[https://doi\.org/10\.1038/s42256\-021\-00302\-5](https://doi.org/10.1038/s42256-021-00302-5)\.
- \[7\]S\. Wang, H\. Wang, P\. Perdikaris, Learning the solution operator of parametric partial differential equations with physics\-informed DeepONets, Sci\. Adv\. 7 \(40\) \(2021\) eabi8605\.[https://doi\.org/10\.1126/sciadv\.abi8605](https://doi.org/10.1126/sciadv.abi8605)\.
- \[8\]T\. Du, Z\. Huang, Y\. Li, Approximation and Generalization of DeepONets for Learning Operators Arising from a Class of Singularly Perturbed Problems, East Asian J\. Appl\. Math\. 14 \(4\) \(2024\) 841–873\.[https://doi\.org/10\.4208/eajam\.2023\-128\.051023](https://doi.org/10.4208/eajam.2023-128.051023)\.
- \[9\]T\. Sun, J\. Zu, PVD\-ONet: A Multi\-scale Neural Operator Method for Singularly Perturbed Boundary Layer Problems, arXiv:2507\.21437 \(2025\)\.[https://doi\.org/10\.48550/arXiv\.2507\.21437](https://doi.org/10.48550/arXiv.2507.21437)\.
- \[10\]J\. Lee, Y\. Hong, S\. Ko, J\.Y\. Lee, Data\-free Asymptotics\-Informed Operator Networks for Singularly Perturbed PDEs, arXiv:2512\.22006 \(2025\)\.[https://doi\.org/10\.48550/arXiv\.2512\.22006](https://doi.org/10.48550/arXiv.2512.22006)\.
- \[11\]C\. Mou, Y\. Zhang, X\. Zhu, Q\. Zhuang, PAS\-Net: Physics\-informed Adaptive Scale Deep Operator Network, arXiv:2511\.14925 \(2025\)\.[https://doi\.org/10\.48550/arXiv\.2511\.14925](https://doi.org/10.48550/arXiv.2511.14925)\.
- \[12\]L\. Lu, X\. Meng, S\. Cai, Z\. Mao, S\. Goswami, Z\. Zhang, G\.E\. Karniadakis, A comprehensive and fair comparison of two neural operators \(with practical extensions\) based on FAIR data, Comput\. Methods Appl\. Mech\. Engrg\. 393 \(2022\) 114778\.[https://doi\.org/10\.1016/j\.cma\.2022\.114778](https://doi.org/10.1016/j.cma.2022.114778)\.
- \[13\]J\. Choi, T\. Yun, N\. Kim, Y\. Hong, Spectral operator learning for parametric PDEs without data reliance, Comput\. Methods Appl\. Mech\. Engrg\. 420 \(2024\) 116678\.[https://doi\.org/10\.1016/j\.cma\.2023\.116678](https://doi.org/10.1016/j.cma.2023.116678)\.
- \[14\]Z\. Liu, H\. Wang, H\. Zhang, K\. Bao, X\. Qian, S\. Song, Render unto Numerics: Orthogonal Polynomial Neural Operator for PDEs with Nonperiodic Boundary Conditions, SIAM J\. Sci\. Comput\. 46 \(4\) \(2024\) C323–C348\.[https://doi\.org/10\.1137/23M1556320](https://doi.org/10.1137/23M1556320)\.
- \[15\]M\. Abid, O\. San, Spectral Embedding via Chebyshev Bases for Robust DeepONet Approximation, arXiv:2512\.09165 \(2025\)\.[https://doi\.org/10\.48550/arXiv\.2512\.09165](https://doi.org/10.48550/arXiv.2512.09165)\.
- \[16\]T\.M\. Squires, R\.J\. Messinger, S\.R\. Manalis, Making it stick: convection, reaction and diffusion in surface\-based biosensors, Nat\. Biotechnol\. 26 \(4\) \(2008\) 417–426\.[https://doi\.org/10\.1038/nbt1388](https://doi.org/10.1038/nbt1388)\.
- \[17\]T\. Gervais, K\.F\. Jensen, Mass transport and surface reactions in microfluidic systems, Chem\. Eng\. Sci\. 61 \(4\) \(2006\) 1102–1121\.[https://doi\.org/10\.1016/j\.ces\.2005\.06\.024](https://doi.org/10.1016/j.ces.2005.06.024)\.
- \[18\]S\. Chevalier, Semianalytical modeling of the mass transfer in microfluidic electrochemical chips, Phys\. Rev\. E 104 \(3\) \(2021\) 035110\.[https://doi\.org/10\.1103/PhysRevE\.104\.035110](https://doi.org/10.1103/PhysRevE.104.035110)\.
- \[19\]A\. Sescu, R\. Alaziz, M\.Z\. Afsar, Effect of wall transpiration and heat transfer on Görtler vortices in high speed flows, AIAA J\. 57 \(3\) \(2019\) 1159–1171\.[https://doi\.org/10\.2514/1\.J057330](https://doi.org/10.2514/1.J057330)\.
- \[20\]S\. Hillcoat, J\.\-P\. Hickey, Pressure–velocity coupling in transpiration cooling, Int\. J\. Heat Mass Transf\. 239 \(2025\) 126532\.[https://doi\.org/10\.1016/j\.ijheatmasstransfer\.2024\.126532](https://doi.org/10.1016/j.ijheatmasstransfer.2024.126532)\.
- \[21\]D\. Xu, J\. Wang, S\. Chen, Skin\-friction and heat\-transfer decompositions in hypersonic transitional and turbulent boundary layers, J\. Fluid Mech\. 941 \(2022\) A4\.[https://doi\.org/10\.1017/jfm\.2022\.269](https://doi.org/10.1017/jfm.2022.269)\.
- \[22\]D\. Passiatore, L\. Sciacovelli, P\. Cinnella, G\. Pascazio, Finite\-rate chemistry effects in turbulent hypersonic boundary layers: A direct numerical simulation study, Phys\. Rev\. Fluids 6 \(5\) \(2021\) 054604\.[https://doi\.org/10\.1103/PhysRevFluids\.6\.054604](https://doi.org/10.1103/PhysRevFluids.6.054604)\.
- \[23\]N\. Perakis, O\.J\. Haidn, M\. Ihme, Heat transfer augmentation by recombination reactions in turbulent reacting boundary layers at elevated pressures, Int\. J\. Heat Mass Transf\. 178 \(2021\) 121628\.[https://doi\.org/10\.1016/j\.ijheatmasstransfer\.2021\.121628](https://doi.org/10.1016/j.ijheatmasstransfer.2021.121628)\.
- \[24\]R\.K\. Shah, A\.L\. London, Laminar Flow Forced Convection in Ducts: A Source Book for Compact Heat Exchanger Analytical Data, Academic Press, New York, 1978\.[https://doi\.org/10\.1016/C2013\-0\-06152\-X](https://doi.org/10.1016/C2013-0-06152-X)\.
- \[25\]A\.S\. Haase, S\.J\. Chapman, P\.A\. Tsai, D\. Lohse, R\.G\.H\. Lammertink, The Graetz–Nusselt problem extended to continuum flows with finite slip, J\. Fluid Mech\. 764 \(2015\) R3\.[https://doi\.org/10\.1017/jfm\.2014\.733](https://doi.org/10.1017/jfm.2014.733)\.
- \[26\]A\.S\. Popel, J\.F\. Gross, Mass transfer in the entrance region of a circular tube, Int\. J\. Heat Mass Transf\. 21 \(8\) \(1978\) 1133–1141\.[https://doi\.org/10\.1016/0017\-9310\(78\)90112\-6](https://doi.org/10.1016/0017-9310(78)90112-6)\.
- \[27\]V\. Debarnot, J\. Fehrenbach, F\. de Gournay, L\. Martire, The case of Neumann, Robin, and periodic lateral conditions for the semi\-infinite generalized Graetz problem and applications, SIAM J\. Appl\. Math\. 78 \(4\) \(2018\) 2227–2251\.[https://doi\.org/10\.1137/17M1157507](https://doi.org/10.1137/17M1157507)\.
- \[28\]T\. Aquino, Equilibrium distributions under advection–diffusion in laminar channel flow with partially absorbing boundaries, J\. Fluid Mech\. 985 \(2024\) A16\.[https://doi\.org/10\.1017/jfm\.2024\.294](https://doi.org/10.1017/jfm.2024.294)\.
- \[29\]Y\. Wang, S\. Chen, X\. Wu, A rational spectral collocation method for solving a class of parameterized singular perturbation problems, J\. Comput\. Appl\. Math\. 233 \(10\) \(2010\) 2652–2660\.[https://doi\.org/10\.1016/j\.cam\.2009\.11\.011](https://doi.org/10.1016/j.cam.2009.11.011)\.
- \[30\]Y\. Nakatsukasa, O\. Sète, L\.N\. Trefethen, The AAA Algorithm for Rational Approximation, SIAM J\. Sci\. Comput\. 40 \(3\) \(2018\) A1494–A1522\.[https://doi\.org/10\.1137/16M1106122](https://doi.org/10.1137/16M1106122)\.
- \[31\]NIST Digital Library of Mathematical Functions, Section 10\.35: Generating Function and Associated Series for modified Bessel functions\.[https://dlmf\.nist\.gov/10\.35](https://dlmf.nist.gov/10.35)\.
- \[32\]K\.B\. Athreya, Modified Bessel function asymptotics via probability, Stat\. Probab\. Lett\. 5 \(5\) \(1987\) 325–327\.[https://doi\.org/10\.1016/0167\-7152\(87\)90004\-6](https://doi.org/10.1016/0167-7152(87)90004-6)\.