Generative Monte Carlo replaces repeated within-cell scattering histories by samples from a learned conditional transmission distribution. The method uses conditional flow matching to generate a particle's exit state from its entry state and optical cell parameters. In the tested regimes, it retains the characteristic statistical convergence of Monte Carlo while reducing the cost of optically thick cell crossings. The work also examines learned bias, training coverage, and diagnostics for use outside the training regime.

@article{farmer2026generative,
  title = {Generative Monte Carlo Sampling for Constant-Cost Particle Transport},
  author = {Farmer, Joseph A. and Murray, Aidan and Krotz, Johannes and McClarren, Ryan G.},
  journal = {Proceedings of the National Academy of Sciences},
  note = {Accepted for publication},
  year = {2026}
}

Monte Carlo transport simulations follow many randomly sampled particle paths. In very opaque materials, simulating every scattering event can be expensive. This work trains a model to sample the overall effect of many collisions at once, then uses those samples within a standard particle simulation. The goal is a faster calculation without losing the uncertainty information that makes Monte Carlo useful.

This work develops a labeled-particle framework for hybrid Monte Carlo--discrete ordinates transport. Particles can switch between explicit-scattering and hybridized descriptions, while an accounting identity ensures that scattering is neither omitted nor counted twice. The formulation supports switching policies based on physical or computational information and allows quasi-Monte Carlo sampling on fixed-dimensional source variables.

@article{krotz2026labeled,
  title = {Labeled-particle hybrid Monte Carlo--discrete ordinates transport with quasi-Monte Carlo sampling},
  author = {Krotz, Johannes and McClarren, Ryan G.},
  journal = {Nuclear Science and Engineering},
  note = {Accepted, to appear},
  year = {2026}
}

The paper develops a flexible way to decide when a particle simulation should explicitly follow scattering events and when it should hand smoother behavior to a deterministic solver. The labels keep track of what each representation is responsible for, so the two pieces can be combined without losing or double-counting physical effects.

We introduce a nonlinear three-temperature reduction for non-grey radiative transfer. The model adds a spectral degree of freedom to a grey-like transport structure, establishes energy conservation and positive steady states, and characterizes local stability and spectral relaxation in an auxiliary infinite-medium problem. In near-Planckian tests, it captures spectral information unavailable to grey transport at a cost closer to low-group calculations.

@article{krotz2026threetemperature,
  title = {A Nonlinear, 3 Temperature Model for Non-Grey Radiative Transfer},
  author = {Krotz, Johannes and McClarren, Ryan G.},
  journal = {SIAM Journal on Applied Mathematics},
  note = {Under review},
  year = {2026}
}

Light and radiation can have different energies, and tracking all of them is costly. This model keeps a small amount of extra spectral information beyond a standard grey approximation. It is designed to be much cheaper than a full multigroup calculation while making clear when its simplified description is reliable and when it is not.

We study sliding magnetic rotor arrays through a reduced nonlinear ODE model. The work decomposes the energy landscape into chambers, classifies critical points and bifurcation boundaries, and derives asymptotic expressions for magnetic friction. Numerical integration verifies the predicted switching mechanisms and connects friction peaks with dissipative jumps between unstable and stable configurations.

@article{krotz2026rotor,
  title = {Energy landscape and dissipation of sliding magnetic rotor arrays},
  author = {Krotz, Johannes and Lueders, Anton},
  journal = {Physical Review E},
  note = {Under review},
  year = {2026}
}

Arrays of small magnets can resist motion in sudden, stick-slip-like jumps. This work uses a reduced mathematical model to identify the stable arrangements, the points where those arrangements change, and how those changes produce dissipation and magnetic friction.

This work develops a neural closure for the M1 radiation-transfer system that preserves hyperbolicity. The closure Jacobian is parameterized using a symmetric matrix network and a strictly convex entropy network whose Hessian supplies a positive-definite symmetrizer. The resulting architecture retains learning flexibility while enforcing real characteristic speeds required for a stable hyperbolic PDE solver.

