Papers

The fission matrix (FM) method recasts neutron transport in matrix form, enabling fast, interpolation-based reactor calculations from a pre-computed database of Monte Carlo-derived coefficients. Propagation of the underlying Monte Carlo statistical uncertainty to the FM eigenvalue and eigenvector is rarely rigorous: prior approaches assume independent coefficients without testing this against their actual covariance. This work derives a full-covariance Delta-method framework for propagating uncertainty to keff and the fission source, and evaluates it with two cheaper alternatives, a diagonal-covariance approximation and direct resampling of the FM coefficients, using an OpenMC model of the JSI TRIGA Mark-II reactor. For the all-rods-out case, these three estimators and an empirical benchmark of 100 independent Monte Carlo runs agree within 5.7% for keff and 1% for the fission source. The sign of the diagonal approximation's bias flips between the two statistics levels, so it must be checked case by case. With a global fixed-source FM generation strategy, the diagonal, resampling, and independent-run estimators agree within about 2% and lie roughly 10% below the full-covariance estimate, a gap attributed to covariance sampling noise. The large off-diagonal covariance is nearly symmetric and mostly cancels when propagated, rather than being small or negative as previously assumed. For the FM control-rod interpolation and combination strategy (FM-CRd), keff matches an exact perturbed reference within 60 pcm (about two standard deviations) and the fission source within 0.7%, while the diagonal and resampling estimators agree closely and are conservative by about a factor of two. Given this and the prohibitive memory cost of the full covariance, the diagonal Delta method is recommended as the practical default, with the full covariance reserved for one-time characterization of a system.

The primary objective of this paper is to investigate sharp constant problems in harmonic Hardy spaces. We first establish a sharp asymptotic constant for a Littlewood-Paley type inequality applicable to weak harmonic $(K, K_0)$-quasiregular mappings, where $K \geq 1$ and $K_0 \geq 0$. Subsequently, we derive sharp Hardy-Littlewood type constants for harmonic mappings.

We study episodic continuous-time control with an unknown vector control gain, scalar state, quadratic action cost, and smooth convex terminal cost. In the scalar Gaussian experiment, let $\Delta(H)$ denote the minimax improvement over zero control and set $\delta=\sqrt2H^2-1$. We prove the critical law $$ \Delta(H)\asymp \delta^4\sqrt{\log(1/\delta)} \qquad (\delta\downarrow0), $$ with a matching lower and upper bound. The lower bound follows from a uniform deficit–energy inequality valid for fully adaptive controls with unbounded amplitudes, while a moving soft-threshold feedback attains the rate. For local parameters $\theta=N^{-1/4}h$, $|h|\le H$, we show that the normalized minimax regret over $N$ episodes is within $O(N^{-1/2})$ of a fixed-horizon Gaussian sequential control problem, without an additional dimension factor. The terminal task enters the limit only through $c_g=\operatorname{Var}(g(Z))$. The Gaussian problem exhibits an exact no-learning phase transition: $$ C_d^T(H)=TH^2/2 \iff H^4T\le d^2/2, $$ yielding the asymptotic task boundary $c_gH^4=d^2/2$.

Muon is emerging as a promising alternative to AdamW for large-scale neural network training, yet theoretical understanding of its practical implementation remains incomplete, as existing analyses often simplify or omit two key components: (i) practical Newton–Schulz iterations with empirically tuned polynomial coefficients $(3.4445,-4.7750,2.0315)$; and (ii) decoupled weight decay for regularization. In this paper, we provide an optimization interpretation and establish convergence for practical Muon, jointly accounting for both components. Specifically, we interpret practical Muon as right-preconditioned optimization of the original loss with a dynamic weighted $\ell_2$ regularizer that vanishes as stationarity is approached, so that the optimization target remains the original objective. We then establish, to our best knowledge, the first convergence guarantee for practical Muon in the stochastic nonconvex setting, with an $\mathcal{O}(T^{-1/4})$ convergence rate in terms of the expected Frobenius norm of the gradient, improving the dimension dependence of the best known AdamW's convergence rate by a factor of $\sqrt{d}$, where $T$ is the iteration horizon and $d$ is the parameter dimension. Experiments further support the theoretical convergence results.

