2) Understanding Defects in Liquid Crystals

Defects in liquid crystals are very interesting phenomena and have been the subject of research for decades, but still, predicting defect patterns accurately are a challenging problem. There are two major difficulties. The first one is how to model defects. Traditional vector models, such as the Oseen-Frank model, are not enough to capture the symmetry breaking behaviors caused by defects. Tensor model, such as the Landau-de Gennes model, has a larger configuration space, and is believed to be able to give a more realistic description. With his colleagues, Zhang showed that different models can lead to quite different defect structures. The second one is how to compute defects. Near defects, order parameters are close to singular. In addition, the stability of defects is very sensitive to model parameters and boundary condition. In order to capture defects, numerical algorithms need to be highly accurate. To this end, they implemented a high accurate spectral method of Q-tensor model and studied defect patterns of liquid crystals in confined space, and discovered many interesting properties of defects that have not been found before. In fact, various defect patterns are obtained, among them disclination lines are more common than point defects; the local structure of disclination lines is studied in detail; the profiles of disclination lines are analytically given. They also found that vector models can only describe point defects, but not disclination lines. See [Hu-Qu-Zhang,Sci. Sin. Math.(2015); Hu-Qu-Zhang,Commun. Comput. Phys.(2016); Qu-Wei-Zhang,Commun. Comput. Phys.(2017)] for more details.

3) Exploring Ordered Structures in Polymers

Block copolymers constructed by linking chemically distinct subchains together can often self-assemble into exquisitely ordered structures. These properties can be used to make novel materials such as microelectronic materials, advanced plastics, nanotemplates. A major challenge in the study of the phase behavior of structured polymer materials is to obtain different stable and metastable phases of the system. Experiments have unveiled an increasing number of interesting structures for block copolymers. In order to search for novel patterns, guidance from numerical prediction can be extremely valuable. Mathematically, predicting phase behavior amounts to minimize a free energy functional for the system (self-consistent field theory), which is a difficult task because of its nonlinearity, multi-solutions, and high dimension. Zhang and his colleagues developed a toolbox for numerically predicting the ordered structure in polymers, consisting of a set of strategies for choosing initial conditions, which they believe is essential in order to obtain certain types of structures, targeted discretization methods that can help finding the correct solution fast and accurately, and self-adapting bounding box which allows them to obtain the optimized periodic structures. They applied the above numerical method to AB diblock copolymers, rod-coil diblock copolymers and ABC star triblock copolymers and obtained a large number of patterns, including new phases that have not been found before. They also investigated the dynamics of phase transitions in polymers. The idea is to find the minimum energy path connecting two distinct phases over the high-dimensional energy landscape. Using the string method, the nucleation process in which a new order emerges from within another order can be obtained. This work demonstrates the power of numerical simulation in understanding complex systems. More details can be found in [Cheng-Lin-E-Zhang- Shi,Phys. Rev. Lett.(2010); Lin-Cheng-E-Shi-Zhang,J. Comput. Phys.(2010); Jiang-Huang-Zhang,J. Comput. Phys.(2010); Xu-Jiang-Zhang,Journal of Physical Chemistry B(2013)].

2. Moment Methods in Kinetic Theory (LI Ruo)

The kinetic equations, including Boltzmann equation, radiative transfer equation and Wigner equation, model physical problems in mesoscopic scaling. On one hand, it is more reliable than classical macroscopic models; on the other hand, it is a high dimensional equation that direct numerical study is prohibitive. A model reduction has to be applied to reduce the dimension of independent variables. The major contributions of Ruo Li and his collaborators’ studies in the past few years on the model reduction and numerical method for Boltzmann equation are in three folds:

1) Further Understanding on the Loss of Hyperbolicity of Grad's 13 Moment System

