Thank you for your comments which have provided an opportunity to clarify key aspects of our work and the critical difference.

Part 1

Weakness 1. Novelty and Fundamental Differences from PINNs: SC-FNO uniquely adds time-integrated parameter sensitivities (∂u/∂p) as a supervisory signal (Lines 368-372). PINNs do not have access to this information, as usual PDEs do not contain ∂u/∂p, leaving it unconstrained. Procedurally, SC-FNO and PINNs differ entirely: SC-FNO utilizes a forward numerical model with differentiable solvers (or finite difference, or adjoint) to prepare data to supervise ∂u/∂p, whereas PINNs rely on collocation points for equation-based loss optimization to supervise (∂u/∂x, ∂u/∂t), but not ∂u/∂p.

The experimental benefits of ∂u/∂p are clear:

• SC-FNO consistently achieves R² > 0.94 for all parameters in PDE1, surpassing FNO and FNO-PINN which remain below 0.82.
• In parameter inversion tasks, SC-FNO reaches near-perfect recovery (R² = 0.998), significantly outperforming FNO (R² = 0.905).
• For complex systems like Navier-Stokes equations, SC-FNO enhances solution accuracy (R² from 0.883 to 0.994) and sensitivity estimation (R² from 0.246 to 0.997).

Given these profound differences in theory, procedure, and performance, the critique of lacking novelty is surprising. Is it fair to dismiss the differentiated outcomes in Figure 1? Does using gradient terms in the loss truly negate novelty?

Weakness 2. Data Efficiency and Computational Considerations: Our experiments confirm that leveraging sensitivity information via our SC-FNO model enhances both efficiency and robustness. Specifically:

• Computational Efficiency: Gradient information, a byproduct of our differentiable solver, increases computation time only moderately. Memory overhead remains low (722MB for SC-FNO vs. 764MB for FNO). Dataset preparation time was in the original Table D.9. For PDE2 (with 1300 training samples), training time per epoch increased only from 21 to 28 seconds. Also, see the new example below for training time and L2 error. We will add training time info to all Tables.

• Reduced Data Requirements: SC-FNO matches the accuracy of FNO with four times fewer data points (500 vs. 2000, original Figure 6). Under limited data conditions (500 samples), SC-FNO (R²>0.9) is significantly better than FNO’s (R²=0.8).

• Robust Generalization: SC-FNO sustains accuracy with parameter perturbations up to 40% beyond the training range and excels in parameter inversion tasks, boosting R² from 0.642 to 0.986 for PDE1 (original Figure 1).

In summary, SC-FNO's sensitivity supervision not only minimizes data needs but also ensures greater generalization, more than offsetting the slight increase in computational demands.