The main purpose of this paper is to address an open problem posed by Huang, Rasila, and Zhu in [Anal. Math. 48(2022), 1069-1080]. We first establish a sharp Lipschitz continuity result for locally univalent harmonic Bloch mappings with respect to the pseudo hyperbolic metric. This result is then applied to demonstrate that harmonic $(K,K_0)$-quasiregular Bloch type mappings are Lipschitz continuous under the same metric, where $K\geq1$ and $K_{0}\geq0$. Moreover, the estimates obtained are asymptotically sharp as $K \to 1^+$ and $K_0 \to 0^+$. In particular, when $K_0 = 0$, our theorem provides a solution to the aforementioned open problem in the asymptotically sharp sense.

Let ${S}_n(321)$ and ${S}_n(\underline{321})$ denote the sets of $n$-permutations avoiding the classical pattern $321$ and the consecutive pattern $\underline{321}$, respectively. Permutations are in bijection with Lehmer codes, a type of inversion sequence. Using Lehmer codes, we create the corresponding weighted posets ${L}_n(321)$ and ${L}_n(\underline{321})$, where the weight of a code is the inversion number of its permutation. We show that there are $2^{n-2}$ maximal elements of ${L}_n(321)$, while the maximal elements of ${L}_n(\underline{321})$ are enumerated by the Padovan numbers. We also show that when $n$ is even, each of these posets has a unique maximum weight element, and that when $n$ is odd, there are two maximum weight elements. These maximum weight Lehmer codes correspond to the pattern avoiding permutations with maximum inversion number.

Black-box Optimization (BO) is a powerful framework for optimizing expensive objective functions or unknown functions with a limited number of evaluations. A central step of standard BO such as Bayesian optimization is the optimization of a surrogate-based acquisition criterion, which is commonly performed using nonlinear optimization or heuristic search. Therefore, conventional black-box optimization generally does not guarantee global optimality in candidate selection. In this study, we propose Integer Linear Programming-based Black-box Optimization (ILP-BO), a quasi-Bayesian optimization framework that transforms kernel-based surrogate optimization over discrete domains into an Integer Linear Programming (ILP) problem. The key idea is to represent nonlinear kernel functions exactly on finite discrete distance levels by introducing binary one-hot auxiliary variables. This transformation converts the nonlinear surrogate into a linear objective with linear constraints and binary variables. To incorporate exploration while preserving the linear structure, we further introduce a Hamming-distance margin that excludes neighborhoods around previously observed points. We derive the proposed formulation for several standard kernels and obtain an analytical upper bound on the Hamming-distance threshold based on the measure in the binary search space. The resulting candidate-selection problem can be solved by integer programming solvers with certificates of optimality. Thus, our methodology has the potential to serve as a highly transparent black-box optimization framework. Experiments on synthetic and discrete optimization benchmarks show that ILP-BO achieves competitive optimization performance compared with practical Bayesian optimization methods.

Hyun, Kim, Wu and Yue constructed optimal and minimal binary linear codes from order ideals of hierarchical posets with two levels. Two different generalizations of the underlying antichain (simplicial complex) setting to odd characteristic are known: down-sets of $\mathbb{F}_p^n$ under the componentwise order, and support-closed subsets of $\mathbb{F}_q^m$. No generalization of the poset setting itself has appeared. We introduce generalized order ideals of a poset of order $p-1$, obtained by attaching multiplicities in ${0,\dots,p-1}$ to the elements of a poset, and study the two natural notions of order ideal that arise for hierarchical posets with two levels. Whenever the ideal meets the upper level, the resulting defining sets are neither down-sets nor support-closed. We determine the weight distributions of the associated complement codes, exhibit a family of Griesmer codes in which the upper element carries an arbitrary multiplicity, and, via the characteristic function of a generalized order ideal, obtain an infinite family of minimal $p$-ary codes of length $p^n-1$ and dimension $n+1$ violating the Ashikhmin–Barg condition.

