††thanks: LLNL-JRNL-835975

# Constraining Model Uncertainty in Plasma Equation-of-State Models with a
Physics-Constrained Gaussian Process

Jim A. Gaffney<EMAIL_ADDRESS>Lin Yang  Suzanne Ali  Lawrence Livermore
National Laboratory, 7000 East Ave, Livermore, CA 94550, USA

###### Abstract

Equation-of-state (EOS) models underpin numerical simulations at the core of
research in high energy density physics, inertial confinement fusion,
laboratory astrophysics, and elsewhere. In these applications EOS models are
needed that span ranges of thermodynamic variables that far exceed the ranges
where data are available, making uncertainty quantification (UQ) of EOS models
a significant concern. Model uncertainty, arising from the choice of
functional form assumed for the EOS, is a major challenge to UQ studies for
EOS that is usually neglected in favor of parameteric and data uncertainties
which are easier to capture without violating the physical constraints on
EOSs. In this work we introduce a new statistical EOS construction that
naturally captures model uncertainty while automatically obeying the
thermodynamic consistency constraint. We apply the model to existing data for
$B_{4}C$ to place an upper bound on the uncertainty in the EOS and Hugoniot,
and show that the neglect of thermodynamic constraints overestimates the
uncertainty by factors of several when data are available and underestimates
when extrapolating to regions where they are not. We discuss extensions to
this approach, and the role of GP-based models in accelerating simulation and
experimental studies, defining portable uncertainty-aware EOS tables, and
enabling uncertainty-aware downstream tasks.

##  I Introduction

Understanding high energy density physics, inertial confinement fusion and
laboratory astrophysics experiments relies heavily on equation-of-state (EOS)
models due to their importance in closing the equations of
hydrodynamicsGaffney  _et al._ (2018). In these applications different
portions of the system pass through a range of thermodynamic states that
greatly exceeds that accessible by focussed EOS measurements. As a result EOS
models routinely extrapolate away from the data used to benchmark them, making
model validation and uncertainty quantification a significant concern.

Typical EOS models are built by collecting a patchwork of experimental data
and first-principles simulations of different types and fitting a parametric
function that can be used to interpolate and extrapolate over a wide range of
input statesMilitzer (2021); Zhang (2020). This parametric approach is
attractive since it allows important physical constraints and limiting
behaviors to be trivially enforced by picking a physically-motivated
functional formWu  _et al._ (2019). Uncertainty in the resulting EOS can be
measured by statistical sampling of the fitted parameters, and propagating to
the quantities of interest. This approach captures the data component of the
total EOS uncertainty but misses another: the uncertainty in the choice of
functional form used, which we call _model uncertainty_. In this paper we
report on work to constrain model uncertainty in equation-of-state models
taking into account the physics-driven constraints on the model form.

Constraining model uncertainty involves exploring the potentially infinite
space of viable functional forms. This can be approximated by considering an
ensemble of parametric functions and fitting each individually, or by fitting
a highly flexible non-parametric model that implicitly represents an entire
class of functions. The first approach has the advantage that every member of
the ensemble can be a physically viable form, but risks underestimating the
uncertainty since it is unlikely to consider all possible functions; the
second approach allows for a much more complete quantification of the
uncertainty but efforts must be made to ensure the space of functions explored
is physically reasonable.

In this work we develop a non-parametric fitting procedure suitable for EOS
applications. Our approach is based on a variation of the well-known Gaussian
process (GP) regression MacKay  _et al._ (1998), where we introduce
modifications to ensure that all functions generated by the GP obey the laws
of thermodynamics and that the known ideal behavior at high temperature is
recovered. We demonstrate our model by applying it to published data for Boron
Carbide ($B_{4}C$), and use it to explore the uncertainty in pressure $P$ and
internal energy $E$ when interpolating first-principles simulation data.

Uncertainty quantification of materials properties is being actively
researched by several groups. Recently, Ali _et al._ Ali  _et al._ (2020)
demonstrated a full quantification of parametric uncertainty in the EOS of
copper by Monte-Carlo sampling data uncertainties and propegating through the
EOS fitting process. Kamga _et al._ Kamga  _et al._ (2014) have quantified the
uncertainty in a single model by exploring the discrepancy with experiment.
Root _et al._ Root (2010), Walters _et al._ Walters  _et al._ (2018), and
Brown _et al._ Brown and Hund (2018) have used Gaussian process based model
calibration Kennedy and O’Hagan (2001); Higdon  _et al._ (2008) to update
physics models to better agree with experiment. While important steps, none of
these approaches fully acount for the range of models that are consistent with
experimental data, or constrain the calibrated model to be consistent with
fundamental physics. This is the main contribution of our work.

