Daily Papers

Bessel functions are critical in scientific computing for applications such as machine learning, protein structure modeling, and robotics. However, currently, available routines lack precision or fail for certain input ranges, such as when the order v is large, and GPU-specific implementations are limited. We address the precision limitations of current numerical implementations while dramatically improving the runtime. We propose two novel algorithms for computing the logarithm of modified Bessel functions of the first and second kinds by computing intermediate values on a logarithmic scale. Our algorithms are robust and never have issues with underflows or overflows while having relative errors on the order of machine precision, even for inputs where existing libraries fail. In C++/CUDA, our algorithms have median and maximum speedups of 45x and 6150x for GPU and 17x and 3403x for CPU, respectively, over the ranges of inputs and third-party libraries tested. Compared to SciPy, the algorithms have median and maximum speedups of 77x and 300x for GPU and 35x and 98x for CPU, respectively, over the tested inputs. The ability to robustly compute a solution and the low relative errors allow us to fit von Mises-Fisher, vMF, distributions to high-dimensional neural network features. This is, e.g., relevant for uncertainty quantification in metric learning. We obtain image feature data by processing CIFAR10 training images with the convolutional layers of a pre-trained ResNet50. We successfully fit vMF distributions to 2048-, 8192-, and 32768-dimensional image feature data using our algorithms. Our approach provides fast and accurate results while existing implementations in SciPy and mpmath fail to fit successfully. Our approach is readily implementable on GPUs, and we provide a fast open-source implementation alongside this paper.

Graph neural networks have recently achieved great successes in predicting quantum mechanical properties of molecules. These models represent a molecule as a graph using only the distance between atoms (nodes). They do not, however, consider the spatial direction from one atom to another, despite directional information playing a central role in empirical potentials for molecules, e.g. in angular potentials. To alleviate this limitation we propose directional message passing, in which we embed the messages passed between atoms instead of the atoms themselves. Each message is associated with a direction in coordinate space. These directional message embeddings are rotationally equivariant since the associated directions rotate with the molecule. We propose a message passing scheme analogous to belief propagation, which uses the directional information by transforming messages based on the angle between them. Additionally, we use spherical Bessel functions and spherical harmonics to construct theoretically well-founded, orthogonal representations that achieve better performance than the currently prevalent Gaussian radial basis representations while using fewer than 1/4 of the parameters. We leverage these innovations to construct the directional message passing neural network (DimeNet). DimeNet outperforms previous GNNs on average by 76% on MD17 and by 31% on QM9. Our implementation is available online.

We construct the quantized electromagnetic field and the corresponding photon propagators in a rigidly rotating cylindrical reference frame. Starting from covariant Maxwell theory, we formulate the gauge-fixed field equations and solve them in cylindrical coordinates. The vector nature of the electromagnetic field requires distinguishing orbital and intrinsic angular momentum, with the total angular momentum as the natural quantum number. We construct the two independent physical photon modes in terms of Bessel functions by solving the coupled differential equations, and provide an alternative construction based on a covariant basis of vector eigenfunctions generated by a covariantly constant tetrad. After normalizing the modes with the Maxwell inner product, we quantize the electromagnetic field and define the corresponding photon states and greater and lesser two-point functions. The physical time-ordered propagator is obtained directly from the quantized modes, with its transverse dependence encoded in polarization kernels that retain the cylindrical spatial structure and its spectral poles exhibiting the rotational energy shift. The full gauge-fixed photon propagator is constructed in two complementary ways: from the scalar Green function combined with the tensor transformation of the inertial Feynman-gauge propagator, and from the covariant vector eigensystem through its spectral decomposition and the Fock--Schwinger proper-time representation. We verify that both constructions satisfy the Green-function equation for the complete gauge-fixed Maxwell operator, providing the propagator required for covariant perturbative calculations with internal photon lines.

Long short-term memory (LSTM) models are a particular type of recurrent neural networks (RNNs) that are central to sequential modeling tasks in domains such as urban telecommunication forecasting, where temporal correlations and nonlinear dependencies dominate. However, conventional LSTMs suffer from high parameter redundancy and limited nonlinear expressivity. In this work, we propose the Quantum-inspired Kolmogorov-Arnold Long Short-Term Memory (QKAN-LSTM), which integrates Data Re-Uploading Activation (DARUAN) modules into the gating structure of LSTMs. Each DARUAN acts as a quantum variational activation function (QVAF), enhancing frequency adaptability and enabling an exponentially enriched spectral representation without multi-qubit entanglement. The resulting architecture preserves quantum-level expressivity while remaining fully executable on classical hardware. Empirical evaluations on three datasets, Damped Simple Harmonic Motion, Bessel Function, and Urban Telecommunication, demonstrate that QKAN-LSTM achieves superior predictive accuracy and generalization with a 79% reduction in trainable parameters compared to classical LSTMs. We extend the framework to the Jiang-Huang-Chen-Goan Network (JHCG Net), which generalizes KAN to encoder-decoder structures, and then further use QKAN to realize the latent KAN, thereby creating a Hybrid QKAN (HQKAN) for hierarchical representation learning. The proposed HQKAN-LSTM thus provides a scalable and interpretable pathway toward quantum-inspired sequential modeling in real-world data environments.

In this paper, we investigate the quantum behaviour near the classical singularity of a black hole. We describe the black-hole interior by the Kantowski-Sachs metric, with a string propagating in this region. This approach is motivated by the fact that, near the classical singularity, spacetime curvature reaches Planck or string scales, necessitating a description in terms of string-theoretic degrees of freedom. The string dynamics is formulated via the Polyakov action, and the combined system is canonically quantized to derive the corresponding Wheeler-DeWitt equation. Although the equation is separable in the radial and azimuthal string sectors, the remaining geometric and string degrees of freedom are nontrivially coupled. To address this, we employ a perturbative approach and obtain the wave function up to first order. At the zeroth-order, the angular part is described by the general Heun function, while the geometric part is expressed in terms of the Bessel function. The first-order wave function is governed by a source term containing the zeroth-order solution, and we obtain its explicit expression in terms of the complementary function and the particular integral. We confirm regularity of the wave function at both zeroth and first order, at the singularity as well as at the horizon. Regularity at the horizon imposes a discretization condition on one eigenvalue, while another eigenvalue is quantized due to the azimuthal periodicity of the wave function. Importantly, both zeroth- and first-order solutions satisfy the DeWitt boundary condition, implying resolution of the black hole singularity within this quantum framework.