We prove an $e$-positive formula for the chromatic symmetric function of proper circular arc graphs solving the $q=1$ case of a conjecture of Ellzey and of Alexandersson and Panova. In doing so, we provide a new proof of the $e$-positivity of unit interval graphs, which alongside Guay-Paquet's reduction gives a new proof of the Stanley–Stembridge conjecture. We define color matrices, which count proper colorings, and tableau matrices, whose entries are nonnegative rational numbers and ratios of elementary symmetric functions. We prove the two matrices are related by a single family of change of basis matrices, which become invertible after restricting to finitely many colors, and we show the chromatic symmetric function of proper circular arc graphs comes from taking the trace of these matrices. Using Hikita's tableaux, this gives an explicit formula for the chromatic symmetric function of a proper circular arc graph as a weighted sum over tableaux whose first and last $k$ vertices lie in the same columns.

We construct maximal subgroups of topological full groups and their commutator subgroups using factor maps and wide open subgroupoids of minimal effective ample groupoids with Cantor unit spaces, assuming almost finiteness or pure infiniteness. For two-to-one coverings and prime factor maps with principal target and fibers of cardinality at most two, surjectivity and injectivity of the same induced map on full-group abelianizations characterize maximality of the factor full group and factor commutator subgroup, respectively. For subgroupoids arising from prime-order cyclic actions free on the unit space, the normalizer of the smaller commutator subgroup in the ambient full group is a semidirect product, maximal exactly when the inclusion induces a surjection on full-group abelianizations. In both constructions, the embedded commutator subgroup has a unique maximal overgroup in the ambient commutator subgroup, namely its normalizer. For Cantor minimal systems, the factor construction yields partition stabilizers of the type in the Grigorchuk-Vorobets conjecture. The subgroupoid normalizers act topologically primitively. Groupoid homology makes the criteria computable for these systems and shifts of finite type. Examples include finitely generated nonsimple and locally finite simple maximal subgroups of infinite index in finitely generated infinite simple amenable groups.

We study the family of functions $V_z:s \mapsto \zeta(s) \exp (zs)$ with numerical experiments. For selected sequences $z_n=x_nu$ on rays, our experiments suggest that suitable fixed points of $V_{z_n}$ approach nontrivial zeros of $\zeta$ in a nearly logarithmic pattern. We formulate this behavior as conjectures with explicit restrictions on the ray. We prove a conditional result identifying finite limits as zeros, and give sufficient hypotheses for the convergence of refinement diagnostics. Substantive AI assistance in the preparation of this draft is disclosed in Section 5.

We prove that every finite loopless signed multigraph $\Sigma=(G,\sigma)$ satisfies $\chi'(\Sigma)\leq\Delta(G)+\mu(G)$, where $\chi'(\Sigma)$ is its chromatic index, and $\Delta(G)$ and $\mu(G)$ are the maximum degree and maximum multiplicity of $G$, respectively. This bound is sharp, even when both positive and negative edges are present, and generalizes both Vizing's theorem for ordinary multigraphs and Behr's theorem for signed simple graphs.