##  II Model Uncertainty in Equations-of-State

A general EOS model takes a pair of natural variables as input and predicts
all other thermodynamic variables of the system as outputs, taking into
account phase changes, molecular and atomic physics, coupling with a radiation
field, etc. In a hydrodynamics simulation, the input variables are those
evolved on the computational mesh, typically the temperature $T$ and mass
density $\rho$ or volume $V$; the outputs are pressure, internal energy,
entropy, specific heat, and any other material property of interest. In this
way the EOS summarizes the microscopic response of the material with the aim
of allowing the complex sub-grid physics to be modeled separately from the
hydrocode itself.

The usual process of building an EOS model is to use as much data as possible
to learn a mapping from the input space
$\mathbf{x}\in\mathbb{R}_{>0}\times\mathbb{R}_{>0}$ into the outputs
$\mathbf{y}\in\mathbb{R}^{N}$. The data are usually a combination of
experimental results from a variety of facilities and platforms, sometimes
combined with results from multiple first-principles simulation studies with
varying fidelities and regions of applicability. The resulting model may be
written

$EOS:\mathbf{x},\mathbf{\alpha}\mapsto\mathbf{y}~{}{\rm,}$ (1)

where $\mathbf{\alpha}$ is a vector of parameters whose values are determined
by fitting to the available data over the entire range of the EOS.

Quantifying our confidence in the EOS in (1) involves admitting that, for a
given input points $x$, there may be a range of outputs which are consistent
with the available data. The variations come from two sources: there may be
multiple possible mappings (model uncertainty); and for a fixed mapping there
may be multiple sets of parameters $\alpha$ (parametric uncertainty). These
two contributions may contribute equally to the total uncertainty in the
predicted quantities, however capturing model uncertainty is more challenging
since it involves enumerating all of the potential mappings (functional forms)
that could be valid EOSs. As a result of this difficulty, most previous work
has considered only parametric uncertainty and can be expected to
underestimate the total uncertainty in the EOS model.

In the remainder of this sections we will present a method for automatically
exploring a continous space of EOS models to generate predictions that capture
all sources of uncertainty. We being by defining the function space in terms
of basic physical constraints, and then present a Gaussian process method for
sampling general EOSs conditional on data.

###  II.1 Thermodynamic Constraints on EOS Functional Forms

The outputs of a physically realistic EOS model are strongly constrained by
fundamental thermodynamics. On such constraint, known as _thermodynamic
consistency_ , arises from the fact that _all_ output quantities - the
elements of $\mathbf{y}$ \- may be derived directly from the thermodynamic
potential, a single scalar function of the inputs $\mathbf{x}$. These
relationships give rise to well-known differential relationships which can be
used to define the space of ‘reasonable’ EOS models. In this section, we will
recap one common form of the consistency constraint.

In hydrocode applications it’s common to use $\mathbf{x}=(T,V)^{T}$, where $T$
is the temperature and $V$ is the volume. In that case, the relevant
thermodynamic potential is the Helmholtz free energy,

$F=U-TS~{}{\rm,}$ (2)

where $U$ is the internal energy and $S$ is the entropy. Combining with the
fundamental thermodynamic relation,

$dU=-SdT-PdV~{}{\rm,}$ (3)

where $P$ is the pressure, allows the EOS quantities $(P,U,S)$ to be defined
in terms of $F$;

$\displaystyle\begin{aligned} P&=-\left.\frac{\partial F}{\partial
V}\right|_{T}~{}{\rm,}\\\ U&=F-T\left.\frac{\partial F}{\partial
T}\right|_{V}~{}{\rm,}\\\ S&=-\left.\frac{\partial F}{\partial
T}\right|_{V}~{}{\rm.}\end{aligned}$ (4)

Finally, differentiating the first two equations by $T$ and $V$ respectively
gives the thermodynamic constistency constraint between EOS ouputs,

$P=T\frac{\partial P}{\partial T}+\frac{\partial U}{\partial V}~{}{\rm.}$ (5)

Obeying the consistency relation (5) is often considered a fundamental
requirement for a useful EOS model, and has been used to test the numerical
convergence of first-principles simualtions. The space of functions that obey
the relation is, therefore, the maximal set of functions that may be used for
EOS work, and may be used to provide an upper limit on model uncertainty. In
the next section, we present a method for sampling from the space of functions
that identically satisfy equation (5).

###  II.2 A Thermodynamically Constrained Gaussian Process