@article{kim2026hyperbolic,
  title = {A Hyperbolic Neural Closure for M1 Radiation Transfer},
  author = {Kim, Bongseok and Zhang, Jiahao and Krotz, Johannes and Balsara, Dinshaw and McClarren, Ryan G. and Lin, Guang},
  journal = {Computer Methods in Applied Mechanics and Engineering},
  note = {Under review},
  year = {2026}
}

A machine-learned shortcut can make a physics simulation fail if it breaks a property that the numerical solver needs. This work builds that property into the neural-network architecture, so the learned closure remains compatible with wave-like radiation equations rather than merely fitting sample data.

This conference paper studies a scatter-limited hybrid Monte Carlo--deterministic transport method together with quasi-Monte Carlo sampling. It investigates how a finite-scatter work-sharing rule can balance particle and deterministic components while retaining accurate transport calculations.

@inproceedings{krotz2026scatterlimited,
  title = {Scatter-Limited Hybrid Monte Carlo, Deterministic Transport with Quasi-Monte Carlo Sampling},
  author = {Krotz, Johannes and McClarren, Ryan G.},
  booktitle = {Proceedings of PHYSOR 2026},
  year = {2026}
}

This work explores how to share a transport calculation between random particle histories and a deterministic solver. Limiting the number of explicitly simulated scattering events can reduce cost, and quasi-Monte Carlo sampling offers a way to make the remaining random inputs more evenly distributed.

This conference manuscript studies statistical models for generative Monte Carlo particle transport, where learned components are used to represent expensive portions of particle histories.

@inproceedings{farmer2027statistical,
  title = {Statistical Models for Generative Monte Carlo Particle Transport},
  author = {Farmer, Joseph A. and Krotz, Johannes and McClarren, Ryan G. and Palmer, Todd S.},
  note = {Submitted to M&C 2027}
}

This manuscript examines the statistical modeling questions that arise when generative models accelerate Monte Carlo particle simulations.

This conference manuscript investigates reversible-computing ideas for Monte Carlo particle transport.

@inproceedings{denico2027reversible,
  title = {Reversible Computing for Monte Carlo Particle Transport},
  author = {Denico, Aubrey M. and Krotz, Johannes and McClarren, Ryan G.},
  note = {Submitted to M&C 2027}
}

This work explores whether reversible-computing ideas can help make particle-transport simulations more efficient.

This conference manuscript studies the recovery of independent Monte Carlo histories in nonlinear Boltzmann transport.

@inproceedings{krotz2027independent,
  title = {Recovering Independent Monte Carlo Histories in Nonlinear Boltzmann Transport},
  author = {Krotz, Johannes and McClarren, Ryan G.},
  note = {Submitted to M&C 2027}
}

This work concerns how a nonlinear particle-transport simulation can recover statistically independent particle histories, a key requirement for reliable Monte Carlo uncertainty estimates.

We derive a general expression for the streaming term in radiative transport equations and other transport problems when formulated in curvilinear coordinates, emphasizing coordinate systems adapted to the geometry of the domain and the directional dependence of particle transport. By parametrizing the angular variable using a local orthonormal frame, we express directional derivatives in terms of curvature-related quantities that reflect the geometry of underlying spatial manifolds. Our formulation highlights how the interaction between coordinate choices and curvature influences the streaming operator, offering geometric interpretations of its components. The resulting framework offers intuitive insight into when and how angular dependence can be simplified and may guide the selection of coordinate systems that balance analytical tractability and computational efficiency.

