Papers

We present a general framework for the high-precision numerical evaluation of multivariate hypergeometric functions defined as solutions of holonomic systems of partial differential equations. Our approach adapts and extends methods originally developed for multi-loop Feynman integrals to the setting of hypergeometric functions of many variables. In particular, we construct Pfaffian systems for arbitrary multivariate hypergeometric functions by applying the Laporta reduction algorithm to suitable systems of differential relations. Next, we construct a numerical scheme based on the Frobenius method, which allows us to compute local power-series solutions with controlled precision and to transport them along prescribed paths in the space of variables. A central part of the paper is devoted to a systematic analysis of multivaluedness and branch structure: we show how the Frobenius method can be used to access different Riemann sheets in a controlled way and to track changes of the solution under analytic continuation around singular loci.

GenEO (`Generalised Eigenvalue problems on the Overlap') is a method for constructing coarse spaces used in the preconditioning of iterative solvers for discrete PDEs. This method combines a (small) number of modes of local PDE eigenproblems to obtain a global coarse space. A coarse solve is then combined with local solves of the global PDE to obtain the preconditioner. A substantial theory for GenEO has been developed for the case when the local elgenproblems are positive semi-definite. This has been applied mostly to positive definite global PDEs, but also recently extended to the case of convection–diffusion–reaction problems, which may be neither self-adjoint, nor positive definite. However, when the global problem is highly indefinite, coarse spaces built from positive semi-definite local eigenproblems fail to be robust in practice. In this paper we consider highly indefinite global PDE problems, characterised by a large parameter $k$ (allowing also highly variable coefficients), and we develop a new spectral coarse space built from solving eigenvalue problems based on local copies of the global problem. We put no constraint on the diameters of the local domains, thus allowing the local eigenvalue problems to be indefinite. The new method (which we call $H_k$-GenEO) is seen to be much more robust as $k$ increases than methods based on positive semi-definite eigenproblems. We provide sufficient conditions for robustness of the preconditioned GMRES iterative method, in terms of the tolerance of the local eigenproblems and the size of the subdomains for the local PDE solves. In practice the method is observed to be robust with respect to $k$ under even weaker conditions on the local eigenproblem tolerance. The experiments also suggest the method can be resilient to high variation in PDE coefficients.

Microwave linear analog computers (MiLACs) have recently emerged to enable high-performance and efficient beamforming in the analog domain. In this paper, we introduce a dual-functionality framework for MiLAC-aided transceivers. Beyond analog-domain precoding/combining (active beamforming), a MiLAC and its antenna array can simultaneously act as a reconfigurable intelligent surface (RIS) (passive beamforming). This allows the MiLAC to execute beamforming for transmission/reception while reflecting external incident signals. We provide an optimal reconfiguration strategy for this dual-functional MiLAC, and characterize the fundamental limits on the trade-off between active and passive rate, namely the capacity region bounds and the sum-rate capacity.

Dynamical systems reconstruction (DSR) aims to learn surrogate models that capture the dynamics underlying time-series data. Reliable scientific inference with these surrogates requires uncertainty estimates consistent with the learned dynamics. We expose a dynamic-probabilistic consistency (DPC) gap: optimizing finite-horizon probabilistic objectives can degrade the learned dynamics or decouple predictive uncertainty from the local tangent dynamics it should reflect. We isolate three mechanisms behind this gap: core collapse, noise masking, and blind uncertainty. Specifically, we show that open-loop rollout objectives can penalize covariance growth driven by the dynamics, encouraging optimization shortcuts that weaken local expansion. To mitigate this gap, we propose KAFFEE (Kalman-Aware Framework For Ergodic Emulation), a differentiable extended Kalman filter-based training framework that jointly learns dynamics and noise covariances by coupling filtered innovation likelihoods with covariance propagation through learned local Jacobians, while differentiating through the Kalman updates. On stochastic hyperchaotic Lorenz-96, KAFFEE mitigates the identified failure modes and improves both reconstruction of dynamical invariants and predictive scores relative to open-loop objectives. The DPC gap also appears when probabilistically adapting a DSR foundation model across 13 chaotic systems, where KAFFEE improves predictive likelihood over open-loop adaptation while preserving the pretrained model's long-term dynamics and enabling in-context Bayesian filtering.

We introduce weak and strong poset Ramsey-Turán numbers for $t$-chains in host poset families, focusing on the Boolean lattice family $\mathcal{B}=\{B_n:n\ge 1\}$. For any poset $P$, we show $\operatorname{RT}(\mathcal{B};n,P,l,t)\le \operatorname{RT}^{\sharp}(\mathcal{B};n,P,l,t)$, with equality when $P$ is a chain. In particular, for $t=1$, $\operatorname{RT}(\mathcal{B};n,C_k,l)=\operatorname{RT}^{\sharp}(\mathcal{B};n,C_k,l)=(k-1)(l-1)$. We also give universal upper bounds for both versions. For fixed $k,l,t$ with $\min\{l-1,k-1\}\ge 1$, we prove $\operatorname{RT}^{\sharp}(\mathcal{B};n,A_k,l,t)=\Theta(n^t)$. More generally, for every non-chain poset $P$, the strong number is $\Theta(n^t)$ for fixed $l,t$. Finally, if $h(P)=r>t$ and $l(n)=\lfloor M_n^\beta\rfloor$ with $0<\beta\le \alpha<1$, then both weak and strong versions admit lower bounds of order $\Omega\!\left(2^{\beta n}n^{-\beta/2}\right)$.

Orbital transfers in multi-body systems are often studied as isolated trajectory design problems, making it difficult to identify the larger transport structure connecting families of periodic orbits, including which families act as hubs, gateways, relays, or persistently difficult-to-access regions. This work introduces a reachable-set-based framework for constructing orbital networks in the circular restricted three-body problem. Finite-$\Delta V$ and finite-time-of-flight reachable-set overlaps are used to infer accessibility relationships between representative periodic orbit families on a common Jacobi energy manifold and to assemble these relationships into a weighted orbital network. Applied to the Earth-Moon system, the resulting network reveals distinct accessibility regimes in which direct reachability, graph connectedness, and feasible multileg closure emerge separately. The analysis identifies multi-orbiter cycler orbits as the dominant hub, gateway, and relay families, with the (3,2)-cycler dominating across much of the sampled budget plane and the short-period (1,1)-cycler dominating in the low-time-of-flight regime, while the stable 2:1 resonant orbit remains persistently difficult to access. Although the maximum-budget network is nearly complete in a binary sense, its weighted accessibility remains strongly non-uniform. Selected proxy-supported connections are refined into concrete trajectories through differential correction, with corrected transfer costs remaining below the proxy estimates in all tested cases. Together, the results demonstrate how reachable-set overlap geometry can expose large-scale transport structure in nonlinear gravitational systems without requiring exhaustive pairwise trajectory optimization.

It has been conjectured that generic smooth solutions of the two-dimensional Euler equation exhibit linear growth of vorticity gradients. We prove an elementary arbitrary-background perturbation principle in the odd symmetric setting. More precisely, for any compactly supported nonnegative function in the half-plane, one can find an arbitrarily small smooth nonnegative perturbation whose associated solution undergoes linear-in-time filamentation in the quadrant. The main ingredients are the lower bound of the center of mass given by Iftimie-Sideris-Gamblin, and the velocity estimate for the sparse part to capture those slowly moving particles.

A hyperbolic structure on a group $G$ is a (not necessarily properly discontinuous) cobounded action of $G$ on a Gromov hyperbolic space, considered up to coarsely $G$-equivariant quasi-isometry. We show that for groups $G$ acting geometrically and positively on a horocyclic product $X\bowtie Y$, all hyperbolic structures on $G$ come from the two actions on the factors. The ingredients of the proof include Malcev rigidity and a new ("vertical") boundary of horocyclic products. We also give the first example of a group acting geometrically on a horocyclic product of two mixed millefeuille spaces.

The polar-angle marginal of a centred bivariate Gaussian distribution, obtained after integrating out the radial coordinate, gives the two-dimensional angular central Gaussian (ACG) distribution of Tyler. While its trigonometric and vector-valued moments have been studied in detail, to our knowledge there are no explicit closed-form expressions for the linear moments $\mathbf{E}[\theta]$ and $\mathbf{E}[\theta^{2}]$ on the natural domain $\theta\in\left]-\pi/2,\pi/2\right[$. Here linear refers to the ordinary moments $\int\theta^{k}f(\theta)\,d\theta$ of the angle regarded as a real-valued variable, in contrast to the circular (trigonometric) moments $\mathbf{E}[e^{ik\theta}]$ customary in directional statistics. We provide such expressions: the mean is a simple arctangent of the parameters, while the second moment is given by the real part of a dilogarithm. The derivation, based on a contour integration around the branch cut of $\arctan z$, is elementary. These quantities naturally arise in physics, where $\theta$ is interpreted as a real-valued phase rather than a circular variable.

Maritime weather routing must optimize route geometry under dynamic wind-wave conditions, obstacle constraints, and fixed-arrival requirements. We present Bézier Evolve and Refine Strategy (BERS), a two-stage framework that combines global evolutionary search (CMA-ES) with local variational refinement (FMS). Routes are parametrized as Bézier curves and evaluated with dense along-path sampling, enabling smooth trajectories while preserving practical feasibility constraints and accounting for mid-segment effects. We evaluate BERS on synthetic benchmarks designed to stress seven operational criteria: continuity, obstacle avoidance, dynamic adaptation, flexible objective design, constant-load feasibility, just-in-time arrival, and local optimality. Across these tests, BERS matches or improves published baselines while maintaining robust convergence under challenging flow fields and land geometries. We then validate the method on real ocean data using hourly ERA5 forcing over 366 daily departures in 2024 for two trans-oceanic corridors (Atlantic and Pacific), with a physics-based model of an 88m cargo vessel with optional rigid wingsails. In real-ocean experiments, route optimization alone reduces mean propulsive energy by 23–59% versus great-circle baselines of the same propulsion mode. Combined with wind-assisted propulsion, total savings reach up to 75%. These results show that BERS provides a practical and scalable foundation for just-in-time, energy-efficient weather routing in maritime decarbonization workflows.

Satellite communications are envisioned as a key enabler for ubiquitous coverage in future 6G networks, yet the broadcast nature renders them vulnerable to eavesdropping, especially given the long-distance transmissions and associated high uncertainties. In this paper, we propose the physical layer security enhancement for multi-beam satellite communications with the assistance of an aerial reconfigurable intelligent surface (ARIS). Considering the high dynamics and uncertainties of channels, we characterize the channel distribution with moment-based ambiguity sets. Accordingly, a distributionally robust secrecy rate optimization is formulated through joint design of transmit and reflection beamforming. We then introduce a conditional value-at-risk-based reformulation to convert the probabilistic constraints into deterministic forms. An alternating optimization framework is subsequently employed to iteratively update the transmit and reflective beamforming vectors until convergence. Simulation results demonstrate that the proposed distributionally robust scheme significantly enhances secrecy performance, and maintains reliable performance across various channel error distributions.

This paper was triggered by recent results on the maximal `average distance between consecutive points' of uniformly distributed sequences (u.f.d.s.). Here we address a generalized version of this question, consider pairwise maximal/minimal/total distances in blocks/segments of $ d \geq 2 $ consecutive points of u.f.d.s., and derive sharp upper bounds for all three aggregations. Our main idea of proof consists in, firstly, adding degrees of freedom, secondly, translating the resulting problem to a solvable optimization problem over the compact family of $ d $-stochastic measures, and, thirdly, showing that the obtained bounds are also sharp bounds for the original problem.

In the finite-dimensional setting, it is known that STFT phase retrieval is always possible when the window's ambiguity function does not vanish. However, it is not known how many zeros are allowed in the ambiguity function for the window still to allow phase retrieval. In order to tackle this problem, we first consider a two-window approach where the second window equals the Fourier transform of the first window. This allows us to apply the uncertainty principle in order to obtain sufficient conditions for phase retrieval. Using the relation between STFT phase retrieval and ambiguity sampling, we can prove sufficient conditions for the single-window phase retrieval problem, showing that only approximately eight ninths of the entries of the window's ambiguity function (and only three quarters in prime dimensions) are required to be nonzero.

We propose the sketched landing method, a randomized variant of the landing method for optimization under orthogonality constraints. Each landing step consists of the sum of a normal component, which reduces infeasibility, and a tangent component, which decreases the objective function. Our main contribution is the introduction of low-dimensional random sketch matrices to reduce the computational cost of these directions. We consider both dense (Gaussian) and sparse (subsampling) sketch matrices, and show how they reduce the per-iteration cost while preserving convergence guarantees in expectation.

We study the Hausdorff dimension spectrum of asymptotic approximation rates of complex Diophantine approximation and that of the asymptotic average excursion time of cusp excursions on the Bianchi orbifold $\mathbb{H}^3/\operatorname{PSL}(2,\mathbb {Z}[i])$ via a unified approach using the Hurwitz map. In particular, we construct a conformal graph directed system (CGDS) for the Hurwitz map and show that the Lyapunov exponent of the Hurwitz CGDS simultaneously captures the asymptotic approximation rate and the the asymptotic average excursion time. Applying the multifractal analysis of Lyapunov exponents for this system, we obtain a formula and real-analyticity for the Hausdorff dimension spectrum functions.

We propose a new family of polynomial approximation kernels for approximating nonnegative polynomials on the hypercube $[-1,1]^n$. Our Kernels produce polynomial sums-of-squares of degree $r$, achieving an $O(\log^3 r/r^2)$ error in the $\ell_1$-norm of the coefficients. This improves on the known error bound $O(1/r)$ from the literature. As a corollary, we obtain an improved convergence rate for the Lasserre hierarchy for polynomial optimization on the hypercube, again improving a known rate by Baldi and Slot from $O(1/r)$ to $O(\log^3 r/r^2)$.

We study the space of conformally covariant $n$-point functions of spinning operators using methods from invariant theory, commutative algebra, and combinatorics. We show that the rational part of any such function can be expressed in terms of the basic building blocks introduced by Costa, Penedones, Poland, Rychkov, thereby providing a rigorous proof of a result that is widely used in the physics literature. We reformulate the problem of enumeration of $n$-point structures in terms of counting lattice points in fractional matching polytopes, and compute these counts using vector partition functions, Hilbert functions, and Kostka numbers. We show that all algebraic relations between the building blocks follow from Gram constraints and compute the number of algebraically independent building blocks. For three-point functions, we derive closed counting formulas for arbitrary integer spins, both with and without Bose symmetry, and discuss a necessary and sufficient condition for the partial conservation operator to lift to a differential operator written in terms of the building blocks. We provide code that generates a basis of three-point structures satisfying these constraints for given values of spins and scaling dimensions.

We show that every pair of 2-torsion elements in the Tate-Shafarevich group of an elliptic curve are visible in the same abelian surface. This was previously only known for a single 2-torsion element. Our result explains some of the original observations on visibility in the paper of Cremona and Mazur.

Architecture diagrams are ubiquitous in deep learning, but they are usually only representational: the tensor-program identities they suggest are still proved by prose and tensor-axis manipulation. We introduce a formal graphical calculus for the structural fragment of tensor programming underlying einops, making such diagrams proof-enabling. Our calculus represents tensor axes as nested graded tubes around a base type. The tube boundary recovers the undirected tensor-network view of axes, while the directed interior retains the operational reading of computation graphs. The key rewrite is grade-naturality: sliding spectacles over tubes. Standard equivariance proofs become short diagrammatic derivations. We additionally demonstrate how our rewrite system may be applied to convert attention masks into pre-processing operations, recovering efficient implementations of sparse attention blocks.

We confirm a compactness conjecture of M. Li. If a complete Riemannian manifold has nonnegative Ricci curvature and uniformly convex boundary in the sense that the second fundamental form satisfies $h\ge1$. Then we prove it is compact, and consequently has finite fundamental group. The proof uses monotone quantities constructed via positive proper harmonic functions with Neumann condition.

We establish intrinsic area–length–topology inequalities for compact free boundary constant mean curvature (CMC) surfaces in three-manifolds with Ricci curvature bounded from below. Our main result is obtained from a conformal upper bound for a constrained first Robin eigenvalue of the Jacobi operator, derived via a balancing argument. This yields a quantitative inequality that does not require stability and captures both interior and boundary contributions. As an application, we obtain explicit topological restrictions for stable free boundary CMC surfaces under a natural curvature pinching condition. In particular, in weakly convex domains, stability forces low topological complexity, with genus at most three and a small number of boundary components. These results show that effective topological control persists even in negatively curved settings, where classical rigidity phenomena are no longer available.

We investigate uniform weak convergence rates for probabilistic numerical methods applied to backward time-fractional diffusion equations whose dynamics are driven by diffusions with possibly unbounded coefficients, such as the Geometric Brownian Motion. The fractional structure is represented through a random time-change by the inverse of a stable subordinator. To approximate the underlying fractional dynamics, we combine discrete Markov chain schemes for the diffusion component with heavy-tailed random walk approximations of the time change. Our analysis builds on Feller semigroup techniques and a high-order sensitivity framework for diffusion semigroups based on the Kunita stochastic flows and tensor fields. We derive uniform bounds for all orders of sensitivities, establish a quasi-contraction property for the associated semigroup, and transfer these estimates to the fractional setting via the convolution representation with the inverse subordinator. As a result, under killing conditions which dominate at least the base-space semigroup growth, we obtain weak convergence rates for the combined continuous-time-random-walk scheme to the time-fractional diffusion, with a logarithmic regime before the discount dominates the stronger smooth-space growth.

A nonparametric variant of the Kiefer-Weiss problem is proposed and solved. The objective is to minimize a weighted sum of the error probabilities of a binary sequential test subject to a constraint on its maximum expected sample size. This maximum is taken over all possible probability distributions on the given sequence space. First, it is shown that the nonparametric Kiefer-Weiss problem can be reduced to an optimal stopping problem. Then, the optimal stopping policy is derived under the assumption that at most k uses of randomization are permitted during any run of the test. The solution to the original problem is then obtained by letting k go to infinity. The optimal cost function is shown to be the solution of a nonlinear Bellman equation. The corresponding optimal stopping policy is shown to be based on a two-dimensional test statistic, with one component tracking the likelihood ratio and the other one tracking the expected remaining sample size. Critically, the stopping policy uses randomization to increase the remaining expected sample size for some runs, while stopping early for others. The optimal randomization rule is shown to be determined by a function that maps the likelihood ratio to an integer-valued sample size. Two approximations of this function are proposed that can be evaluated easily in practice. The results are illustrated with two numerical examples of nonparametric Kiefer-Weiss tests, one for a shift in the success probability of a Bernoulli distribution, and one for a shift in the mean of a normal distribution.

This paper studies the local well-posedness and the behavior at potential finite blow-up times of mild solutions to the three-dimensional fractional Boussinesq equations in Lei-Lin and Lei-Lin-Gevrey spaces $\mathcal{X}_{a,\sigma}^s(\mathbb{R}^3)$. By combining fixed point arguments with Fourier estimates adapted to these spaces, we obtain local existence and uniqueness results for data in $\mathcal{X}_{a,\sigma}^s(\mathbb{R}^3)$, including the usual Lei-Lin case $a=0$. We also establish several blow-up criteria for maximal mild solutions and derive lower bounds for the growth of the corresponding norms as the maximal time is approached. In the strict Lei-Lin-Gevrey regime, and when the dissipative exponents coincide, these estimates yield an exponential type blow-up criterion.

The underdamped harmonic oscillator is analyzed through the complex mapping $\zeta = e^{-i\varphi}we^{-w}$ with $w = \beta t + i\Omega t$, which transforms the dynamics into a logarithmic spiral. Within this framework, the displacement extrema correspond to crossings of the imaginary axis by $\zeta(t)$, yielding the explicit times $t_n = (\theta - \varphi - \pi/2 + n\pi)/\Omega$, where $\theta = \arctan(\Omega/\beta)$. The Lambert $W$ function provides closed-form solutions $t = -\beta^{-1}W_k(-\beta A/\omega_0)$ for the times at which the spiral radius attains a given threshold $A$, covering both the rising and decaying branches. The quality factor $Q = \omega_0/(2\beta) = \tfrac{1}{2}\sec\theta$ is directly encoded in the ray angle $\theta$ of the $(u,v)$-plane. Key geometric invariants are derived: the winding number $N_\varepsilon \approx (Q/\pi)\ln(2Q/\varepsilon)$ for large $Q$, the enclosed area $A = \omega_0^2\Omega/(8\beta^3) \approx Q^3$ in the lightly damped limit, and the energy decay $E(t) = E_0 e^{-\omega_0 t/Q}$. Three methods for determining $Q$ from experimental data are compared: logarithmic decrement, ray-angle measurement, and spiral turn counting. The turn-counting method proves particularly robust for high-$Q$ systems, where successive amplitude peaks differ by tiny fractions. The framework unifies classical damped oscillations with complex analysis and special functions.

Free-fermion solvability in quantum spin systems is increasingly understood to be governed by a graph Clifford algebra defined from the frustration graph of the Hamiltonian. When the frustration graph belongs to certain classes, such as the even-hole-free and claw-free (ECF) class, the Hamiltonian is solvable by hidden free fermions: it admits a free-fermion solution although it does not reduce to a Majorana bilinear under the Jordan-Wigner transformation. However, unlike in the Jordan-Wigner case, where each mode is a linear combination of single Majorana fermions, the explicit operator structure of the hidden free-fermion modes—and that of the local conserved charges—has remained obscure. In this work, we derive a path-product expansion that expresses each free-fermion mode as a linear combination of products along induced paths in the extended frustration graph. The expansion follows from the Krylov generating function and yields the modes directly, without using the transfer matrix or nonlocal conserved charges; the resulting mode decomposition also computes infinite-temperature dynamical correlation functions. We further obtain explicit expressions for local conserved charges as linear combinations of path products along induced paths; these charges apply beyond the free-fermion (ECF) class to more general claw-free frustration graphs. We also identify a family of generalized conserved charges containing both the known nonlocal conserved charges and these local charges. For the homogeneous periodic Fendley model, the local conserved charges exhibit the same Catalan-tree pattern as those of the spin-$1/2$ XXX chain.

We investigate the relationship between the global convergence of Cesàro averages and the pointwise statistical behavior of dynamical systems. First, we prove that if the Cesàro averages accumulate on at least two different measures (a property we call the minimal oscillation property) then the system is non-statistical. Second, we show that a system possesses a natural measure in the strong sense if and only if it is uniquely ergodic. As a consequence, every minimal homeomorphism on the circle possesses a natural measure in the strong sense which is physical and whose basin is the entire circle.

We present a regression-adjustment framework designed for the estimation of longitudinal treatment effects in randomized experiments under static regimes. While regression-adjustment methods are useful for variance reduction in randomized experiments by using pre-treatment covariates, they usually focus only on average effects, from which we cannot obtain valuable insights into when the effects appear and how long they continue. To address this issue, we consider intermediate outcomes and evolving post-treatment covariates over time, and we represent such dynamic trajectories using transition kernels. Furthermore, we establish the asymptotic normality and the semiparametric efficiency bound for our estimator, enabling more powerful statistical inference. Simulation studies and empirical analysis using A/B test data from a streaming platform in Japan show the practical advantages of our method.

The spectral density estimation problem asks for an algorithm that, given an $n\times n$ matrix $A$, outputs a probability measure that is a good approximation to the uniform distribution on the eigenvalues of $A$, called the spectral density of $A$. This paper considers the setting where $A$ is a large normal matrix that is accessible only through matrix-vector product queries. We provide an algorithm that makes just $m$ matrix-vector queries to $A$ and returns, with high probability, a measure within earth mover's distance $O(1/m+\log m/{\sqrt n})$ of the true spectral density of $A$. We provide a complementary lower bound that any algorithm producing an $\varepsilon$-approximation to the true spectral density for large matrices must make $\Omega(1/\varepsilon)$ matrix-vector queries. The lower bound holds even for the more restricted case of real symmetric input matrices. In combination with our upper bound, it shows that spectral density estimation is essentially no harder for complex normal matrices than for real symmetric matrices.

Image Scanning Microscopy (ISM) is a fluorescence imaging technique that combines detector-array acquisition and computational reconstruction to achieve the theoretical resolution of an ideal confocal microscope, i.e., one operating with an infinitesimally small pinhole, while maintaining high signal-to-noise ratio. Among the reconstruction methods for obtaining the super-resolved image, multi-image deconvolution (MID) and its extension aimed at preserving the optical sectioning capability of confocal microscopy, known as super-resolution sectioning ISM (s$^2$ISM), are among the most widely used approaches. Both methods rely on Richardson–Lucy-type iterative schemes, whose semi-convergent behavior requires early stopping and often leads to noise amplification and reconstruction artifacts. In this work, we introduce a self-tuning explicit regularization framework for both MID and s$^2$ISM reconstruction. Within a Bayesian maximum a posteriori formulation, we combine a multi-frame Poisson data fidelity term with explicit regularization, considering $\ell_1$ and smoothed total variation penalties as representative examples. We further develop an automatic and ground-truth-free strategy for regularization parameter selection by adapting the residual whiteness principle to the multi-frame Poisson setting and introducing a spectral high-pass extension tailored to s$^2$ISM. The resulting framework enables stable reconstructions without empirical stopping rules. To demonstrate the proposed framework, we consider first-order optimization schemes based on proximal gradient and mirror descent methods with adaptive backtracking strategies. Experiments on simulated and real fluorescence ISM datasets demonstrate improved reconstruction stability and image quality with respect to unregularized approaches, while enabling robust super-resolution and optical sectioning in low-photon conditions.