Gaussian processes (GPs) are a popular method for uncertainty quantification
since since they are lightweight and have a clear statistical interpretation,
resulting in a natural description of prediction uncertainties. They are also
highly flexible meaning that they can represent a wide range of functional
forms in a single closed form. This interpretability and flexibility makes GPs
a popular choice when quantifying model uncertainty.

Another property of GPs of particular relevance here is that a linear
transformation of a GP is another GP. This is useful result when working with
constrained functions, since a weakly constrained GP of standard form can be
transformed into a strongly constrained one by applying suitable operators.
This approach has recently been used to build GPs that automatically solve
differential equationsLange-Hegermann (2018) and to enforce a curl-free vector
output suitable for fitting magnetic fieldsJidling  _et al._ (2017), and is
the approach we take in this work.

A Gaussian process models the joint distribution of function outputs as a
multivariate normal distribution. The correlation between output values is
determined by their relative input values - typically, points closer in input
space will be more strongly correlated in output space. The details of this
correlation structure is controlled by the GP kernel, a scalar function which
determines the correlation between two input points $\mathbf{x}$ and
$\mathbf{x}^{\prime}$. At the core of this work is a GP kernel that enforces
correlations between EOS quantities and their derivatives in such a way as to
guarentee thermodynamic consistency.

We start by introducing a general GP for the Helmholtz free energy,

$F\sim GP(F_{R}(\mathbf{x}),k_{F}(\mathbf{x},\mathbf{x}^{\prime}))~{}{\rm,}$
(6)

where we have introduced a reference free energy $F_{R}$ as the prior mean,
and $k_{F}$ is the free energy kernel. For further details on GPs and their
implementation, we refer the reader to any textbook or tutorial on the
subjectSchulz, Speekenbrink, and Krause (2018); uqp ; sta . For our purposes
here, it suffices to note that this model for $F$ is unconstrained except for
a weak prior on the smoothness of the functional form introduced by the
kernel, and the reference free energy that determines the long-range
extrapolation behavior of the free energy.

The next stage in our analysis is to transform from the unconstrained free
energy to the constrained space of EOS outputs. Using equations (4) to define
a linear operator $\mathcal{\hat{Y}}_{\mathbf{x}}$ so that
$(P,U,S)^{T}=\mathcal{\hat{Y}}_{\mathbf{x}}F(\mathbf{x})$, we can use a
standard identity for the action of a linear operator on a multivariate normal
to define the multi-output EOS GP (see appendix A for some mathematical
details),

$\displaystyle\begin{pmatrix}P\\\ U\\\ S\end{pmatrix}\sim
GP(\mathcal{\hat{Y}}_{\mathbf{x}}F_{R},\mathcal{\hat{Y}}_{\mathbf{x}^{\prime}}k_{F}\mathcal{\hat{Y}}_{\mathbf{x}}^{T})~{}{\rm,}$
(7)

wheere we have supressed the dependance on input $\mathbf{x}$ except where
necessary to clarify the argument of the kernel function on which
$\mathcal{\hat{Y}}$ acts. The transformed GP (7) is constrained to sample
functions with the correct thermodynamic consistency relations. Given a choice
of reference and kernel functions for the free energy, the mean and kernel for
the EOS model can be found through straightforward manipulations, easily
handled by a symbolic algebra code. Finally, since any thermodynamic quantity
can be written as a linear transformation of the free energy, our model is
easily generalized by adding the relevant rows to the operator
$\mathcal{\hat{Y}}$.

Equation (7) provides a general approach for assessing uncertainty and
confidence in EOS models and tables. The trained GP efficiently represents a
continous space of EOS models that are consistent with the data; generating
multiple tables for use in hydrocodes can be done simply by sampling from a
single multivariate normal distribution. Finally, since GPs are a popular tool
in all areas of uncertainty quantification, this picture opens the door for a
multitude of new developments including uncertainty propegation, acceleration
of simulation studies, and principled design and analysis of new experiments.

##  III Application to high-temperature Boron Carbide

The development of new theoretical EOS models typically proceeds by running a
large set of first-principles simulations and using them to train a wide
ranging model. One example is the recent work for $B_{4}C$ by Zhang _et al._
Zhang (2020), which combined simulations using path integral Monte Carlo
(PIMC), quantum molecular dynamics (QMD), Green’s function Korringa-Kohn-
Rostoker (KKR) theory, activity expansion (ACTEX) and other methods to build a
wide ranging interpolation which was cross-validated by experiments at the
National Ignition Facility. For the purposes of this paper, we will use the
results of the PIMC and QMD simulations which form a set of 231 samples of the
mapping $(T,\rho)\rightarrow(P,U)$, with statistical uncertainties on
$(T,P,U)$. Note that for consistency with Zhang _et al._ we will plot the mass
density $\rho$, which is related to the volume used in section II.2 by
$\rho({\rm g/cc})=619.8/V(a^{3}_{0}/{\rm atom})$, and we do not include the
entropy output in equation (7) since those values were not reported. We use
these data to find optimal hyperparameters for our EOS GP regression model as
described previously, and to make predictions of the uncertainties in the EOS.

