Search

An algorithm for robust initial orbit determination (IOD) under perturbed orbital dynamics is presented. By leveraging map inversion techniques defined in the algebra of Taylor polynomials, this tool returns a highly accurate solution to the IOD problem and estimates a range centered on the aforementioned solution in which the true orbit should lie. To meet the specified accuracy requirements, automatic domain splitting is used to wrap the IOD routines and ensure that the local truncation error, introduced by a polynomial representation of the state estimate, remains below a predefined threshold. The algorithm is presented for three types of ground-based sensors, namely range radars, Doppler-only radars, and optical telescopes, by considering their different constraints in terms of available measurements and sensor noise. Finally, the improvement in performance with respect to a Keplerian-based IOD solution is demonstrated using large-scale numerical simulations over a subset of tracked objects in low Earth orbit., Comment: accepted for publication on Astrodynamics

The relativistic hydrodynamics (RHD) equations have three crucial intrinsic physical constraints on the primitive variables: positivity of pressure and density, and subluminal fluid velocity. However, numerical simulations can violate these constraints, leading to nonphysical results or even simulation failure. Designing genuinely physical-constraint-preserving (PCP) schemes is very difficult, as the primitive variables cannot be explicitly reformulated using conservative variables due to relativistic effects. In this paper, we propose three efficient Newton--Raphson (NR) methods for robustly recovering primitive variables from conservative variables. Importantly, we rigorously prove that these NR methods are always convergent and PCP, meaning they preserve the physical constraints throughout the NR iterations. The discovery of these robust NR methods and their PCP convergence analyses are highly nontrivial and technical. As an application, we apply the proposed NR methods to design PCP finite volume Hermite weighted essentially non-oscillatory (HWENO) schemes for solving the RHD equations. Our PCP HWENO schemes incorporate high-order HWENO reconstruction, a PCP limiter, and strong-stability-preserving time discretization. We rigorously prove the PCP property of the fully discrete schemes using convex decomposition techniques. Moreover, we suggest the characteristic decomposition with rescaled eigenvectors and scale-invariant nonlinear weights to enhance the performance of the HWENO schemes in simulating large-scale RHD problems. Several demanding numerical tests are conducted to demonstrate the robustness, accuracy, and high resolution of the proposed PCP HWENO schemes and to validate the efficiency of our NR methods., Comment: 49 pages

This paper proposes and analyzes a novel efficient high-order finite volume method for the ideal magnetohydrodynamics (MHD). As a distinctive feature, the method simultaneously preserves a discretely divergence-free (DDF) constraint on the magnetic field and the positivity-preserving (PP) property, which ensures the positivity of density, pressure, and internal energy. To enforce the DDF condition, we design a new discrete projection approach that projects the reconstructed point values at the cell interface into a DDF space, without using any approximation polynomials. This projection method is highly efficient, easy to implement, and particularly suitable for standard high-order finite volume WENO methods, which typically return only the point values in the reconstruction. Moreover, we also develop a new finite volume framework for constructing provably PP schemes for the ideal MHD system. The framework comprises the discrete projection technique, a suitable approximation to the Godunov--Powell source terms, and a simple PP limiter. We provide rigorous analysis of the PP property of the proposed finite volume method, demonstrating that the DDF condition and the proper approximation to the source terms eliminate the impact of magnetic divergence terms on the PP property. The analysis is challenging due to the internal energy function's nonlinearity and the intricate relationship between the DDF and PP properties. To address these challenges, the recently developed geometric quasilinearization approach is adopted, which transforms a nonlinear constraint into a family of linear constraints. Finally, we validate the effectiveness of the proposed method through several benchmark and demanding numerical examples. The results demonstrate that the proposed method is robust, accurate, and highly effective, confirming the significance of the proposed DDF projection and PP techniques., Comment: 26 pages

This paper deals with the problem of constructing a flight scheme to Venus, in which a spacecraft flying to the planet after a gravity assist maneuver and transition to a resonant orbit in order to re-encounter with Venus, makes a passage of a minor celestial body. The 117 candidate asteroids from the NASA JPL catalogue, whose diameter exceeds 1 km, were selected. The flight trajectories which meet the criteria of impulse-free both flyby Venus and asteroid, and the subsequent landing on the surface of Venus were found within the interval of launch dates from 2029 to 2050. The trajectory of the spacecraft flight from the Earth to Venus including flyby of Venus and asteroids with a subsequent landing on the surface of Venus was analyzed.

We present a new divergence-free and well-balanced hybrid FV/FE scheme for the incompressible viscous and resistive MHD equations on unstructured mixed-element meshes in 2 and 3 space dimensions. The equations are split into subsystems. The pressure is defined on the vertices of the primary mesh, while the velocity field and the normal components of the magnetic field are defined on an edge-based/face-based dual mesh in two and three space dimensions, respectively. This allows to account for the divergence-free conditions of the velocity field and of the magnetic field in a rather natural manner. The non-linear convective and the viscous terms are solved at the aid of an explicit FV scheme, while the magnetic field is evolved in a divergence-free manner via an explicit FV method based on a discrete form of the Stokes law in the edges/faces of each primary element. To achieve higher order of accuracy, a pw-linear polynomial is reconstructed for the magnetic field, which is guaranteed to be divergence-free via a constrained L2 projection. The pressure subsystem is solved implicitly at the aid of a classical continuous FE method in the vertices of the primary mesh. In order to maintain non-trivial stationary equilibrium solutions of the governing PDE system exactly, which are assumed to be known a priori, each step of the new algorithm takes the known equilibrium solution explicitly into account so that the method becomes exactly well-balanced. This paper includes a very thorough study of the lid-driven MHD cavity problem in the presence of different magnetic fields. We finally present long-time simulations of Soloviev equilibrium solutions in several simplified 3D tokamak configurations even on very coarse unstructured meshes that, in general, do not need to be aligned with the magnetic field lines., Comment: 57 pages, 33 figures, 13 tables, reference-data (supplementary electronic material) will be available after publication on the Journal web-page