After given the 13 moment system in 1949, Grad himself should have realized that the 13 moment system, together with else modes in his paper, are not globally hyperbolic. It was explicitly pointed out by Muller et al. in 1998 that the 13 moment system in 1D is not globally hyperbolic, while hyperbolic around the local equilibrium. Intuitively, one may believe that the 13 moment system in 3D is hyperbolic around the local equilibrium. Li et al. proved in [Cai-Fan-Li, Kinetic and Related Models (2014)] that this is NOT the case. In fact, the local equilibrium is on the boundary of the hyperbolicity region of the 13 moment system in 3D. Consequently, an arbitrary small perturbation on the local equilibrium may lead to the non-existence of the solution of the system.

2) Globally Hyperbolic Regularization of Grad's Moment System

Li et al. proposed a regularization to the original Grad's moment system, for arbitrary order and for any dimension [Cai-Fan-Li, Comm. Pure Appl. Math. (2014)]. The regularized systems enjoy the following advantages: (1) the systems are globally hyperbolic and symmetric; (2) the spectral approximation accuracy of Grad's expansion is perfectly preserved; (3) the structure of the characteristic fields, including the wave speeds and the convectors, are fully clarified. The regularized moment system, as a symmetric hyperbolic quasi-linear PDE system, is then locally well posed in time. For 1D case, it can be shown that the regularization they proposed is unique under certain well-understood constraints [Cai-Fan-Li, Commun. Math. Sci. (2012)]. This essentially solved the problem on the loss of hyperbolicity of Grad's moment system. In [Cai-Fan-Li, SIAM J. Appl. Math. (2015)], Li et al. further revealed why Grad's moment system may lose its hyperbolicity, and developed a framework to carry out hyperbolic model reduction for generic kinetic equations.

3) Numerical Regularized Moment System for Arbitrary OrderNumerical Regularized Moment System for Arbitrary Order

Li et al. developed the numerical scheme to systematically solve the regularized moment system. After deriving the reduced model with local blessedness in time, it is a formidable task to develop corresponding numerical schemes and codes, due to the possible huge number of equations in the system. In their new method, the order of the numerical discretisation and the moment expansion in the canonical procedures is exchanged [Cai-Li, SIAM J. Sci. Comput. (2010)]. By resolving several technical difficulties, Li et al. successfully gave the so-called NRxx method, which may solve the Boltzmann equation using moment method for arbitrary order, with the same cost on code developing. This makes it possible to apply the regularized moment method to different applications, including polyatomic gases [Cai-Li, J. Comput. Phys. (2014)], plasma [Cai-Li-Wang, SIAM J. Sci. Comput. (2013)], electron transportation [Li-Lu-Wang-Yao, Commun. Comput. Phys. (2014)], micro-flows [Cai-Li-Qiao, Comput. & Fluids (2013); Cai-Li-Qiao, SIAM J. Sci. Comput. (2012)] and so on.

3. Moment Model of Arbitrary Order for Special Relativistic Boltzmann Equation (Huazhong Tang)

It is difficult to derive the relativistic moment system of higher order since the family of orthogonal polynomials can not be found easily. Several authors have tried to construct the family of orthogonal polynomials analogous to the Hermite polynomials. Unfortunately, there is no explicit expression of the moment systems if the order of the moment system is larger than 3. Moreover, the hyperbolicity of existing general moment systems is not proved, even for the second order moment system (e.g. the general Israel and Stewart system). Huazhong Tang and his colleagues derived arbitrary order globally hyperbolic moment system for the one-dimensional special relativistic Boltzmann equation and studied some properties of the moment system: the eigenvalues and their bound as well as eigenvectors, hyperbolicity, characteristic fields, linear stability, and Lorentz covariance [Kuang-Tang, J. Stat. Phys. (2017)]. Their work was built on their careful study of two families of the complicate Grad type orthogonal polynomials depending on a parameter.

4. Mixed Finite Element Method of Elasticity Problems and Fast Iterative Methods (HU Jun)