Figure 1: First-principals data used to train our thermodynamically consistent
Gaussian Process, consisting of previously-published path integral Monte Carlo
(PIMC) and quantum molecular dynamics (QMD) simulations for boron carbideZhang
(2020). We use the published pressure $P$ and internal energy $U$, as a
function of temperature $T$ and mass density $\rho$, along with statistical
uncertainties on $T$, $P$ and $U$. The two simulation methods overlap in a
small region around $T=100{\rm eV}$, where both are pushed to the limits of
their applicability.

We plot the pressure $P$ and internal energy $U$, along with their quoted
uncertainties, in figure 1. While the EOS quantities themselves are in good
agreement between the PIMC and DFT methods, the uncertainties are quite
different and show a discontinuous change where the two methods overlap at
$T\sim 100\rm{eV}$. In the PIMC simulations statistical uncertainties grow to
almost 10% at lower temperatures while the QMD uncertainties are much lower.
One goal of our GP regression technique is to combine the results from both
simultions and provide a measure of the true uncertainty in this region.

The data in figure 1 also identify a practical challenge to fitting GPs, and
machine learning models in general, to EOS data: the inputs and outputs must
span many orders of magnitude. For the inputs $T$ and $V$ we can take the
standard approach and log-scale, however our method of enforcing constraints
requires that any transformation of the outputs $P$ and $U$ is linear, making
log-scaling inappropriate for those quantities. We instead rely on the
reference model to capture the large-scale changes in output and use the GP to
model deviations around the reference, known as the _excess_ pressure $\Delta
P$ and internal energy $\Delta U$. For this dataset, the excess EOS quantities
vary smoothly between 0.5 and 3, a range which is suitable for reliable
training of the GP model.

###  III.1 Specifying and Training the EOS GP

Our model is specified by the reference free energy $F_{R}$ and free energy
kernel function $k_{F}$. The choice of reference function does not effect the
behavior of the GP close to training data, but does determine the functional
form on extrapolation since the posterior mean reverts to the reference in
this limit. Here, we will use the ideal gas free energy. The choice of kernel
function is also relatively free and is intuitively related to our prior
belief on the form of $F$. In this work, we use a Gaussian or radial basis
function (RBF) kernel,

$\displaystyle k_{F}(lnT,ln\rho;lnT^{\prime},ln\rho^{\prime})=\sigma
e^{-\frac{(lnT-lnT^{\prime})^{2}}{2\lambda_{lnT}^{2}}}e^{-\frac{(ln\rho-
ln\rho^{\prime})^{2}}{2\lambda_{ln\rho}^{2}}}~{}{\rm,}$ (8)

which introduces an unknown correlation length for each input direction,
$\lambda_{lnT}$ and $\lambda_{ln\rho}$, and an overall correlation strength
$\sigma>0$. These parameters are learned from the EOS data based on the
probability of the data defined by equation (7). For this work, the optimized
correlation lengths are of particular interest since they describe the
lengthscale over which prediction uncertainty increases as the GP is
extrapolated away from available data.

Given the choice of free energy kernel function (8) we use the computer
algebra package sympy to determine the thermodynamically consistency kernel
function. Mathematical details and a description of the resulting kernel are
given in appendix A. With the kernel function in hand we ‘train’ the EOS GP by
optimizing the parameters $\sigma$, $\lambda_{lnT}$ and $\lambda_{lnV}$ by
maximizing the log-probability of the data. It is well known that the training
process can run into numerical issues when data are close together and have
small uncertainties, as is the case here; we solve this problem by randomly
perturbing the training inputs according to the statistical uncertainties in
the simulated temperatures, and using a stochastic gradient descent algorithm
that is robust to these variations. We have found that this is a reliable
regularization that does not require _ad-hoc_ modifications to the GP model,
and has the added bonus of providing a relatively simply way of accounting for
uncertainties in the simulated temperature. At the end of this process we find
optimal parameters $\lambda_{lnT}=4.78$, $\lambda_{ln\rho}=4.03$, for which
the $R^{2}$ correlation coefficient between simulation data and GP prediction
is better that $99.5\%$.

###  III.2 Uncertainties in the $B_{4}C$ EOS

The optimized EOS GP, for the first time, provides us with a wide ranging EOS
model that captures both data and model uncertainty. We visualize the
uncertain model by plotting the mean and standard deviation of the predicted
pressure and energy as a function of input density $\rho$ and temperature $T$
in figure 2.