We derive the tree-level interaction vertices of the linear sigma model with quarks in the presence of global rotation. Starting from the scalar and fermionic field solutions in a rotating background, we construct the corresponding spectral representations and use them to evaluate the mesonic self-interactions and Yukawa couplings of the model. We show that rotation modifies the usual momentum-space structure of the interaction vertices: conservation of longitudinal momentum is accompanied by angular-momentum selection rules, while the transverse dynamics is encoded in rotation-induced radial form factors involving products of Bessel functions. The resulting expressions provide a set of interaction rules adapted to the symmetries and finite geometry of the rotating system. Although derived within the linear sigma model with quarks, the structures obtained for scalar self-interactions and scalar-fermion couplings can be extended to a broader class of quantum field theories under global rotation. These results provide the building blocks for perturbative calculations and contribute toward a microscopic quantum-field-theoretical description of interacting matter in rotating environments.

Formulas involving fundamental mathematical constants had a great impact on various fields of science and mathematics, for example aiding in proofs of irrationality of constants. However, the discovery of such formulas has historically remained scarce, often perceived as an act of mathematical genius by great mathematicians such as Ramanujan, Euler, and Gauss. Recent efforts to automate the discovery of formulas for mathematical constants, such as the Ramanujan Machine project, relied on exhaustive search. Despite several successful discoveries, exhaustive search remains limited by the space of options that can be covered and by the need for vast amounts of computational resources. Here we propose a fundamentally different method to search for conjectures on mathematical constants: through analysis of integer sequences. We introduce the Enumerated Signed-continued-fraction Massey Approve (ESMA) algorithm, which builds on the Berlekamp-Massey algorithm to identify patterns in integer sequences that represent mathematical constants. The ESMA algorithm found various known formulas for e, e^2, tan(1), and ratios of values of Bessel functions. The algorithm further discovered a large number of new conjectures for these constants, some providing simpler representations and some providing faster numerical convergence than the corresponding simple continued fractions. Along with the algorithm, we present mathematical tools for manipulating continued fractions. These connections enable us to characterize what space of constants can be found by ESMA and quantify its algorithmic advantage in certain scenarios. Altogether, this work continues in the development of augmenting mathematical intuition by computer algorithms, to help reveal mathematical structures and accelerate mathematical research.

We prove a single-logarithm bound for the eigenvalues in the plunge region of Fourier concentration operators, removing the extra logarithmic factor in the Kulikov--Dam Larsen bound for two general admissible sets throughout each fixed exponential window. For bounded spatial and frequency sets A,BsubsetR^d of positive measure with finite upper codimension-one Minkowski boundary content, the number Λ_varepsilon(cA,B) of eigenvalues in (varepsilon,1-varepsilon) satisfies Λ_varepsilon(cA,B)le C(d,A,B,β)c^{d-1}D(1+log_+(c/D)), where D=log(1/(varepsilon(1-varepsilon))), for every fixed β>0, cge2, 0<varepsilon<1/2, and Dleβc, with log_+t=max(0,log t). The argument uses positive band-limited smoothing and requires no symmetry, convexity, or connectedness. In dimension two, the corresponding exponential-window estimate has a constant uniform in the window parameter. We also determine the optimal coordinate-Lipschitz constant d-1 for recursive decompositions of balls and standard simplices in every dimension, and obtain box-frequency estimates with explicit geometric constants. For C^1 convex bodies in three dimensions, the optimal infimum of recursive coordinate slopes is 2. Uniform real-order Bessel estimates give independent disk and ball bounds at every spectral threshold, including an explicit deep-tail term. For smooth strictly convex planar frequency bodies with positive curvature, we prove correlation estimates for aligned Gaussian packets at arbitrary relative scales and quantitative lower bounds on the cost of estimating spatial blocks separately. Together these results connect the global spectral bound with sharp geometric constructions and the cancellation retained by grouping boundary interactions.

We study the degree of reliability of extrapolation of complex electromagnetic permittivity functions based on their analyticity properties. Given two analytic functions, representing extrapolants of the same experimental data, we examine how much they can differ at an extrapolation point outside of the experimentally accessible frequency band. We give a sharp upper bound on the worst case extrapolation error, in terms of a solution of an integral equation of Fredholm type. We conjecture and give numerical evidence that this bound exhibits a power law precision deterioration as one moves further away from the frequency band containing measurement data.

In this work, we focus on a fractional differential equation in Riesz form discretized by a polynomial B-spline collocation method. For an arbitrary polynomial degree p, we show that the resulting coefficient matrices possess a Toeplitz-like structure. We investigate their spectral properties via their symbol and we prove that, like for second order differential problems, also in this case the given matrices are ill-conditioned both in the low and high frequencies for large p. More precisely, in the fractional scenario the symbol has a single zero at 0 of order α, with α the fractional derivative order that ranges from 1 to 2, and it presents an exponential decay to zero at π for increasing p that becomes faster as α approaches 1. This translates in a mitigated conditioning in the low frequencies and in a deterioration in the high frequencies when compared to second order problems. Furthermore, the derivation of the symbol reveals another similarity of our problem with a classical diffusion problem. Since the entries of the coefficient matrices are defined as evaluations of fractional derivatives of the B-spline basis at the collocation points, we are able to express the central entries of the coefficient matrix as inner products of two fractional derivatives of cardinal B-splines. Finally, we perform a numerical study of the approximation behavior of polynomial B-spline collocation. This study suggests that, in line with non-fractional diffusion problems, the approximation order for smooth solutions in the fractional case is p+2-α for even p, and p+1-α for odd p.

We study a space-time Brownian motion with drift B(t)=(t_0+t,y_0+W(t)+t) killed at the moving boundary of the cone {(t,x):0<x<t}. This article determines the parabolic Martin boundary and all harmonic functions associated with this process. To that end, the asymptotics of Green's functions are determined along all directions. We also find the exit probabilities at the edges, the probability of remaining in the cone forever and the laws of the exit point and exit time. From this, we derive an explicit formula for the transition kernel of the process. These results arise from two different methods initially introduced to study random walks. An analytical approach, developed in the 1970s by Malyshev and based on the steepest descent method on a Riemann surface, is used to determine the asymptotics of the Green's functions. A recursive compensation approach, inspired by the method developed in the 1990s by Adan, Wessels and Zijm, is used to determine the harmonic functions.