Let $p$ be an odd prime, $q=p^k$, and $r=p^s$ with $s\mid k$. We study $\mathbb{F}_r$-valued generalized Maiorana-McFarland functions $H(x,y,z)=P^{(z)}(x)+\operatorname{Tr}_s^k(yz)$ on $V_n\times\mathbb{F}_q\times\mathbb{F}_q$. The ingredients $P^{(z)}$ are vectorial dual-bent functions with weakly regular components, constant on the cosets of $\mathbb{F}_r^*$ in $\mathbb{F}_q^*$, and, for $z\ne 0$, of inversion type with respect to $t\mapsto t^{-1}$. Under a mild condition on $P^{(0)}$, we prove that the partition into the $\mathbb{F}_r$-valued level sets of $H$, refined according to the $\mathbb{F}_r$-lines of the relevant coordinates, is Fourier-reflexive. It therefore induces a translation association scheme with $\frac{q-1}{r-1}(r+1)+r$ classes (one fewer in a degenerate case), symmetric when the ingredients are even. For $s=1$ this recovers recent scalar constructions, while for $s=k$ it yields $(2q+1)$-class schemes from vectorial bent functions whose components may all be non-weakly regular. We also construct fusions along $\mathbb{F}_r$-subspaces, obtaining $(2r+1)$- and $(3r+2)$-class schemes and connections between different levels. Counterexamples show that the inversion, coset-constancy, subspace, and weak-regularity hypotheses cannot in general be omitted. To our knowledge, these are the first association schemes obtained from vectorial bent functions with non-weakly regular components.

This paper investigates the fractional Laplacian equation in half-spaces, focusing on cases with nonlocal Neumann boundary conditions. We consider both single equations and systems, and establish two key principles: the decay at infinity principle and the narrow region principle. These foundational results, combined with the direct method of moving planes, allow us to derive symmetry and monotonicity properties for decaying solutions. To the best of our knowledge, these are the first results concerning the symmetry of solutions to fractional Laplacian equations with nonlocal Neumann conditions.

For complex Hessian equations on compact Kähler manifolds, we give an alternative proof of the gradient estimate without using blow-up argument, and clarify its dependence on the right-hand side. In addition, using the similar approach, the case of Kähler manifolds with boundary is also investigated.

Plug-and-Play (PnP) methods replace proximal operators with learned denoisers, which produce state-of-the-art reconstruction quality, but sacrifice the variational interpretation and convergence guarantees of proximal methods. Learned Proximal Networks (LPNs) address this problem by designing the denoiser architecture so that it is exactly the proximal operator of a regularizer. In this work, we clarify the theoretical foundations of LPNs and extend the framework to a broader class of activation functions. We then study two practical mechanisms for controlling the learned prior: (i) averaging the LPN with the identity, a common heuristic in PnP methods, and (ii) directly scaling the implicit regularizer induced by the proximal operator. In particular, we characterize the regularizer induced by averaging and develop a convergent method for evaluating the scaled proximal operator.

We study local eigenvalue statistics and dynamical (de)-localization properties for long-range Anderson models, including the fractional Anderson model. These systems are random Anderson-type perturbations of operators with long-range (non random) hopping terms of the form $|T(n,m)| \sim \|n-m\|^{-(d+2\beta)}$ for $\beta>0$, on the $d$-dimensional lattice. In the presence of a strong enough random potential, these models exhibit dense pure point spectrum with polynomially decaying eigenvectors, almost surely. We show that in this strong disorder regime, and for $\beta>d/2$, the local eigenvalue statistics centered at any $E$ in the deterministic spectrum is a Poisson point process with intensity given by the density of states function $n(E)$. Moreover, we prove that, at strong disorder, there is no dynamical localization for these models, verifying a conjecture of Disertori et al. We achieve this by establishing explicit, sharp lower bounds on the Green's functions fractional moments which imply that large moments of the position operator are infinite. Hence, these models exhibit Poisson eigenvalue statistics and dynamical delocalization, that is, the absence of dynamical localization. In particular, this shows that in dimension $d=1$, the fractional Anderson model, a random perturbation of the fractional Laplacian $(-\Delta)^\alpha$ with $\alpha\in (0,1)$, exhibits Poisson eigenvalue statistics and dynamical delocalization if the disorder is strong and the exponent $\frac{1}{2}<\alpha<1$. This is in stark contrast with what is known for the usual Anderson model that has only nearest-neighbor hopping.