Figure 2: Predictions from our constrained Gaussian process EOS model for
$B_{4}C$. Points show the simulation values (top row) and uncertainties
(bottom row) used to train the GP, and contours show the GP interpolation, on
a common colorscale. The thermodynamically consistent GP faithfully reproduces
the training data ($R^{2}=0.995$) with a significantly reduced uncertainty.

The prediction means in the top row of 2 compare the mean of the predicted
excess pressure and energy (contours) with the simulations (points) on a
common color scale. At the locations of the simulation points the GP mean is
in excellent agreement with the simulations, with an average error of $1\%$,
and the GP prediction uncertainty is typically a factor of $10$ lower than the
uncertainty in the simulation themselves. Away from the training data the EOS
GP picks up important features from the simulations and smoothly interpolates
over the whole space. In principal the EOS GP model can predict the pressure
and energy at _any_ value of $lnT$ and $ln\rho$, however prediction
uncertainties rapidly increase away from the available data as shown in plots
of the prediction uncertainty in the bottom row of figure 2. The rate of
increase is a combination of the inferred smoothness of the free energy,
described by the optimized values of $\lambda_{lnT}$ and $\lambda_{ln\rho}$,
and the effect of the constraints. We find that while the inferred length
scales suggest that the free energy can be extrapolated for a decade or more
in $T$ and $\rho$, uncertainty in the EOS qantities grows much faster as a
result of the thermodynamic constraints.

Our GP model provides a principled and efficient way of sampling from the
space of EOS models for downstream tasks, for example in hydrodynamics
simulations or the analysis of experimental data. One important example is the
calculation of the _Hugoniot_ , the set of thermodynamic states that can be
accessed through shock compression of a material, along which the majority of
experimental data are constrained to lie. We show the distribution of shock
Hugoniot curves for $B_{4}C$ in figure 3, generated by sampling from our GP
model and postprocessing each sample to generate a set of Hugoniot curves. The
uncertainty in the Hugoniot, measured here by the standard deviation in
predicted density at given temperature, is lower than the uncertainty in the
underlying EOS quantities as a consequence of the dense sampling of
simulations in the Hugoniot region and the physically correlated uncertainties
in EOS quantities that our model provides.

Figure 3: Predicted Hugoniot, with uncertainty, for $B_{4}C$. We use our GP
model to sample from the distribution of equation-of-state models that are
consistent with our first-principles simulations. For each EOS sample, we
calculated the Hugoniot curve and plot the standard devition in density at
constant temperature.

###  III.3 Importance of the thermodynamic consistency constraint

Our EOS GP model ensures thermodynamic consistency by limiting the space of
functions that it can sample. A result of this is that our GP model can be
expected to have lower prediction uncertainties than a standard unconstrained
on. In effect, the laws of thermodynamics have been used as an additional
source of information resulting in lower uncertainty.

Figure 4: Including the thermodynamic consistency constraint on model form
significantly reduces the uncertainty in the resulting EOS. We compare the
uncertainty in the predicted pressure for $B_{4}C$ from a Gaussian process
models that neglect (left) and include (right) the constraint. Inclusion of
the constraint reduces the uncertainty by a factor of 10 over most of the
space where training simulations are available.

We compare uncertainties with and without the thermodynamic constraint by
plotting the fractional uncertainty in predicted pressure from two models in
figure 4. In the figure, the left panel shows the uncertainty in pressure for
a GP with a standard RBF kernel which samples from a space of smooth functions
but does not enforce any other constraints; the right panel shows the results
from our EOS GP model. The effect of the constraint is clear; for the standard
GP the prediction uncertainty falls below 1% only where simulations are
densely sampled and the temperature is low enough that statistical noise on
the DFT simulations is low, while for the EOS GP the uncertainty is less than
1% over almost the entire range of the simulated dataset. In the extrapolation
region at the bottom left of each plot the effect of the constraint is also
clear, in increasing uncertainty as discussed in the previous section.
Clearly, an approach to model uncertainty that neglects the constraints would
dramatically overestimate uncertainty where data are available, and
underestimate when extrapolating into regions where no data can be collected.

##  IV Summary & Further Work

In this work we have developed a thermodynamically constrained Gaussian
process for use in equation-of-state (EOS) modeling. Our EOS GP is constrained
to explore the space of thermodynamically consistent functions by explicitly
relating all thermodynamic quantities to an underlying free energy model. When
trained on simulation or experimental data, the EOS GP naturally estimates
model uncertainty in the predicted thermodynamic quantities which, until now,
have been very difficult to capture. This is the main contribution of this
work so far.

