Automated Material Model Discovery

• Automated material model discovery is a data-driven framework that identifies constitutive laws and thermodynamic potentials from experimental or simulated data while ensuring physical admissibility.
• It leverages techniques such as sparse regression, symbolic generation, and physics-augmented neural networks to balance predictive accuracy, sparsity, and interpretability.
• The approach integrates model selection, parameter identification, and uncertainty quantification to address complex material behaviors and diverse loading conditions.

Searching arXiv for papers on automated material model discovery and related constitutive discovery frameworks. Automated material model discovery denotes the data-driven identification of constitutive laws, free-energy densities, dissipation potentials, or related predictive material relations directly from experiments or simulations while preserving interpretability and, in many formulations, embedding mechanical admissibility by design. In the recent literature, the topic spans sparse regression over curated feature libraries, constitutive neural networks that output thermodynamic potentials, grammar-based and genetic symbolic regression, unsupervised full-field inverse methods based on equilibrium, and differentiable finite-element model updating for history-dependent materials. Across these variants, the common objective is simultaneous model selection and parameter identification under explicit trade-offs between predictive accuracy, sparsity, physical consistency, and computational tractability (McCulloch et al., 2023, Flaschel et al., 2022, Ferreira et al., 12 May 2025).

1. Emergence of the field and its problem classes

Automated constitutive discovery has been framed as an inverse problem in which one seeks a constitutive law that maps deformation or strain histories to stresses in an interpretable and physically admissible form. For hyperelasticity, the target is typically a strain-energy density $W$ or $\psi$ from which stresses are obtained by differentiation; for inelasticity, the target extends to elastic and inelastic potentials together with evolution laws for internal variables (Linka et al., 2022, Ji et al., 19 Feb 2026).

Several distinct problem classes now coexist. A first class is local/direct discovery, in which strain–stress pairs or stress–deformation curves are available, as in supervised hyperelastic characterization of brain tissue from uniaxial and torsional data (Flaschel et al., 2023) or local history-dependent fitting in ADiMU (Ferreira et al., 12 May 2025). A second class is global/indirect discovery, in which stresses are not observed and the constitutive law is inferred from full-field displacement measurements and net reaction forces by enforcing equilibrium in weak form, as in EUCLID for plasticity and generalized standard materials (Flaschel et al., 2022, Flaschel et al., 2022). A third class is symbolic or grammar-based discovery, where the search variable is the mathematical expression itself rather than only coefficients in a fixed library (Kissas et al., 2024, Hou et al., 2024). A fourth class is physics-augmented neural discovery, where trainable architectures encode kinematics, thermodynamics, convexity, or monotonicity and can subsequently be sparsified or symbolified (Linka et al., 2022, Jadoon et al., 2024, Ji et al., 19 Feb 2026).

This heterogeneity reflects the fact that the constitutive search space is strongly problem-dependent. Hyperelastic discovery often reduces to selecting functions of invariants or principal stretches, whereas elastoplastic and viscoelastic discovery must also identify internal variables, rate equations, and dissipation mechanisms. A plausible implication is that “automated material model discovery” is best understood not as a single algorithmic family but as a family of inverse constitutive identification strategies organized by data modality, physics assumptions, and representational bias.

2. Physical structure and admissibility constraints

A defining feature of the field is that successful methods do not treat constitutive modeling as unconstrained function approximation. They instead encode objectivity, symmetry, incompressibility or compressibility structure, and thermodynamic consistency either a priori in the representation or a posteriori through admissibility checks.

For hyperelasticity, the usual kinematic basis is the deformation gradient $F$, the right Cauchy–Green tensor $C = F^\top F$, the Jacobian $J = \det F$, and the invariants

$$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$

Energy-based formulations then define stresses by differentiation, for example

$$P = \frac{\partial W}{\partial F}, \qquad \sigma = \frac{1}{J} F \left(\frac{\partial W}{\partial F}\right)^\top,$$

or, in incompressible settings,

$P = \frac{\partial W}{\partial F} - p F^{-T},$

with $p$ acting as a Lagrange multiplier enforcing $J=1$ (McCulloch et al., 2023, Linka et al., 2022).

These constructions are not merely formal. They are used to guarantee that learned stresses remain derivable from a potential and therefore satisfy the Clausius–Duhem inequality for purely elastic processes when the potential is well posed. In constitutive neural networks, this is achieved by making the network output a scalar free energy rather than stress components directly, with inputs chosen as invariants or principal stretches so that frame indifference and isotropy are enforced by construction (Linka et al., 2022). In grammar-based hyperelastic discovery, objectivity and isotropy are embedded by generating only expressions in the invariants of $\psi$0 and by appending normalization corrections so that $\psi$1 and stress is near zero at the identity (Kissas et al., 2024). In symbolic regression for brain cortex, invariant-based models are constructed as sums of convex functions in invariants, while stretch- and strain-based models undergo posterior Hessian or ellipticity checks (Hou et al., 2024).