We study in this paper an unsupervised, two-step, model-informed deep learning framework for solving the phaseless inverse scattering problem. The objective is to reconstruct a compactly supported function that characterizes a scatterer from boundary measurements of the modulus of the total wave corresponding to multiple incident waves. In the first step, the phaseless data are transformed into an imaging function that is directly related to the underlying scatterer. Motivated by the contrast source method, we derive a coupled system of model equations using the Lippmann-Schwinger integral equation and a spectral-based equation involving the imaging function. In the second step, the unknown scatterer and contrast source function are each parameterized as independent feedforward neural networks, which are trained simultaneously at each iteration using the derived system of model equations. The Lippmann-Schwinger integral equation enforces the underlying physics as a model constraint, while the spectral-based equation incorporating the imaging function supplies geometric information on the shape and location of the unknown scatterer. The resulting framework achieves accurate and efficient reconstructions while maintaining robustness to measurement noise. Furthermore, transforming some aspects of the problem to the spectral domain enables an effective reduction in the computational cost, allowing us to recover the unknown scatterer quickly. Numerical experiments in both two- and three-dimensional settings demonstrate the effectiveness, stability, and practical potential of the proposed approach.

By eliminating the win condition in the game of Connect Four and extending the board to infinite height, a rich state space of positions is obtained. We investigate the number of positions reachable on an $n$-column board after $k$ color-alternating moves. For fixed $k$ we demonstrate polynomiality, derive a partial formula for the polynomial coefficients, and precisely characterize the asymptotic behavior as $n \to \infty$. We then turn our attention to the fixed-$n$ case and show that, under a natural addition operation, positions reachable in an even number of moves form a monoid with a highly symmetric finite generating set; by examining certain free submonoids, we bound the exponential growth rate as $k \to \infty$.

Let $u=(u_n)_{n\ge0}$ be a binary sequence, and define \[ X_0=1,\qquad X_{n+1}=2X_n+u_n \qquad(n\ge0). \] Let $A_u$ be the set of all nonempty finite subset sums of the sequence $(X_n)$, and let $L_u(N)$ denote the maximum length of an arithmetic progression contained in $A_u\cap[1,N]$. We prove that there are absolute constants $c>0$ and $N_0\geq 1$ such that, for every binary sequence $u$, \[ L_u(N)\ge \exp\!\left(c\sqrt{\frac{\log N}{\log\log N}}\right) \] for all $N\geq N_0$. Moreover, if the random variables $u_n$ are independent and uniformly distributed on $\{0,1\}$, then, almost surely, \[ L_u(N)\ll_u N^{2/3}\exp\!\left(C\sqrt{\log N\log\log N}\right), \] where $C>0$ is an absolute constant. Furthermore, every eventually periodic binary sequence satisfies $L_u(N)\gg_u N^{1/2}$.

Chu proved that a symplectic structure $\omega$ on an even-dimensional Lie algebra induces a left-symmetric structure which is defined by $\omega(x\scriptstyle{\Delta}_{\omega} y,z)=-\omega(y,[x,z])$. El Bourkadi and Mansouri proved that a cosymplectic structure on an odd-dimensional Lie algebra $\mathfrak{g}$ induces the left-symmetric structure, which is defined by using a linear isomorphism from $\mathfrak{g}$ to $\mathfrak{g}^*$ associated with the cosymplectic structure. In this paper, for a non-degenerate bilinear form $\phi$ on a Lie algebra $\mathfrak{g}$, we give a necessary and sufficient condition for the product $\scriptstyle{\Delta}_{\phi}$ on $\mathfrak{g}$ defined by $\phi(x\scriptstyle{\Delta}_{\phi} y,z)=-\phi(y,[x,z])$ to be a left-symmetric structure. We also prove that the left-symmetric structure $\scriptstyle{\Delta}_{\phi}$ is complete if and only if the Lie algebra $\mathfrak{g}$ is unimodular. Moreover, we formulate the notion of double extension for non-degenerate bilinear forms and prove that a certain class of non-degenerate bilinear forms is obtained by double extension.