We have demonstrated the importance of model uncertainty, and the effect of
the thermodynamic constraints, by applying our model to existing first-
principals simulations for $B_{4}C$. We show that our EOS GP provides a high
quality representaion of the EOS with uncertainties of less than $1\%$ over a
wide range of temperatures and densities, and that prediction uncertainties
are generally much lower than the statistical uncertainties in the simulations
themselves. Comparison with an unconstrained GP model shows that neglecting
the constraint results in an _overestimate_ of the uncertainty by factors of
several in regions where data are available, and an _underestimate_ when
extrapolating to new parts of the parameter space. As a result, an attempt to
measure model uncertainty without accounting for the thermodynamic constraints
is likely to be unreliable.

The thermodynamic consistency constraint considered here is not the only
constraint on EOS models, and as a result the uncertainties we have presented
are an upper bound. Perhaps the most important additional constraint is
thermodynamic stability resulting in positivity of the specific heat and bulk
modulus. Since these are inequality constraints enforcing them requires a
different approach to the one adopted hereSwiler  _et al._ (2020), and will be
the subject of further work.

Known limiting behaviour is also an essential component of a meaningful EOS
model. Our EOS GP controls long-range behaviour through the introduction of a
reference free energy as the prior mean which ensures that the mean prediction
reverts to the reference far from any available data. However, in the current
approach uncertainty will (incorrectly) continue to grow while the model
reverts to the known behaviour. Controlling the increase in uncertainty, as
well as implementing multiple material phases, in order to provide a
sufficiently wide-ranging EOS to be of use in hydrodynamics simulations is
ongoing.

Despite the further work needed, the use of Gaussian processes has already
provided important insights, and promises to significantly change the way the
EOS tables are generated and stored. Prediction uncertainties can be used to
adaptively choose points for new simulations or experiments to focus on
uncertain regions, potentially reducing the number of data needed and
accelerating the development of new EOSs. GPs also provide a natural and
memory-efficient representation of the full distribution of EOS models in a
single object, and do not require stuctured grids in $T$ and $\rho$, which
suggests a useful extension to the current EOS table formats to include
uncertainty. As a result we believe that GPs could form the core of a new
uncertainty-aware EOS modeling framework.

###### Acknowledgements.

Prepared by LLNL under Contract DE-AC52-07NA27344.

## References

  * Gaffney  _et al._ (2018) J. Gaffney, S. Hu, P. Arnault, A. Becker, L. Benedict, T. Boehly, P. Celliers, D. Ceperley, O. Čertík, J. Clérouin, G. Collins, L. Collins, J.-F. Danel, N. Desbiens, M. Dharma-wardana, Y. Ding, A. Fernandez-Pañella, M. Gregor, P. Grabowski, S. Hamel, S. Hansen, L. Harbour, X. He, D. Johnson, W. Kang, V. Karasiev, L. Kazandjian, M. Knudson, T. Ogitsu, C. Pierleoni, R. Piron, R. Redmer, G. Robert, D. Saumon, A. Shamp, T. Sjostrom, A. Smirnov, C. Starrett, P. Sterne, A. Wardlow, H. Whitley, B. Wilson, P. Zhang, and E. Zurek, High Energy Density Physics 28, 7 (2018).
  * Militzer (2021) B. Militzer, Physical Review E 103 (2021), 10.1103/PhysRevE.103.013203.
  * Zhang (2020) S. Zhang, Physical Review E 102 (2020), 10.1103/PhysRevE.102.053203.
  * Wu  _et al._ (2019) C. J. Wu, D. A. Young, P. A. Sterne, and P. C. Myint, The Journal of Chemical Physics 151, 224505 (2019).
  * MacKay  _et al._ (1998) D. J. MacKay _et al._ , NATO ASI series F computer and systems sciences 168, 133 (1998).
  * Ali  _et al._ (2020) S. J. Ali, D. C. Swift, C. J. Wu, and R. G. Kraus, Journal of Applied Physics 128, 185902 (2020).
  * Kamga  _et al._ (2014) P.-H. Kamga, B. Li, M. McKerns, L. Nguyen, M. Ortiz, H. Owhadi, and T. Sullivan, Journal of the Mechanics and Physics of Solids 72, 1 (2014).
  * Root (2010) S. Root, Physical Review Letters 105 (2010), 10.1103/PhysRevLett.105.085501.
  * Walters  _et al._ (2018) D. J. Walters, A. Biswas, E. C. Lawrence, D. C. Francom, D. J. Luscher, D. A. Fredenburg, K. R. Moran, C. M. Sweeney, R. L. Sandberg, J. P. Ahrens, and C. A. Bolme, Journal of Applied Physics 124, 205105 (2018).
  * Brown and Hund (2018) J. L. Brown and L. B. Hund, Journal of the Royal Statistical Society: Series C (Applied Statistics) 67, 1023 (2018).
  * Kennedy and O’Hagan (2001) M. C. Kennedy and A. O’Hagan, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63, 425 (2001).
  * Higdon  _et al._ (2008) D. Higdon, J. Gattiker, B. Williams, and M. Rightley, Journal of the American Statistical Association 103, 570 (2008).
  * Lange-Hegermann (2018) M. Lange-Hegermann, Advances in Neural Information Processing Systems 31 (2018).
  * Jidling  _et al._ (2017) C. Jidling, N. Wahlström, A. Wills, and T. B. Schön, Advances in Neural Information Processing Systems 30 (2017).
  * Schulz, Speekenbrink, and Krause (2018) E. Schulz, M. Speekenbrink, and A. Krause, Journal of Mathematical Psychology 85, 1 (2018).
  * (16) “Uqpy documentation: Gaussian processes,” .
  * (17) “Stan documentaion: Gaussian processes,” .
  * Swiler  _et al._ (2020) L. P. Swiler, M. Gulian, A. L. Frankel, C. Safta, and J. D. Jakeman, Journal of Machine Learning for Modeling and Computing 1 (2020).