For inelastic materials, admissibility extends to dissipation. iCKAN formulates finite-strain inelasticity with the multiplicative decomposition $\psi$2, dual elastic and inelastic potentials, and a convexification operator applied to the learned inelastic potential so that the Clausius–Plank inequality reduces to a non-negative form $\psi$3 (Ji et al., 19 Feb 2026). Finite-strain elastoplastic PANNs similarly define free-energy and yield potentials so that each term entering the reduced dissipation is non-negative by construction (Jadoon et al., 2024). EUCLID generalizes this logic to generalized standard materials by representing behavior through a convex Helmholtz free energy $\psi$4 and a convex dissipation potential $\psi$5 or $\psi$6, using convexity to guarantee stability and thermodynamic consistency (Flaschel et al., 2022).

This emphasis on hard constraints narrows the admissible hypothesis space. The literature repeatedly argues that such restriction is not a limitation but a precondition for robust extrapolation, sparse term selection, and physically interpretable parameter recovery (McCulloch et al., 2023, Linka et al., 2022).

3. Representations used for discovery

The main representational choices can be organized by the object being discovered: coefficients in a predefined library, symbolic expressions generated from a grammar or genetic program, or trainable neural potentials that are later sparsified or symbolified.

Representation	Core mechanism	Representative papers
Sparse library regression	Select a small subset of predefined constitutive features and fit coefficients	(McCulloch et al., 2023, Flaschel et al., 2023, Flaschel et al., 14 Jul 2025)
Symbolic generation	Generate candidate energy expressions by grammar or genetic programming	(Kissas et al., 2024, Hou et al., 2024)
Physics-augmented neural potentials	Learn $\psi$7, $\psi$8, $\psi$9, or hardening potentials with architectural constraints	(Linka et al., 2022, Jadoon et al., 2024, Ji et al., 19 Feb 2026, Ferreira et al., 12 May 2025)
Full-field unsupervised inverse discovery	Infer constitutive law from equilibrium residuals without stress labels	(Flaschel et al., 2022, Flaschel et al., 2022)

In sparse library approaches, the dictionary is usually built from classical constitutive motifs. Constitutive neural networks in one influential formulation use eight functional building blocks in $F$0 and $F$1, including linear, quadratic, and exponential terms, so that neo-Hookean, Blatz–Ko, Mooney–Rivlin, Yeoh, and Demiray-type models become special cases of the same architecture (Linka et al., 2022). Best-in-class modeling extends this idea to a library of sixteen building blocks, adding anisotropic terms in $F$2 and $F$3 and then performing bottom-up densification rather than top-down pruning (Linka et al., 2024). Supervised EUCLID for brain tissue constructs an invariant library with generalized Mooney–Rivlin monomials and a $F$4 feature, together with a dense generalized Ogden stretch library using fixed exponents $F$5 (Flaschel et al., 2023).

Symbolic approaches broaden the hypothesis class. Formal grammars generate large libraries of admissible hyperelastic laws in Polish notation while biasing toward objectivity, isotropy, normalization, and growth properties (Kissas et al., 2024). Symbolic regression with genetic programming searches over invariant-based, principal-stretch-based, and normal-strain-based formulas and selects compact expressions by combining normalized fit error with expression-complexity control (Hou et al., 2024). iCKAN occupies an intermediate position: the learned functions are represented by trainable B-splines inside a Kolmogorov–Arnold architecture, then converted to closed-form symbolic expressions by a post-training symbolification step (Ji et al., 19 Feb 2026).

A further representational distinction concerns linear versus nonlinear dependence on parameters. Principal-stretch Ogden-type networks with fixed exponents are linear in the weights and therefore reduce to convex least-squares problems with unique global minima for the chosen feature set (McCulloch et al., 2023). In contrast, invariant-based networks with exponentials or free exponents induce non-convex optimization landscapes with multiple local minima (McCulloch et al., 2023, Linka et al., 2024). Much of the recent optimization literature is a response to this specific distinction.

4. Optimization, sparsity, and model selection

Sparsity control is central because constitutive models are expected to remain interpretable, robust in extrapolation, and computationally deployable. A standard objective is

$F$6

with $F$7 for ridge, $F$8 for lasso, $F$9 for explicit cardinality, and $C = F^\top F$0 for fractional penalties (McCulloch et al., 2023).

