Research
Our group works on the mathematical analysis and numerical approximation of multiscale problems arising from cell biology, neuroscience, and quantum mechanics. Below is an overview of our five main directions.
1. Collective Cell Dynamics
What mathematical laws govern the shape, size, and internal structure of a growing tumor?
Tissues are not passive fluids; they are active, living materials where individual cells consume nutrients, proliferate, die, and respond to mechanical stress. At the macroscopic scale, this produces striking patterns: fingering invasion fronts, necrotic cores, and morphological phase transitions. Understanding these patterns requires bridging cell-level biochemistry with tissue-level mechanics.
Mathematical challenges. Classical reaction–diffusion models ignore cell–cell adhesion and volume exclusion. Free-boundary formulations are more faithful but introduce geometric nonlinearities and singular limits (e.g., the Hele–Shaw or porous-media limit). We ask: Can one rigorously derive these macroscopic laws from cell-density models? And how do nutrient supply and mechanical feedback determine the final morphology?
Our approach. Together with collaborators, we analyzed the incompressible limit of cell-density models with pressure laws, rigorously justifying the emergence of free-boundary dynamics. We proved existence of nonsymmetric traveling waves in Hele–Shaw type tumor models and showed how nutrient consumption can destabilize the tumor boundary. A recent line of work treats tumor growth with a necrotic core as an obstacle problem in pressure.
- J.-G. Liu, M. Tang, L. Wang & Z. Zhou, An accurate front capturing scheme for tumor growth models with a free boundary limit, J. Comput. Phys. 364 (2018). arXiv
- Y. Feng, M. Tang, X. Xu & Z. Zhou, Tumor boundary instability induced by nutrient consumption and supply, Z. Angew. Math. Phys. 74 (2023). arXiv
- X. Dou, C. Shen & Z. Zhou, Tumor growth with a necrotic core as an obstacle problem in pressure, Acta Appl. Math. 191 (2024). arXiv
2. Neurodynamics & Kinetic Theory
Can a population of ten thousand noisy neurons be described by a single partial differential equation?
The brain is a multiscale tower: ion channels open in milliseconds, single neurons spike in milliseconds, and networks synchronize over seconds. Rather than simulating every neuron, kinetic theory coarse-grains the population into probability densities over voltage and conductance variables. The resulting PDEs are reminiscent of statistical physics—but with discontinuities (reset and threshold) that create fascinating blow-up and synchronization phenomena.
Mathematical challenges. The Fokker–Planck equations for integrate-and-fire networks contain nonlocal reset terms and can develop singularities in finite time. Pulse-coupled oscillators lead to mean-field systems with discontinuous velocities. We need to prove global well-posedness, long-time convergence to equilibrium, and the validity of the mean-field limit itself.
Our approach. We introduced a generalized solution framework for the NNLIF neuron model that continues past blow-up via a dilation of time, proving global well-posedness. For the voltage-conductance kinetic system, we established exponential ergodicity through a probabilistic reformulation. We also designed multiscale solvers that capture synchronization in noisy integrate-and-fire networks without resolving every spike.
- Z. Du, Y. Xie & Z. Zhou, A synchronization-capturing multi-scale solver to the noisy integrate-and-fire neuron networks, Multiscale Model. Simul. 22 (2024). arXiv:2305.05915
- J. Carrillo, X. Dou & Z. Zhou, A simplified voltage-conductance kinetic model for interacting neurons and its asymptotic limit, SIAM J. Math. Anal. 56 (2024). arXiv
- J. Hu, J.-G. Liu, Y. Xie & Z. Zhou*, A structure preserving numerical scheme for Fokker-Planck equations of neuron networks: numerical analysis and exploration, J. Comput. Phys. 433 (2021), 110195. arXiv
3. Quantum Dynamics & Nonadiabatic Computation
How do electrons and atomic nuclei move together when their time scales differ by orders of magnitude?
In molecular quantum mechanics, the Born–Oppenheimer approximation decouples fast electrons from slow nuclei. Yet many chemical reactions—photosynthesis, vision, photovoltaics—rely on nonadiabatic transitions where this separation breaks down. The wave function lives in a high-dimensional configuration space, making direct simulation prohibitively expensive.
Mathematical challenges. We need asymptotic methods that capture quantum interference without resolving every oscillation, and stochastic algorithms whose sampling error can be controlled rigorously.
Our approach. We developed Frozen Gaussian Sampling—a mesh-free Monte Carlo method based on the Gaussian wave-packet transform—and proved its convergence for semiclassical Schrödinger equations with vector potentials. For nonadiabatic dynamics, we gave the first rigorous justification of surface-hopping algorithms in the mixed quantum-classical scaling, and designed efficient samplers for metal-surface dynamics. Recent work also addresses Bloch electrons with Weyl nodes using multiscale analysis.
- Y. Xie & Z. Zhou, Frozen Gaussian Sampling: A Mesh-free Monte Carlo Method For Approximating Semiclassical Schrödinger Equations, Comm. Math. Sci. 22 (2024). arXiv:2112.05405
- Z. Huang, L. Xu & Z. Zhou, Efficient Frozen Gaussian Sampling Algorithms for Nonadiabatic Quantum Dynamics at Metal Surfaces, J. Comput. Phys. 474 (2023). arXiv:2206.02173
- J. Lu & Z. Zhou, Frozen Gaussian approximation with surface hopping for mixed quantum-classical dynamics, Math. Comp. 87 (2018). arXiv:1602.06459
4. Statistical Sampling & Molecular Dynamics
How can we compute thermal averages of quantum systems without paying the exponential price of dimension?
Many properties of materials—heat capacity, reaction rates, spectral densities—are thermal averages. For quantum systems, Feynman’s path-integral formulation converts the quantum partition function into a classical-looking integral over polymer ring configurations, but the dimension scales with the number of particles and imaginary-time slices.
Mathematical challenges. Standard molecular dynamics mixes physical time and imaginary time, leading to resonances and slow convergence. Interacting particle systems require O(N²) force evaluations per step. We need algorithms that are simultaneously ergodic, dimension-free, and cheap per step.
Our approach. We proved dimension-free ergodicity for path-integral molecular dynamics (PIMD) and gave quantitative convergence rates for the Random Batch Method (RBM) in interacting particle systems. By combining RBM with PIMD, we achieved efficient thermal-average sampling for quantum particle systems. We also developed multi-level Monte Carlo strategies for nonadiabatic regimes and infinite-swapping algorithms to accelerate rare-event sampling.
- X. Ye & Z. Zhou, Dimension-free Ergodicity of Path Integral Molecular Dynamics, Comm. Comput. Phys. 38 (2025). arXiv:2307.06510
- S. Jin, L. Li, X. Ye & Z. Zhou, Ergodicity and long-time behavior of the Random Batch Method for interacting particle systems, Math. Models Methods Appl. Sci. 33 (2023). arXiv:2202.04952
- J. Lu, Y. Lu & Z. Zhou, Continuum limit and preconditioned Langevin sampling of the path integral molecular dynamics, J. Comput. Phys. 423 (2020). arXiv:1811.10995
5. Learning Theory & Optimal Control
Can classical mathematical tools make complex agent systems more transparent and principled?
Large language models and multi-agent systems are, at their core, high-dimensional dynamical systems with nonlinear feedback. Classical control theory and game theory offer rigorous languages—Hamiltonians, adjoint equations, mean-field limits, Pontryagin principles—that can illuminate their training dynamics and decision-making structures. Our interest is not in building bigger models, but in understanding them mathematically: what do they optimize, where do they fail, and how can we certify their behavior?