The problems that are most frequently solved in scientific and engineering computing may probably be the elasticity equations. The finite element method (FEM) was invented in analyzing the stress of the elastic structures in the 1950s. The mixed FEM within the Hellinger-Reissner (H-R) principle for elasticity yields a direct stress approximation since it takes both the stress and displacement as an independent variable; while the displacement FEM only gives an indirect stress approximation.  Hu, a founder of the celebrated Hu-Washizu principle for elasticity, pointed out that, the H-R principle is more general than both the minimum potential and complementary energy principles, and is much fitter for numerical solutions [Hu,Acta Phys. Sin.(1954)]. Indeed, the mixed FEM can be free of locking for nearly incompressible materials, and be applied to plastic materials, and approximate both the equilibrium and traction boundary conditions more accurate. However, the symmetry of the stress plus the stability conditions make the design of the mixed FEM for elasticity surprisingly hard.

Since the 1960s, many mathematicians have worked on this problem, which leads to weakly symmetric elements [Arnold-Brezzi-Douglas,Japan J. Appl. Math.(1984); Arnold-Falk-Winther,Math. Comp.(2007); Boffi-Brezzi-Fortin,Commun. Pure Appl. Anal.(2009)], and composite elements [Johnson-Mercier,Numer. Math.(1978)]. In 2002, using the elasticity complexes, Arnold and Winther designed the first family of symmetric mixed elements with polynomial shape functions in 2D which was extended to 3D [Arnold, Awanou and Winther,Math. Comp.2008].

1) Mixed Finite Element Methods of Elasticity Problems

Jun Hu developed a new framework to design and analyze the mixed FEM of elasticity problems by establishing the following three main results:

A structure of the discrete stress space: on simplicial grids, the discrete stress space can be selected as the symmetric matrix-valued Lagrange element space, enriched by a symmetric matrix-valued polynomial H(div) bubble function space on each simplex; a corresponding choice applies to product grids.

Two basic algebraic results: (1) on each simplex, the symmetric matrices of rank one produced by the tensor products of the unit tangent vectors of the n(n+1)/2 edges of the simplex, form a basis of the space of the symmetric matrices in n dimensions;   (2) on each simplex, the divergence space of the above H(div) bubble function space is equal to the orthogonal complement space of the rigid motion space with respect to the corresponding discrete displacement space (A similar result holds on a macroelement for the product grids).

These define a two-step stability analysis. As a result, on both simplicial and product grids, Hu et al. were able to define the first families of both symmetric and optimal mixed elements with polynomial shape functions in any space dimension. Further more, the discrete stress space has a simple basis which essentially consists of symmetric matrix-valued Lagrange element basis functions. See [Hu-Man-Zhang,J. Sci. Comput.(2014); Hu-Zhang,Sci. China Math.(2015); Hu,J. Comput. Math.(2015); Hu,SIAM J. Numer. Anal.(2015); Hu-Zhang,Math. Models Methods Appl. Sci.(2016); Hu-Man-Wang-Zhang,Comput. Math. Appl.(2016); Chen-Hu-Huang,Comput. Methods Appl. Math.(2017)] for more details.

2) Fast Auxiliary Space Preconditioner for Linear Elasticity in Mixed Form

For most of conforming mixed finite elements of linear elasticity, they use vertex degrees of freedom. Hence, it is difficult to employ the usual hybridization technique for them. To avoid such a difficulty, by combing the mesh dependent norm  technique (namely, equipping  the discrete stress with the L2norm, the discrete displacement with the discrete H1norm), the block diagonal (triangular)  preconditioner technique,  the auxiliary space preconditioner  technique, and the (generalized) minimal residual method, Hu et al. proposed a fast auxiliary space  preconditioner for mixed finite elements of linear elasticity and proved  the uniform convergence with respect to both the meshsize and crucial Lame constant. As a result, one only needs in some sense to solve a symmetric and positive definite system of the H1conforming linear element of the linear elasticity problem [Chen-Hu-Huang,arXiv(2016)].