A recurring conclusion is that the choice of penalty qualitatively changes the discovery outcome. In constitutive neural network studies, $C = F^\top F$1 or ridge regularization is reported as unsuitable for model discovery because it stabilizes coefficients without inducing exact zeros; $C = F^\top F$2 or lasso promotes sparsity but introduces strong bias; and $C = F^\top F$3 provides the most transparent handle on the trade-off between interpretability and predictability, simplicity and accuracy, and bias and variance (McCulloch et al., 2023). The same emphasis on non-smooth sparsity motivates later work devoted specifically to algorithms for minimizing objectives of the form $C = F^\top F$4 in constitutive discovery (Flaschel et al., 14 Jul 2025).

For quadratic mismatch functions, coordinate descent implements the classical LASSO, while LARS computes the entire piecewise-linear regularization path and identifies the critical values of $C = F^\top F$5 at which coefficients enter or leave the active set (Flaschel et al., 14 Jul 2025). For non-quadratic or non-convex constitutive mappings, ISTA and pathwise ISTA are proposed as practical proximal-gradient schemes for obtaining approximate regularization paths (Flaschel et al., 14 Jul 2025). These methods are particularly relevant when the model outputs depend nonlinearly on weights through exponentials or chain-rule stress maps.

Non-convexity has also led to a shift from top-down sparsification to bottom-up model growth. Best-in-class modeling starts from the best single-term model, then iteratively adds the term that most decreases the objective, thereby converting a combinatorial search such as $C = F^\top F$6 possible term combinations into a sequence of smaller nonlinear fits (Linka et al., 2024). Closely related guidance appears in the constitutive-network study, which recommends exact subset enumeration for one-term and two-term models and a bottom-up “densify” strategy for larger libraries (McCulloch et al., 2023).

Normalization is another technical theme. In multi-test hyperelastic fitting, stress residuals are normalized by the maximum recorded stress in each loading modality to prevent one mode from dominating the objective (McCulloch et al., 2023). Supervised EUCLID for brain tissue similarly scales and concatenates uniaxial and torsional regression systems with empirically chosen weights $C = F^\top F$7 and $C = F^\top F$8 to equalize signal magnitudes before applying non-negative Lasso and Pareto-based model selection (Flaschel et al., 2023).

For full-field and history-dependent problems, optimization is inseparable from differentiable state updates and PDE solves. ADiMU backpropagates through return mapping, Newton iterations, and vectorized finite-element assembly, thereby enabling model updating for conventional, hybrid, and neural models without introducing extra hyperparameters beyond those intrinsic to the selected architecture and optimizer (Ferreira et al., 12 May 2025). This suggests that automated material model discovery is increasingly converging with differentiable scientific computing rather than remaining a standalone sparse-regression problem.

5. Data modalities, benchmark domains, and representative discoveries

The field now covers a broad spectrum of datasets. Synthetic stress–deformation data remain the standard vehicle for verification because they allow exact benchmarking against known ground-truth models and parameters (McCulloch et al., 2023, Hou et al., 2024). Experimental applications include human brain gray and white matter, human brain cortex, VHB viscoelastic polymers, rubber, porcine skin, arteries, mild steel under cyclic loading, and soft-matter systems such as artificial meat (Flaschel et al., 2023, Ji et al., 19 Feb 2026, Linka et al., 2024, Jadoon et al., 2024).

Several representative discoveries recur across studies. In human brain tissue, multiple approaches converge on strong $C = F^\top F$9 dependence. Constitutive-network experiments on real brain data report that the best unregularized fit in the screened window often selects the Blatz–Ko term with $J = \det F$0 and $J = \det F$1, while the best one-term and two-term $J = \det F$2 models are quadratic or exponential functions of $J = \det F$3 (McCulloch et al., 2023). Best-in-class modeling likewise identifies one-term and two-term models centered on $J = \det F$4 and $J = \det F$5 for gray and white matter (Linka et al., 2024). Supervised EUCLID applied to 81 human-brain specimens finds that the most frequently discovered model is one-term Ogden, followed by two-term Ogden and mixed Ogden-plus-$J = \det F$6 forms, with average standardized mean squared error about $J = \det F$7 across specimens (Flaschel et al., 2023). Symbolic regression on human brain cortex reaches a related conclusion in the invariant representation: all discovered invariant-based models are functions of $J = \det F$8 only, including an optimal compact form

$J = \det F$9

(Hou et al., 2024).

For viscoelastic and thermo-viscoelastic polymers, iCKAN demonstrates that symbolic inelastic discovery can recover explicit elastic and inelastic potentials from sequential finite-strain data. On VHB 4910, a two-branch Maxwell-like architecture discovers exponential-plus-polynomial elastic potentials and convexified inelastic potentials with training and test NMSE of approximately $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$0 and $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$1, while the symbolified versions retain similar accuracy (Ji et al., 19 Feb 2026). On VHB 4905, temperature dependence is learned as explicit cubic or piecewise-quadratic functions $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$2 inside the elastic potentials, with training and test NMSE approximately $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$3 and $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$4 (Ji et al., 19 Feb 2026).