Zipf's law appears in many application areas but does not have a closed form expression, which may make its use cumbersome. Since it coincides with the truncated version of the Zeta distribution, in this paper we propose three approximate closed form expressions for the truncated Zeta distribution, which may be employed for Zipf's law as well. The three approximations are based on the replacement of the sum occurring in Zipf's law with an integral, and are named respectively the integral approximation, the average integral approximation, and the trapezoidal approximation. While the first one is shown to be of little use, the trapezoidal approximation exhibits an error which is typically lower than 1\%, but is as low as 0.1\% for the range of values of the Zipf parameter below 1.

A comparative analysis of the special shapes (patterns, profiles) of the eclipses applied for the phenomenological modeling of the light curves of eclipsing binary stars is conducted. Families of functions are considered, generalizing local approximations (Andronov, 2010, 2012) and the functions theoretically unlimited in a width, based on a Gaussian (Mikulasek, 2015). For an analysis, the light curve of the star V0882 Car = 2MASS J11080308 - 6145589 of the classic Algol - subtype (\beta Persei) is used. By analyzing dozens of modified functions with additional parameters, it was chosen the 14 best ones according to the criterion of the least sum of squares of deviations. The best are the functions with an additional parameter, describing profiles, which are limited in phase.

Wentzel, Kramers, Brillouin (WKB) approximation for fractional systems is investigated in this paper using the fractional calculus. In the fractional case the wave function is constructed such that the phase factor is the same as the Hamilton's principle function "S". To demonstrate our proposed approach two examples are investigated in details.

We discuss methods for modeling eclipsing binary stars using the "physical", "simplified" and "phenomenological" models. There are few realizations of the "physical" Wilson-Devinney (1971) code and its improvements, e.g. Binary Maker, Phoebe. A parameter search using the Monte-Carlo method was realized by Zola et al. (2010), which is efficient in expense of too many evaluations of the test function. We compare existing algorithms of minimization of multi-parametric functions and propose to use a "combined" algorithm, depending on if the Hessian matrix is positively determined. To study methods, a simply fast-computed function resembling the "complete" test function for the physical model. Also we adopt a simplified model of an eclipsing binary at a circular orbit assuming spherical components with an uniform brightness distribution. This model resembles more advanced models in a sense of correlated parameter estimates due to a similar topology of the test function. Such a model may be applied to detached Algol-type systems, where the tidal distortion of components is negligible.

Earlier, in the CINT program at Los Alamos National Laboratory, we focused ultrafast mode-locked lasers on the tip-sample junction of a scanning tunneling microscope to generate currents at hundreds of harmonics of the laser pulse repetition frequency. Each harmonic has a signal-to-noise ratio of 20 dB with a 10-dB linewidth of only 3 Hz. Now we model closed quantum nanocircuits with rectangular, triangular, or delta-function barrier, shunted by a beryllium filament for quasi-coherent electron transport over mean-free paths as great as 68 nm. The time-independent Schrödinger equation is solved with the boundary conditions that the wavefunction and its derivative are continuous at both connections. These four boundary conditions are used to form a four-by-four complex matrix equation with only zeros in the right-hand column vector which is required to have a non-trivial solution with each of the closed nanocircuits. Each model has four parameters: (1) the barrier length, (2) the height and shape of the barrier, (3) the length of the pre-barrier, and (4) the electron energy. Any three of these may be specified and then the fourth is varied to bring the determinant to zero to find the solutions on lines or surfaces in the space defined by the four parameters. First, we use a simplistic model having a rectangular barrier. The second model has a triangular barrier as a first approximation to field emission, and we are considering applying this approach for a self-contained nanoscale extension of our earlier effort to generate the harmonics at Los Alamos. The third model has a delta-function barrier, and the fourth model is an extension of the first one where the width of the rectangular barrier is varied inversely with its height.

We consider N classical particles interacting via the Coulomb potential in spatial dimension d and in the presence of an external trap, at equilibrium at inverse temperature beta. In the large N limit, the particles are confined within a droplet of finite size. We study smooth linear statistics, i.e. the fluctuations of sums of the form {cal L}_N = sum_{i=1}^N f({bf x}_i), where {bf x}_i's are the positions of the particles and where f({bf x}_i) is a sufficiently regular function. There exists at present standard results for the first and second moments of {cal L}_N in the large N limit, as well as associated Central Limit Theorems in general dimension and for a wide class of confining potentials. Here we obtain explicit expressions for the higher order cumulants of {cal L}_N at large N, when the function f({bf x})=f(|{bf x}|) and the confining potential are both rotationnally invariant. A remarkable feature of our results is that these higher cumulants depend only on the value of f'(|{bf x}|) and its higher order derivatives evaluated exactly at the boundary of the droplet, which in this case is a d-dimensional sphere. In the particular two-dimensional case d=2 at the special value beta=2, a connection to the Ginibre ensemble allows us to derive these results in an alternative way using the tools of determinantal point processes. Finally we also obtain the large deviation form of the full probability distribution function of {cal L}_N.

We establish a unified framework for understanding the results on the Weil quadratic form obtained by Yoshida (1992), Bombieri (2001, 2003), Connes--Consani (2023), and Connes--Consani--Moscovici (2025+) from the perspective of the screw function introduced in Suzuki (2023). An advantage of the approach via the screw function is that it provides a method to study the Weil quadratic form, which is originally defined in terms of distributions, by means of continuous functions. Based on this framework, we formulate a conjecture stating that a self-adjoint operator whose eigenvalues are the imaginary parts of the nontrivial zeros of the Riemann zeta function can be obtained as the limit, as a to infty, of self-adjoint operators arising from nonlocal realizations of the first-order differential operator on the finite interval [-a,a]. All these results are obtained without assuming the Riemann Hypothesis. This conjecture may be compared with the limit formula for the Riemann zeta function expressed in terms of zeta-regularized products proposed by Connes, Consani, and Moscovici, and it sheds new light on the spectral-theoretic interpretation of the nontrivial zeros of the Riemann zeta function.