@misc{krotz2025curvilinearcoordinatescurvatureradiative,
      title={Curvilinear coordinates and curvature in radiative transport}, 
      author={Johannes Krotz and Ryan G. McClarren},
      year={2025},
      eprint={2508.20852},
      archivePrefix={arXiv},
      primaryClass={math.NA},
      url={https://arxiv.org/abs/2508.20852}, 
}

Transport equations are used to model how light, radiation, or particles move through a medium. The part of the equation that describes particles moving through space and changing direction is the streaming term. In flat Cartesian coordinates this term has a familiar form. In curved or geometry-adapted coordinates, the same term becomes more delicate.

Many transport problems are not naturally box-shaped. They may involve spherical shells, cylindrical devices, curved boundaries, or coordinate systems chosen to follow the geometry of the domain. In those settings, the coordinates affect how particle directions are represented. A coordinate system can simplify the equations, but it can also introduce extra terms that have to be understood carefully.

What the Paper Does

This paper derives the streaming term for transport equations written in curvilinear coordinates. Instead of treating every coordinate system as a separate calculation, it gives a unified formula based on a local orthonormal frame and curvature-related quantities.

The point is that the streaming operator depends not only on position and direction, but also on how the chosen coordinate surfaces bend and how the local directions rotate from point to point. Written this way, the extra terms in the transport equation have a direct geometric meaning.

• The coordinate geometry determines which additional streaming terms appear.
• Flat, curved, and twisting directions contribute in different ways.
• Curvature gives a compact way to see how geometry enters the transport equation.

Why This Is Useful

Deriving the streaming term in a curved coordinate system is often done case by case. Spherical, cylindrical, elliptical, and other coordinates each lead to their own formulas. The curvature-based formulation gives a common structure behind these calculations.

• Simplification: it shows when terms disappear because of the geometry.
• Interpretation: it replaces some long coordinate calculations with geometric quantities.
• Guidance: it helps identify coordinate systems that make the equations easier to analyze or compute.

Examples

The paper works through standard settings such as cylindrical and spherical coordinates, and also considers more involved examples such as ellipsoids and translating graphs. The same framework explains why some geometries produce simple formulas while others require additional curvature terms.

• Cylinders: several terms vanish because parts of the geometry are flat or do not twist.
• Spheres: curvature enters in a symmetric and recognizable way.
• Ellipsoids and translating surfaces: more terms remain, but the curvature formulation still shows where they come from.

Takeaway

The paper does not introduce a new numerical solver. It gives a geometric way to derive and interpret the streaming term. This can make derivations cleaner, help remove unnecessary terms, and support better coordinate choices for transport simulations in curved geometries.

A Bayesian data assimilation scheme is formulated for advection-dominated advective and diffusive evolutionary problems, based upon the Dynamic Likelihood (DLF) approach to filtering. The DLF was developed specifically for hyperbolic problems –waves–, and in this paper, it is extended via a split step formulation, to handle advection-diffusion problems. In the dynamic likelihood approach, observations and their statistics are used to propagate probabilities along characteristics, evolving the likelihood in time. The estimate posterior thus inherits phase information. For advection-diffusion the advective part of the time evolution is handled on the basis of observations alone, while the diffusive part is informed through the model as well as observations. We expect, and indeed show here, that in advection-dominated problems, the DLF approach produces better estimates than other assimilation approaches, particularly when the observations are sparse and have low uncertainty. The added computational expense of the method is cubic in the total number of observations over time, which is on the same order of magnitude as a standard Kalman filter and can be mitigated by bounding the number of forward propagated observations, discarding the least informative data.

@article{KROTZ_2025_DynamicLikelihoodFilter,
title = {A dynamic likelihood approach to filtering transport processes: advection-diffusion dynamics},
journal = {Journal of Computational Physics},
volume = {536},
pages = {114089},
year = {2025},
issn = {0021-9991},
doi = {https://doi.org/10.1016/j.jcp.2025.114089},
url = {https://www.sciencedirect.com/science/article/pii/S0021999125003729},
author = {Johannes Krotz and Juan M. Restrepo and Jorge Ramirez},
keywords = {Data assimilation, Bayesian estimation, Dynamic likelihood filter, 
Advection-diffusion, Transport, Kalman filter},
abstract = {A Bayesian data assimilation scheme is formulated for advection-dominated advective 
        and diffusive evolutionary problems, based upon the Dynamic Likelihood (DLF) approach to 
        filtering. The DLF was developed specifically for hyperbolic problems –waves–, and in this
        paper, it is extended via a split step formulation, to handle advection-diffusion problems. 
        In the dynamic likelihood approach, observations and their statistics are used to propagate 
        probabilities along characteristics, evolving the likelihood in time. The estimate posterior
        thus inherits phase information. For advection-diffusion the advective part of the time 
        evolution is handled on the basis of observations alone, while the diffusive part is informed 
        through the model as well as observations. We expect, and indeed show here, that in 
        advection-dominated problems, the DLF approach produces better estimates than other 
        assimilation approaches, particularly when the observations are sparse and have low 
        uncertainty. The added computational expense of the method is cubic in the total number of 
        observations over time, which is on the same order of magnitude as a standard Kalman filter 
        and can be mitigated by bounding the number of forward propagated observations, discarding 
        the least informative data.}}}
        

Advection-diffusion equations model quantities that are carried by a flow and also spread out over time. Examples include smoke in air, pollutants in a river, temperature, or concentration fields. Predictions usually combine a physical model with sensor data, a process known as data assimilation.

A standard approach is the Kalman filter, which updates a model whenever new observations arrive. It is effective in many settings, but it can struggle when observations are sparse. In transport problems, sparse data are especially difficult because the main error is often not just the amount of material, but where that material has moved.

What the Paper Does

This paper extends the Dynamic Likelihood Filter (DLF) to advection-diffusion problems. The DLF uses observations not only at the time they are taken, but also propagates their information forward along the flow of the system. This creates additional likelihood information at later times and locations, even where no sensor is present.

For the advective part of the dynamics, the method follows the motion suggested by the observations. For the diffusive part, it combines the model with the observed information. This split lets the filter retain useful phase and location information while still accounting for spreading and smoothing.

Where It Helps

The method is designed for cases where transport dominates diffusion and the data are limited but reliable.

• Sparse observations: only a few sensors are available.
• Advection-dominated dynamics: flow or wind moves the quantity more strongly than diffusion spreads it.
• Imperfect models: the governing equations are useful, but they do not capture every detail.

These conditions appear in weather forecasting, environmental monitoring, contaminant tracking, and engineering transport problems.

What the Results Show

The numerical tests show that the DLF gives better estimates than a standard Kalman filter when advection is dominant and observations are sparse with low uncertainty. It improves both the magnitude and location of the transported quantity, corrects initial-condition errors more effectively, and remains useful when the model is not exact.

The method has an added computational cost, scaling cubically with the number of observations used over time. This is comparable in order to a standard Kalman filter, and the cost can be reduced by limiting how many propagated observations are kept.

Takeaway

The paper shows how to extend the Dynamic Likelihood idea from wave problems to advection-diffusion. By carrying observational information forward with the flow, the method makes sparse measurements more useful for tracking transported quantities.

We present a hybrid method for time-dependent particle transport problems that combines Monte Carlo (MC) estimation with deterministic solutions based on discrete ordinates. For spatial discretizations, the MC algorithm computes a piecewise constant solution and the discrete ordinates use bilinear discontinuous finite elements. From the hybridization of the problem, the resulting problem solved by Monte Carlo is scattering free, resulting in a simple, efficient solution procedure. Between time steps, we use a projection approach to “relabel” collided particles as uncollided particles. From a series of standard 2-D Cartesian test problems we observe that our hybrid method has improved accuracy and reduction in computational complexity of approximately an order of magnitude relative to standard discrete ordinates solutions.

@article{KROTZ_2024_HybridMCDG,
title = {A hybrid Monte Carlo, discontinuous Galerkin method for linear kinetic transport equations},
journal = {Journal of Computational Physics},
volume = {514},
pages = {113253},
year = {2024},
issn = {0021-9991},
doi = {https://doi.org/10.1016/j.jcp.2024.113253},
url = {https://www.sciencedirect.com/science/article/pii/S0021999124005011},
author = {Johannes Krotz and Cory D. Hauck and Ryan G. McClarren},
keywords = {Hybrid stochastic-deterministic method, Monte Carlo, Kinetic equations, Particle transport},
abstract = {We present a hybrid method for time-dependent particle transport problems that combines
         Monte Carlo (MC) estimation with deterministic solutions based on discrete ordinates. For 
         spatial discretizations, the MC algorithm computes a piecewise constant solution and the 
         discrete ordinates use bilinear discontinuous finite elements. From the hybridization of 
         the problem, the resulting problem solved by Monte Carlo is scattering free, resulting in 
         a simple, efficient solution procedure. Between time steps, we use a projection approach to 
         “relabel” collided particles as uncollided particles. From a series of standard 2-D Cartesian
         test problems we observe that our hybrid method has improved accuracy and reduction in 
         computational complexity of approximately an order of magnitude relative to standard discrete
         ordinates solutions.}}
    

Particle transport problems appear in radiation shielding, nuclear engineering, atmospheric radiative transfer, medical physics, plasma physics, and astrophysics. The goal is to predict how particles move through a material, scatter, and sometimes get absorbed.

Two common approaches are deterministic solvers and Monte Carlo methods. Deterministic solvers turn the physics into large systems of equations. They can be accurate, but the cost grows quickly at high resolution. Monte Carlo methods follow randomly sampled particle paths. They are flexible, but the results contain sampling noise unless many samples are used.

What the Paper Does

This paper combines the two approaches by separating particles into two groups: particles that have not collided yet and particles that have already scattered.

• Uncollided particles: these move along relatively simple paths and are handled efficiently with Monte Carlo.
• Collided particles: after scattering, their distribution is smoother and is handled with a deterministic discrete ordinates/discontinuous Galerkin solver.

The split makes the Monte Carlo part scattering-free, which simplifies that portion of the computation. Between time steps, a projection or relabeling step moves information back into the uncollided category so the method can continue efficiently over time.

Results

On standard two-dimensional Cartesian test problems, the hybrid method improves accuracy compared with standard discrete ordinates solutions and reduces computational complexity by about an order of magnitude. It also avoids some of the typical drawbacks of the separate methods: less sampling noise than a purely Monte Carlo calculation and fewer ray-effect artifacts than a purely deterministic angular discretization.

Why It Matters

The method is relevant anywhere particles stream and scatter through matter, including nuclear systems, radiation transport, atmospheric modeling, medical physics, and fusion-related simulations. The main advantage is not that it replaces either Monte Carlo or deterministic methods, but that it uses each where it is better suited.

What Comes Next

The paper focuses on single-energy problems. Natural extensions include multi-energy neutron transport, nonlinear radiative transfer, and adaptive strategies that decide how to balance the stochastic and deterministic parts during a simulation.

We present the near-Maximal Algorithm for Poisson-disk Sampling (nMAPS) to generate point distributions for variable resolution Delaunay triangular and tetrahedral meshes in two and three-dimensions, respectively. nMAPS consists of two principal stages. In the first stage, an initial point distribution is produced using a cell-based rejection algorithm. In the second stage, holes in the sample are detected using an efficient background grid and filled in to obtain a near-maximal covering. Extensive testing shows that nMAPS generates a variable resolution mesh in linear run time with the number of accepted points. We demonstrate nMAPS capabilities by meshing three-dimensional discrete fracture networks (DFN) and the surrounding volume. The discretized boundaries of the fractures, which are represented as planar polygons, are used as the seed of 2D-nMAPS to produce a conforming Delaunay triangulation. The combined mesh of the DFN is used as the seed for 3D-nMAPS, which produces conforming Delaunay tetrahedra surrounding the network. Under a set of conditions that naturally arise in maximal Poisson-disk samples and are satisfied by nMAPS, the two-dimensional Delaunay triangulations are guaranteed to only have well-behaved triangular faces. While nMAPS does not provide triangulation quality bounds in more than two dimensions, we found that low-quality tetrahedra in 3D are infrequent, can be readily detected and removed, and a high quality balanced mesh is produced.

@article{KROTZ_2022_PoissonDiskDFN,
title = {Variable resolution Poisson-disk sampling for meshing discrete fracture networks},
journal = {Journal of Computational and Applied Mathematics},
volume = {407},
pages = {114094},
year = {2022},
issn = {0377-0427},
doi = {https://doi.org/10.1016/j.cam.2022.114094},
url = {https://www.sciencedirect.com/science/article/pii/S0377042722000073},
author = {Johannes Krotz and Matthew R. Sweeney and Carl W. Gable and Jeffrey D. Hyman and Juan M. Restrepo},
keywords = {Discrete fracture network, Maximal Poisson-disk sampling, Mesh generation, Conforming Delaunay triangulation},
abstract = {We present the near-Maximal Algorithm for Poisson-disk Sampling (nMAPS) to generate point 
        distributions for variable resolution Delaunay triangular and tetrahedral meshes in two and 
        three-dimensions, respectively. nMAPS consists of two principal stages. In the first stage, 
        an initial point distribution is produced using a cell-based rejection algorithm. In the 
        second stage, holes in the sample are detected using an efficient background grid and filled
        in to obtain a near-maximal covering. Extensive testing shows that nMAPS generates a variable 
        resolution mesh in linear run time with the number of accepted points. We demonstrate nMAPS 
        capabilities by meshing three-dimensional discrete fracture networks (DFN) and the surrounding 
        volume. The discretized boundaries of the fractures, which are represented as planar polygons, are 
        used as the seed of 2D-nMAPS to produce a conforming Delaunay triangulation. The combined mesh
        of the DFN is used as the seed for 3D-nMAPS, which produces conforming Delaunay tetrahedra 
        surrounding the network. Under a set of conditions that naturally arise in maximal Poisson-disk 
        samples and are satisfied by nMAPS, the two-dimensional Delaunay triangulations are guaranteed
        to only have well-behaved triangular faces. While nMAPS does not provide triangulation quality 
        bounds in more than two dimensions, we found that low-quality tetrahedra in 3D are infrequent, 
        can be readily detected and removed, and a high quality balanced mesh is produced.}}
   

Simulations of groundwater, oil, gas, or contaminant transport in fractured rock need a mesh that represents both the fractures and the surrounding volume. The mesh is made of triangles in two dimensions or tetrahedra in three dimensions. Its quality matters because poorly shaped elements can make numerical solvers inaccurate or unstable.

There is also a cost issue. A mesh that is too coarse misses important geometric detail, while a mesh that is too fine is expensive to use. This paper introduces nMAPS, the near-Maximal Algorithm for Poisson-disk Sampling, to generate variable-resolution point sets for meshes of discrete fracture networks.

What the Method Does

nMAPS uses Poisson-disk sampling: points are placed so that they are not too close to one another, but still cover the domain well. The method allows the spacing to vary, so more points can be placed near intersections, boundaries, or other important features, while smoother regions can use fewer points.

After the points are placed, they are connected using a Delaunay triangulation. In two dimensions, the method gives well-behaved triangular faces under conditions satisfied by the samples produced by nMAPS. In three dimensions, it does not give the same formal quality bounds, but low-quality tetrahedra are uncommon and can be detected and removed in practice.

Why It Matters

The method gives a direct way to build meshes that are fine where they need to be and coarser elsewhere. Tests show that nMAPS runs in linear time with respect to the number of accepted points and produces balanced meshes for three-dimensional fracture networks and their surrounding volume.

This is useful for simulations of fractured media, including groundwater flow, energy extraction, and contaminant transport. The same idea may also be useful in other applications where variable-resolution meshes are needed around complicated geometry.

Summary

nMAPS provides an efficient way to generate point distributions for Delaunay meshes. It is not a complete solution to every meshing problem, but it gives a practical route to meshes that balance resolution, cost, and element quality for discrete fracture networks.