Wavelets are widely used in various disciplines to analyse signals both in space and scale. Whilst many fields measure data on manifolds (i.e., the sphere), often data are only observed on a partial region of the manifold. Wavelets are a typical approach to data of this form, but the wavelet coefficients that overlap with the boundary become contaminated and must be removed for accurate analysis. Another approach is to estimate the region of missing data and to use existing whole-manifold methods for analysis. However, both approaches introduce uncertainty into any analysis. Slepian wavelets enable one to work directly with only the data present, thus avoiding the problems discussed above. Applications of Slepian wavelets to areas of research measuring data on the partial sphere include gravitational/magnetic fields in geodesy, ground-based measurements in astronomy, measurements of whole-planet properties in planetary science, geomagnetism of the Earth, and cosmic microwave background analyses., Comment: 4 pages

Multi-scale deconvolution is an ill-posed inverse problem in imaging, with applications ranging from microscopy, through medical imaging, to astronomical remote sensing. In the case of high-energy space telescopes, multi-scale deconvolution algorithms need to account for the peculiar property of native measurements, which are sparse samples of the Fourier transform of the incoming radiation. The present paper proposes a multi-scale version of CLEAN, which is the most popular iterative deconvolution method in Fourier space imaging. Using synthetic data generated according to a simulated but realistic source configuration, we show that this multi-scale version of CLEAN performs better than the original one in terms of accuracy, photometry, and regularization. Further, the application to a data set measured by the NASA Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) shows the ability of multi-scale CLEAN to reconstruct rather complex topographies, characteristic of a real flaring event.

The continuously growing number of objects orbiting around the Earth is expected to be accompanied by an increasing frequency of objects re-entering the Earth's atmosphere. Many of these re-entries will be uncontrolled, making their prediction challenging and subject to several uncertainties. Traditionally, re-entry predictions are based on the propagation of the object's dynamics using state-of-the-art modelling techniques for the forces acting on the object. However, modelling errors, particularly related to the prediction of atmospheric drag may result in poor prediction accuracies. In this context, we explore the possibility to perform a paradigm shift, from a physics-based approach to a data-driven approach. To this aim, we present the development of a deep learning model for the re-entry prediction of uncontrolled objects in Low Earth Orbit (LEO). The model is based on a modified version of the Sequence-to-Sequence architecture and is trained on the average altitude profile as derived from a set of Two-Line Element (TLE) data of over 400 bodies. The novelty of the work consists in introducing in the deep learning model, alongside the average altitude, three new input features: a drag-like coefficient (B*), the average solar index, and the area-to-mass ratio of the object. The developed model is tested on a set of objects studied in the Inter-Agency Space Debris Coordination Committee (IADC) campaigns. The results show that the best performances are obtained on bodies characterised by the same drag-like coefficient and eccentricity distribution as the training set.