5. Numerical Methodsin Fluid Mechanics (TANG Huazhong, LI Tiejun)

We studied monotone schemes of hyperbolic conservation laws, and physical- constraints-preserving schemes of both relativistic hydrodynamics and relativistic magnetohydrodynamics equations.

1) Local Oscillations in Finite Difference Solutions of Hyperbolic Conservation Laws

It was generally expected that monotone schemes are oscillation-free for hyperbolic conservation laws. However, recently local oscillations were observed and usually understood to be caused by relative phase errors. In order to explain this, they used the discrete Fourier analysis and the modified equation analysis to distinguish the dissipative and dispersive effects of numerical schemes for low frequency and high frequency modes, and revealed that in order to avoid numerical oscillations, the relative phase errors should be offset by numerical dissipation of at least the same order. Numerical damping, i.e. the zero order term in the corresponding modified equation, is important to dissipate the oscillations caused by the relative phase errors of high frequency modes. This is in contrast to the role of numerical viscosity, the second order term, which is the lowest order term usually present to suppress the relative phase errors of low frequency modes. More details can be found in in [Li-Tang-Warnecke-Zhang, Math. Comp.(2009)].

2) Physical-Constraints-Preserving Schemes of both Relativistic Hydrodynamics and Relativistic Magnetohydrodynamics Equations

Relativistic hydrodynamics (RHD) is important in investigating numerous astrophysical phenomena, from stellar to galactic scales; e.g., gamma-ray bursts, astrophysical jets, core collapse super-novae, coalescing neutron stars, formation of black holes, etc. The RHD equations are highly nonlinear, making their analytical treatment extremely difficult. Tang et al. proved the convexity and other properties of the admissible state set and discover a concave function with respect to the conservative vector replacing the pressure which is an important ingredient to enforce the positivity-preserving property for the non-relativistic case [Wu-Tang, J. Comput. Phys.(2015); Wu-Tang,Astrophys. J. Suppl. Ser.(2017)]. They also studied the admissible state setGof relativistic magnetohydrodynamics (RMHD). To overcome the difficulties arising from the extremely strong nonlinearities and no explicit formulas of the primitive variables and flux vectors with respect to the conservative vector, two equivalent forms ofGwith explicit constraints on the conservative vector were skillfully discovered. The first was derived by analyzing roots of several polynomials and transferring successively them, and further used to prove the convexity ofGwith the aid of semi-positive definiteness of the second fundamental form of a hypersurface. While the second was derived based on the convexity and then used to show the orthogonal invariance ofG. The Lax-Friedrichs (LxF) splitting property does not hold generally for nonzero magnetic field, but by a constructive inequality and pivotal techniques, we discovered the generalized LxF splitting properties, combining the convex combination of some LxF splitting terms with a discrete divergence-free condition of magnetic field [Accepted byMath. Model. Meth. Appl. Sci.(2017)]. Based on these theoretical results, they developed physical-constraints-preserving (PCP) schemes for both special RHD equations with an ideal or a general equation of state and one- and two-dimensional RMHD equations. Our analysis revealed in theory for the first time that the discrete divergence-free condition was closely connected with the PCP property.

3) Rare Event Study for Chemical Reaction Kinetics

In developmental biology, it is believed that the change of phenotype of cells can be described through a dynamics driven by a “potential energy landscape”. The famous concept “epigenetic landscape” was even proposed by the British biologist Waddington in the 1950s. Though intuitive and influential, the rationale of the epigenetic landscape is missing and a mathematical theory about the energy landscape is desired in biophysics community. Based on the quasi-potential concept from Freidlin-Wentzell large deviation theory (LDT), Li et al. propose to rationalize the Waddington landscape by constructing the global quasi-potential energy landscape for a biological system. Furthermore, the realization of this construction to typical problems motivates some interesting theory. For genetic switching models involving positive feedback in prokaryotic cells, we encounter the problem to establish the two- scale large deviations for Gillespie-type jump process. In this model, the DNA switches between active and inactive states very fast, and once it is in active state, the transcription and translation process will be started. The classical large volume limit does not hold in this case and we need to combine the Freidlin-Wentzell type LDT for the translation process and the Donsker-Varadhan type LDT for the DNA switching process. With this LDT, they get an explicit formula of the Hamiltonian and find that the strict correspondence between the drift and diffusion terms in the classical chemical Langevin approximation is lost. This enriches the classical theory about the chemical reaction kinetics [Lv-Li-Li-Li, PLoS ONE(2014); Lv-Li-Li-Li, PLoS Comp. Biol.(2015); Li-Lin, Commun. Math Sci.(2017)].