Structural survival of offshore structures is crucial for the growing marine economy. Calculating the added mass, radiation damping, and excitation coefficients to quantify wave loads with the traditional boundary element method (BEM) presents a computational bottleneck. The matched eigenfunction expansion method (MEEM), a long-known but rarely-used alternative, offers computational benefits due to its semi-analytical nature. However, previous work fails to directly compare its accuracy and computational performance with BEM, leaving the extent of its utility unknown. Furthermore, the geometry-dependent convergence for cylindrical and slanted geometries has not yet been documented, making the method's practicality for general geometries unclear. This paper presents a unifying MEEM framework for modeling an arbitrary number of fixed or heaving surface-piercing annular cylinders with continuous and radially-monotonic body profiles, and explores the method's block matrix structure, convergence behavior, ability to accurately approximate slanted geometries, and computational advantages over the BEM solver Capytaine. The numerical experiments show that MEEM can compute hydrodynamic coefficients of slanted geometries within 5% of Capytaine, even for angles as steep as 15 degrees from vertical. Finally, MEEM can achieve 2% convergence of its hydrodynamic coefficients an order of magnitude faster than Capytaine with a matrix size two orders of magnitude smaller, making it a computationally effective alternative to traditional BEM solvers. These contributions enable hydrodynamic analysis of a broad range of shapes with increased speed and confidence, paving the way for future optimization studies to yield improved designs.

Schrödinger operators with nonlocal δ-interaction are studied with the use of the Lax-Phillips scattering theory methods. The condition of applicability of the Lax-Phillips approach in terms of non-cyclic functions is established. Two formulas for the S-matrix are obtained. The first one deals with the Krein-Naimark resolvent formula and the Weyl-Titchmarsh function, whereas the second one is based on modified reflection and transmission coefficients. The S-matrix S(z) is analytical in the lower half-plane C_- when the Schrödinger operator with nonlocal δ-interaction is positive self-adjoint. Otherwise, S(z) is a meromorphic matrix-valued function in C_- and its properties are closely related to the properties of the corresponding Schrödinger operator. Examples of S-matrices are given.

We study a paradigm of coding for compression of the natural numbers via the zeta distribution and develop a statistical-mechanical interpretation, both in terms of Hagedorn systems and a Bose gas with energy levels given by logarithms of prime numbers. We also propose a simple coding scheme for the zeta distribution that nearly achieves the ideal code length. For block coding of vectors of natural numbers, we derive the micro-canonical entropy function and demonstrate its asymptotic linearity implying that its behavior is analogous to that of a Hagedorn system. We also derive the large deviations rate function, and provide a formula for the best coding parameter in the large deviations sense. We show that due the Hagedorn-type phase transition there is only partial equivalence of ensembles, due to the degeneration of the domain of the partition function.

We introduce featom, an open source code that implements a high-order finite element solver for the radial Schr\"odinger, Dirac, and Kohn-Sham equations. The formulation accommodates various mesh types, such as uniform or exponential, and the convergence can be systematically controlled by increasing the number and/or polynomial order of the finite element basis functions. The Dirac equation is solved using a squared Hamiltonian approach to eliminate spurious states. To address the slow convergence of the kappa=pm1 states due to divergent derivatives at the origin, we incorporate known asymptotic forms into the solutions. We achieve a high level of accuracy (10^{-8} Hartree) for total energies and eigenvalues of heavy atoms such as uranium in both Schr\"odinger and Dirac Kohn-Sham solutions. We provide detailed convergence studies and computational parameters required to attain commonly required accuracies. Finally, we compare our results with known analytic results as well as the results of other methods. In particular, we calculate benchmark results for atomic numbers (Z) from 1 to 92, verifying current benchmarks. We demonstrate significant speedup compared to the state-of-the-art shooting solver dftatom. An efficient, modular Fortran 2008 implementation, is provided under an open source, permissive license, including examples and tests, wherein particular emphasis is placed on the independence (no global variables), reusability, and generality of the individual routines.

Studies of density matrices for random quantum states lead naturally to the fixed trace Laguerre ensemble in random matrix theory. Previous studies have uncovered explicit rational function formulas for moments of purity statistic (trace of the squared density matrix), and also a third order linear differential equation satisfied by the eigenvalue density. We further probe the origin of these results from the viewpoint of integrability, which is taken here to mean wider classes of recursions and differential equations, and give extensions. Prominent in our study are first order linear matrix differential equations. One application given is to the derivation of the third order scalar equation for the density. Another is to obtain the explicit rational function formula for the variance of the purity statistic in the β generalised fixed trace Laguerre ensemble. In the original case (β= 2), the purity cumulants are expressed in terms of the large argument expansion of a particular σ-Painlevé IV transcendent. In a different but related direction, the exact computation of the two-point correlation for the fixed determinant circular unitary ensemble SU(N) is given the Appendix.

We prove new large value estimates for the Riemann zeta function on the critical line. For instance, we improve the upper bound on the measure of tin [T, 2T] with |ζ(1/2 + it)| geq T^{1/8} for the first time since Hardy-Littlewood (1923). Our results imply improved upper bounds on the k-th moment of zeta for every 4 < k leq 12. We show in particular that equation* \int_{0}^{T} |ζ(1{2} + it)|^6 \,d t \ll_\varepsilon T^{5{4} - 1{60} + \varepsilon}. equation* The main ingredient is a new large value estimate for exponential sums with square-root phases, and certain perturbations. Those arise in the expansion of the short second moment of zeta.

In this paper, we prove the existence of stationary pulse solutions for a bulk-surface reaction-diffusion system modeling cell polarization in two-dimensional domains. Using a singular perturbation method, we construct pulse solutions whose centers lie near non-degenerate critical points of a geometry-dependent potential expressed in terms of the Neumann Green's function, as formally predicted in watanabe2026. The main analytical difficulty arises from the half-Laplacian induced by the logarithmic singularity of the Green's function, resulting in algebraically rather than exponentially decaying inner solutions. Our analysis may also serve as a useful guide for applying singular perturbation methods to other bulk-surface type systems.

In this article, we continue the development of the Riemann-Hilbert formalism for studying the asymptotics of Toeplitz+Hankel determinants with non-identical symbols, which we initiated in GI. In GI, we showed that the Riemann-Hilbert problem we formulated admits the Deift-Zhou nonlinear steepest descent analysis, but with a special restriction on the winding numbers of the associated symbols. In particular, the most natural case, namely zero winding numbers, is not allowed. A principal goal of this paper is to develop a framework that extends the asymptotic analysis of Toeplitz+Hankel determinants to a broader range of winding-number configurations. As an application, we consider the case in which the winding numbers of the Szego-type Toeplitz and Hankel symbols are zero and one, respectively, and compute the asymptotics of the norms of the corresponding system of orthogonal polynomials.