This work outlines a fast, high-precision time-domain solver for scalar, electromagnetic and gravitational perturbations on hyperboloidal foliations of Kerr space-times. Time-domain Teukolsky equation solvers have typically used explicit methods, which numerically violate Noether symmetries and are Courant-limited. These restrictions can limit the performance of explicit schemes when simulating long-time extreme mass ratio inspirals, expected to appear in LISA band for 2-5 years. We thus explore symmetric (exponential, Pad\'e or Hermite) integrators, which are unconditionally stable and known to preserve certain Noether symmetries and phase-space volume. For linear hyperbolic equations, these implicit integrators can be cast in explicit form, making them well-suited for long-time evolution of black hole perturbations. The 1+1 modal Teukolsky equation is discretized in space using polynomial collocation methods and reduced to a linear system of ordinary differential equations, coupled via mode-coupling arrays and discretized (matrix) differential operators. We use a matricization technique to cast the mode-coupled system in a form amenable to a method-of-lines framework, which simplifies numerical implementation and enables efficient parallelization on CPU and GPU architectures. We test our numerical code by studying late-time tails of Kerr spacetime perturbations in the sub-extremal and extremal cases., Comment: 28 pages, 9 figures

Inspired by recent interest in geometric deep learning, this work generalises the recently developed Slepian scale-discretised wavelets on the sphere to Riemannian manifolds. Through the sifting convolution, one may define translations and, thus, convolutions on manifolds - which are otherwise not well-defined in general. Slepian wavelets are constructed on a region of a manifold and are therefore suited to problems where data only exists in a particular region. The Slepian functions, on which Slepian wavelets are built, are the basis functions of the Slepian spatial-spectral concentration problem on the manifold. A tiling of the Slepian harmonic line with smoothly decreasing generating functions defines the scale-discretised wavelets; allowing one to probe spatially localised, scale-dependent features of a signal. By discretising manifolds as graphs, the Slepian functions and wavelets of a triangular mesh are presented. Through a wavelet transform, the wavelet coefficients of a field defined on the mesh are found and used in a straightforward thresholding denoising scheme., Comment: 12 pages, 12 figures

Global seismicity on all three solar system's bodies with in situ measurements (Earth, Moon, and Mars) is due mainly to mechanical Rieger resonance (RR) of the solar wind's macroscopic flapping, driven by the well-known PRg=~154-day Rieger period and detected commonly in most heliophysical data types and the interplanetary magnetic field (IMF). Thus, InSight mission marsquakes rates are periodic with PRg as characterized by a very high (>>12) fidelity {\Phi}=2.8 10^6 and by being the only >99%-significant spectral peak in the 385.8-64.3-nHz (1-180-day) band of highest planetary energies; the longest-span (v.9) release of raw data revealed the entire RR, excluding a tectonically active Mars. For check, I analyze rates of Oct 2015-Feb 2019, Mw5.6+ earthquakes, and all (1969-1977) Apollo mission moonquakes. To decouple magnetosphere and IMF effects, I study Earth and Moon seismicity during traversals of Earth magnetotail vs. IMF. The analysis showed with >99-67% confidence and {\Phi}>>12 fidelity that (an unspecified majority of) moonquakes and Mw5.6+ earthquakes also recur at Rieger periods. About half of spectral peaks split but also into clusters that average to usual Rieger periodicities, where magnetotail reconnecting clears the signal. Moonquakes are mostly forced at times of solar-wind resonance and not just during tides, as previously and simplistically believed. Earlier claims that solar plasma dynamics could be seismogenic are confirmed. This result calls for reinterpreting the seismicity phenomenon and for reliance on global magnitude scales. Predictability of solar-wind macroscopic dynamics is now within reach, which paves the way for long-term physics-based seismic and space weather prediction and the safety of space missions. Gauss-Vanicek Spectral Analysis revolutionizes geophysics by computing nonlinear global dynamics directly (renders approximating of dynamics obsolete)., Comment: 46 printed pages, 8 figures, 5 tables; includes time-series of all Apollo natural moonquakes (1969-1977), all Mw5.6+ natural earthquakes between Oct 2015-Feb 2019, and all v.9 InSight natural marsquakes (Jan 2019-Sep 2021); published version

We present an open-source CUDA-based package that consists of a compilation of exponential integrators where the action of the matrix exponential or the $\varphi_l$ functions on a vector is approximated using the method of polynomial interpolation at Leja points. Using a couple of test examples on an NVIDIA A100 GPU, we show that one can achieve significant speedups using CUDA over the corresponding CPU code. LeXInt, written in a modular format, facilitates easy integration into any existing software package, and can be used for temporal integration of any differential equation., Comment: Code available at https://github.com/Pranab-JD/LeXInt; comments welcome

As astronomical data grows in volume and complexity, the scalability of analysis software becomes increasingly important. At the same time, astrophysics analysis software relies heavily on open-source contributions, so languages and tools that prioritize both performance and readability are especially valuable. Julia, with its just-in-time compiler and high level syntax, offers a compelling alternative to traditional languages like Python or C. In this paper, we outline ShOpt.jl, a new software package for point spread function (PSF) characterization written in Julia. ShOpt.jl features a number of performance optimizations, such as multithreading, the use of preconditioners, and the implementation of the memory-limited Broyden-Fletcher-Goldfarb-Shanno algorithm, as well as the flexibility to choose between principal component analysis, an autoencoder, and analytic profiles for PSF characterization. As observatories like the James Webb Space Telescope bring astrophysics into a new era of wide-field, high-resolution imaging, the challenges of PSF modeling become more pronounced. Tools like ShOpt.jl provide the community with a scalable, efficient, and accurate solution to these challenges, while also demonstrating the potential of Julia as a language that meets the demands of modern astrophysical research., Comment: 7 pages, 2 figures; submitted to the Journal of Open Source Software

In this work, we explore how to classify asteroids in co-orbital motion with a given planet using Machine Learning. We consider four different kinds of motion in mean motion resonance with the planet, nominally Tadpole, Horseshoe and Quasi-satellite, building 3 datasets defined as Real (taking the ephemerides of real asteroids from the JPL Horizons system), Ideal and Perturbed (both simulated, obtained by propagating initial conditions considering two different dynamical systems) for training and testing the Machine Learning algorithms in different conditions. The time series of the variable theta (angle related to the resonance) are studied with a data analysis pipeline defined ad hoc for the problem and composed by: data creation and annotation, time series features extraction thanks to the tsfresh package (potentially followed by selection and standardization) and the application of Machine Learning algorithms for Dimensionality Reduction and Classification. Such approach, based on features extracted from the time series, allows to work with a smaller number of data with respect to Deep Learning algorithms, also allowing to define a ranking of the importance of the features. Physical Interpretability of the features is another key point of this approach. In addition, we introduce the SHapley Additive exPlanations for Explainability technique. Different training and test sets are used, in order to understand the power and the limits of our approach. The results show how the algorithms are able to identify and classify correctly the time series, with a high degree of performance.

We consider particle systems described by moments of a phase-space density and propose a realizability-preserving numerical method to evolve a spectral two-moment model for particles interacting with a background fluid moving with nonrelativistic velocities. The system of nonlinear moment equations, with special relativistic corrections to $\mathcal{O}(v/c)$, expresses a balance between phase-space advection and collisions and includes velocity-dependent terms that account for spatial advection, Doppler shift, and angular aberration. This model is closely related to the one promoted by Lowrie et al. (2001; JQSRT, 69, 291-304) and similar to models currently used to study transport phenomena in large-scale simulations of astrophysical environments. The method is designed to preserve moment realizability, which guarantees that the moments correspond to a nonnegative phase-space density. The realizability-preserving scheme consists of the following key components: (i) a strong stability-preserving implicit-explicit (IMEX) time-integration method; (ii) a discontinuous Galerkin (DG) phase-space discretization with carefully constructed numerical fluxes; (iii) a realizability-preserving implicit collision update; and (iv) a realizability-enforcing limiter. In time integration, nonlinearity of the moment model necessitates solution of nonlinear equations, which we formulate as fixed-point problems and solve with tailored iterative solvers that preserve moment realizability with guaranteed convergence. We also analyze the simultaneous Eulerian-frame number and energy conservation properties of the semi-discrete DG scheme and propose an "energy limiter" that promotes Eulerian-frame energy conservation. Through numerical experiments, we demonstrate the accuracy and robustness of this DG-IMEX method and investigate its Eulerian-frame energy conservation properties.

This work presents a large-scale simulation study investigating the deployment and operation of distributed swarms of CubeSats for interplanetary missions to small celestial bodies. Utilizing Taylor numerical integration and advanced collision detection techniques, we explore the potential of large CubeSat swarms in capturing gravity signals and reconstructing the internal mass distribution of a small celestial body while minimizing risks and Delta V budget. Our results offer insight into the applicability of this approach for future deep space exploration missions.

This work outlines a time-domain numerical integration technique for linear hyperbolic partial differential equations sourced by distributions (Dirac $\delta$-functions and their derivatives). Such problems arise when studying binary black hole systems in the extreme mass ratio limit. We demonstrate that such source terms may be converted to effective domain-wide sources when discretized, and we introduce a class of time-steppers that directly account for these discontinuities in time integration. Moreover, our time-steppers are constructed to respect time reversal symmetry, a property that has been connected to conservation of physical quantities like energy and momentum in numerical simulations. To illustrate the utility of our method, we numerically study a distributionally-sourced wave equation that shares many features with the equations governing linear perturbations to black holes sourced by a point mass., Comment: 29 pages, 4 figures

This study explores the implications of Dark Matter in Neutron Stars (DMANS) by focusing on two specific astronomical objects: HESS J1731-347 and PSR J0952-0607. Varying the Fermi momentum k$f^{\rm DM}$ of DM, the study analyzes the EOS for the INRS model with and without DM. Results show the robustness of the model, with most EOS curves within chiral Effective Field Theory bounds. Our model predicts a maximum mass of $2.343 \ M_\odot$ for PSR J0952-0607, satisfying NICER bounds. The analysis suggests HESS J1731-347 could be a DMANS. Constraints on DM within NSs are established, and tidal deformability lies within GW event limits. Nonradial $f$-mode oscillations increase with DM, concluding low mass stars pulsate at higher frequencies.

Anisotropic diffusion is imperative in understanding cosmic ray diffusion across the Galaxy, the heliosphere, and the interplay of cosmic rays with the Galactic magnetic field. This diffusion term contributes to the highly stiff nature of the cosmic ray transport equation. To conduct numerical simulations of time-dependent cosmic ray transport, implicit integrators (namely, Crank-Nicolson (CN)) have been traditionally favoured over the CFL-bound explicit integrators in order to be able to take large step sizes. We propose exponential methods to treat the linear anisotropc diffusion equation in the presence of advection and time-independent and time-dependent sources. These methods allow us to take even larger step sizes that can substantially speed-up the simulations whilst generating highly accurate solutions. In or subsequent work, we will use these exponential solvers in the Picard code to study anisotropic cosmic ray diffusion and we will consider additional physical processes such as continuous momentum losses and reacceleration., Comment: The 38th International Cosmic Ray Conference (ICRC2023)

Aims: To formulate, implement, and validate a user-independent release of CLEAN for Fourier-based image reconstruction of hard X-rays flaring sources. Methods: CLEAN is an iterative deconvolution method for radio and hard X-ray solar imaging. In a specific step of its pipeline, CLEAN requires the convolution between an idealized version of the instrumental Point Spread Function (PSF), and a map collecting point sources located at positions on the solar disk from where most of the flaring radiation is emitted. This convolution step has highly heuristic motivations and the shape of the idealized PSF, which depends on the user's choice, impacts the shape of the overall reconstruction. Here we propose the use of an interpolation/extrapolation process to avoid this user-dependent step, and to realize a completely unbiased version of CLEAN. Results: Applications to observations recorded by the Spectrometer/Telescope for Imaging X-rays (STIX) on-board Solar Orbiter show that this unbiased release of CLEAN outperforms the standard version of the algorithm in terms of both automation and reconstruction reliability, with reconstructions whose accuracy is in line with the one offered by other imaging methods developed in the STIX framework. Conclusions: This unbiased version of CLEAN proposes a feasible solution to a well-known open issue concerning CLEAN, i.e., its low degree of automation. Further, this study provided the first application of an interpolation/extrapolation approach to image reconstruction from STIX experimental visibilities.

Estimating the trace of the inverse of a large matrix is an important problem in lattice quantum chromodynamics. A multilevel Monte Carlo method is proposed for this problem that uses different degree polynomials for the levels. The polynomials are developed from the GMRES algorithm for solving linear equations. To reduce orthogonalization expense, the highest degree polynomial is a composite or double polynomial found with a polynomial preconditioned GMRES iteration. Added to some of the Monte Carlo pieces is deflation of eigenvalues that reduces the variance. Deflation is also used for finding a reduced degree deflated polynomial. The new Multipolynomial Monte Carlo method can significantly improve the trace computation for matrices that have a difficult spectrum due to small eigenvalues.

Taking insight from the theory of general relativity, where space and time are treated on the same footing, we develop a novel geometric variational discretization for second order initial value problems (IVPs). By discretizing the dynamics along a world-line parameter, instead of physical time directly, we retain manifest translation symmetry and conservation of the associated continuum Noether charge. A non-equidistant time discretization emerges dynamically, realizing a form of automatic adaptive mesh refinement (AMR), guided by the system symmetries. Using appropriately regularized summation by parts finite difference operators, the continuum Noether charge, defined via the Killing vector associated with translation symmetry, is shown to be exactly preserved in the interior of the simulated time interval. The convergence properties of the approach are demonstrated with two explicit examples., Comment: 35 pages, 12 figures

We demonstrate that gauge-equivariant pooling and unpooling layers can perform as well as traditional restriction and prolongation layers in multigrid preconditioner models for lattice QCD. These layers introduce a gauge degree of freedom on the coarse grid, allowing for the use of explicitly gauge-equivariant layers on the coarse grid. We investigate the construction of coarse-grid gauge fields and study their efficiency in the preconditioner model. We show that a combined multigrid neural network using a Galerkin construction for the coarse-grid gauge field eliminates critical slowing down., Comment: 8 pages, 5 figures

We present an algorithm to approximate the real trilogarithm for a real argument with IEEE 754-1985 double precision accuracy. The approximation is structured such that it can make use of instruction-level parallelism when executed on appropriate CPUs., Comment: 4 pages, 3 tables, attached source code

Using the blockwise matrix inversion, inversions of large matrices with different ways of memory handling are presented in this article. Algorithm for performing inversion of matrix which is partitioned into large number of blocks is presented in which inversions and multiplications involving the blocks can be carried out with parallel processing. Optimized memory handling and efficient methods for intermediate multiplications among the partitioned blocks are implemented in this algorithm., Comment: 22 pages, 7 figures

The identification of interesting substructures within jets is an important tool for searching for new physics and probing the Standard Model at colliders. Many of these substructure tools have previously been shown to take the form of optimal transport problems, in particular the Energy Mover's Distance (EMD). In this work, we show that the EMD is in fact the natural structure for comparing collider events, which accounts for its recent success in understanding event and jet substructure. We then present a Shape Hunting Algorithm using Parameterized Energy Reconstruction (SHAPER), which is a general framework for defining and computing shape-based observables. SHAPER generalizes N-jettiness from point clusters to any extended, parametrizable shape. This is accomplished by efficiently minimizing the EMD between events and parameterized manifolds of energy flows representing idealized shapes, implemented using the dual-potential Sinkhorn approximation of the Wasserstein metric. We show how the geometric language of observables as manifolds can be used to define novel observables with built-in infrared-and-collinear safety. We demonstrate the efficacy of the SHAPER framework by performing empirical jet substructure studies using several examples of new shape-based observables., Comment: 45+15 pages, 18 figures, 5 tables. Code available at https://github.com/rikab/SHAPER; v2: Minor fixes; v3: Updated to match JHEP version; v4: Minor formatting fixes

Using the blockwise matrix inversion, inversions of large matrices with different ways of memory handling are presented in this article. Algorithm for performing inversion of matrix which is partitioned into large number of blocks is presented in which inversions and multiplications involving the blocks can be carried out with parallel processing. Optimized memory handling and efficient methods for intermediate multiplications among the partitioned blocks are implemented in this algorithm., Comment: 22 pages, 7 figures

The identification of interesting substructures within jets is an important tool for searching for new physics and probing the Standard Model at colliders. Many of these substructure tools have previously been shown to take the form of optimal transport problems, in particular the Energy Mover's Distance (EMD). In this work, we show that the EMD is in fact the natural structure for comparing collider events, which accounts for its recent success in understanding event and jet substructure. We then present a Shape Hunting Algorithm using Parameterized Energy Reconstruction (SHAPER), which is a general framework for defining and computing shape-based observables. SHAPER generalizes N-jettiness from point clusters to any extended, parametrizable shape. This is accomplished by efficiently minimizing the EMD between events and parameterized manifolds of energy flows representing idealized shapes, implemented using the dual-potential Sinkhorn approximation of the Wasserstein metric. We show how the geometric language of observables as manifolds can be used to define novel observables with built-in infrared-and-collinear safety. We demonstrate the efficacy of the SHAPER framework by performing empirical jet substructure studies using several examples of new shape-based observables., Comment: 45+15 pages, 18 figures, 5 tables. Code available at https://github.com/rikab/SHAPER; v2: Minor fixes; v3: Updated to match JHEP version; v4: Minor formatting fixes

Trajectory optimization of a robot manipulator consists of both optimization of the robot movement as well as optimization of the robot end-effector path. This paper aims to find optimum movement parameters including movement type, speed, and acceleration to minimize robot energy. Trajectory optimization by minimizing the energy would increase the longevity of robotic manipulators. We utilized the particle swarm optimization method to find the movement parameters leading to minimum energy consumption. The effectiveness of the proposed method is demonstrated on different trajectories. Experimental results show that 49% efficiency was obtained using a UR5 robotic arm.

Identifying parameters of computational models from experimental data, or model calibration, is fundamental for assessing and improving the predictability and reliability of computer simulations. In this work, we propose a method for Bayesian calibration of models that predict morphological patterns of diblock copolymer (Di-BCP) thin film self-assembly while accounting for various sources of uncertainties in pattern formation and data acquisition. This method extracts the azimuthally-averaged power spectrum (AAPS) of the top-down microscopy characterization of Di-BCP thin film patterns as summary statistics for Bayesian inference of model parameters via the pseudo-marginal method. We derive the analytical and approximate form of a conditional likelihood for the AAPS of image data. We demonstrate that AAPS-based image data reduction retains the mutual information, particularly on important length scales, between image data and model parameters while being relatively agnostic to the aleatoric uncertainties associated with the random long-range disorder of Di-BCP patterns. Additionally, we propose a phase-informed prior distribution for Bayesian model calibration. Furthermore, reducing image data to AAPS enables us to efficiently build surrogate models to accelerate the proposed Bayesian model calibration procedure. We present the formulation and training of two multi-layer perceptrons for approximating the parameter-to-spectrum map, which enables fast integrated likelihood evaluations. We validate the proposed Bayesian model calibration method through numerical examples, for which the neural network surrogate delivers a fivefold reduction of the number of model simulations performed for a single calibration task.

The objective of this study is to address the difficulty of simplifying the geometric model in which a differential problem is formulated, also called defeaturing, while simultaneously ensuring that the accuracy of the solution is maintained under control. This enables faster and more efficient simulations, without sacrificing accuracy. More precisely, we consider an isogeometric discretisation of an elliptic model problem defined on a two-dimensional hierarchical B-spline computational domain with a complex boundary. Starting with an oversimplification of the geometry, we build a goal-oriented adaptive strategy that adaptively reintroduces continuous geometrical features in regions where the analysis suggests a large impact on the quantity of interest. This strategy is driven by an a posteriori estimator of the defeaturing error based on first-order shape sensitivity analysis, and it profits from the local refinement properties of hierarchical B-splines. The adaptive algorithm is described together with a procedure to generate (partially) simplified hierarchical B-spline geometrical domains. Numerical experiments are presented to illustrate the proposed strategy and its limitations.

Cardiac digital twins provide a physics and physiology informed framework to deliver predictive and personalized medicine. However, high-fidelity multi-scale cardiac models remain a barrier to adoption due to their extensive computational costs and the high number of model evaluations needed for patient-specific personalization. Artificial Intelligence-based methods can make the creation of fast and accurate whole-heart digital twins feasible. In this work, we use Latent Neural Ordinary Differential Equations (LNODEs) to learn the temporal pressure-volume dynamics of a heart failure patient. Our surrogate model based on LNODEs is trained from 400 3D-0D whole-heart closed-loop electromechanical simulations while accounting for 43 model parameters, describing single cell through to whole organ and cardiovascular hemodynamics. The trained LNODEs provides a compact and efficient representation of the 3D-0D model in a latent space by means of a feedforward fully-connected Artificial Neural Network that retains 3 hidden layers with 13 neurons per layer and allows for 300x real-time numerical simulations of the cardiac function on a single processor of a standard laptop. This surrogate model is employed to perform global sensitivity analysis and robust parameter estimation with uncertainty quantification in 3 hours of computations, still on a single processor. We match pressure and volume time traces unseen by the LNODEs during the training phase and we calibrate 4 to 11 model parameters while also providing their posterior distribution. This paper introduces the most advanced surrogate model of cardiac function available in the literature and opens new important venues for parameter calibration in cardiac digital twins.

A novel approach for the stabilization of the Discontinuous Galerkin method based on the Dafermos entropy rate crition is presented. First, estimates for the maximal possible entropy dissipation rate of a weak solution are derived. Second, families of conservative Hilbert-Schmidt operators are identified to dissipate entropy. Steering these operators using the bounds on the entropy dissipation results in high-order accurate shock-capturing DG schemes for the Euler equations, satisfying the entropy rate criterion and an entropy inequality.

We present and analyse a hybridized discontinuous Galerkin method for incompressible flow problems using non-affine cells, proving that it preserves a key invariance property that illudes most methods, namely that any irrotational component of the prescribed force is exactly balanced by the pressure gradient and does not influence the velocity field. This invariance property can be preserved in the discrete problem if the incompressibility constraint is satisfied in a sufficiently strong sense. We derive sufficient conditions to guarantee discretely divergence-free functions are exactly divergence-free, and give examples of divergence-free finite elements on meshes containing triangular, quadrilateral, tetrahedral, or hexahedral cells generated by a (possibly non-affine) map from their respective reference cells. In the case of quadrilateral cells, we prove an optimal error estimate for the velocity field that does not depend on the pressure approximation. Our theoretical analysis is supported by numerical results.

We introduce two hybridizable discontinuous Galerkin (HDG) methods for numerically solving the Monge-Ampere equation. The first HDG method is devised to solve the nonlinear elliptic Monge-Ampere equation by using Newton's method. The second HDG method is devised to solve a sequence of the Poisson equation until convergence to a fixed-point solution of the Monge-Ampere equation is reached. Numerical examples are presented to demonstrate the convergence and accuracy of the HDG methods. Furthermore, the HDG methods are applied to r-adaptive mesh generation by redistributing a given scalar density function via the optimal transport theory. This r-adaptivity methodology leads to the Monge-Ampere equation with a nonlinear Neumann boundary condition arising from the optimal transport of the density function to conform the resulting high-order mesh to the boundary. Hence, we extend the HDG methods to treat the nonlinear Neumann boundary condition. Numerical experiments are presented to illustrate the generation of r-adaptive high-order meshes on planar and curved domains., Comment: 25 pages, 9 figures, and 8 tables

Finite element methods and kinematically coupled schemes that decouple the fluid velocity and structure's displacement have been extensively studied for incompressible fluid-structure interaction (FSI) over the past decade. While these methods are known to be stable and easy to implement, optimal error analysis has remained challenging. Previous work has primarily relied on the classical elliptic projection technique, which is only suitable for parabolic problems and does not lead to optimal convergence of numerical solutions to the FSI problems in the standard $L^2$ norm. In this article, we propose a new kinematically coupled scheme for incompressible FSI thin-structure model and establish a new framework for the numerical analysis of FSI problems in terms of a newly introduced coupled non-stationary Ritz projection, which allows us to prove the optimal-order convergence of the proposed method in the $L^2$ norm. The methodology presented in this article is also applicable to numerous other FSI models and serves as a fundamental tool for advancing research in this field.

This paper complements the empirical justification of the revised scheme in Part I of this work with a mathematical justification leveraging a semi-discrete analysis framework for assessing the splitting error of process coupling methods. The novelty of the framework is that splitting error is distinguished from the process time integration errors, i.e., the errors caused by discrete time integration of individual processes, leading to expressions that are more easily interpreted utilizing existing physical understanding of the processes that the terms represent. This application of this framework to dust life cycle in EAMv1 showcases such an interpretation, using the leading-order splitting error that results from the framework to confirm (i) that the original EAMv1 scheme artificially strengthens the effect of dry removal processes, and (ii) that the revised splitting reduces that artificial strengthening. While the error analysis framework is presented in the context of the dust life cycle in EAMv1, the framework can be broadly leveraged to evaluate process coupling schemes, both in other physical problems and for any number of processes. This framework will be particularly powerful when the various process implementations support a variety of time integration approaches. Whereas traditional local truncation error approaches require separate consideration of each combination of time integration methods, this framework enables evaluation of coupling schemes independent of particular time integration approaches for each process while still allowing for the incorporation of these specific time integration errors if so desired. The framework also explains how the splitting error terms result from (i) the integration of individual processes in isolation from other processes, and (ii) the choices of input state and timestep size for the isolated integration of processes.

This paper presents a new discretization error quantification method for the numerical integration of ordinary differential equations. The error is modelled by using the Wishart distribution, which enables us to capture the correlation between variables. Error quantification is achieved by solving an optimization problem under the order constraints for the covariance matrices. An algorithm for the optimization problem is also established in a slightly broader context.

Neural networks are universal approximators and are studied for their use in solving differential equations. However, a major criticism is the lack of error bounds for obtained solutions. This paper proposes a technique to rigorously evaluate the error bound of Physics-Informed Neural Networks (PINNs) on most linear ordinary differential equations (ODEs), certain nonlinear ODEs, and first-order linear partial differential equations (PDEs). The error bound is based purely on equation structure and residual information and does not depend on assumptions of how well the networks are trained. We propose algorithms that bound the error efficiently. Some proposed algorithms provide tighter bounds than others at the cost of longer run time., Comment: 10 page main artichle + 5 page supplementary material

Image classification is one of the most fundamental tasks in Computer Vision. In practical applications, the datasets are usually not as abundant as those in the laboratory and simulation, which is always called as Data Hungry. How to extract the information of data more completely and effectively is very important. Therefore, an Adaptive Data Augmentation Framework based on the tensor T-product Operator is proposed in this paper, to triple one image data to be trained and gain the result from all these three images together with only less than 0.1% increase in the number of parameters. At the same time, this framework serves the functions of column image embedding and global feature intersection, enabling the model to obtain information in not only spatial but frequency domain, and thus improving the prediction accuracy of the model. Numerical experiments have been designed for several models, and the results demonstrate the effectiveness of this adaptive framework. Numerical experiments show that our data augmentation framework can improve the performance of original neural network model by 2%, which provides competitive results to state-of-the-art methods.

We propose, analyze and implement a virtual element discretization for an interfacial poroelasticity-elasticity consolidation problem. The formulation of the time-dependent poroelasticity equations uses displacement, fluid pressure, and total pressure, and the elasticity equations are written in the displacement-pressure formulation. The construction of the virtual element scheme does not require Lagrange multipliers to impose the transmission conditions (continuity of displacement and total traction, and no-flux for the fluid) on the interface. We show the stability and convergence of the virtual element method for different polynomial degrees, and the error bounds are robust with respect to delicate model parameters (such as Lame constants, permeability, and storativity coefficient). Finally, we provide numerical examples that illustrate the properties of the scheme.

In this work, we prove rigorous error estimates for a hybrid method introduced in [15] for solving the time-dependent radiation transport equation (RTE). The method relies on a splitting of the kinetic distribution function for the radiation into uncollided and collided components. A high-resolution method (in angle) is used to approximate the uncollided components and a low-resolution method is used to approximate the the collided component. After each time step, the kinetic distribution is reinitialized to be entirely uncollided. For this analysis, we consider a mono-energetic problem on a periodic domains, with constant material cross-sections of arbitrary size. To focus the analysis, we assume the uncollided equation is solved exactly and the collided part is approximated in angle via a spherical harmonic expansion ($\text{P}_N$ method). Using a non-standard set of semi-norms, we obtain estimates of the form $C(\varepsilon,\sigma,\Delta t)N^{-s}$ where $s\geq 1$ denotes the regularity of the solution in angle, $\varepsilon$ and $\sigma$ are scattering parameters, $\Delta t$ is the time-step before reinitialization, and $C$ is a complicated function of $\varepsilon$, $\sigma$, and $\Delta t$. These estimates involve analysis of the multiscale RTE that includes, but necessarily goes beyond, usual spectral analysis. We also compute error estimates for the monolithic $\text{P}_N$ method with the same resolution as the collided part in the hybrid. Our results highlight the benefits of the hybrid approach over the monolithic discretization in both highly scattering and streaming regimes.

Motivated by various computational applications, we investigate the problem of estimating nested expectations. Building upon recent work by the authors, we propose a novel Monte Carlo estimator for nested expectations, inspired by sparse grid quadrature, that does not require sampling from inner conditional distributions. Theoretical analysis establishes an upper bound on the mean squared error of our estimator under mild assumptions on the problem, demonstrating its efficiency for cases with low-dimensional outer variables. We illustrate the effectiveness of our estimator through its application to problems related to value of information analysis, with moderate dimensionality. Overall, our method presents a promising approach to efficiently estimate nested expectations in practical computational settings.

In this work, we examine a numerical phase-field fracture framework in which the crack irreversibility constraint is treated with a primal-dual active set method and a linearization is used in the degradation function to enhance the numerical stability. The first goal is to carefully derive from a complementarity system our primal-dual active set formulation, which has been used in the literature in numerous studies, but for phase-field fracture without its detailed mathematical derivation yet. Based on the latter, we formulate a modified combined active-set Newton approach that significantly reduces the computational cost in comparison to comparable prior algorithms for quasi-monolithic settings. For many practical problems, Newton converges fast, but active set needs many iterations, for which three different efficiency improvements are suggested in this paper. Afterwards, we design an iteration on the linearization in order to iterate the problem to the monolithic limit. Our new algorithms are implemented in the programming framework pfm-cracks [T. Heister, T. Wick; pfm-cracks: A parallel-adaptive framework for phase-field fracture propagation, Software Impacts, Vol. 6 (2020), 100045]. In the numerical examples, we conduct performance studies and investigate efficiency enhancements. The main emphasis is on the cost complexity by keeping the accuracy of numerical solutions and goal functionals. Our algorithmic suggestions are substantiated with the help of several benchmarks in two and three spatial dimensions. Therein, predictor-corrector adaptivity and parallel performance studies are explored as well., Comment: 49 pages, 45 figures, 9 tables

For the two-layer shallow water equations, a high-order compact gas-kinetic scheme (GKS) on triangular mesh is proposed. The two-layer shallow water equations have complex source terms in comparison with the single layer equations. The main focus of this study is to construct a time-accurate evolution solution at a cell interface and to design a well-balanced scheme. The evolution model at a cell interface provides not only the numerical fluxes, but also the flow variables. The time-dependent flow variables at the closed cell interfaces can be used to update the cell-averaged gradients for the discretization of the the source terms inside each control volume in the development of the well-balanced scheme. Based on the cell-averaged flow variable and their gradients, high-order initial data reconstruction can be achieved with compact stencils. The compact high-order GKS has advantages to simulate the flow evolution in complex domain covered by unstructured mesh. Many test cases are used to validate the accuracy and robustness of the scheme for the two-layer shallow water equations.

Many integral equation-based methods are available for problems of time-harmonic electromagnetic scattering from perfect electric conductors. Moreover, there are numerous ways in which the geometry can be represented, numerous ways to represent the relevant surface current and/or charge densities, numerous quadrature methods that can be deployed, and numerous fast methods that can be used to accelerate the solution of the large linear systems which arise from discretization. Among the many issues that arise in such scattering calculations are the avoidance of spurious resonances, the applicability of the chosen method to scatterers of non-trivial topology, the robustness of the method when applied to objects with multiscale features, the stability of the method under mesh refinement, the ease of implementation with high-order basis functions, and the behavior of the method as the frequency tends to zero. Since three-dimensional scattering is a challenging, large-scale problem, many of these issues have been historically difficult to investigate. It is only with the advent of fast algorithms and modern iterative methods that a careful study of these issues can be carried out effectively. In this paper, we use GMRES as our iterative solver and the fast multipole method as our acceleration scheme in order to investigate some of these questions. In particular, we compare the behavior of the following integral equation formulations with regard to the issues noted above: the standard electric, magnetic, and combined field integral equations with standard RWG basis functions, the non-resonant charge-current integral equation, the electric charge-current integral equation, the augmented regularized combined source integral equation and the decoupled potential integral equation DPIE. Various numerical results are provided to demonstrate the behavior of each of these schemes., Comment: 27 pages, 9 figures

Solving inverse problems, where we find the input values that result in desired values of outputs, can be challenging. The solution process is often computationally expensive and it can be difficult to interpret the solution in high-dimensional input spaces. In this paper, we use a problem from additive manufacturing to address these two issues with the intent of making it easier to solve inverse problems and exploit their results. First, focusing on Gaussian process surrogates that are used to solve inverse problems, we describe how a simple modification to the idea of tapering can substantially speed up the surrogate without losing accuracy in prediction. Second, we demonstrate that Kohonen self-organizing maps can be used to visualize and interpret the solution to the inverse problem in the high-dimensional input space. For our data set, as not all input dimensions are equally important, we show that using weighted distances results in a better organized map that makes the relationships among the inputs obvious., Comment: 16 figures and 4 tables

We investigate the propagation of acoustic singular surfaces, specifically, linear shock waves and nonlinear acceleration waves, in a class of inhomogeneous gases whose ambient mass density varies exponentially. Employing the mathematical tools of singular surface theory, we first determine the evolution of both the jump amplitudes and the locations/velocities of their associated wave-fronts, along with a variety of related analytical results. We then turn to what have become known as Krylov subspace spectral (KSS) methods to numerically simulate the evolution of the full waveforms under consideration. These are not only performed quite efficiently, since KSS allows the use of `large' CFL numbers, but also quite accurately, in the sense of capturing theoretically-predicted features of the solution profiles more faithfully than other time-stepping methods, since KSS customizes the computation of the components of the solution corresponding to the different frequencies involved. The presentation concludes with a listing of possible, acoustics-related, follow-on studies.

In this paper, we develop a fast and accurate pseudospectral method to approximate numerically the fractional Laplacian $(-\Delta)^{\alpha/2}$ of a function on $\mathbb{R}$ for $\alpha=1$; this case, commonly referred to as the half Laplacian, is equivalent to the Hilbert transform of the derivative of the function. The main ideas are as follows. Given a twice continuously differentiable bounded function $u\in \mathit{C}_b^2(\mathbb{R})$, we apply the change of variable $x=L\cot(s)$, with $L>0$ and $s\in[0,\pi]$, which maps $\mathbb{R}$ into $[0,\pi]$, and denote $(-\Delta)_s^{1/2}u(x(s)) \equiv (-\Delta)^{1/2}u(x)$. Therefore, by performing a Fourier series expansion of $u(x(s))$, the problem is reduced to computing $(-\Delta)_s^{1/2}e^{iks} \equiv (-\Delta)^{1/2}(x + i)^k/(1+x^2)^{k/2}$. On a previous work, we considered the case with $k$ even and $\alpha\in(0,2)$, so we focus now on the case with $k$ odd. More precisely, we express $(-\Delta)_s^{1/2}e^{iks}$ for $k$ odd in terms of the Gaussian hypergeometric function ${}_2F_1$, and also as a well-conditioned finite sum. Then, we use a fast convolution result, that enable us to compute very efficiently $\sum_{l = 0}^Ma_l(-\Delta)_s^{1/2}e^{i(2l+1)s}$, for extremely large values of $M$. This enables us to approximate $(-\Delta)_s^{1/2}u(x(s))$ in a fast and accurate way, especially when $u(x(s))$ is not periodic of period $\pi$., Comment: 22 pages, 5 figures

In this paper, a class of smoothing modulus-based iterative method was presented for solving implicit complementarity problems. The main idea was to transform the implicit complementarity problem into an equivalent implicit fixed-point equation, then introduces a smoothing function to obtain its approximation solutions. The convergence analysis of the algorithm was given, and the efficiency of the algorithms was verified by numerical experiment