In this paper, we propose an unsupervised, two-step model-informed deep learning framework for solving the inverse source problem. In the first step, boundary measurement data are transformed into an imaging function that encodes information about the shape and location of the unknown source. Leveraging this information, we derive a Fourier-based model equation in the spectral domain that relates the imaging function to the Fourier coefficients of the source function. In the second step, this model equation is incorporated into the training of a model-informed neural network to recover the parameters of interest. The proposed approach enables fast and accurate reconstruction of the source function while maintaining robustness to noise. The effectiveness of the method is demonstrated through numerical experiments in both two- and three-dimensional settings. For the two-dimensional case, we further compare our approach with the traditional least-squares method to validate its computational efficiency and reconstruction accuracy.

The Berger-Charpin conjecture predicts that the permutation automorphism group of a cyclic code is generally the $q$-affine group, which is generated by the shift and the Frobenius multiplier on the index set. The word generally is not made precise in the literature, and the goal of this series of papers is to study when the Berger-Charpin conjecture or its variants hold. In this first part, we establish a natural connection between the theory of cyclic codes and the theory of cyclotomic association schemes (equivalently, Schur rings over cyclic groups). Using it, we show that the $q$-affine group is not the correct group to expect for certain code lengths. We characterize these lengths by an arithmetic condition, and show that they have density zero.

The lower asymptotic proportion of distinct nontrivial zeros of the Riemann zeta function, relative to the total number counted with multiplicity, is at least $0.8369928814\ldots$. Earlier work proves $0.83625\ldots$, and a report we have not verified claims $0.83672\ldots$. As in the proof of the bound $0.83625$, an unconditional version of Montgomery's pair-correlation theorem gives an asymptotic energy estimate. The new ingredient is a short matrix inequality with a free clipping parameter. It strengthens the lower bound for this energy in terms of the number of distinct zeros. The gain is a nonnegative correction from overlaps between different nearby zeros on the critical line, which is retained even when some of these zeros are double. The matrix inequality, the threshold lemma, the block dichotomy, the counting assembly and the exact arithmetic are proved formally in Lean 4. The constant relies on a computer-assisted local inequality from recent work that has not yet been refereed. That computation was re-run independently, and every imported input is listed. This paper is primarily an experiment in AI-assisted mathematical research (Section 4).

A graph G is uniquely C_t-saturated if G contains no cycle of length t and, for every edge e of the complement, G+e contains exactly one cycle of length t; it is nontrivial if it has at least t vertices. Wenger and West proved that no nontrivial uniquely C_6- or C_7-saturated graphs exist, and conjectured the same for every t >= 6; the case t=8 was verified but never published, and t >= 9 has remained open. We prove the conjecture for t=9: there is no nontrivial uniquely C_9-saturated graph. The proof is a case analysis on the length L of a longest even cycle of length at most 12 (L is 4, 6, 8, 10, or 12). The cases L=12, L=10, and L=4 are settled by hand; for L=6 and L=8, hand classifications of the components outside the cycle reduce each case to a bounded finite family of configurations, eliminated by a short, replayable computer enumeration with known-answer controls and completeness certified by zero cap hits. Along the way we prove, in sharpened form, the t=9 instance of a structural lemma Wenger and West stated without proof.