We study simplified bootstrap problems for probability distributions on the infinite line and the circle. We show that the rapid convergence of the bootstrap method for problems on the infinite line is related to the fact that the smallest eigenvalue of the positive matrices in the exact solution becomes exponentially small for large matrices, while the moments grow factorially. As a result, the positivity condition is very finely tuned. For problems on the circle we show instead that the entries of the positive matrix of Fourier modes of the distribution depend linearly on the initial data of the recursion, with factorially growing coefficients. By positivity, these matrix elements are bounded in absolute value by one, so the initial data must also be fine-tuned. Additionally, we find that we can largely bypass the semi-definite program (SDP) nature of the problem on a circle by recognizing that these Fourier modes must be asymptotically exponentially small. With a simple ansatz, which we call the shoestring bootstrap, we can efficiently identify an interior point of the set of allowed matrices with much higher precision than conventional SDP bounds permit. We apply this method to solving unitary matrix model integrals by numerically constructing the orthogonal polynomials associated with the circle distribution.

This paper presents the first numerical realization of Suzuki's Weil-Quadratic-Form operator, a candidate for the Hilbert--Pólya program linking spectral positivity to the Riemann Hypothesis (RH). Suzuki's 2026 construction was purely theoretical; here, the operator is instantiated via P1 finite-element discretization and Richardson extrapolation. Key results include: (R1) In the prime-free regime, the spectrum follows a closed Archimedean law A_k(a) = log(1/a) + log(k-2) + B_0 + O(a), with B_0 = log q - 2log 2, confirmed to 30-digit precision. (R2) A Mellin double-pole argument proves the head coefficient B(ν) and shows B_0 depends only on the conductor q, independent of the Archimedean parameter. (R2b) The degree d of an L-function appears directly as the logarithmic slope of the spectrum. (R3) Total spectral intensity follows the prime number theorem, S(a) sim (2a)^3/6. (R4) Nontrivial zeros are not eigenvalues but occur in the explicit-formula error term of the prime symbol. (R5) The best-match line σ^*(a) descends toward the critical line. (R6) Weil's positivity criterion is realized in operator form: bounded residual growth corresponds to all zeros on the line, while an injected off-line zero causes exponential blow-up. (R7) The lowest eigenvalue λ_1(a) is strictly positive, decays superexponentially, and passes smoothly through the first prime threshold. (R8) The characteristic function W(a,0;z) is computed for the first time, with all zeros confirmed real. (R9) Indirect traces of GUE statistics appear in the moment structure, even where direct detection is blocked. The authors emphasize that this work does not prove RH. All results are Archimedean and universal, with significance lying in the faithful numerical realization of classical identities rather than new arithmetic.

We have developed a novel activation function, named the Cauchy Activation Function. This function is derived from the Cauchy Integral Theorem in complex analysis and is specifically tailored for problems requiring high precision. This innovation has led to the creation of a new class of neural networks, which we call (Comple)XNet, or simply XNet. We will demonstrate that XNet is particularly effective for high-dimensional challenges such as image classification and solving Partial Differential Equations (PDEs). Our evaluations show that XNet significantly outperforms established benchmarks like MNIST and CIFAR-10 in computer vision, and offers substantial advantages over Physics-Informed Neural Networks (PINNs) in both low-dimensional and high-dimensional PDE scenarios.

A single two-input gate suffices for all of Boolean logic in digital hardware. No comparable primitive has been known for continuous mathematics: computing elementary functions such as sin, cos, sqrt, and log has always required multiple distinct operations. Here I show that a single binary operator, eml(x,y)=exp(x)-ln(y), together with the constant 1, generates the standard repertoire of a scientific calculator. This includes constants such as e, pi, and i; arithmetic operations including addition, subtraction, multiplication, division, and exponentiation as well as the usual transcendental and algebraic functions. For example, exp(x)=eml(x,1), ln(x)=eml(1,eml(eml(1,x),1)), and likewise for all other operations. That such an operator exists was not anticipated; I found it by systematic exhaustive search and established constructively that it suffices for the concrete scientific-calculator basis. In EML (Exp-Minus-Log) form, every such expression becomes a binary tree of identical nodes, yielding a grammar as simple as S -> 1 | eml(S,S). This uniform structure also enables gradient-based symbolic regression: using EML trees as trainable circuits with standard optimizers (Adam), I demonstrate the feasibility of exact recovery of closed-form elementary functions from numerical data at shallow tree depths up to 4. The same architecture can fit arbitrary data, but when the generating law is elementary, it may recover the exact formula.

In recent decades, a growing number of discoveries in fields of mathematics have been assisted by computer algorithms, primarily for exploring large parameter spaces that humans would take too long to investigate. As computers and algorithms become more powerful, an intriguing possibility arises - the interplay between human intuition and computer algorithms can lead to discoveries of novel mathematical concepts that would otherwise remain elusive. To realize this perspective, we have developed a massively parallel computer algorithm that discovers an unprecedented number of continued fraction formulas for fundamental mathematical constants. The sheer number of formulas discovered by the algorithm unveils a novel mathematical structure that we call the conservative matrix field. Such matrix fields (1) unify thousands of existing formulas, (2) generate infinitely many new formulas, and most importantly, (3) lead to unexpected relations between different mathematical constants, including multiple integer values of the Riemann zeta function. Conservative matrix fields also enable new mathematical proofs of irrationality. In particular, we can use them to generalize the celebrated proof by Ap\'ery for the irrationality of zeta(3). Utilizing thousands of personal computers worldwide, our computer-supported research strategy demonstrates the power of experimental mathematics, highlighting the prospects of large-scale computational approaches to tackle longstanding open problems and discover unexpected connections across diverse fields of science.

According to the generalized Polya theorem, the Gaussian distribution on the real line is characterized by the property of equidistribution of a monomial and a linear form of independent identically distributed random variables. We give a complete description of a-adic solenoids for which an analog of this theorem is true. The proof of the main theorem is reduced to solving some functional equation in the class of continuous positive definite functions on the character group of an a-adic solenoid