##  Appendix A Mathematical Details

Our goal is to derive a Gaussian process kernel that results in a statistical
model that samples the space of thermodynamic functions. In this section we
will provide some mathematical details and explicit details on the approach
taken to derive the kernel used in this work. Much of the details here follow
the work of Jidling _et al._ Jidling  _et al._ (2017) with some
specializations for the the application to EOS modeling. Finally we will make
some general observations about the structure of the EOS GP kernel.

###  A.1 General Approach

We start by defining a free energy

$F:\mathbb{R}^{2}\mapsto\mathbb{R};~{}\mathbf{x}\mapsto F(\mathbf{x})$ (9)

and an equation of state

$EOS:\mathbb{R}^{2}\mapsto\mathbb{R}^{N};~{}\mathbf{x}\mapsto\mathbf{y}(\mathbf{x})~{}{\rm,}$
(10)

both of which act on a pair of natural thermodynamic variables $\mathbf{x}$.
As discussed in the main text, the choice of natural variables defines the
form of $F$ and the laws of thermodynamics give operators relating free energy
to EOS outputs, which we write

$\mathcal{\hat{Y}}_{\mathbf{x}}:\mathbb{R}\mapsto\mathbb{R}^{N};~{}F(\mathbf{x})\mapsto\mathcal{\hat{Y}}_{\mathbf{x}}F(\mathbf{x})=\begin{pmatrix}\mathcal{\hat{Y}}_{\mathbf{x}}^{(1)}F(x)\\\
\vdots\\\ \mathcal{\hat{Y}}_{\mathbf{x}}^{(N)}F(x)\end{pmatrix}~{}{\rm.}$ (11)

A key observation for this work is that each of the
$\mathcal{\hat{Y}}_{\mathbf{x}}^{(i)}$ is a linear operator, in which case
introducing a Gaussian process for the free energy $F(\mathbf{x})\sim
GP(F_{R}(\mathbf{x}),\,K_{F}(\mathbf{x},\mathbf{x}^{\prime}))$ results in
another GP for the EOS quantities with a transformed mean and covariance
function,

$y(\mathbf{x})\sim
GP(\mathcal{\hat{Y}}_{\mathbf{x}}F_{R}(\mathbf{x}),\,\mathcal{\hat{Y}}_{\mathbf{x}^{\prime}}K_{F}(\mathbf{x},\mathbf{x}^{\prime})\mathcal{\hat{Y}}_{\mathbf{x}}^{T})~{}{\rm.}$
(12)

As a result of the transformation to EOS outputs, the free energy kernel has
become an $N\times N$ matrix of kernel functions with elements

$(\mathcal{\hat{Y}}_{\mathbf{x}^{\prime}}K_{F}(\mathbf{x},\mathbf{x}^{\prime})\mathcal{\hat{Y}}_{\mathbf{x}}^{T})_{ij}=\mathcal{\hat{Y}}_{\mathbf{x}^{\prime}}^{(i)}\,\mathcal{\hat{Y}}_{\mathbf{x}}^{(j)}\,K_{F}(\mathbf{x},\mathbf{x}^{\prime})~{}{\rm.}$
(13)

The off-diagonal kernel functions, $i\neq j$, play a critical role in our
model by enforcing correlations such that the thermodynamic constraints are
satisfied.