For cyclic metal plasticity at finite strain, physics-augmented neural networks discover hardening potentials that outperform classical Armstrong–Frederick and Ohno–Wang parameter fitting on both synthetic and experimental data. In the reported experimental mild-steel case, the 4NN model achieves best, mean, and standard-deviation losses of $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$5, $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$6, and $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$7, compared with $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$8, $$I_1 = \operatorname{tr}(C), \qquad I_2 = \tfrac{1}{2}\big[(\operatorname{tr} C)^2 - \operatorname{tr}(C^2)\big], \qquad I_3 = \det(C) = J^2.$$9, and $$P = \frac{\partial W}{\partial F}, \qquad \sigma = \frac{1}{J} F \left(\frac{\partial W}{\partial F}\right)^\top,$$0 for Armstrong–Frederick parameter fitting (Jadoon et al., 2024). ADiMU further generalizes this differentiable-update paradigm to local stress–strain discovery and global full-field discovery for conventional, hybrid, and neural constitutive models, including elasto-plasticity and GRU-based surrogates with millions of parameters (Ferreira et al., 12 May 2025).

At the laboratory interface, not all automation concerns constitutive laws themselves. The AMT GUI exemplifies a neighboring workflow in which predictive surrogates and particle swarm optimization are used to recommend the next experiments from materials datasets without requiring programming expertise (Shakeel et al., 2023). This is not constitutive discovery in the narrow mechanics sense, but it illustrates a broader trend: materials-model automation increasingly links model fitting, experiment design, and closed-loop data acquisition.

6. Limitations, controversies, and likely directions

The literature is explicit about unresolved limitations. One is combinatorial complexity. Exact $$P = \frac{\partial W}{\partial F}, \qquad \sigma = \frac{1}{J} F \left(\frac{\partial W}{\partial F}\right)^\top,$$1 subset search is tractable only for small libraries or very low cardinalities; larger spaces require heuristics such as forward densification, clustering, or latent-space search (McCulloch et al., 2023, Linka et al., 2024). A second is non-convexity: invariant-based exponential libraries, fractional penalties, neural hardening laws, and symbolic-generation pipelines all admit multiple local minima, so global optimality is typically not guaranteed (McCulloch et al., 2023, Flaschel et al., 14 Jul 2025, Kissas et al., 2024).

A third limitation is identifiability under restricted loading. Multiple papers note that limited loading modes cannot identify all parameters; lack of shear can preclude inference of shear-related terms, and moderate-strain training data may be insufficient to excite the nonlinear structure of Ogden-type models, causing the discovery algorithm to settle on accurate surrogates rather than the literal ground-truth form (McCulloch et al., 2023, Kissas et al., 2024). This also underlies some apparent controversies in the literature, such as whether brain tissue is best represented by invariant-based or stretch-based laws. The published evidence does not establish a universal winner; rather, it shows that different representations can fit the same dataset well while emphasizing different aspects of asymmetry, convexity, or extrapolation (Flaschel et al., 2023, Hou et al., 2024).

A fourth issue concerns convexity, ellipticity, and extrapolation. Several methods bias toward polyconvexity or convexity but do not guarantee it globally (McCulloch et al., 2023, Kissas et al., 2024). Stretch-based symbolic models for brain cortex can lose convexity under large deformations even when they pass training-range Hessian checks (Hou et al., 2024). Best-in-class modeling explicitly notes that convexity and polyconvexity are not enforced (Linka et al., 2024). This suggests that physical admissibility remains a spectrum: some frameworks guarantee it by architecture, while others rely on curated libraries and posterior checks.

Future directions described across the papers are comparatively aligned. They include extending hyperelastic discovery to viscoelasticity, plasticity, damage, growth, remodeling, anisotropy, orthotropy, and temperature-dependent behavior (Ji et al., 19 Feb 2026, Jadoon et al., 2024, Flaschel et al., 2022). They also include uncertainty quantification through Bayesian or probabilistic regularization, stronger convexity guarantees through attribute grammars or convex architectures, and tighter integration with differentiable finite-element solvers and inverse design loops (McCulloch et al., 2023, Kissas et al., 2024, Ferreira et al., 12 May 2025). A plausible implication is that the field is moving toward unified pipelines in which constitutive structure, experiment design, and downstream simulation are co-optimized rather than treated as separate stages.

Automated material model discovery therefore sits at the intersection of continuum mechanics, sparse optimization, symbolic computation, and differentiable programming. Its distinctive claim is not merely that models can be fit automatically, but that the form of the constitutive law itself can be selected, constrained, and interpreted from data in a way that remains compatible with the mathematical and thermodynamic structure of materials theory (Linka et al., 2022, Ji et al., 19 Feb 2026).