We establish lower bounds on the number of resonances for three-dimensional semiclassical Dirac operators with mass-type potentials. For a pair Dl_j=Dl_0+mbfβv_j, j=1,2, satisfying the exterior analyticity, decay and non-degeneracy assumptions of the paper, we associate a relative distribution μ to the corresponding mass functions. We prove that if E_0>1 is a boundary point of suppμ, then every complex neighbourhood of E_0, and likewise of -E_0, contains in total at least Chbar^{-3} resonances of the two operators, counted with multiplicity. The proof transfers an analytic wave-front singularity of μ to the relative classical density of states and combines this with Khochman's local trace formula and a relative semiclassical trace formula. As a consequence, if v is a single bounded mass-type potential, analytic with power decay in an exterior sector, with inf(1+v)>0 and sup v>0, then Dl=Dl_0+mbfβv has at least Chbar^{-3} resonances in every complex neighbourhood of E_0=1+sup v and of -E_0, for all sufficiently small hbar. The exponent agrees with Khochman's upper bound, the two resonance families are related by spectral symmetry, and E_0 is necessarily a trapping energy. We also prove the relative semiclassical trace formula needed in the argument for general Hermitian matrix-valued perturbations.

Motivated by W. P. Thurston, we ask: What is the shape of a meromorphic function on a simply connected Riemann surface Ω_z? We consider Speiser functions, i.e. meromorphic functions on a simply connected Riemann surface, that have a finite number q at least 2 of singular (critical or asymptotic) values. As a first result, we make precise the correspondence between: Speiser functions w(z), Speiser Riemann surfaces R_w(z), Speiser q-tessellation, and analytic Speiser graphs of index q. As the second main result, we characterize tessellations with alternating colors (equivalently abstract pre-Speiser graphs) that are realized by Speiser functions on Ω_z. The characterization is in terms of the q-regular extension problem of bipartite planar graphs. As third main results, the Speiser Riemann surface R_w(z) can be constructed by isometric glueing of a finite number of types of sheets, where each sheet is a maximal domain of single-valuedness of the inverse of w(z). Furthermore, a unique decomposition of R_w(z) into maximal logarithmic towers and a soul is provided. Using vector fields we recognize that logarithmic towers come in two flavors: exponential or h-tangent blocks, directly related to the exponential or the hyperbolic tangent functions on the upper half plane. The surface R_w(z) of a finite Speiser function is characterized by surgery of a rational block and a finite number of exponential or h-tangent blocks.

The numerical solution of implicit and stiff differential equations by implicit numerical integrators has been largely investigated and there exist many excellent efficient codes available in the scientific community, as Radau5 (based on a Runge-Kutta collocation method at Radau points) and Dassl, based on backward differentiation formulas, among the others. When solving fractional ordinary differential equations (ODEs), the derivative operator is replaced by a non-local one and the fractional ODE is reformulated as a Volterra integral equation, to which these codes cannot be directly applied. This article is a follow-up of the article by the authors (Guglielmi and Hairer, SISC, 2025) for differential equations with distributed delays. The main idea is to approximate the fractional kernel t^{α-1}/ Γ(α) (α>0) by a sum of exponential functions or by a sum of exponential functions multiplied by a monomial, and then to transform the fractional integral (of convolution type) into a set of ordinary differential equations. The augmented system is typically stiff and thus requires the use of an implicit method. It can have a very large dimension and requires a special treatment of the arising linear systems. The present work presents an algorithm for the construction of an approximation of the fractional kernel by a sum of exponential functions, and it shows how the arising linear systems in a stiff time integrator can be solved efficiently. It is explained how the code Radau5 can be used for solving fractional differential equations. Numerical experiments illustrate the accuracy and the efficiency of the proposed method. Driver examples are publicly available from the homepages of the authors.

We advocate a strategy of bootstrapping Feynman integrals from just knowledge of their singular behavior. This approach is complementary to other bootstrap programs, which exploit non-perturbative constraints such as unitarity, or amplitude-level constraints such as gauge invariance. We begin by studying where a Feynman integral can become singular, and the behavior it exhibits near these singularities. We then characterize the space of functions that we expect the integral to evaluate to, in order to formulate an appropriate ansatz. Finally, we derive constraints on where each singularity can appear in this ansatz, and use information about the expansion of the integral around singular points in order to determine the value of all remaining free coefficients. Throughout, we highlight how constraints that have previously only been derived for integrals with generic masses can be extended to integrals involving particles of equal or vanishing mass. We illustrate the effectiveness of this approach by bootstrapping a number of examples, including the four-point double box with a massive internal loop.

We propose a penalized nonparametric approach to estimating the quantile regression process (QRP) in a nonseparable model using rectifier quadratic unit (ReQU) activated deep neural networks and introduce a novel penalty function to enforce non-crossing of quantile regression curves. We establish the non-asymptotic excess risk bounds for the estimated QRP and derive the mean integrated squared error for the estimated QRP under mild smoothness and regularity conditions. To establish these non-asymptotic risk and estimation error bounds, we also develop a new error bound for approximating C^s smooth functions with s >0 and their derivatives using ReQU activated neural networks. This is a new approximation result for ReQU networks and is of independent interest and may be useful in other problems. Our numerical experiments demonstrate that the proposed method is competitive with or outperforms two existing methods, including methods using reproducing kernels and random forests, for nonparametric quantile regression.

The analysis of the spectral features of a Toeplitz matrix-sequence left{T_{n}(f)right}_{ninmathbb N}, generated by a symbol fin L^1([-π,π]), real-valued almost everywhere (a.e.), has been provided in great detail in the last century, as well as the study of the conditioning, when f is nonnegative a.e. Here we consider a novel type of problem arising in the numerical approximation of distributed-order fractional differential equations (FDEs), where the matrices under consideration take the form \[ T_{n}=c_0T_{n}(f_0)+c_{1} h^h T_{n}(f_{1})+c_{2} h^{2h} T_{n}(f_{2})+\cdots+c_{n-1} h^{(n-1)h}T_{n}(f_{n-1}), \] c_0,c_{1},ldots, c_{n-1} in [c_*,c^*], c^*ge c_*>0, independent of n, h=1{n}, f_jsim g_j, g_j=|θ|^{2-jh}, j=0,ldots,n-1. Since the resulting generating function depends on n, the standard theory cannot be applied and the analysis has to be performed using new ideas. Few selected numerical experiments are presented, also in connection with matrices that come from distributed-order FDE problems, and the adherence with the theoretical analysis is discussed together with open questions and future investigations.