###  A.2 Specific Example

In this work we take the natural variables $\mathbf{x}=(lnT,lnV)^{T}$ and
limit our EOS outputs to the pressure and internal energy,
$\mathbf{y}=(P,U)^{T}$. The transformation operator comes directly from the
thermodynamic relations (4), which following log-transforming the input points
become

$\mathcal{\hat{Y}}_{\mathbf{x}}\equiv\mathcal{\hat{Y}}_{(lnT,lnV)^{T}}=\begin{pmatrix}-\frac{1}{V}\frac{\partial}{\partial
lnV}\\\ 1-\frac{\partial}{\partial lnT}\end{pmatrix}{\rm.}$ (14)

The log transformation of inputs allows for a more numerically stable handling
of the wide range of thermodynamic conditions for whith the EOS model must be
applied. The outputs also vary over a wide range, however log scaling of the
outputs results in a non-linear operator $\mathcal{\hat{Y}}_{\mathbf{x}}$,
precluding our simple use of GPs. Instead we use a GP to model the excess free
energy, $\Delta F(\mathbf{x})=F(\mathbf{x})/F_{R}(\mathbf{x})\sim
GP(1,\,K_{\Delta F}(\mathbf{x},\mathbf{x}^{\prime}))$ and the EOS model
becomes

$\frac{\mathbf{y}(\mathbf{x})}{\mathbf{y}_{R}(\mathbf{x})}\sim
GP(\mathcal{\hat{Y}}_{\mathbf{x}}\,1,\,\mathcal{\hat{Y}}_{\mathbf{x}^{\prime}}K_{\Delta
F}(\mathbf{x},\mathbf{x}^{\prime})\mathcal{\hat{Y}}_{\mathbf{x}}^{T})~{}{\rm,}$
(15)

where $\mathbf{y}_{R}(\mathbf{x})$ are the EOS outputs for the reference free
energy model.

Using equations (12) and (13) results in a two dimensional matrix of kernel
functions giving the physics-driven correlations between EOS quantities,

$\displaystyle\begin{aligned}
K_{\Delta\mathbf{y}}(\mathbf{x},\mathbf{x}^{\prime})&\equiv\mathcal{\hat{Y}}_{\mathbf{x}^{\prime}}K_{\Delta
F}(\mathbf{x},\mathbf{x}^{\prime})\mathcal{\hat{Y}}_{\mathbf{x}}^{T}\\\
&=\begin{pmatrix}k_{PP}(\mathbf{x},\mathbf{x}^{\prime})&k_{PU}(\mathbf{x},\mathbf{x}^{\prime})\\\
k_{UP}(\mathbf{x},\mathbf{x}^{\prime})&k_{UU}(\mathbf{x},\mathbf{x}^{\prime})\\\
\end{pmatrix}{\rm,}\end{aligned}$ (16)

where, for example, the kernel for correlations between pressure values is
given by

$\displaystyle k_{PP}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle\equiv\mathcal{\hat{Y}}_{\mathbf{x}^{\prime}}^{(P)}\,\mathcal{\hat{Y}}_{\mathbf{x}}^{(P)}\,K_{\Delta
F}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=\frac{1}{V^{\prime}}\frac{\partial}{\partial
lnV^{\prime}}\frac{1}{V}\frac{\partial}{\partial lnV}K_{\Delta
F}(lnT,lnV;lnT^{\prime},lnV^{\prime})~{}{\rm.}$

For our choice of inner kernel function $K_{\Delta
F}(lnT,lnV;lnT^{\prime},lnV^{\prime})$ given in equation (8), we evaluate the
4 resulting EOS kernels using the symbolic algebra package sympy. This allows
an easy derivation of all the required terms to implement our multi-output GP,
as well as gradients to accelerate the optimization of hyperparameters.

We plot the resulting kernel functions in figure 5. This plot shows each of
the elements of equation (A.2) as a function of $\mathbf{x}$ and with
$\mathbf{x}^{\prime}=(0,0)^{T}$ shown in red. The thermodynamic constraints
result in asymmetric correlations between EOS quantities with important
implications for the mutual information between different EOS quantities at
the same and different points.

Figure 5: The correlation structure in the outputs of our Gaussian Process,
induced by the thermodynamic constraints. We plot the degree of correlation
between outputs at the red point and a test point, as a function of test point
location. Panels are labeled with the output quantity at the test point first,
i.e., the "Pressure-Energy" panel shows the correlation between the pressure
at the test point and the energy at the red point. The asymmetric correlation
structure, and correlations between different EOS quantities, are a direct
result of the thermodynamic consistency of our approach.