6. Numerical methods in Quantum Mechanics (Sihong Shao, Tiao Lu)

Sihong Shao and Tiao Lu developed numerical methods for the Wigner equation and the Dirac equation both of which serve as the fundamental mathematical models for quantum mechanics.

1) Efficient Algorithms for Many-body Wigner Quantum Dynamics

Distinct from the Schrödinger wave function, the Wigner function shares many analogies to the classical mechanism and simply reduces to the classical counterpart as the Planck constant vanishes. However, numerical resolutions for the Wigner equation remain one of the most challenging problems in computational physics due to the high dimensionality and the highly oscillatory structure. Shao and Lu completed the design and implementation of an accurate numerical scheme for the Wigner quantum dynamics in 4D phase space, combining an efficient conservative semi-Lagrangian scheme in the temporal-spatial space with an accurate spectral element method in the momentum space [Shao-Lu-Cai, Commun. Comput. Phys. (2011); Xiong-Chen-Shao, SIAM J. Sci. Comput. (2016)]. With an auxiliary function, they cast the Wigner equation into a renewal-type integral equation and proved that its solution is equivalent to the first moment of a stochastic branching random walk and thus revealed the inherent relation between the Wigner equation and a stochastic branching random walk model [Shao-Sellier, J. Comput. Phys. (2015); Shao-Xiong, arXiv (2016, 2017)]. A device-adaptive inflow boundary condition in simulating a resonant tunneling diode was proposed to save the computational cost [Jiang-Lu-Cai, J. Comput. Phys. (2014)]. They proved without any additional prerequisite conditions that the solution of the Wigner equation with inflow boundary conditions will be symmetric only if the potential is symmetric [Li-Lu-Sun, SIAM J. Appl. Math. (2014)]. This improves the result in [Taj-Genovese-Rossi, Europhys. Lett. (2006)]. By numerical studies, they present the convergence of the numerical solution to the symmetric profile for three different numerical schemes. This implies that the upwind schemes can also yield a symmetric numerical solution, contrary to the argument given in Taj-Genovese-Rossi’s work.

2) Mathematical Analysis and Numerical Methods for Dirac Equation

The spectrum of 1-electron Dirac operator is well understood but little is known on the spectrum of N-electron Dirac operator. Even for the simplest 2-electron Dirac-Coulomb operator, the question whether it has bound electronic states is still open. Shao and his colleagues took the first step in this direction and derived the relativistic electron-electron coalescence conditions by virtue of matrix rings and Tracy-Singh products. By making use of these asymptotic behaviors, a further asymptotic analysis implies strongly that the Dirac-Coulomb Hamiltonian has no bound electronic states [Li-Shao-Liu, J. Chem. Phys. (2012)]. For the nonlinear Dirac solitary waves, the stability properties have remained a complete mystery. Shao and his colleagues studied their multi-hump structure and explored its relation with the stability in both theory and numeric. Their results inferred that both the multi-hump profile and high order nonlinearity could undermine the stability during the scattering [Shao-Tang, Phys. Lett. A (2005); Xu-Shao-Tang-Wei, Commun. Math. Sci. (2015); Xu-Shao-Tang, J. Comput. Phys. (2013); Shao et al., Phys. Rev. E (2014)].