We review Bohr's atomic model and its extension by Sommerfeld from a mathematical perspective of wave mechanics. The derivation of quantization rules and energy levels is revisited using semiclassical methods. Sommerfeld-type integrals are evaluated by elementary techniques, and connections with the Schr\"{o}dinger and Dirac equations are established. Historical developments and key transitions from classical to quantum theory are discussed to clarify the structure and significance of the old quantum mechanics.

We prove a density theorem for the auxiliar function mathcal R(s) found by Siegel in Riemann papers. Let α be a real number with frac12< αle 1, and let N(α,T) be the number of zeros ρ=β+iγ of mathcal R(s) with 1ge βgeα and 0<γle T. Then we prove \[N(α,T)\ll T^{\frac32-α}(\log T)^3.\] Therefore, most of the zeros of mathcal R(s) are near the critical line or to the left of that line. The imaginary line for π^{-s/2}Γ(s/2)mathcal R(s) passing through a zero of mathcal R(s) near the critical line frequently will cut the critical line, producing two zeros of ζ(s) in the critical line.

Straight line equation y=mx with slope m, when singularly perturbed as ay^3+y=mx with a positive parameter a, results in S-shaped curves or S-curves on a real plane. As arightarrow 0, we get back y=mx which is a cumulative distribution function of a continuous uniform distribution that describes the occurrence of every event in an interval to be equally probable. As arightarrowinfty, the derivative of y has finite support only at y=0 resembling a degenerate distribution. Based on these arguments, in this work, we propose that these S-curves can represent maximum entropy uniform distribution to a zero entropy single value. We also argue that these S-curves are superposable as they are only parametrically nonlinear but fundamentally linear. So far, the superposed forms have been used to capture the patterns of natural systems such as nonlinear dynamics of biological growth and kinetics of enzyme reactions. Here, we attempt to use the S-curve and its superposed form as statistical models. We fit the models on a classical dataset containing flower measurements of iris plants and analyze their usefulness in pattern recognition. Based on these models, we claim that any non-uniform pattern can be represented as a singular perturbation to uniform distribution. However, our parametric estimation procedure have some limitations such as sensitivity to initial conditions depending on the data at hand.

This paper deals with a non-self-adjoint differential operator which is associated with a diffusion process with random jumps from the boundary. Our main result is that the algebraic multiplicity of an eigenvalue is equal to its order as a zero of the characteristic function Delta(lambda) . This can be used to determine the multiplicities of eigenvalues for concrete operators.

We provide conditional and unconditional asymptotic formulae for the exponential sums sum_γ,γ^{-iτ}, where the summation is over the ordinates of the nontrivial zeros ρ=β+iγ of the Riemann zeta-function. In particular, the obtained results are related to the Lindelöf Hypothesis for these ordinates (in the sense of Gonek et al. [10]).

Phenomenological modeling of variable stars allows determination of a set of the parameters, which are needed for classification in the "General Catalogue of Variable Stars" and similar catalogs. We apply a recent method NAV ("New Algol Variable") to eclipsing binary stars of different types. Although all periodic functions may be represented as Fourier series with an infinite number of coefficients, this is impossible for a finite number of the observations. Thus one may use a restricted Fourier series, i.e. a trigonometric polynomial (TP) of order s either for fitting the light curve, or to make a periodogram analysis. However, the number of parameters needed drastically increases with decreasing width of minimum. In the NAV algorithm, the special shape of minimum is used, so the number of parameters is limited to 10 (if the period and initial epoch are fixed) or 12 (not fixed). We illustrate the NAV method by application to a recently discovered Algol-type eclipsing variable 2MASS J11080308-6145589 (in the field of previously known variable star RS Car) and compare results to that obtained using the TP fits. For this system, the statistically optimal number of parameters is 44, but the fit is still worse than that of the NAV fit. Application to the system GSC 3692-00624 argues that the NAV fit is better than the TP one even for the case of EW-type stars with much wider eclipses. Model parameters are listed.

We performed the optical spectroscopy of 16 classical Be stars in 11 open clusters older than 100 Myr. Ours is the first spectroscopic study of classical Be stars in open clusters older than 100 Myr. We found that the H alpha emission strength of most of the stars is less than 40 Angstrom, in agreement with previous studies. Our analysis further suggests that one of the stars, KW97 35 12, might be a weak H alpha emitter in nature, showing H alpha equivalent width of negative 0.5 Angstrom. Interestingly, we also found that the newly detected classical Be star LS III 47 37b might be a component of the possible visual binary system LS III 47 37, where the other companion is also a classical Be star. Hence, the present study indicates the possible detection of a binary Be system. Moreover, it is observed that all 16 stars exhibit a lesser number of emission lines compared to classical Be stars younger than 100 Myr. Furthermore, the spectral type distribution analysis of B type and classical Be stars for the selected clusters points out that the existence of CBe stars can depend on the spectral type distribution of B type stars present in these clusters.

This thesis addresses a fundamental problem in deformation quantization: the difficulty of calculating the star-exponential, the symbol of the evolution operator, due to convergence issues. Inspired by the formalism that connects the star-exponential with the quantum propagator for bosonic systems, this work develops the analogous extension for the fermionic case. A rigorous method, based on Grassmann variables and coherent states, is constructed to obtain a closed-form expression for the fermionic star-exponential from its associated propagator. As a primary application, a fermionic version of the Feynman-Kac formula is derived within this formalism, allowing for the calculation of the ground state energy directly in phase space. Finally, the method is validated by successfully applying it to the simple and driven harmonic oscillators, where it is demonstrated that a simplified ("naive") approach (with an ad-hoc "remediation") is a valid weak-coupling limit of the rigorous ("meticulous") formalism, thereby providing a new and powerful computational tool for the study of fermionic systems.

This course, intended for undergraduates familiar with elementary calculus and linear algebra, introduces the extension of differential calculus to functions on more general vector spaces, such as functions that take as input a matrix and return a matrix inverse or factorization, derivatives of ODE solutions, and even stochastic derivatives of random functions. It emphasizes practical computational applications, such as large-scale optimization and machine learning, where derivatives must be re-imagined in order to be propagated through complicated calculations. The class also discusses efficiency concerns leading to "adjoint" or "reverse-mode" differentiation (a.k.a. "backpropagation"), and gives a gentle introduction to modern automatic differentiation (AD) techniques.