A global classification of the asymptotic oscillatory dynamics is established for a symmetric zero-divergence class of three-dimensional piecewise-linear Filippov systems, up to a set of initial conditions of zero Lebesgue measure. The affine fields are related by an involution, and the switching plane contains a visible–visible two-fold. In canonical coordinates, the eigenvalues are \(\mu\pm i\) and \(-2\mu\), while \(H\) measures the focal-line inclination. For every \(\mu>0\), a simple-period crossing cycle exists if and only if \(H\in\mathcal I_\mu\), and is unique, symmetric, hyperbolic, and orbitally asymptotically stable. Its half-period parametrizes \(\mathcal I_\mu\) and determines the crossing points, period, and Floquet multipliers. Global dissipation excludes crossing cycles with any higher number of crossings and makes the classified cycle the \(\omega\)-limit set of every crossing-only trajectory. If attractive sliding has no interior pseudo-equilibria, sliding is transient unless the trajectory reaches the two-fold. The initial conditions leading to the two-fold lie in a countable union of analytic surfaces and curves. Under sufficiently small perturbations in the Whitney \(C^1\) topology, the cycle persists as the unique simple-period crossing limit cycle and remains hyperbolic and orbitally asymptotically stable.

We formalize in Lean 4 the Galerkin construction of global Leray–Hopf weak solutions for the unforced two-dimensional Navier–Stokes equations on arbitrary bounded open domains with no-slip boundary conditions. For every initial datum in the solenoidal \(L^2\) velocity space, the solution has a weakly continuous velocity path, satisfies the weak equation locally in time, and obeys the energy inequality at every time. For rectangles, we also formalize finite-horizon solutions with continuous forcing in the dual energy space. The development includes divergence-free graph spaces, compact spectral coordinates, transport cancellation, Ladyzhenskaya estimates, and simultaneous space–time compactness. The abstract Hilbert-space results apply to evolution equations with a compact energy-to-velocity embedding and a skew trilinear nonlinearity.

In this note we prove $\mathrm{forb}(n,\mathbf{1}_{r-2}\times C_4)=\Theta(n^{r-1/2})$, thus providing a series of counterexamples of a 2005 conjecture of the first two authors. Here $\mathbf{1}_{r-2}\times C_4$ is the $r$-extension of the 4-cycle, that is the $r$-uniform hypergraph of four edges obtained by joining a set of $r-2$ vertices to each edge of $C_4$. $\mathrm{forb}(n,F)$ is the largest number of edges of a simple hypergraph $G$ on $n$ vertices without trace (induced subhypergraph) $F$.

We prove a conjecture proposed by Robert J. Lang over 20 years ago, stating that in a circle-packed origami design, any axial polygon with 5 or more sides can be reduced by properly adding a new circle. The proof is based on the structure of the medial axis of the region enclosed by the circles, and yields a concrete algorithm for finding the added circle. The result extends to configurations involving rivers.

This study develops an energy formulation for linear anisotropic Maxwell–Cole–Cole media with multiple fractional polarization channels, passive conductivity, prescribed currents, and finite-terminal Caputo memory with zero prehistory. A weighted solution class captures the integrable fractional initial layer generated by quiescent polarization. Fractional polarization power then separates into a nonnegative history storage and a nonnegative constitutive dissipation rate. Although the instantaneous rate may diverge at initialization, it remains integrable, while stored history energy and accumulated loss vanish initially. Combined with the Poynting theorem, this decomposition yields a global field-plus-memory identity that establishes electromagnetic passivity, closed-system energy decay, and conditional uniqueness under zero-supply boundary conditions. The framework also recovers the Debye, instantaneous-dielectric, vanishing-susceptibility, and vacuum limits under convergence assumptions. A positive diffusive representation provides a second decomposition whose storage never exceeds the history storage. Their difference is an exact nonnegative weighted variance that vanishes only for a zero polarization history. In agreement with these time-domain results, the finite-time harmonic response approaches the anisotropic multichannel Cole–Cole permittivity, whose Hermitian imaginary part ensures nonnegative cycle-averaged absorption. A reproducible single-channel study shows that fractional order broadens the dielectric-loss band, redistributes wave attenuation, strengthens the integrable initial power layer, and prolongs the step-response tail. An exact electric-step response verifies both energy representations. The formulation connects fractional initial layers, hereditary storage, passivity, diffusive relaxation, dispersion, absorption, and attenuation in a local constitutive model.