We also refined our experiment with PDE2 (Burger's equation), segmenting the spatial domain into 40 zones. Each zone had specific advection (α) and forcing amplitude (δ), alongside global parameters viscosity (γ) and frequency (ω), totaling 82 parameters (2*M+2). The results are:

These new results (to be included in revision) show that, with higher input dimensions, the default FNO had a much larger error for the state variable while FNO stayed highly accurate.

Weakness 3. “SC-FNO Should Be Harder to Train”:

The FNO's direct work in Fourier space actually helps with SC-FNO training, not hinders it. Our experiments show this: training remains stable with learnable loss coefficients and adds minimal overhead (only 10% computation time, 5.8% memory increase from 722MB to 764MB). Instead of making training harder, the sensitivity information actually helps guide the model to better solutions - that's why we see such improved performance across different equations. The complete process detailed in Algorithms 1-3 shows how this straightforward approach brings major benefits without significant complexity.

Dear reviewer,

We greatly appreciate your thorough review, which has helped improve the quality of our manuscript. We have carefully considered all your comments and have updated our manuscript accordingly to address your concerns.

We would be most grateful for your input if any additional revisions are needed during the remaining time of the discussion period.

We sincerely appreciate the thoughtful and detailed feedback provided in your review. Your comments have identified key areas where further clarification could significantly strengthen our manuscript.

Part 1

Question 1 (Discussion on Memory Overhead in AD and FD): Our discussion in Section 3.4 primarily focuses on demonstrating the feasibility of computing Jacobians using both traditional finite difference methods and automatic differentiation approaches. The experimental results show that both methods are practical, with AD being more efficient. As detailed in Table D.9, generating 2000 samples with gradients takes 252.54 seconds using AD versus 1907.32 seconds using finite differences, while achieving comparable accuracy (R² > 0.96). The memory overhead remains modest - training FNO on PDE1 used 722MB while SC-FNO required 764MB, representing only a 5.8% increase. We would be happy to provide a more detailed profiling of memory usage in the revised manuscript for other cases, but generally, the requirements are quite modest.

Dear reviewer,

We greatly appreciate your thorough review, which has helped improve the quality of our manuscript. We have carefully considered all your comments and have updated our manuscript accordingly to address your concerns.

We would be most grateful for your input if any additional revisions are needed during the remaining time of the discussion period.

Part 2

Question4 (Assessing SC-FNO's Performance Under Bifurcation-Rich Systems):

Excellent question and thanks for asking! While our original submission included cases with bifurcation characteristics like Navier-Stokes, we've conducted additional experiments with two well-known bifurcation-rich systems with varying input dimensions: First is the Allen-Cahn equation with five parameters, represented as ∂u/∂t = ϵ(x)∂²u/∂x² + α(x)u - β(x)u³, with initial condition u(x,0) = Atanh(kx). Parameters include diffusion coefficient ϵ, linear term α, cubic term β, and initial condition parameters A, k, with ranges: ϵ ∈ [0.01,0.1], α ∈ [0.01,1.0], β ∈ [0.01,1.0], A ∈ [0.1,0.9], k ∈ [5.0,10.0]. This equation exhibits phase transitions with sharp solution changes near critical parameter values. With 100 training samples, SC-FNO has relative L2: 0.015 vs FNO’s 0.021, only a noticeable benefit for u (however, keep watching as the second case gets more interesting), while also showing much improved Jacobian accuracy (relative L2: 0.049 vs 0.583).

The second case is the Ginzburg-Landau equation: ∂u/∂t = α(x)u + β(x)∂²u/∂x² + γ(x)u³ + δu⁵, with initial condition u(x,0) = Acos(2πx/L). Here, α(x), β(x), γ(x) are divided into M=40 zones, plus global parameters δ and A, resulting in 122 parameters (3M + 2).

We used the following ranges for the parameters:

• α(x) ∈ [0.05,0.5] per zone (linear term)
• β(x) ∈ [0.01,0.2] per zone (diffusion)
• γ(x) ∈ [-1.0,-0.5] per zone (cubic term) Plus global parameters:
• δ ∈ [0.001,0.01] (quintic term)
• A ∈ [0.01,0.05] (initial amplitude)

Training with both N=500 and N=100 samples shows SC-FNO maintains excellent state prediction even with small data. With N=100, it achieves nearly 4x better state prediction accuracy (relative L2: 0.007 vs FNO’s 0.025) and even 3 times smaller error than FNO with N=500 and more training time. It also shows improved Jacobian accuracy (relative L2: 0.162 vs 1.132), demonstrating robustness for high-dimensional parameters with limited data.

These highly nonlinear, bifurcation-rich cases demonstrate our method's ability to handle these challenges. This, to our knowledge, maybe the only known surrogate model scheme that can handle these many input parameters, and actually, we do not see the limit anywhere in sight! Future papers will explore much higher degrees of input..

We greatly appreciate the reviewers' thorough examination. Our apologies for posting the reply to you later than to others, as we were diligently running the requested bifurcation cases, which was a great suggestion.

Part 1

Questions 1 (Model Comparisons and Computational Details): Yes! In all FNO models (FNO, SC-FNO, FNO-PINN, SC-FNO-PINN), we used identical FNO architectures for each case, with differences only in loss configurations. We have updated Table C.5 to include the requested information and added a new Table C.6 showing computational times across all models and cases reported in Tables 1-3.

These additions will provide complete transparency regarding model complexity and computational requirements across all comparisons.

Question 2 (Potential advantages or limitations of including higher order terms): Thanks for the chance to clarify. Rather than incorporating higher-order derivatives in our loss function, we focus specifically on first-order parameter sensitivities (du/dp). This choice is motivated by our goal of improving parameter-dependent tasks, and our results demonstrate that these first-order sensitivities alone provide significant improvements in accuracy and robustness without the need for higher-order terms. This leads to significant improvements in accuracy and robustness particularly for out-of-sample predictions and parameter inference and reduced training data demand. In terms of disadvantages, we note the moderately larger effort in data preparation, the need for either a differentiable solver, finite difference, or adjoint, and one might be tempted to think of unstable gradients (although we saw this in none of the test cases). Paying an acceptance cost to substantially raise the performance ceiling should be worth it in many cases, especially when the input dimension gets higher, as we show below.

Question 3 (Extension to Complex Geometries): Indeed an interesting direction would be to handle complex geometries beyond rectangular domains. This will take more time (for many other methods in the literature, irregular geometry is tackled in the 2nd or 3rd papers). Since our parameter-solution Jacobian supervision is fundamentally geometry-independent, such integration should be natural and could indeed be highly valuable for complex engineering applications. For example, architectures like the Geometry-Informed Neural Operator (Li et al., 2024) seems plausible; There is also no fundamental barriers to using traditional irregular geometry methods to solve for stencil weights as a nondifferentiable operations. We will write this as a future direction.

Dear reviewer,

We greatly appreciate your thorough review, which has helped improve the quality of our manuscript. We have carefully considered all your comments and have updated our manuscript accordingly to address your concerns.

We would be most grateful for your input if any additional revisions are needed during the remaining time of the discussion period.

Thank you for your detailed and thoughtful response, particularly regarding the mode, number of learnable parameters, training time, and other aspects of FNO training. Your clarifications have been highly informative and are sincerely appreciated.

I fully understand that addressing my earlier questions about incorporating a loss function involving higher derivatives (Question 2) and conducting simulations on complex domains (Question 3) would require substantial time and effort. It is entirely reasonable that these issues are not addressed in the current version of the manuscript.

At this stage, I would like to kindly request clarification on the following two points:

1. Could you provide further details on how $\alpha(x)$ and $\beta(x)$ were generated for the Allen-Cahn equation or Ginzburg-Landau equations? I was unable to find this information in the appendix. While I understand that further revisions may no longer be feasible at this stage, any clarification on how the parameter set was constructed would be greatly appreciated.

2. Could you elaborate on the computation of $\partial \mathbf{u} / \partial\mathbf{p}$ in equation (6) for high-dimensional parameter cases? Given that $\mathbf{p}$ represents a function in the context of the Allen Cahn equation and Ginzburg-Landau equation, computing this derivative appears to be nontrivial. Any additional explanation regarding this computation would be extremely helpful.

I would be deeply grateful if the authors could provide clarifications, should it still be possible at this stage. Thank you again for your time and effort in addressing these inquiries.

We greatly appreciate the reviewers' detailed examination of our work. Despite the constraints of the word limit, we have addressed all feedback thoroughly. Additionally, we have responded separately to the points you identified as weaknesses in our work.

Part 1

Question 1. Runtime & Memory Requirements: We apologize for this oversight. The overhead is small compared to the benefits. Our experiments with PDE1 show that using 500 training points with gradients achieves similar accuracy to 2000 points without gradients while the gradient computation added an acceptable amount of time. As shown in Table D.9, in the original manuscript Appendix, generating 2000 samples with gradients takes just 252.54 seconds using our AD solver. Training SC-FNO in this case, on the other hand, required 20%-40 more time compared to the original FNO depending on the case. In terms of memory consumption, the difference is also modest—training FNO on PDE1 used 722 MB, while SC-FNO required 764 MB. This modest computational and memory overhead is justified by the significant improvements in model performance and generalization capabilities, demonstrating how sensitivity regularisation and time-step-free methods like FNO complement each other exceptionally well. In the next response, you will further see a new experiment that contains 84 parameters where FNO’s accuracy even for u dropped noticeably while SC-FNO maintained accuracy. For this case (to be added to the revised version), the added training time was ~35% extra, but the errors were reduced by 4 times.

Question 2. Handling Functional Inputs (du/du0): Great question. Yes! We ran a new experiment with PDE2 (Burger's equation) using zoned parameters - dividing the spatial domain into M=40 segments with different advection α and forcing amplitude δ in each zone, along with global parameters viscosity γ and forcing frequency ω (total 2M+2=82 parameters).

The results demonstrate SC-FNO's clear advantages over FNO (Table below):

• With N=100 samples, SC-FNO got less than 1/3 of the relative L2 error as N=500 for FNO (L2=2.8e-2) in less computational time.
• With 100 training samples, SC-FNO achieves ~4x less error in state prediction (relative L2: 0.00729 vs 0.02822) and 11x better Jacobian accuracy (relative L2: 0.17698 vs 1.96270)
• Even with only 50 training samples, SC-FNO maintains high accuracy while FNO degrades significantly.
• This dramatic improvement comes with only a modest runtime increase (11.23s vs 8.09s per epoch for N=100)

These results clearly demonstrate SC-FNO's superior capability in handling high-dimensional functional parameters while maintaining accuracy and stability, even with limited training data. In fact, we see clearly we can allow much larger amounts of inputs, but we believe that would require careful benchmarking in a future paper.

Thank you for your thoughtful review. Below, you will find our detailed responses to the points identified as weaknesses.

Writing Clarity and Parameter Definitions: The comprehensive parameter definitions, ranges, and equations that the reviewer notes are missing are detailed in Appendix B (Table B.4) and Section 2. For example, Table B.4 explicitly lists parameter ranges for all cases: ODE1 (α,β ∈ [1,3], γ ∈ [0,1]), PDE1 (c ∈ [0,0.25], α,β,γ,ω ∈ [0,0.25]), etc. Due to page limitations, we focused the main text on key findings and methodology while preserving technical details in appendices.

Response to Parameter Implementation in FNO:

As shown in Figure A.7, parameters (P) are incorporated through the lifting layer along with spatial coordinates (X) and initial conditions (a(x)). The architecture processes parameters as additional input channels that are combined with other inputs through neural network layers. Full parameter definitions and ranges are provided in Table B.4. Computational considerations regarding this implementation are addressed in our main rebuttal under memory and runtime discussions.

Response to Accuracy and Performance Analysis (Methodology):

We began our study acknowledging FNO's strong accuracy with large training samples, deliberately starting from conditions where original FNO excels. This approach allowed us to demonstrate SC-FNO's significant advantages in both state values (u) and sensitivities (∂u/∂p) under more challenging, real-world conditions. For parameter perturbations up to 40% beyond training ranges, SC-FNO maintains R² > 0.91 while FNO drops to 0.53 (Table 1). With limited data (500 samples), SC-FNO achieves R² > 0.9 while FNO falls to 0.8 (Figure 6). These results underscore our method's practical value, showing improved robustness and efficiency precisely where standard FNO struggles most.

Dear reviewer,

Thank you for your reply and insightful feedback.

If we understand your question correctly, of course, p is used as input to both FNO and SC-FNO and it's a fair comparison. Their input data are identical. We wouldn't make a mistake like not providing p to the default FNO. FNO has access to the same input information but it does not seem to fully understand how to use it like SC-FNO does.

The SC-FNO architecture processes parameters (P) alongside spatial coordinates and initial conditions through the lifting layer (fc0) (possibly as function inputs). This layer reshapes and repeats parameters to match the problem's spatial-temporal dimensions, then concatenates them with other inputs before neural network processing.

Data preparation time for SC-FNO was previously in Table D.9. We now amend it with the time for FNO without gradients as shown in the table below for PDE2 (Burger's equation) and PDE3 (Navier-Stokes equation), measured per 10 samples using 1 GPU Tesla P100-PCIE-12GB.

This represents a one-time cost per equation, and since input impacts are resolved upfront, no repeated preparation is needed for different parameters. Furthermore, when considering the supervision of high-dimensional inputs, where default FNO requires far more than 5x more data (and thus more preparation time), we argue this additional computational cost is well justified. Based on simple logic, the data preparation time of the default FNO will grow exponentially more to achieve similar accuracy as the input dimension grows.

We hope by now we have presented something truly unique and powerful here. If there is anything else we can do to convince you that is doable within the remaining time and reasonable scope, please let us know.

Dear reviewer,

Thank you for your reply and insightful feedback.

If we understand your question correctly, of course, p is used as input to both FNO and SC-FNO and it's a fair comparison. Their input data are identical. We wouldn't make a mistake like not providing p to the default FNO. FNO has access to the same input information but it does not seem to fully understand how to use it like SC-FNO does.

The SC-FNO architecture processes parameters (P) alongside spatial coordinates and initial conditions through the lifting layer (fc0) (possibly as function inputs). This layer reshapes and repeats parameters to match the problem's spatial-temporal dimensions, then concatenates them with other inputs before neural network processing.

Data preparation time for SC-FNO was previously in Table D.9. We now amend it with the time for FNO without gradients as shown in the table below for PDE2 (Burger's equation) and PDE3 (Navier-Stokes equation), measured per 10 samples using 1 GPU Tesla P100-PCIE-12GB.

This represents a one-time cost per equation, and since input impacts are resolved upfront, no repeated preparation is needed for different parameters. Furthermore, when considering the supervision of high-dimensional inputs, where default FNO requires far more than 5x more data (and thus more preparation time), we argue this additional computational cost is well justified. Based on simple logic, the data preparation time of the default FNO will grow exponentially more to achieve similar accuracy as the input dimension grows.

We hope by now we have presented something truly unique and powerful here. If there is anything else we can do to convince you that is doable within the remaining time and reasonable scope, please let us know.

Dear reviewer,

We greatly appreciate your thorough review, which has helped improve the quality of our manuscript. We have carefully considered all your comments and have updated our manuscript accordingly to address your concerns.

We would be most grateful for your input if any additional revisions are needed during the remaining time of the discussion period.

Many thanks for your thoughtful feedback and for raising your score. Your concerns about scalability and flexibility highlight important considerations for the method's practical applications.

Please trust us that SC-FNO can scale up to handle high-dimensional parameter spaces - we have preliminarily tested it with 2500 parameters with promising results. However, properly validating and presenting these findings requires additional time and rigorous analysis to meet publication standards. Given the current paper's focus on establishing the fundamental benefits of sensitivity supervision for parameter inversion, whose importance is showcased even with as few as one or a handful of parameters, we felt including preliminary high-dimensional results might dilute this core message. Meanwhile, the 82-parameter case has already been included in the revised manuscript.

Furthermore, we should highlight two key features of our framework:

• Neural operators have resolution-independent capabilities - FNO can be trained on lower resolution grids (e.g., 32×32) yet make predictions on much finer resolutions (128×128 or higher). This means parameter dimensionality can be managed efficiently without sacrificing accuracy at inference time (Please, see Figure 1, Page 2 at https://arxiv.org/pdf/2010.08895 ).
• Our training procedure employs strategic sampling - instead of computing gradients at all points, we randomly select a subset of spatial-temporal points in each epoch (n < N spatial points × t < T time points) where n < N and t < T. The neural operator's predicted gradients at these points are computed using AD and compared with the pre-computed true gradients. This sampling varies between epochs to eventually cover the full solution space. We explained this implementation in our revised manuscript.

We are committed to demonstrating these strong capabilities comprehensively in our next paper. This is a solemn promise. Given the time constraints and scope, we cannot include such results in the current work. We appreciate your understanding of the need to progress methodically in establishing this new approach. We also want to thank you again for your constructive comments that helped us shape the manuscript.