Bezigons, i.e., closed paths composed of Bézier curves, have been widely employed to describe shapes in image vectorization results. However, most existing vectorization techniques infer the bezigons by simply approximating an intermediate vector representation (such as polygons). Consequently, the resultant bezigons are sometimes imperfect due to accumulated errors, fitting ambiguities, and a lack of curve priors, especially for low-resolution images. In this paper, we describe a novel method for vectorizing clipart images. In contrast to previous methods, we directly optimize the bezigons rather than using other intermediate representations; therefore, the resultant bezigons are not only of higher fidelity compared with the original raster image but also more reasonable because they were traced by a proficient expert. To enable such optimization, we have overcome several challenges and have devised a differentiable data energy as well as several curve-based prior terms. To improve the efficiency of the optimization, we also take advantage of the local control property of bezigons and adopt an overlapped piecewise optimization strategy. The experimental results show that our method outperforms both the current state-of-the-art method and commonly used commercial software in terms of bezigon quality.

We investigate singularities of radial solutions of the classical two-phase and one-phase Stefan problem. In a remarkable result, M.A. Herrero and J.J.L. Velazquez show existence of a solution such that the free boundary described by g(t) satisfies in dimension ngeq 3 $g(t)=c_n t |log t|^{-1{n - 2}} (1 + o(1)), tto 0,$ and they conjecture that these asymptotics hold in general. Here we prove the conjecture for all radial solutions, which means complete asymptotic rigidity in that class. In order to obtain these precise asymptotics for the two-phase Stefan problem, for which frequency formulas and monotonicity formulas are unknown, we apply the Potential Reduction to ODE method introduced in arXiv:2608.18913 for the Neumann Bernoulli problem. The philosophy of that method is that the asymptotic shape of the free surface is governed by a precise ordinary differential equation. The fact that we can apply it in this paper to an equation of different type with quasilinear principal part, which is in our regime not the perturbation of a linear operator, confirms the versatility of the method.

Solutions of nonlinear functional equations are generally not expressed as a finite number of combinations and compositions of elementary and known special functions. One of the approaches to study them is, firstly, to find formal solutions (that is, series whose terms are described and ordered in some way but which do not converge apriori) and, secondly, to study the convergence or summability of these formal solutions (the existence and uniqueness of actual solutions with the given asymptotic expansion in a certain domain). In this paper we deal only with the convergence of formal functional series having the form of an infinite sum of power functions with (complex, in general) power exponents and satisfying analytical functional equations of the following three types: a differential, q-difference or Mahler equation.

Ap\'{e}ry-type (inverse) binomial series have appeared prominently in the calculations of Feynman integrals in recent years. In our previous work, we showed that a few large classes of the non-alternating Ap\'ery-type (inverse) central binomial series can be evaluated using colored multiple zeta values of level four (i.e., special values of multiple polylogarithms at fourth roots of unity) by expressing them in terms of iterated integrals. In this sequel, we shall prove that for several classes of the alternating versions we need to raise the level to eight. Our main idea is to adopt hyperbolic trigonometric 1-forms to replace the ordinary trigonometric ones used in the non-alternating setting.

We study the parity-perturbed Hofstadter Q-recursion $ widetilde Q(1)=widetilde Q(2)=1,qquad widetilde Q(n)=widetilde Q(n-widetilde Q(n-1)) +widetilde Q(n-widetilde Q(n-2))+(-1)^n, and the associated two-variable Dirichlet series Z_{widetilde Q}(s,t)=sum_{nge1}n^{-s}widetilde Q(n)^{-t}. The estimate \widetilde Q(n)=n/2+O(n/\log n) gives the exact domain of absolute convergence Re(s+t)>1. With w=s+t, we separate the universal term 2^tζ(w) and derive exact transport, frequency-position, and dyadic renormalization identities. The main result concerns t=-1. For E(n)=2\widetilde Q(n)-n and A(X)=\sum_{n\le X}E(n), the binary-arch clock yields A(X)=Xlog_2X+XΩ!left(log_23X{32}right) +O!left(X{log X}right), where Ω is an explicit continuous periodic function. This continues the normalized correction to Rew>0 and yields a boundary resonance lattice: a double resonance at w=0 and simple resonances at 2πi m/\log2. After subtracting the full-slice order-X skeleton, we analyze the negative-even arch channel. Its companion-forest layers have a weak Gaussian limit, and a canonical subsequence realizes the optimal n/\log n pointwise scale with an explicit signed constant. The negative-arch mass satisfies A_r=512{92π}16^r{sqrt r} left(1-13{16r}+O(r^{-2})right). We do not claim a full-slice continuation across Rew=0$.

To model high-risk traffic densities a cellular automaton model is constructed, exhibiting Ising model-like properties. The attraction-repulsion forces between vehicles are evaluated as coupling, and the coupling function p(1-p)=g/8, where p is the initial distribution of cells with state 1, and g is the Moore neighbor count, which has roots cos^2(π/8) and sin^2(π/8) for g=1. This is achieved without the use of any trigonometric functions in the code. From the roots, the tangent polynomial tan^2(x) + tan(x) emerges. Kramers-Wannier duality is recovered and it is conjectured that the 1/2 difference between the roots serves as a fixed point for a projection mechanism from the 2-dimensional Ising model onto a 1-dimensional Ising chain through the Gudermannian function. The sigmoid evaluated at the proposed fixed point is substituted into the derivative of the logistic coupling function, σ(1/2)(1-σ(1/2)), yielding a numerical approximation to the three-dimensional Ising inverse critical coupling. Finally, the results are linked to risk densities in traffic and vehicle types, accounting for the amplification of fatal accidents.

We made a comparative analysis of numerical methods for multidimensional optimization. The main parameter is a number of computations of the test function to reach necessary accuracy, as it is computationally "slow". For complex functions, analytic differentiation by many parameters can cause problems associated with a significant complication of the program and thus slowing its operation. For comparison, we used the methods: "brute force" (or minimization on a regular grid), Monte Carlo, steepest descent, conjugate gradients, Brent's method (golden section search), parabolic interpolation etc. The Monte-Carlo method was applied to the eclipsing binary system AM Leo.