text:
\section{Introduction}

The formation of complex structures during solidification often results from
a subtle interplay of non-equilibrium, nonlinear processes, for which
seemingly small changes in molecular dynamics at the nanoscale can produce
large morphological changes at all scales. One popular example of this
phenomenon is the formation of snow crystals, which are ice crystals that
grow from water vapor in an inert background gas. Although this is a
relatively simple physical system, snow crystals display a remarkable
variety of columnar and plate-like forms, and much of the phenomenology of
their growth remains poorly understood. (For a review of the physics of snow
crystal growth, see \cite{libbrechtreview}.)

Observations of snow crystal growth dating back to the 1930s \cite{nakaka, libbrechtreview} reveal a complex and puzzling
dependence on temperature. At pressures near one bar and water vapor supersaturation levels often found
in clouds, for example, ice crystals typically grow into thin plate-like
forms near -2 C, slender columns and needles near -5 C, very thin plates
again near -15 C, and columns again below -30 C. The observed variations in
snow crystal structure with temperature and supersaturation are often
displayed in a snow crystal morphology diagram, as shown in Figure \re
{morphdiagram}. Despite more than a half-century of study, we still do not
understand the basic physical mechanisms that are responsible for the
unusual temperature-dependent morphologies of growing ice crystals \cit
{libbrechtreview, nelson, pruppacher}.

\begin{figure}[htb]
  \centering

  \includegraphics[width=5in,height=3.64in,keepaspectratio]{SnowMorph}
  \caption{A snow crystal morphology
diagram showing the variations in growth forms as a function of temperature
and supersaturation for crystals grown in air at a pressure of one bar (from
\protect\cite{libbrechtreview}). The water saturation line shows the supersaturation of
supercooled liquid water with respect to ice. Note that thin plate-like
crystals appear only in narrow temperature ranges near -2 C and -15 C.}
  \label{morphdiagram}
\end{figure}

Our focus here will be on the growth of thin plates near -15 C, as the
largest and thinnest crystals form in this temperature region, making this a
dominant feature in the morphology diagram. Aspect ratios of
L_{a}/L_{c}\approx 100$ or larger are commonly found with atmospheric ice
crystals, where $L_{a}$ is the crystal size along the a-axis of the crystal
(roughly equal to the diameter of the plate), and $L_{c}$ is the size of the
plate in the c-direction (equal to the plate thickness). Large stellar
plates typically have plate diameters of 2-4 mm with thicknesses of 10-20
\mu $m \cite{nakaka}.

In modeling the growth of ice crystals from water vapor, we write the
perpendicular growth velocity of a surface in the usual Hertz-Knudsen form
\begin{eqnarray}
v &=&\alpha \frac{\Omega (p_{surf}-p_{sat})}{\sqrt{2\pi mkT}}  \label{hertz}
\\
&=&\alpha v_{kin}\sigma _{surf}  \nonumber
\end{eqnarray
where $p_{surf}$ is the water vapor pressure at the growing surface, $p_{sat}
$ is the water vapor pressure in equilibrium with ice, $m$ is the mass of a
water molecule, and $\Omega $ is the molecular volume \cite{libbrechtreview
. This relation defines the condensation coefficient, $\alpha ,$ that
parameterizes the attachment kinetics at the ice surface, as well as the
kinetic velocity $v_{kin}$, where $\sigma _{surf}=(p_{surf}-p_{sat})/p_{sat}$
is the water vapor supersaturation at the surface. In general $\alpha $ will
be different for the two principal facets of ice, and for fast kinetics on
non-faceted surfaces we expect $\alpha \rightarrow 1.$

If the supersaturation had a constant value independent of location near the
ice surface, then the growth of a thin plate would imply $\alpha
_{prism}/\alpha _{basal}\approx v_{prism}/v_{basal}\approx L_{a}/L_{c}\gg 1.$
Diffusion modeling shows that the supersaturation at the edges of a thin
plate is not substantially higher than at the middle of the basal facets,
and for a large fast-growing plate the minimum supersaturation is actually
found at the edges of the plate \cite{libbrechtreview}. Thus
diffusion-limited growth will generally not promote the growth of thin
plates in the absence of a large $\alpha _{prism}/\alpha _{basal}$ ratio,
and in many cases diffusion requires that we must have $\alpha
_{prism}/\alpha _{basal}>L_{a}/L_{c},$ which further accentuates the
difference in attachment kinetics on the two facet surfaces.

The problem with this simple inference is that measurements of ice crystal
growth rates for faceted surfaces have yielded values of $\alpha
_{prism}/\alpha _{basal}$ that are far smaller than what is needed to form
thin plate-like crystals. For example, numerous researchers have found that
ice crystals grow into roughly isometric shapes at low background pressures,
where effects from particle diffusion are reduced \cite{gonda, beckman, lamb
, suggesting that $\alpha _{prism}/\alpha _{basal}$ is near unity for the
conditions in these experiments. Since the growth of thin plates requires
\alpha _{prism}/\alpha _{basal}\gg 1,$ the discrepancy between low-pressure
and high-pressure measurements necessarily requires a nontrivial explanation.

In an earlier paper we proposed a mechanism to solve this problem that we
called \textit{structure-dependent attachment kinetics} (abbreviated \textit
SDAK}) \cite{sdak}. In essence, we proposed that $\alpha _{prism}$ depends
strongly on the structure of the prism facet surface, such that $\alpha
_{prism}$ increases when the width of the prism facet decreases, becoming
near unity when the width approaches atomic dimensions. For the case of a
thin plate, the width of the facet surface on the edge of the plate is
approximately $w\approx \sqrt{aR},$ where $a$ is the molecular step height
and $R$ is the effective radius of curvature of the edge of the growing
plate \cite{sdak}. Taking $R\approx 0.5$ $\mu $m for a typical thin plate
yields $w\approx 40a$ in the case of ice. We proposed that the SDAK
mechanism leads to an edge-enhancing growth instability that promotes the
formation of thin plates.

Below we present observations of a hysteresis behavior in the growth of thin
plates that strongly supports our proposed SDAK mechanism. Quantitative
measurements reveal that $\alpha _{prism}$ is indeed substantially larger on
the edge of a thin plate than on a flat prism facet. This behavior then
leads to an edge-enhancing growth instability that is responsible for the
formation thin plate-like crystals near -15 C.

\section{Observations}

\subsection{Intrinsic Growth Rates of Ice Facets}

We recently published improved measurements of the growth rates of faceted
basal surfaces as a function of temperature \cite{basalgrowth}, indicating
\alpha _{basal}\approx \exp (-\sigma _{basal}/\sigma _{surf})$ over a broad
temperature range, and specifically $\sigma _{basal}=\sigma
_{basal,0}\approx 2.0$ percent at -15 C. Similar measurements of prism
growth rates \cite{unpub} give the same functional form for $\alpha _{prism}$
at -15 C, and for flat prism facets we have $\sigma _{prism}=\sigma
_{prism,0}\approx 4.2$ percent. These measurements were made in a background
of air at a pressure of approximately 0.03 bar, since lower pressures reduce
effects from diffusion-limited growth. Similar measurements at 1 bar yield
\alpha $ values that are consistent with the low-pressure data \cit
{basalgrowth}. These measurements are roughly consistent with previous data,
although several improvements in our experimental techniques have yielded
\sigma _{prism,0}$ values that are higher than quoted in \cit
{libbrechtreview}.

These data indicate that $\alpha _{prism}/\alpha _{basal}\approx \exp
(-\Delta \sigma /\sigma )$ with $\Delta \sigma =2.2$ percent at -15 C. Since
this ratio is less than unity for all supersaturations, it suggests (in the
absence of the SDAK mechanism) that the intrinsic growth form at -15 C
should be a columnar crystal. Demonstrating this intrinsic columnar growth
directly in low-pressure observations is difficult, however, since any ice
surface contacting a substrate may be affected by substrate interactions.
Generally we have found that these interactions reduce the nucleation
barrier and yield untrustworthy growth measurements, especially at low
supersaturations when the intrinsic 2D nucleation rate is low. Nevertheless,
when one basal facet is flat against the substrate (as in \cite{basalgrowth
), we do occasionally observe the growth of tall columns. Since substrate
interactions cannot reduce the nucleation barrier all along a tall column,
our interpretation of these observations is that sometimes the substrate
interactions are greatly reduced, due to some impurity coating on the
substrate that we have not yet determined. When this happens we witness the
intrinsic growth rates on all facets, giving columnar crystals. More
typically substrate interactions increase the growth of the prism facets via
enhanced nucleation where the crystal contacts the substrate, resulting in
more isometric ice prisms.

Our overall conclusion from this series of measurements is that the
intrinsic attachment coefficients for large facet surfaces are adequately
described by the functional form above with $\sigma _{basal,0}\approx 2.0$
percent and $\sigma _{prism,0}\approx 4.2$ percent at -15 C. Growth at low
background pressures yields mainly large faceted surfaces, and the available
evidence suggests that the intrinsic growth form at -15 C may indeed be a
columnar crystal. Additional low-pressure measurements, preferably with free
falling or levitated crystals \cite{levitated} to avoid substrate
interactions, are needed to confirm this surprising conclusion. The fact
that thin plates are usually observed at -15 C is then entirely because of
the SDAK mechanism. The resulting edge-enhancing instability appears only at
higher background pressures, since the instability depends in part on the
effects of diffusion-limited growth. We now examine additional measurements
supporting this claim.

\subsection{Plate-on-Pedestal Growth}

We begin with a detailed description of the growth of a representative
crystal at -15 C in air at a pressure of one bar. Using the apparatus
described in \cite{vigapparatus}, we began the experiment by dropping a
small hexagonal plate crystal onto a temperature-controlled sapphire
substrate. The initial radius of this crystal (taken to be half the distance
from one prism facet to the opposite facet) was 13.5 $\mu $m, and the
initial thickness was 2.5 $\mu $m. The radius was determined by direct
optical microscopy while the thickness was measured by broad-band
interferometry, as described in \cite{vigapparatus}. Initially the
temperature $T_{IR}$ of an ice reservoir above the crystal was equal to the
substrate temperature $T_{subst}$, so the supersaturation was equal to zero.
We then slowly increased the ice reservoir temperature and thus the
supersaturation above the test crystal, causing it to grow.

\begin{figure}[p]
  \centering

  \includegraphics[bb=0 0 993 1237,width=5in,height=6.23in,keepaspectratio]{TimeSeries}
  \caption{An example of \textquotedblleft
plate-on-pedestal\textquotedblright\ ice crystal growth. Top: The radius and
thickness of the test crystal as a function of time. One basal surface of
the crystal was resting on a substrate, and the test chamber was filled with
air at a pressure of one bar. Middle: The basal $v_{h}$ and prism $v_{r}$
growth velocities. Bottom: The supersaturation far above the crystal
\protect\sigma _{\infty }$ and the supersaturation at the basal surface
\protect\sigma _{surf}.$ The latter quantity was determined from the basal
growth velocity using Equation \protect\ref{hertz}, assuming the measured
condensation coefficient $\protect\alpha _{basal}=\exp (-\protect\sigma
_{basal,0}/\protect\sigma _{surf})$ described in the text.}
  \label{timeseries}
\end{figure}

\begin{figure}[p]
  \centering

  \includegraphics[width=5in,height=6.45in,keepaspectratio]{0528-1-composite1}
  \caption{An example of \textquotedblleft
plate-on-pedestal\textquotedblright\ ice crystal growth. Each image is
labeled with the time in seconds, and the crystal sizes and growth
velocities at each time are given in Figure \protect\ref{timeseries}.
Between $t=130$ and $t=180$ a thin plate began growing from the top edge of
the ice prism. Throughout its subsequent growth this plate did not contact
the substrate except via its supporting pedestal. By $t=301$ the outer edge
of the plate had begun to thicken as the radial velocity diminished.}
  \label{composite}
\end{figure}

The crystal radius, its thickness, and the growth velocities as a function
of time are shown in Figure \ref{timeseries}, along with the supersaturation
above the crystal $\sigma _{\infty }$. Note that $\sigma _{\infty }$ was
changed with time during the experiment by changing $T_{IR},$ giving $\sigma
_{\infty }(t)$ as shown in the Figure, and this was the only external
experimental parameter that was changed during the run. Since the distance
from the substrate to the ice reservoir was much greater than the size of
the test crystal, $\sigma _{\infty }$ can be assumed (for modeling purposes)
to be the supersaturation at a hemispherical boundary far away from the
crystal.

As $\sigma _{\infty }$ was increased, the ice prism initially maintained its
simple shape, namely that of a thick hexagonal plate. At roughly $t\approx
160$ seconds, the morphology of the crystal changed as a thin plate-like
crystal began to grow out from the upper edge of the prism. We call this
\textquotedblleft plate-on-pedestal\textquotedblright\ growth because the
thin plate grew outward above the substrate, supported by a central
pedestal, like one half of a capped column crystal \cite{libbrechtreview}.

Figure \ref{composite} shows images of the plate-on-pedestal crystal at
various stages during its growth. The morphology of the thin plate during
its rapid growth phase is essentially that of a sectored plate crystal \cit
{libbrechtreview}, and from additional observations we determined that the
ridges along the $a$-axes of the crystal were on the underside of the plate,
on the basal surface nearest the substrate. The top of the plate appeared to
be a flat basal facet surface without ridges. In Figure \ref{timeseries} the
radius refers to that of the top plate, while the thickness gives the
overall thickness of the crystal from the substrate to the top basal
surface. Once $\sigma _{\infty }$ was reduced, the edge of the thin plate
thickened, and eventually a thick \textquotedblleft rib\textquotedblright\
appeared around the circumference of the plate, as seen in the last image in
Figure \ref{composite}.

\begin{figure}[t]
  \centering

  \includegraphics[bb=0 0 1239 984,width=4.5in, keepaspectratio]{vhvrplot}
  \caption{A plot of the basal growth
velocity $v_{h}$ and prism growth velocity $v_{r}$ of the test crystal (see
also Figure \protect\ref{timeseries}). By $t=166$ the thin plate had formed and
the prism velocity quickly increased. After $t=178$ the supersaturation was
reduced and the edge of the plate thickened as both $v_{h}$ and $v_{r}$
became smaller.}
  \label{vhvr}
\end{figure}

It is instructive to examine a $v_{h}$-$v_{r}$ plot of the crystal growth
shown in Figure \ref{vhvr}, as this exhibits a characteristic hysteresis
behavior. In the early stages of the experiment, both $v_{h}$ and $v_{r}$
increased slowly as $\sigma _{\infty }$ was increased, in this case with
v_{r}\approx 1.5v_{h}.$ Since the prism facets were in contact with the
substrate during this phase of the growth, $v_{r}$ was likely influenced by
substrate interactions that increased nucleation on the prism facets and
thus increased $v_{r}$ compared to that for a free-standing prism facet
surface. Had substrate interactions been absent, we would expect $v_{r}<v_{h}
$ because $\sigma _{prism,0}>\sigma _{basal,0}.$ By $t=166$ the thin plate
had begun to grow out from the top edge of the ice prism, which changed the
structure of the prism surface from a broad facet surface to a sharp edge.
According to our interpretation of the subsequent growth behavior, the sharp
edge of the prism surface resulted in a higher $\alpha _{prism}$ (and thus a
higher $v_{r})$ because of the SDAK effect. This brought into play a growth
instability that sharpened the edge of the plate further, leading to a rapid
increase in $v_{r}$ while $v_{h}$ (and thus the supersaturation near the
crystal surface) remained essentially constant. Around $t=178$ the
supersaturation $\sigma _{\infty }$ was reduced as shown in Figure \re
{timeseries}, which caused $v_{h}$ to drop. Because the edge of the plate
was still sharp, however, $\alpha _{prism}$ (and thus $v_{r}$) remained
high. As $\sigma _{\infty }$ was reduced further, eventually the edge of the
plate thickened and $\alpha _{prism}$ dropped back down as it approached the
intrinsic value for a flat facet.

The supersaturation field around the growing crystal was determined by
diffusion, and we had no direct measure of $\sigma _{surf}.$ However we can
use the flat basal facet as a \textquotedblleft witness
surface\textquotedblright\ to estimate $\sigma _{surf}$ as follows. Since we
know $\alpha _{basal}(\sigma _{surf})=\exp (-\sigma _{basal,0}/\sigma
_{surf})$ from other measurements of the growth of large basal facets, we
can determine $\sigma _{surf}$ from this and the measured basal velocity
v_{h}$ using Equation \ref{hertz}. The result is plotted in the bottom panel
of Figure \ref{timeseries}. We then further assume that $\alpha
_{prism}=\exp (-\sigma _{prism}/\sigma _{surf})$ and use $v_{r}$ and
Equation \ref{hertz} to estimate $\sigma _{prism},$ with the result shown in
Figure \ref{sigmacalcs}. Although somewhat crude, this quantitative estimate
clearly confirms our hypothesis that the rapid increase in $v_{r}$ when the
thin plate formed is indicative of a rapid increase in $a_{prism}.$ Modeling
the attachment coefficient as $\alpha _{prism}=\exp (-\sigma _{prism}/\sigma
_{surf}),$ this indicates that $\sigma _{prism}\ll \sigma _{prism,0}$ on the
edge of the thin plate, as expected from the SDAK mechanism.

\begin{figure}[t]
  \centering

  \includegraphics[bb=0 0 1232 970,width=4.5in, keepaspectratio]{SigmaCalcs}
  \caption{Estimates of $\protect\sigma
_{surf}$ and $\protect\sigma _{prism}$ as the test crystal grew, using the
basal facet as a \textquotedblleft witness surface\textquotedblright\ as
described in the text.}
  \label{sigmacalcs}
\end{figure}

\begin{figure}[p]
  \centering

  \includegraphics[bb=0 0 1075 1336,width=5.5in,height=6.84in,keepaspectratio]{timeseries1}
  \caption{A second example crystal grown
at -15 C, showing a more gradual increase in $\protect\sigma _{\infty }$
with time during the early stages of growth.}
  \label{timeseries1}
\end{figure}

\begin{figure}[p]
  \centering

  \includegraphics[width=5in,height=7.25in,keepaspectratio]{hysteresis2a}
  \caption{Top: A plot of the basal growth
velocity $v_{h}$ and prism growth velocity $v_{r}$ of the test crystal shown
in Figure \protect\ref{timeseries1}. Bottom:\ Three images of the crystal at
different times during its growth.}
  \label{hysteresis2}
\end{figure}

A key point in this analysis is an examination of what happened to this test
crystal around $t=166.$ What caused $v_{r}$ to suddenly increase while
v_{h} $ remained essentially constant, as seen in Figure \ref{vhvr}?
Essentially the only change that occurred was that a plate-like appendage
began growing on the top edge of the prism facet, as seen in Figure \re
{composite}. This alone would not substantially alter the diffusion field
around this small crystal, as is evidenced by the fact that $v_{h}$ remains
roughly unchanged. Since the diffusion field did not change abruptly near
t=166$, the rapid increase in $v_{r}$ means that $\alpha _{prism}$ increased
threefold while $\sigma _{surf}$ remained essentially constant. Our
conclusion, therefore, is that $\alpha _{prism}$ depends on the structure of
the prism facet, which confirms our original hypothesis of
structure-dependent attachment kinetics.

Figures \ref{timeseries1} and \ref{hysteresis2} show another example where
\sigma _{\infty }$ was increased more gradually as the plate developed. In
this case the plate initially formed at a lower $v_{h}$ compared to the
previous example, and the maximum $v_{r}$ was lower. Again the $v_{h}$-$v_{r}
$ plot shows significant hysteresis, indicative of the SDAK mechanism. Here
we see that $\alpha _{prism}$ increased by a factor of four as the thin
plate developed, while during this same time $v_{h}$ (and thus $\sigma
_{surf})$ declined slightly. Modeling $\sigma _{prism}$ as before yields a
broad dip with a minimum value of $\sigma _{prism}\approx 0.8$ at $t\approx
450$ seconds.

We observed similar plate-on-pedestal growth behaviors at other temperatures
ranging from -12 C to -17 C. The hysteresis in the $v_{h}$-$v_{r}$ plots was
most pronounced at -15 C, and we observed that the minimum $v_{h}$ needed to
produce a transition to plate-like growth was approximately two times higher
at -17 C than at the other temperatures. Overall these observations are
consistent with what we know from the morphology diagram, namely that the
growth of thin, plate-like crystals is most pronounced near -15 C, as shown
in Figure \ref{morphdiagram}.

\begin{figure}[p]
  \centering

  \includegraphics[width=4.3in, keepaspectratio]{popcomposite2a}
  \caption{Models of plate-on-pedestal
growth as $\protect\sigma _{prism}$ is varied. Shown are cross-sections of
the cylindrically symmetric crystals after 20 seconds of growth. Brightness
around the crystal is proportional to supersaturation. For each of these
models the initial crystal radius and diameter were both 5 $\protect\mu $m,
\protect\sigma _{\infty }=5$ percent, $\protect\alpha _{basal}=\exp (
\protect\sigma _{basal,0}/\protect\sigma _{surf}),$ and $\protect\alpha
_{prism}=\exp (-\protect\sigma _{prism}/\protect\sigma _{surf}).$ Each model
has the same scale, and the $\protect\sigma _{prism}=0$ plate has a final
radius of 26 $\protect\mu $m. The prism growth velocities for the final
crystals shown were $v_{r}=(250,300,450,600,650,900)$ nm/sec, from top to
bottom.}
  \label{popmodels}
\end{figure}

\section{Modeling}

It is also instructive to examine computer models of growing crystals as
\sigma _{prism}$ and other parameters are varied, in order to compare with
our observations. To this end we used a cylindrically symmetric cellular
automata model \cite{ca} for which the functions $\alpha _{basal}(\sigma
_{surf})$ and $\alpha _{prism}(\sigma _{surf})$ were inputs in the code. For
these computations we used $\alpha _{basal}(\sigma _{surf})=\exp (-\sigma
_{basal,0}/\sigma _{surf})$ with $\sigma _{basal,0}=2.0$ percent and $\alpha
_{prism}(\sigma _{surf})=\exp (-\sigma _{prism}/\sigma _{surf})$ with a
constant $\sigma _{prism}$ for each model. These models cannot reproduce the
full SDAK growth instability, or the observed hysteresis, because $\sigma
_{prism}$ was assumed to have a constant value independent of the structure
of the crystal. Therefore we used the models mainly to examine the
transition to plate-like growth as a function of $\sigma _{prism}$ and
\sigma _{\infty }.$

Figure \ref{popmodels} shows one example of a set of models as $\sigma
_{prism}$ was varied. The figure shows that the models produced thin
plate-like crystals growing from the initial pedestal only when $\sigma
_{prism}$ was quite low, confirming the low $\sigma _{prism}$ values we
inferred from our measurements. From this and additional model calculations
we reached two conclusions:\ 1) both the observed plate-on-pedestal
morphology and the measured $v_{r}$ values were approximately reproduced
with $\sigma _{prism}$ values that were consistent with those inferred from
our experimental data, and 2) the initial growth phase with no thin plate
was reproduced only with larger $\sigma _{prism}$ values, again roughly
consistent with our measured values. Although these constant-$\sigma _{prism}
$ models cannot reproduce the full behavior of our observed crystals, they
nevertheless support our overall conclusions regarding the SDAK mechanism
and the resulting edge-enhancing growth instability.

To produce a full model of the SDAK instability, it would be necessary for
\alpha _{prism}$ to vary with the structure of the crystal, as this is an
inherent feature of structure-dependent attachment kinetics. One possibility
would be to let $\alpha _{prism}(\sigma _{surf})=\exp (-\sigma
_{prism}(R)/\sigma _{surf})$ where $R$ is the radius of curvature of the
edge of the crystal plate. The model would assume $\sigma _{prism}(R=\infty
)=\sigma _{prism,0}$, where $\sigma _{prism,0}$ is the measured value for
flat prism facets. At smaller $R,$ $\sigma _{prism}$ would become smaller,
as seen in Figure \ref{sigmacalcs}. Realizing such a dynamical model would
require some guesses for the correct functional form for $\sigma _{prism}(R)
, as this function is not determined from our measurements or from
theoretical considerations. Without more constraints on this and other model
inputs, we believe a full dynamical model would be of limited use at this
time, other than as a proof-of-principle that the SDAK instability can
exhibit a hysteresis behavior that approximates what we observed. However,
additional experiments like those described here may yield a useful
measurement of $\sigma _{prism}(R)$, so additional progress in quantitative
dynamical modeling is certainly possible.

\section{Discussion}

To summarize, the data and modeling presented above provide strong evidence
for structure-dependent attachment kinetics and an edge-enhancing growth
instability at -15 C. Our ice crystal growth measurements and modeling
indicate that:

1) For the growth of flat facet surfaces near -15 C, the condensation
coefficients for both the basal and prism facets are well described by the
functional form $\alpha (\sigma )=\exp (-\sigma _{facet,0}/\sigma ),$ where
\sigma $ is the supersaturation at the surface and $\sigma _{facet,0}$ is a
constant.

2) For flat facet surfaces we have measured $\sigma _{basal,0}\approx 2.0$
percent and $\sigma _{prism,0}\approx 4.2$ percent at -15 C. Because $\alpha
_{basal,0}<\alpha _{prism,0}$, this suggests that the intrinsic growth
morphology at -15 C is a columnar prism.

3) The measurements suggest that the values of $\sigma _{basal,0}$ and
\sigma _{prism,0}$ are independent of background air pressure, at least for
the range 0.01 -- 1 bar. We also expect there would be no intrinsic pressure
dependence in the attachment kinetics from theoretical considerations, since
these pressures are too low to substantially affect the ice surface.

4) While chemical contamination of the ice surface is always present at some
level in experiments, we have found that contamination levels must be quite
high to affect ice growth \cite{contamination}. Furthermore, chemical
contamination typically suppresses the growth of thin plates at -15 C \cit
{contamination}, suggesting that contamination is not an important factor in
the observations described here.

5) Our observations of plate-on-pedestal growth show a rather abrupt
transition to plate-like growth. Using the basal facet as a
\textquotedblleft witness surface\textquotedblright\ allowed us to estimate
\sigma _{prism}$ as described above. During plate growth we found $\sigma
_{prism}\ll \sigma _{prism,0}.$

6) Computer modeling confirms that $\sigma _{prism}\ll \sigma _{prism,0}$ is
necessary for plate-like growth.

7) Hysteresis in the $v_{h}$-$v_{r}$ plots requires that $\alpha
_{prism}(\sigma _{surf})$ is not a simple single-valued function independent
of crystal morphology. Instead, $\alpha _{prism}$ must be large at the edge
of a thin plate and smaller on the edges of thick plates or on flat facet
surfaces.

Together these considerations strongly support the SDAK mechanism we
proposed earlier, and indicate the presence of an edge-enhancing growth
instability. Beginning with a crystal prism containing only large facets, we
initially have $\sigma _{prism}\approx \sigma _{prism,0}$ and $\sigma
_{basal}\approx \sigma _{basal,0}.$ At low pressures the prism and basal
facets remain large and the SDAK mechanism is absent, giving the roughly
isometric growth that has been observed. But at higher pressures, diffusion
causes the top edge of the crystal to grow a bit faster than other parts of
the crystal, so the edge sharpens slightly. This sharper edge leads to a
reduced $\sigma _{prism},$ which causes the edge growth to increase further.
A positive feedback results: as $\sigma _{prism}$ reduces, the edge grows
faster and becomes sharper, thus reducing $\sigma _{prism}$ more, and so on.
This positive feedback is the hallmark of a growth instability, and the end
result is the formation of a thin plate.

An edge-enhancing instability driven by the SDAK mechanism is a relatively
simple physical explanation that nicely explains all the observations. The
existence of a growth instability of this nature may also explain why ice
crystal morphologies change so abruptly with temperature in the morphology
diagram. If the instability tips only slightly in the direction of fast
prism growth, then thin plates will form. But if the same instability tips
toward fast basal growth, then hollow columns (which, like plates, also have
thin edges) become the preferred morphology. Even small changes in the
underlying nature of the instability can result in large changes in the
final crystal morphologies.

We do not yet understand the molecular mechanisms that produce the SDAK
effect, although structure-dependent surface premelting may be playing a
role. We also do not yet understand the temperature dependence of this
instability, and why it is especially prevalent at -15 C. Nevertheless, it
appears that this instability is responsible for the formation of thin
plates at -15 C, which is a dominant feature in the morphology diagram.
Additional observations may yield more insights into the underlying
molecular processes responsible for shaping the growth of snow crystals
under different conditions.


text:
\section{Introduction}

Physical models of many complex natural systems are, at best, ``partially'' known as conservation laws do not provide a closed system of equations. Accurate theoretical models for closing the system of conservation equations are available for homogeneous systems exhibiting time and length scale separation. Examples of accurate closures include Newtonian stress for homogeneous (Newtonian) fluids, Fick's law for mass flux in diffusion processes, and the Darcy law for fluid flux in porous media. For more complex systems, including non-homogeneous turbulence, non-Newtonian fluid flow, multiphase flow and transport in porous media, and granular materials, accurate theoretical closures are not available. Instead,  phenomenological constitutive relationships are used, which are usually accurate for a narrow range of conditions. Even when sufficiently accurate closed-form partial differential equations (PDE) models are available, (space-dependent) parameters are typically unknown.

Computational approaches for parameter and constitutive law estimation cast this inherently ill-posed problem as an optimization or a statistical inference task, usually requiring  repeated evaluation of expensive forward solvers until the parameters that minimize a given error metric are found or until space of parameter configurations satisfying available data is explored. This
  results in significant and often intractable computational cost.
  Furthermore, minimization via gradient-based methods requires computing gradients from said expensive forward models, which either requires additional computational cost or careful formulation of the adjoint problems.
  Additionally, forward modeling requires the knowledge of initial and boundary conditions, which usually are not fully known and, therefore, also must be estimated from data together with unknown parameters and constitutive relations, significantly complicating parameter estimation.

Although significant progress has been made over the last two decades involving high-order schemes for PDEs, automatic differentiation of computer code, PDE-constrained optimization, and optimization under uncertainty, parameter estimation in large-scale problems remains a significant challenge \cite{lieberman2010parameter}.

While there are a number of established methods for parameter estimation in (closed-form) PDE models,
  such as Bayesian inference and maximum a posteriori probability (MAP) estimation~\cite{stuart2010inverse,carrera2005inverse}, %
existing approaches for learning unknown physics from partially known models and data are few and not fully mature.

Unknown forms of equations at a given scale can either be found by upscaling (coarsening) known equations governing the same process at a smaller scale (e.g., \cite{whitaker2013method}), or learned from data. In this work, we are interested in the latter.

In the most general case, system dynamics can be described as
\begin{equation*}
  \frac{\partial u}{\partial t} = F(u),
\end{equation*}
where $F$ is a function or differential operator.
A recent review of methods for learning $F$ can be found in \cite{rudy2018data}. These methods include NARMAX \cite{chen1989representations}, equation-free methods \cite{kevrekidis2003equation}, Laplacian spectral analysis \cite{giannakis2012nonlinear}, and neural networks \cite{gonzalez1998identification}.
Recently, data-driven approximations of Koopman operators, including dynamic mode decomposition \cite{williams2015data}, diffusion maps \cite{williams2015data}, delay coordinates \cite{giannakis2017data,brunton2017chaos,askham2018variable}, and neural networks \cite{yeung2017learning,wehmeyer2018time, mardt2018vampnets,lusch2017deep,raissi2018deep}, have gained significant attention.
This work concerns with problems where conservation laws and other physical knowledge provide significantly more information
about the structure of the operator $F(u)$.
Specifically, we are interested in a steady state of the problem
\begin{equation}
  \label{eq:pde-problem}
\frac{\partial u(\mathbf{x},t)}{\partial t} = -\nabla \cdot (K(\mathbf{x},u) \nabla u(\mathbf{x},t)),
\end{equation}
where $K(\mathbf{x},u)$
is an unknown constitutive relationship, which is a function of space and the PDE state. The boundary conditions may or may not be known.
Among other physical phenomena, this equation describes flow in porous media~\cite{bear2013dynamics}.
  Two special cases of this problem are $K(\mathbf{x}, u) = K(\mathbf{x})$, where Eq (\ref{eq:pde-problem}) becomes a linear diffusion equation with heterogeneous diffusion coefficient $K(\mathbf{x})$; and $K(\mathbf{x}, u) = K(u)$, where Eq (\ref{eq:pde-problem}) becomes a nonlinear diffusion equation with a state-dependent coefficient $K(u)$.
We are interested in learning $K(\mathbf{x},u)=K(\mathbf{x})$ when measurements of both $u$ and $K$ are available, and $K(\mathbf{x},u)=K(u)$ when only $u$ measurements are available.

To achieve this objective, we propose a physics informed DNN method, where PDEs
and data are used to train DNN representations of the PDE states and the unknown parameters constitutive relations.
With application to Eq (\ref{eq:pde-problem}), this approach consists of defining two DNNs, one for $K(\mathbf{x},u)$ and another for $u(\mathbf{x})$, together with auxiliary DNNs obtained by substituting $K$ and $u$ into Eq (\ref{eq:pde-problem}) and using automatic differentiation to evaluate the right-hand side expression and boundary conditions.
These networks are trained simultaneously using available data.
Our work extends the physics informed DNN method proposed in~\cite{raissi2017physicsPart1,raissi2017physicsPart2} for finding unknown constants and solutions to PDEs given states observations. It is important to note that the physics informed DNN method was originally developed for time-dependent PDEs, and the extension of this method proposed herein also can be applied to time-dependent problems.

The work is organized as follows: in Section \ref{philms}, we introduce the PDE model and the physics-informed DNN approach.
In Section \ref{linearPDE}, we provide a detailed study of the accuracy and performance of the proposed approach for a linear parameter estimation problem with $K(\mathbf{x},u)=K(\mathbf{x})$.
In Section \ref{nonlinearPDE}, we demonstrate the method's accuracy for learning the non-linear constitutive relationship $K(\mathbf{x},u) = K(u)$.
Discussion and conclusions are presented in Section \ref{conclusions}.

\section{Physics-informed Deep Neural Network Approach}\label{philms}

Consider the steady-state PDE
\begin{equation}
  \label{eq:steady-pde}
  \mathcal{L}[u,K(\mathbf{x},u)] = 0, \quad \mathbf{x} \in \Omega
\end{equation}
subject to the Dirichlet and Neumann boundary conditions
\begin{align}
  \label{eq:steady-pde-bc-D}
  u(\mathbf{x}) &=g(\mathbf{x}), && \mathbf{x}\in \partial_D \Omega,\\
  \label{eq:steady-pde-bc-N}
  \mathbf{n} \cdot K(\mathbf{x},u) \nabla u(\mathbf{x}) &= q(\mathbf{x}), && \mathbf{x}\in \partial_N \Omega.
\end{align}
Here, $\mathcal{L}$ is the known differential operator, %
  $\Omega \subset \mathbb{R}^d$, $d \in [1, 3]$ is the simulation domain with the boundary $\partial \Omega$, $\partial_D \Omega$ and $\partial_N \Omega$ are the ``Dirichlet'' and ``Neumann'' parts of the boundary satisfying $\partial \Omega = \partial_D \Omega \cup \partial_N \Omega$ and $\partial_D \Omega \cap \partial_N \Omega \equiv \varnothing$, and $\mathbf{n}$ is the outward unit vector normal to $\partial \Omega$.
  We denote the state of the model by $u : \Omega \to \mathcal{U} \subseteq \mathbb{R}$ and the unknown constitutive relation by $K : \Omega \times \mathcal{U} \to \mathcal{K} \subseteq \mathbb{R}$.

We assume that $N_K$ measurements of $K$, $N_u$ measurements of $u$, $N_D$ measurements of $g$, and $N_N$ measurements of $q$ are
  collected at the locations $\{ \mathbf{x}^K_i \}^{N_K}_{i = 1}$, $\{ \mathbf{x}^u_i \}^{N_u}_{i = 1}$, $\{ \mathbf{x}^D_i \}^{N_D}_{i = 1}$, and $\{ \mathbf{x}^N_i \}^{N_N}_{i = 1}$, respectively.
  The observations are denoted by $K^{*}_i \equiv K(\mathbf{x}^K_i, u(\mathbf{x}^K_i))$ ($i = 1, \dots, N_K$), $u^{*} \equiv u(\mathbf{x}^u_i)$ ($i = 1, \dots, N_u$), $g^{*}_i \equiv g(\mathbf{x}^D_i)$ ($i = 1, \dots, N_D$), and $q^{*}_i \equiv q(\mathbf{x}^N_i)$ ($i = 1, \dots, N_N$).

To learn $K(\mathbf{x},u)$, we define the following DNNs for $u(\mathbf{x})$ and $K(\mathbf{x},u)$:
\begin{equation}
  \hat{u}(\mathbf{x};\theta) = \mathcal{NN}_u(\mathbf{x};\theta), \quad \hat{K}(\mathbf{x}, u; \gamma) = \mathcal{NN}_K(\mathbf{x}, u; \gamma),
  \label{DNNs}
\end{equation}
where $\theta$ and $\gamma$ are the DNN parameters.
  Substituting these two DNNs into the governing equation (\ref{eq:steady-pde}) and the Neumann boundary condition (\ref{eq:steady-pde-bc-N}), and evaluating the corresponding spatial derivatives via automatic differentiation, %
yields two additional ``auxiliary'' DNNs:
\begin{gather*}
  f(\mathbf{x}; \gamma, \theta) = \mathcal{L} \left [ \mathcal{NN}_u(\mathbf{x}; \theta),  \mathcal{N}\mathcal{N}_K \left ( \mathbf{x}, \mathcal{NN}_u \left (\mathbf{x}; \theta \right ); \gamma \right )  \right ] = \mathcal{NN}_f (\mathbf{x};\theta,\gamma),\\
  f_N(\mathbf{x}; \gamma, \theta) = \mathbf{n} \cdot \mathcal{NN}_K \left ( \mathbf{x},   \mathcal{NN}_u(\mathbf{x}; \theta); \gamma \right ) \nabla \mathcal{NN}_u(\mathbf{x};\theta) = \mathcal{N} \mathcal{N}_N (\mathbf{x};\theta,\gamma).
\end{gather*}
Next, we define the following loss function to train these four networks simultaneously:
\begin{equation}
  \label{loss_fn_with_BC}
  \begin{split}
    L(\theta, \gamma) &= \frac{1}{N_K}\sum \limits_{i=1}^{N_K} \left [ \hat{K}(\mathbf{x}^K_i, \hat{u}(\mathbf{x}^K_i;\theta); \gamma) - K^*_i \right]^2 + \frac{1}{N_u}\sum \limits_{i=1}^{N_u} \left [\hat{u}(\mathbf{x}^u_i;\theta) - u^*_i \right]^2 \\
    &+ \frac{1}{N_D}\sum \limits_{i=1}^{N_D} \left [ \hat{u}(\mathbf{x}^D_i; \theta) - g^{*}_i \right ]^2 +\frac{1}{N_N}\sum \limits_{i=1}^{N_N} \left [ f_N(\mathbf{x}^N_i;\gamma,\theta) - q^{*}_i \right ]^2\\
    &+ \frac{1}{N_c} \sum \limits_{i=1}^{N_c} f(\mathbf{x}^c_i;\gamma,\theta)^2.
  \end{split}
\end{equation}
The first and second terms in $L$ force the $K$ and $u$ DNNs to match the $K$ and $u$ measurements.
  The third and fourth terms enforce Dirichlet and Neumann boundary conditions.
  Finally, the fifth term enforces the PDE at $N_c$ ``collocation'' points $\{ \mathbf{x}^c_i \}^{N_c}_{i = 1}$ that can be chosen uniformly or non-uniformly over $\Omega$ depending on the problem.

The DNNs are trained, i.e., $\theta$ and $\gamma$ are found, by minimizing the loss function:
\begin{equation}
  \label{eq:min_loss}
  (\theta, \gamma) = \argmin_{\theta,\gamma} L(\theta, \gamma).
\end{equation}
The minimization is carried out using the L-BFGS-B method~\cite{byrd1995limited} together with Xavier's normal initialization scheme~\cite{glorot2010understanding}.
  We use a quasi-Newton optimizer such as L-BFGS-B instead of stochastic gradient descent \cite{zhang2013asynchronous} (a more common optimizer for DNNs) because of its superior rate of convergence and more favorable computational cost for problems with a relatively small number of observed data.

Note that the proposed method does not make any assumptions about the measurement noise.  Also, our method can be easily extended to time-dependent PDEs by defining DNNs in Eq (\ref{DNNs}) as functions of both $\mathbf{x}$ and $t$ \cite{raissi2017physicsPart1,raissi2017physicsPart2}.

In the following two sections, we apply the physics-informed DNNs to learn parameters and constitutive relationships in PDE models
of the form (\ref{eq:steady-pde})--(\ref{eq:steady-pde-bc-N}).

\section{Parameter estimation in a linear diffusion equation}
\label{linearPDE}

In this section, we consider a linear diffusion equation with unknown diffusion coefficient $K(\mathbf{x})$,
\begin{equation}
  \label{eq:linearPDE}
  \nabla \cdot (K(\mathbf{x}) \nabla u(\mathbf{x})) = 0,  \quad \mathbf{x} \equiv (x_1,x_2)^T \in (0,1) \times (0,1)
\end{equation}
subject to the Dirichlet boundary conditions
\begin{equation}
  \label{eq:linearPDE-bc-D}
  u(\mathbf{x}) = 1, \quad x_2 = 0 \quad \textnormal{and}\quad
  u(\mathbf{x}) = 0, \quad x_2 = 1
\end{equation}
and the Neumann boundary conditions
\begin{equation}
  \label{eq:linearPDE-bc-N}
  \frac{\partial u(\mathbf{x})}{\partial x_1} = 0 \quad x_1 = \{ 0, 1 \}.
\end{equation}
Among other problems, this equation describes saturated flow in heterogeneous porous media with hydraulic conductivity $K(\mathbf{x})$~\cite{bear2013dynamics}.
We assume that $N_K$ measurements of $K(\mathbf{x})$ and $N_u$ measurements of $u(\mathbf{x})$ are available:
$K^{*}_i \equiv K(\mathbf{x}^K_i)$ $(i = 1, \dots, N_K)$ and $u^{*}_i \equiv u(\mathbf{x}^u_i)$ $(i = 1, \dots, N_u)$.
We define DNNs for $K(\mathbf{x})$ and $u(\mathbf{x})$, $ \hat{K}(\mathbf{x};\gamma)=\mathcal{N} \mathcal{N}_K(\mathbf{x};\gamma) $, and $ \hat{u}(\mathbf{x};\theta)=\mathcal{N} \mathcal{N}_u(\mathbf{x};\theta), $ together with two auxiliary DNNs $ f(\mathbf{x};\gamma,\theta)=\nabla \cdot [\mathcal{N} \mathcal{N}_K(\mathbf{x};\gamma) \nabla \mathcal{N} \mathcal{N}_u(\mathbf{x};\theta)] = \mathcal{N} \mathcal{N}_f (\mathbf{x};\theta,\gamma) $ and $f_N(\mathbf{x};\theta)= \partial \mathcal{N} \mathcal{N}_u(\mathbf{x};\theta)/\partial x_2 = \mathcal{N} \mathcal{N}_N (\mathbf{x};\theta)$.
For this problem, the loss function takes the form:
\begin{equation}
  \label{eq:loss_fn_linear}
  \begin{split}
    L(\theta, \gamma) &= \frac{1}{N_K} \sum \limits_{i=1}^{N_K} \left [ \hat{K}(\mathbf{x}^K_i;\gamma) - K^{*}_i \right ]^2 + \frac{1}{N_u} \sum \limits_{i=1}^{N_u} \left [ \hat{u}(\mathbf{x}^u_i;\theta) - u^{*}_i \right ]^2 \\
        &+ \frac{1}{N_D} \sum \limits_{i=1}^{N_D} \left [ \hat{u}(\mathbf{x}_i^D;\theta) - g^{*}_i \right ]^2 +\frac{1}{N_N}\sum \limits_{i=1}^{N_N} f_N(\mathbf{x}_i^N;\gamma,\theta)^2\\
    &+ \frac{1}{N_c} \sum \limits_{i=1}^{N_c} f(\mathbf{x}^c_i; \gamma, \theta) ^2.
  \end{split}
\end{equation}
The DNNs are trained by minimizing the loss function~(\ref{eq:loss_fn_linear}) as described in Section \ref{philms}. Throughout this work, we use feed-forward networks with two hidden layers and 50 units per layer.

To demonstrate the proposed approach, we generate a reference $\ln K(\mathbf{x})$ field as a realization of the Gaussian process with zero mean and covariance function $C(\mathbf{x}, \mathbf{x}')=\sigma^2 \exp(-\| \mathbf{x}-\mathbf{x}' \|^2 / 2\lambda^2)$, with $\sigma = 1$ and $\lambda = 0.15$.
The reference $u$ is generated by solving Eqs (\ref{eq:linearPDE})--(\ref{eq:linearPDE-bc-N}) using the finite volume (FV) method with the two-point flux approximation and a cell-centered regular mesh with 1024 cells.
Figure~\ref{fig:references} presents the reference $K$ and $u$ fields.
We randomly choose $N_K$ and $N_u$ FV cell centroids as measurement locations for $K$ and $u$, respectively.
These measurement locations are shown in Figure~\ref{fig:est_with_BC}.
For evaluating the loss function and training the DNNs, we use $N_c = 1024$ uniformly distributed collocation points.

\begin{figure}[htbp]
  \centering%
  \includegraphics[width=0.45\textwidth]{k_ref_fields.png}%
  \includegraphics[width=0.45\textwidth]{u_ref_fields.png}%
  \caption{Reference $K$ (left) and $u$ (right) fields.}
  \label{fig:references}
\end{figure}

\begin{figure}[htbp]
  \centering
   \begin{subfigure}[t]{0.49\textwidth}
  \includegraphics[width=0.99\textwidth]{est_diff_eq_BC_a.png}
    \caption{ $\hat{K}(\mathbf{x})$}
 \end{subfigure}
    \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=0.99\textwidth]{est_diff_eq_BC_b.png}
      \caption{ $\hat{u}(\mathbf{x})$}
 \end{subfigure}
    \begin{subfigure}[t]{0.49\textwidth}
      \includegraphics[width=0.99\textwidth]{est_diff_eq_BC_c.png}
        \caption{ $|{K}(\mathbf{x}) - \hat{K}(\mathbf{x})|$}
 \end{subfigure}
    \begin{subfigure}[t]{0.49\textwidth}
        \includegraphics[width=0.99\textwidth]{est_diff_eq_BC_d.png}
         \caption{ $|{u}(\mathbf{x}) - \hat{u}(\mathbf{x})|$}
 \end{subfigure}
  \caption{Estimated (a) $\hat{K}(\mathbf{x})$ and (b) $\hat{u}(\mathbf{x})$ fields. The red dots indicate the observation locations. Absolute point errors in (c) $\hat{K}(\mathbf{x})$ and (d) $\hat{u}(\mathbf{x})$. $N_K = 250$, $N_u = 100$, and $N_c = 1024$.}
  \label{fig:est_with_BC}
\end{figure}

  We quantify the error between estimated and reference $K$ and $u$ fields in terms of the relative $L_2$ errors, defined as
  %
  \begin{equation*}
    \varepsilon_u = \frac{\int_\Omega \left [ u(\mathbf{x})-\hat{u}(\mathbf{x}) \right ]^2 \mathrm{d} \mathbf{x}}{\int_\Omega u^2(\mathbf{x}) \, \mathrm{d} \mathbf{x}}, \quad \varepsilon_K = \frac{\int_\Omega \left [ K(\mathbf{x})-\hat{K}(\mathbf{x}) \right ]^2 \, \mathrm{d} \mathbf{x}}{\int_\Omega K^2(\mathbf{x}) \, \mathrm{d} \mathbf{x}}.
  \end{equation*}

Figure~\ref{fig:est_with_BC} shows the estimated $\hat{K}$ and $\hat{u}$ with $N_K = 250$, $N_u = 100$, and $N_c = 1024$. %
The relative $L_2$ errors %
are $\varepsilon_u \approx 0.5$\% and $\varepsilon_K \approx 1.7$ \%.
Figure~\ref{fig:est_with_BC} also depicts the point-wise absolute error in estimates of $K$ and $u$.
The point errors in $K$ are concentrated in the upper left corner where no $K$ measurements are available with a maximum point error of approximately $30\%$.
The point errors in $u$ are much smaller (maximum error is approximately $1\%$) and more uniformly distributed throughout the domain.

\begin{figure}[htbp]
  \centering%
  \includegraphics[width=0.45\textwidth]{initialization_k_std_error.png}%
  \includegraphics[width=0.45\textwidth]{initialization_u_std_error.png}
  \caption{Mean and standard deviation of $\hat{K}$ and $\hat{u}$ obtained with 11 different network initializations using Xavier's initialization scheme as a function of $N=N_K=N_u$. The $u$ and $K$ measurement locations are fixed, and $N_c = 1024$ collocation points are used.}
  \label{fig:initialization}
\end{figure}

Next, we study the effect of DNN initialization on the estimated $K$ and $u$.
  For this purpose, we draw multiple initializations of the DNNs employing Xavier's scheme (see Section~\ref{philms}), and, for each initialization, we train the DNNs.
  We quantify the effect of initialization in terms of the mean and standard deviation of the relative $L_2$ errors of the estimated fields obtained for each initialization, i.e.,
  %
  \begin{equation*}
    \overline{\varepsilon}_{(\cdot)}  = \frac1{N_s} \sum_{i=1}^{N_s} \varepsilon_{(\cdot), i}, \quad \sigma_{\varepsilon_{(\cdot)}} =  \sqrt{\frac1{N_s} \sum_{i=1}^{N_s} \left ( \varepsilon_{(\cdot), i} - \overline{\varepsilon}_{(\cdot)} \right )^2},
  \end{equation*}
  where $\varepsilon_{(\cdot), i}$ is the relative $L_2$ error for either $K$ or $u$ for the $i$th initialization, and $N_s$ is the number of network initializations.
Figure \ref{fig:initialization} shows the mean and standard deviation of $\varepsilon_u$ and $\varepsilon_K$
as a function of $N=N_K=N_u$ obtained from $N_s=11$ different network initializations.
 The $u$ and $K$ measurement locations are the same in these simulations, and $N_c = 1024$ collocation points are used.
As before, we see that uncertainty in $\hat{u}$ is much smaller than in $\hat{K}$, i.e., $\sigma_{\varepsilon_K} >> \sigma_{\varepsilon_u}$.
For both $\hat{u}$ and $\hat{K}$, the standard deviation associated with the initialization is approximately 10 times smaller than the mean value (the coefficient of variation is $\approx 0.1$), which indicates that the initialization of the DNNs does not have a significant effect on DNN predictions.

\begin{figure}[htbp]
  \centering%
 \begin{subfigure}[t]{0.49\textwidth}
   \includegraphics[width=0.99\textwidth]{collocation_k_std_error_obs_20.png}%
  \caption{ $\overline{\varepsilon}_{K}$ and $\sigma_{\varepsilon_{K}}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
  \includegraphics[width=0.99\textwidth]{collocation_u_std_error_obs_20.png
    \caption{$\overline{\varepsilon}_{u} $ and $\sigma_{\varepsilon_{u}}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
  \includegraphics[width=0.99\textwidth]{collocation_k_std_error_obs_all.png}%
    \caption{ $\overline{\varepsilon}_{K}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
  \includegraphics[width=0.99\textwidth]{collocation_u_std_error_obs_all.png}%
    \caption{ $\overline{\varepsilon}_{u}$}
  \end{subfigure}
  \caption{(a)-(b): Mean and standard deviation of the estimated $u$ and $K$ relative errors as a function of the number of collocation points $N_c$. For a given $N_c$, the DNNs are trained 11 times for different configurations of collocation points to compute the mean and variance of the relative errors. The $u$ and $K$ measurement locations are fixed, and $N=N_K=N_u=20$. Shaded area width is equal to two standard deviations.  (c)-(d): Mean $K$ and $u$ relative errors as a function of $N_c$ and $N$.}
  \label{fig:collocation}
\end{figure}

For the DNNs training, the governing PDEs are enforced at $N_c$ collocation points. To study the effect of the number and location of collocation points, we compute the relative errors $\varepsilon_u$ and $\varepsilon_K$ as a function of the number and location of the collocation points. Figures \ref{fig:collocation}(a) and (b) show the mean and standard deviation of the relative errors versus $N_c$ for $N=N_u=N_K=20$.
For a given $N_c$, the DNNs are trained $N_s=11$ times for different locations chosen via
Latin hypercube sampling \cite{mckay1979comparison} to compute the mean and variance of the relative errors.
The error in $\hat{K}$ is about eight times larger than in $\hat{u}$.
As expected, the mean and standard deviation of $\varepsilon_u$ and $\varepsilon_K$ decrease with increasing $N_c$ until they reach asymptotic values at approximately $N_c = 300$, which is approximately 33\% of the number of grid points in these simulations.
The location of collocation points has a notable  effect on the errors, especially for a relatively small $N_c$, which is evident from relatively large coefficients of variation
$\sigma_{\varepsilon_{K}} / \overline{\varepsilon}_{K}$
and
$\sigma_{\varepsilon_{u}} / \overline{\varepsilon}_{u}$.

Figures \ref{fig:collocation}(c) and (d) show that  $\overline{\varepsilon}_u$ and $\overline{\varepsilon}_K$ asymptotically decrease for four considered values of $N$. The asymptotic values of $\overline{\varepsilon}_u$ and $\overline{\varepsilon}_K$ also decrease with increasing $N$.
For all considered $N$, imposing PDE constraints reduces the mean error in $K$ and $u$ by close to 50\%.
Here, as in Figures \ref{fig:collocation}(a) and (b), $\varepsilon_u$ is significantly smaller than $\varepsilon_K$.

\begin{figure}[htbp]
  \centering%
   \begin{subfigure}[t]{0.49\textwidth}
  \includegraphics[width=0.99\textwidth]{fixed_k_20_std_error_kfield.png}%
   \caption{$\overline{\varepsilon}_{K}$ and $\sigma_{\varepsilon_{K}}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
  \includegraphics[width=0.99\textwidth]{fixed_k_20_std_error_ufield.png}
   \caption{$\overline{\varepsilon}_{u}$ and $\sigma_{\varepsilon_{u}}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
  \includegraphics[width=0.99\textwidth]{fixed_k_all_mean_error_kfield.png}%
    \caption{$\overline{\varepsilon}_{K}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
  \includegraphics[width=0.99\textwidth]{fixed_k_all_mean_error_ufield.png}%
    \caption{$\overline{\varepsilon}_{u}$}
 \end{subfigure}
  \caption{ (a) and (b): Mean relative $L_2$ error of the predicted $K$ and  $u$ as a function of the number of $u$ observations. The number of $K$ observations is 20. Shaded area width is equal to two standard deviations computed for 11 different configurations of $u$ observations. (c) and (d):  Mean relative $L_2$ error of the predicted $K$ and  $u$ as a function of the number of $u$ and $K$ observations.}
  \label{fig:fixed_k}
\end{figure}

In Figures \ref{fig:fixed_k} and \ref{fig:fixed_u}, we study how the number of $K$ observations versus the number of $u$ observations affects $\varepsilon_u$ and $\varepsilon_K$. Here, the number of collocation points is $N_c = 1024$. Figures \ref{fig:fixed_k}(a) and (b) depict the effect of $N_u$ on the mean and standard deviation of $\varepsilon_u$ and $\varepsilon_K$ for $N_K = 20$. For each $N_u$ value, we estimate $\hat{K}$ and $\hat{u}$ with $N_s=11$ different distributions of the $u$ measurement  locations generated with the Latin hypercube sampler. Then, we compute $\overline{\varepsilon}_u$, $\overline{\varepsilon}_K$, $\sigma_{\varepsilon_u}$, and $\sigma_{\varepsilon_K}$. We see that $N_u$ does not significantly affect $\varepsilon_K$ with $\overline{\varepsilon}_K \approx24$\% and $ \sigma_{\varepsilon_K} \approx6$\%. For $\hat{u}$, $\overline{\varepsilon}_u$ deceases from more than 3\% for $N_u=20$ to $\approx 2$\% for $N_u>100$. The standard deviation $ \sigma_{\varepsilon_u}$ is
1.5\% for $N_u=20$ and less than 1\% for $N_u>50$.

Figures \ref{fig:fixed_k} (c) and (d) show $\overline{\varepsilon}_K$ and $\overline{\varepsilon}_u$ as a function of $N_u$ for different $N_K$.
For all considered $N_K$, $\overline{\varepsilon}_K$ is practically independent of $N_u$ and decreases with increasing $N_K$. On the other hand, $\overline{\varepsilon}_u$ decreases with increasing $N_u$ and/or $N_K$.

\begin{figure}[htbp]
  \centering
     \begin{subfigure}[t]{0.49\textwidth}
 \includegraphics[width=0.99\textwidth]{fixed_u_20_std_error_kfield.png}
    \caption{$\overline{\varepsilon}_{K}$ and $\sigma_{\varepsilon_{K}}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
 \includegraphics[width=0.99\textwidth]{fixed_u_20_std_error_ufield.png}
    \caption{$\overline{\varepsilon}_{u}$ and $\sigma_{\varepsilon_{u}}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
 \includegraphics[width=0.99\textwidth]{fixed_u_all_mean_error_kfield.png}
    \caption{$\overline{\varepsilon}_{K}$}
 \end{subfigure}
 \begin{subfigure}[t]{0.49\textwidth}
 \includegraphics[width=0.99\textwidth]{fixed_u_all_mean_error_ufield.png}
     \caption{$\overline{\varepsilon}_{u}$}
 \end{subfigure}
 \caption{(a) and (b): Mean relative $L_2$ error of the predicted $K$ and  $u$ as a function of the number of $K$ observations. The number of $u$ observations is 20. Shaded area width is equal to two standard deviations computed for 11 different configurations of $K$ observations. (c) and (d):  Mean relative $L_2$ error of the predicted $K$ and  $u$ as a function of the number of $u$ and $K$ observations.}
  \label{fig:fixed_u}
\end{figure}

Figures \ref{fig:fixed_u}(a) and (b) reveal that all  $\overline{\varepsilon}_u$, $\overline{\varepsilon}_K$, $\sigma_{\varepsilon_u}$, and $\sigma_{\varepsilon_K}$ decrease with increasing $N_K$ for fixed $N_u$. Alternatively, Figures \ref{fig:fixed_u}(c) and (d) demonstrate that $N_u$ has a relatively minor effect on  $\overline{\varepsilon}_u$ and $\overline{\varepsilon}_K$. In all considered cases, $\overline{\varepsilon}_u$ is almost an order of magnitude smaller than $\overline{\varepsilon}_K$.

The main conclusion to be drawn from  Figures \ref{fig:fixed_k} and \ref{fig:fixed_u} is that $K$ measurements are more important than $u$ measurements for reducing error in $\hat{K}$ and $\hat{u}$. This may be explained partially by the fact that $u(\mathbf{x})$ is much smoother than $K(\mathbf{x})$, and a relatively small number of $u$ measurements are needed to describe the $u$ field. Beyond this number (approximately 50 for this example), additional $u$ measurements do not have a significant effect on $\overline{\varepsilon}_u$ and $\overline{\varepsilon}_K$.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.75\textwidth]{mapvsnn}
  \caption{Comparison between $L_2$ error in the MAP estimate of $K$ and mean $L_2$ error $\overline{\varepsilon}_K$  in the DNN $K$ estimate.}
  \label{fig:mapvsnn}
\end{figure}

Finally, we compare our approach with the MAP method.
To regularize the $K$ estimate, we penalize the norm of the $K$ gradient, which prefers a smoother estimate \cite{barajassolano-2014-linear}.
For this regularizer and the regular FV discretization with $M = 1024$ cells, the MAP estimate $\hat{\mathbf{k}} \in \mathbb{R}^M$ is defined as the vector of cell-centered values of $K$ computed as the solution to the minimization problem

\begin{equation}
  \begin{aligned}
    & \hat{\mathbf{k}} = \argmin_{\mathbf{k}} & & \| \mathbf{u}^{*} - \mathbf{H}_u \mathbf{u} ||^2_2 + ||\ln \mathbf{k}^{*} - \mathbf{H}_K \ln \mathbf{k} ||_2^2+ \gamma \| \mathbf{L} \ln \mathbf{k} \|^2_2\\
    & \text{subject to } & & \mathbf{l}(\mathbf{u}, \mathbf{k}) = \bm{0},
  \end{aligned}
  \label{MAP}
\end{equation}

where $\mathbf{k}^{*} \equiv (K^{*}_1, \dots, K^{*}_{N_K})^{\top}$ and $\mathbf{u}^{*} \equiv (u^{*}_1, \dots, u^{*}_{N_u})^{\top}$ are the vectors of $K$ and $u$ observations, respectively; $\mathbf{H}_K \in \mathbb{R}^{N_K \times M}$ and $\mathbf{H}_u \in \mathbb{R}^{N_u \times M}$ are observation operators; $\mathbf{u}$ is the discretized state; $\mathbf{l}(\mathbf{u}, \mathbf{k}) = \bm{0}$ is the discretized problem~(\ref{eq:linearPDE})--(\ref{eq:linearPDE-bc-N}); $\mathbf{L}$ is the discrete gradient operator; and $\gamma<1$ is a problem-dependent coefficient.
The minimization problem (\ref{MAP}) is solved via the Levenberg-Marquardt (LM) algorithm~\cite{barajassolano-2014-linear}.

For the considered problems, we find that the smallest $L_2$ error is obtained with $\gamma=10^{-6}$.
Figures \ref{fig:mapvsnn} and \ref{fig:visualmap} compare the two methods. Figure \ref{fig:mapvsnn} presents $\varepsilon_K$ as a function of $N=N_K=N_u$ for $\hat{K}$ found from the two methods. The physics informed DNNs produce $\hat{K}$ with smaller errors for all considered $N$.

\begin{figure}[htbp]
  \centering
 \includegraphics[width=0.30\textwidth]{nn5050}
 \includegraphics[width=0.30\textwidth]{map5050}
 \caption{Comparison between the DNN prediction of $K$ (left) and the MAP prediction of $K$ (right) using 50 observations of both $u$ and $K$.}
  \label{fig:visualmap}
\end{figure}
Figure \ref{fig:visualmap} depicts $\hat{K}$ obtained from the two methods for $N=50$. The $\hat{K}$ field estimated from the physics informed DNNs is significantly smoother (and closer to the reference $K$ shown in Figure \ref{fig:references}) than one estimated from MAP,  even though the $\varepsilon_K$ errors for the two methods are relatively  similar: $19\%$ for NN versus $22\%$ for MAP. The tent-like character  of the MAP prediction in Figure \ref{fig:visualmap} stems from the discrepancy with respect to observations is penalized more than the smoothness of the estimate (result of a relatively small $\gamma$). We find that larger $\gamma$ results in a smoother field, but also in a larger prediction error.

Both the MAP and physics-informed DNN methods involve an objective-function minimization.
In addition, MAP estimation via gradient-based optimization algorithms (such as the LM algorithm) requires computing the gradient of the predicted observations, $\mathbf{H}_u \mathbf{u}$ and $\mathbf{H}_K \ln \mathbf{k}$ with respect to $\mathbf{k}$.
For an FV discretization, this is done via the discrete adjoint method.
The total cost for each iteration of the LM algorithm is one forward solution of the PDE problem to evaluate the objective function and one adjoint solution to compute the gradient.
Therefore, MAP requires careful discretization of the PDE problem and formulation and solution of the corresponding adjoint problem.
In contrast, in the physics-informed DNN method, both spatial derivatives and the gradients with respect to DNN parameters are computed via automatic differentiation, and the methodology does not require solving the PDE problem or formulating an adjoint problem.

Finally, significant gains can be achieved in the physics informed DNN method performance by employing GPU accelerators.

\section{Nonlinear diffusion equation}\label{nonlinearPDE}

In this section, we consider a nonlinear diffusion equation with unknown state-dependent diffusion coefficient $K(u)$,
\begin{equation}
  \label{eq:nonlinearPDE}
  \nabla \cdot [K(u)\nabla u(\mathbf{x})] =0, \quad (x_1,x_2) \in (0,L_1) \times (0,L_2)
\end{equation}
subject to the boundary conditions
\begin{align}
  \label{eq:nonlinearPDE-bc-D}
  u(\mathbf{x}) &= u_0, &&  x_1=L_1, \\
  \label{eq:nonlinearPDE-bc-N-1}
  -K(u) \frac{\partial u(\mathbf{x})}{\partial x_1} &= q, && x_1 = 0\\
  \label{eq:nonlinearPDE-bc-N-2}
  \frac{\partial u(\mathbf{x})}{\partial x_2} &= 0, && x_2 = \{ 0, L_2 \}.
\end{align}

This equation describes a two-dimensional horizontal unsaturated flow (flow of water and air) in a homogeneous porous medium, where $u(\mathbf{x})$ is the water pressure and $K(u)$ is the pressure-dependent partial conductivity of the porous medium~\cite{bear2013dynamics}.
In practice, $K(u)$ is difficult to measure directly.
Therefore, this work assumes that no measurements of $K(u)$ are available, and only $N_u$ measurements of $u$ are given.

We define two DNNs for unknown $K(u)$ and $u(\mathbf{x})$,
\begin{equation}
  \hat{u}(\mathbf{x};\theta) = \mathcal{NN}_u(\mathbf{x}; \theta), \quad \hat{K}(u; \gamma) = \mathcal{NN}_K(u; \phi)
\end{equation}
  and two auxiliary DNNs obtained by substituting the DNNs for $K$ and $u$ into (\ref{eq:nonlinearPDE}), (\ref{eq:nonlinearPDE-bc-N-1}), and~(\ref{eq:nonlinearPDE-bc-N-2}),
\begin{gather*}
  f(\mathbf{x}; \theta, \gamma) = \nabla \cdot \left [ \mathcal{NN}_K \left ( \mathcal{NN}_u(\mathbf{x}; \theta); \gamma \right ) \nabla  \mathcal{NN}_u(\mathbf{x}; \theta)  \right ] = \mathcal{NN}_f (\mathbf{x}; \theta, \gamma), \\
  \mathbf{f}_N(\mathbf{x}; \phi, \gamma) = -\mathcal{NN}_K \left ( \mathcal{NN}_u(\mathbf{x}; \theta); \gamma \right ) \nabla \mathcal{NN}_u(\mathbf{x}; \theta) = \mathcal{N} \mathcal{N}_N (\mathbf{x}; \theta, \gamma),
\end{gather*}
  where $\mathbf{f}_N$ is a vector DNN with two DNN components $f^{(x_1)}_N$ and $f^{(x_2)}_N$.
Then, the loss function becomes
\begin{equation*}
  \begin{split}
    \mathcal{L}(\theta, \gamma) &=  \frac{1}{N_u} \sum \limits_{i=1}^{N_u} \left [ \hat{u}(\mathbf{x}^u_i; \theta) - u^*_i \right ]^2 \\
    &+ \frac{1}{N_c} \sum \limits_{i=1}^{N_c} f(\mathbf{x}_i; \theta, \gamma)^2 + \frac{1}{N_D} \sum \limits_{i=1}^{N_D} \left [ \hat{u}(\mathbf{x}_i^D; \theta) - u_0 \right ]^2\\
    &+ \frac{1}{N_N^{(x_1)}} \sum \limits_{i=1}^{N_N} \left [ f_N^{(x_1)}(\mathbf{x}_i^{N, x_1}; \phi, \gamma) - q \right ]^2 + \frac{1}{N_N^{(x_2)}} \sum \limits_{i=1}^{N_N} \left [f_N^{(x_2)}(\mathbf{x}_i^{N,x_2}; \phi, \gamma) \right]^2,
  \end{split}
  \label{loss_fn_Richards}
\end{equation*}
where
$\mathbf{x}_i^{N,x_1}$ $(i=1,...,N_N^{(x_1)})$ are the collocation points on the Neumann boundary $(x_1=0, x_2)$ and $\mathbf{x}_i^{N,x_2}$ $(i=1,...,N_N^{(x_2)})$ are the collocation points on the Neumann boundaries $(x_1, x_2=0)$ and $(x_1, x_2=L_2)$.

This model is tested with data generated using the Subsurface Transport Over Multiple Phases (STOMP) code \cite{white1995modeling} with the van Genuchten model \cite{van1980closed} for the $K(u)$ function
\begin{gather}
  \label{vGK}
  K(s(u)) = K_s s^{\frac12} \left[1- \left(1-s^{\frac{1}{m}} \right)^m \right]^2,\\
  \label{vGs}
  s(u) = \left \{ 1+\left [ \alpha (u_g-u ) \right]^{\frac{1}{1-m}} \right \}^{-m}.
\end{gather}
Here, $K_s$ is the saturated hydraulic conductivity, $u_g = \frac{P_g}{\rho g}$,  $P_g$ is the air pressure, $\rho$ is the density, $g$ is gravity, and $\alpha$ and $m $ are the van Genuchten parameters. The following parameter values are used in the STOMP simulation: $u_0 = -10$ m, $\alpha = 0.1$, $m = 0.469$, $q=8.25 \times 10^{-5}$ m/s, $u_g = 0$, and $K_s = 8.25 \times 10^{-4}$ m/s.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.75\textwidth]{nonlinear2d_u.png}
  \caption{(Left) Referenced $u(\mathbf{x})$ field generated using STOMP with the van Genuchten model for $K(u)$ and the locations of $u$ observations. (Right) Absolute error in the  $u(\mathbf{x})$ field estimated with the physics informed DNN.}
  \label{fig:unsaturated}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.75\textwidth]{vanGenuchten.png}
  \caption{Comparison of the estimated $\hat{K}(u)$ and the reference $K(u)$ given by the van Genuchten model.}
  \label{fig:vanGenuchten}
\end{figure}

Figure \ref{fig:unsaturated} shows the reference $u(\mathbf{x})$ field generated with STOMP and the assumed locations of $u$ measurements. Figure \ref{fig:vanGenuchten} presents the estimated $\hat{K}(u)$ function and the reference $K(u)$ function given by Eqs  (\ref{vGK}) and (\ref{vGs}). It is evident that the physics informed DNN method provides an accurate estimate of unknown $K(u)$ without any direct measurements of $K$ as a function of $u$.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.4\textwidth]{nonlinear2d_u_noisy.png}
   \includegraphics[width=0.5\textwidth]{nonlinear2d_k_noisy.png}
   \caption{ (Left) Absolute error in $\hat{u}(\mathbf{x})$ in the presence of measurement noise. (Right) Comparison of the reference $K(u)$ given by the van Genuchten model and estimated $\hat{K}(u)$ in the presence of measurement noise.
     }
  \label{fig:unsaturated_noise}
\end{figure}

Finally, we examine the robustness of the physics informed DNN in the presence of observation noise. We use the exact same setup as before except we add 1\% random noise to the values of observed $u$. Figure \ref{fig:unsaturated_noise} shows the difference between predicted and referenced $u(x)$ and $K(u)$.  The added noise increases maximum error in the reconstructed $u$ from 0.002 to 0.012, but the accuracy of the reconstructed $K(u)$ practically does not change. Note that the relative $L_2$ prediction errors for the ``noisy'' case are quite small, including $7.4\times10^{-4}$ for $u$ and $6.4\times10^{-3}$ for $K$. For comparison, the relative $L_2$ prediction errors for the ``noiseless'' case are $5.8\times10^{-5}$ for $u$ and $5.9\times10^{-3}$ for $K$.

\section{Conclusions}\label{conclusions}

In this work, we have presented a physics informed DNN method for estimating parameters and unknown physics (constitutive relationships) in PDE models. The proposed method uses both PDEs and measurements to train DNNs to approximate unknown parameters and constitutive relationships, as well as states (the PDE solution). Physical knowledge increases the accuracy of DNN training with small data sets and affords the ability to train DNNs when no direct measurements of the functions of interest are available.

We have tested this method for estimating an unknown space-dependent diffusion coefficient in a linear diffusion equation and an unknown constitutive relationship in a non-linear diffusion equation. For the parameter estimation problem, we assume that partial measurements of the coefficient and state are available and have demonstrated that the proposed method is more accurate than the state-of-the-art MAP parameter estimation method. For the non-linear diffusion PDE model with unknown constitutive relationship (state-dependent diffusion coefficient), the proposed method has proven able to accurately estimate the non-linear diffusion coefficient without any measurements of the diffusion coefficient and with measurements of the state only.
We also have demonstrated that adding physics constraints in training DNNs could increase the accuracy of DNN parameter estimation by as much as 50\%.

Parameter estimation is an ill-posed problem, and standard parameter estimation methods, including MAP, rely on regularization. In this work, we have trained DNNs for unknown parameters without regularizing  the estimated parameter field or unknown function. In the absence of regularization, we have found that the estimates of the parameter and state mildly depend on the DNN Xavier initialization scheme. For the considered problem, the uncertainty (standard deviation) and mean error decreased with an increasing number of parameter measurements. The coefficient of variation of the relative error (the ratio of the relative error standard deviation to the mean value) was found to be approximately $0.1$.
In future research, we will investigate the effect of regularization in the physics informed DNN method on  parameter estimation accuracy.

\section{Acknowledgments}
This research was partially supported by the U.S. Department of Energy (DOE) Advanced Scientific Computing (ASCR) and Biological Environmental Research (BER) offices and the  Pacific
Northwest National Laboratory (PNNL) "Deep Learning for Scientific Discovery Agile Investment program. PNNL is operated by Battelle for the DOE under Contract DE-AC05-76RL01830.

\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

Let $G\subset\CC^n$ be a domain. For $a\in G$ let
\begin{align*}
\CM_G(a):&=\{|f|: f\in\CO (G,\DD),\;f(a) =0\},\\
\CS^{(p)}_G(a):&=\{\sqrt[p]u:\;\; u:G\too[0,1): \log u\in\PSH(G),\\
&\qquad u\in\Cal C^p(\{a\}),\;\exists_{C>0}: u(z)\leq C\|z-a\|^p,\;z\in G\},\quad p\in\NN,\\
\CK_G(a):&=\{u:\;\;u:G\too [0,1): \log u\in\PSH(G),\\
&\qquad\exists_{C>0}: u(z)\leq C\|z-a\|,\;z\in G\},
\end{align*}
where $\DD\subset\CC$ stands for the unit disc, $\CO(G,\DD)$, resp.~$\PSH(G)$ denote the set of all holomorphic functions on $G$ having values in $\DD$, resp.~the set of all plurisubharmonic functions on $G$, and ``$u\in\Cal C^p(\{a\})$'' means that $u$ is of class $\Cal C^p$ in a neighborhood of $a$  (cf.~\cite{JarPfl2013}, \S\;4.2). Note that $\CS^{(1)}_G(a)$ is different from $\CK_G(a)$ (see Remark \ref{Rem1}\eqref{Rem1bb}). Put
\begin{multline*}
\CS_G(a):=\CS^{(2)}_G(a)=\{\sqrt{u}:\;\; u:G\too[0,1): \log u\in\PSH(G),\\
u\in\Cal C^2(\{a\}),\;u(0)=0\}.
\end{multline*}
Obviously, $\CM_G(a)\subset\CS_G(a)\subset\CK_G(a)$ and $\CS^{(p)}_G(a)\subset\CK_G(a)$, $p\in\NN$.
If $\CF\in\{\CM, \CS^{(p)}, \CK\}$, then we define:
\begin{align*}
d_G^\CF(a,z):&=\sup\{v(z):\; v\in\CF_G(a)\},\quad a, z\in G,\\
\delta_G^\CF(a;X):&=\sup\Big\{\limsup_{\lambda\to0}\frac{v(a+\lambda X)}{|\lambda|}: v\in\CF_G(a)\Big\},\quad a\in G,\;X\in\CC^n.
\end{align*}
For $\CF\in\{\CM, \CS, \CK\}$ the families $(d^\CF_G)_G$ and $(\delta^\CF_G)_G$ are \emph{holomorphically contractible}, i.e.

\bu $d^\CF_{\DD}(0,t)=t$, $t\in[0,1)$, \quad $\delta^\CF_{\DD}(0;1)=1$;

\bu for any domains $G\subset\CC^n$, $D\subset\CC^m$ and for any holomorphic mapping $F:G\too D$ we have
\begin{align}
d_D^\CF(F(a),F(z))&\leq d_G^\CF(a,z),\quad a,z\in G, \label{C21}\\
\delta_D^\CF(F(a);F'(a)(X))&\leq\delta_G^\CF(a;X),\quad a\in G,\;X\in\CC^n. \label{C22}
\end{align}

In particular, the families $(d^\CF_G)_G$ and $(\delta^\CF_G)_G$ are invariant under biholomorphic mappings.

If $\CF=\CM$, then we get the \emph{M\"obius pseudodistance} $\bm_G:=d_G^\CM$ and the \emph{Ca\-ra\-th\'e\-odo\-ry--Reiffen pseudometric}
$\bgamma_G:=\delta_G^\CM$. It is known that
\begin{gather}
\bgamma_G(a,z)=\lim_{\lambda\to0}\frac{\bm_G(a,a+\lambda X)}{|\lambda|}=\max\{|f'(z)(X)|: f\in\CO (G,\DD),\;f(a) =0\}.\label{eq1}
\end{gather}

If $\CF=\CS$, then we get the \emph{Sibony function} $\bbs_G:=d_G^\CS$ and the \emph{Sibony pseudometric} $\bS_G:=\delta_G^\CS$.
It is known that
$$
\bS_G(a;X)=\sup\{\sqrt{\CL u(a;X)}: u\in\CS_G(a)\},
$$
where $\CL u(a;X):=\sum_{j,k=1}^n\pder{^2u}{z_j\partial\overline z_k}(a)X_j\overline X_k$ is the \emph{Levi form} (cf.~\cite{JarPfl2013}, Proposition 4.2.16).
In particular, $\bS_G(a;\cdot)$ is a $\CC$-seminorm.

If $\CF=\CK$, then we get the \emph{pluricomplex Green function} $\bg_G:=d_G^\CK$ and the \emph{Azukawa pseudometric}
$\bA_G:=\delta_G^\CK$. It is known that $\bg_G(a,\cdot)\in\CK_G(a)$, $\log\bA_G(a;\cdot)\in\PSH(\CC^n)$, and
\begin{gather}
\bA_G(a;X)=\limsup_{\lambda\to0}\frac{\bg_G(a,a+\lambda X)}{|\lambda|}\quad\text{(cf.~\cite{JarPfl2013}, Lemma 4.2.3).}
\label{eq2}
\end{gather}

If $\CF=\CS^{(p)}$, $p\neq2$, then we get the \emph{higher order Sibony function} $\bbs^{(p)}_G:=d^{\CS^{(p)}}_G$
and the \emph{higher order Sibony pseudometric} $\bS^{(p)}_G:=\delta^{\CS^{(p)}}_G$.

The basic properties of $\bm_G$, $\bgamma_G$, $\bg_G$,
and $\bA_G$ are well understood. In contrast to that, little is known on properties of $\bS_G$ and almost nothing on $\bbs^{(p)}_G$, $p\in\NN$, and
$\bS^{(p)}_G$, $p\neq2$.

The main aim of this note is to show that the basic properties of $\bbs^{(p)}_G$ and $\bS^{(p)}_G$ differ essentially from the corresponding properties of $\bm_G$, $\bg_G$, $\bgamma_G$, and $\bA_G$.

\section{Holomorphic contractibility}

\begin{Remark}\label{Rem1}
\begin{enumerate}[(a)]
\item\label{Rem1a} $\CS^{(p)}_G(a)=\{\sqrt[p]u:\;u:G\too[0,1): \log u\in\PSH(G),\;u\in\Cal C^p(\{a\}),\;\ord_au\geq p\}$, where $\ord_au$ denotes the order of zero of $u$ at $a$.

\item\label{Rem1b} In view of the Taylor formula we have
$$
\bS^{(p)}_G(a;X)=\sup\big\{(\tfrac1{p!}|u^{(p)}(a)(X)|)^{1/p}: \sqrt[p]{u}\in\CS^{(p)}_G(a)\big\},\quad a\in G,\;X\in\CC^n,
$$
where $u^{(p)}(a):\CC^n\too\RR$ stands for the $p$-th Fr\'echet differential of $u$ at $a$.

\item\label{Rem1bb} In view of \eqref{Rem1b} we get $\bS^{(p)}_G(a;\cdot)\equiv0$ for $p$ odd. In particular, $\bS^{(1)}_\DD(0;1)=0<1=\bA_\DD(0;1)$.

\item\label{Rem1c} $\bbs^{(p)}_G\leq\bg_G$, $\bS^{(p)}_G\leq\bA_G$. In particular,
$\bbs^{(p)}_{\DD}(0,\lambda)\leq\bg_{\DD}(0,\lambda)=|\lambda|$, $\bS^{(p)}_{\DD}(0;1)\leq\bA_{\DD}(0;1)=1$.

\item\label{Rem1d} If $\bg^{p+\eps}_G(a,\cdot)\in\Cal C^p(\{a\})$ for $0<\eps\ll1$, then $\bg^{1+\eps/p}_G(a,\cdot)\in\CS^{(p)}_G(a)$. Consequently, $\bbs^{(p)}_G(a,\cdot)=\bg_G(a,\cdot)$.
In particular, $\bbs^{(p)}_{\DD}(0,\lambda)=|\lambda|$, $\lambda\in\DD$.

\item\label{Rem1e} If $\bg^{2p}_G(a,\cdot)\in\Cal C^{2p}(\{a\})$,  then $\bS^{(2p)}_G(a;\cdot)=\bA_G(a;\cdot)$.
In particular,  $\bS^{(2p)}_\DD(0;1)=1$.

\item\label{Rem1f} If $F:G\too D$ is holomorphic, then $v\circ F\in\CS^{(p)}_G(a)$ for every $v\in\CS^{(p)}_D(F(a))$. Consequently,
the family $(\bbs^{(p)}_G)_G$ (resp.~$(\bS^{(p)}_G)_G$) satisfies \eqref{C21} (resp.~\eqref{C22}).

\item\label{Rem1g} The families $(\bbs^{(p)}_G)_G$ and $(\bS^{(2p)}_G)_G$ are holomorphically contractible. They will be the main objects of our investigation in the sequel.

\item\label{Rem1h} $\bm_G\leq\bbs^{(p)}_G\leq\bg_G$ and $\bgamma_G\leq\bS^{(2p)}_G\leq\bA_G$.
\end{enumerate}
\end{Remark}

\section{Upper semicontinuity}

It is known that for $\CF\in\{\CM, \CK\}$ the functions $G\times G\in(z,w)\tuu d_G^\CF(z,w)$ and $G\times\CC^n\ni(z,X)\tuu\delta_G^\CF(z;X)$
are upper semicontinuous (cf.~\cite{JarPfl2013}, Propositions 2.6.1, 2.7.1(c), 4.2.10(g,k)).
We will prove that in general the functions $\bbs^{(p)}_G(\cdot,z^0)$ and $\bS^{(2p)}_G(\cdot;X^0)$ are not upper semicontinuous
(Examples \ref{Exam1}, \ref{Exam2}).

Recall that $\bS_G(a;\cdot)$ is a seminorm and therefore it is continuous.
We do not know whether the functions $\bbs_G(a,\cdot)$, $p\in\NN$, and $\bS^{(2p)}_G(a;\cdot)$, $p\geq2$, are upper semicontinuous.

\begin{Example}[cf.~\cite{JarPfl2013}, Example 4.2.18]\label{Exam1}
Let
$$
G:=\{(z_1,z_2,z_3)\in\CC^3: |z_1|e^{\phi(z_2,z_3)}<1\}
$$
with
$$
\phi(\xi,\eta):=\sum_{k=1}^\infty\lambda_k\log\Big(\frac{|\xi-a_k|^2+|\eta|}k\Big), \quad(\xi,\eta)\in\CC^2,
$$
where $(a_k)^\infty_{k=1}\subset\DD\setminus\{0\}$ is a dense subset of $\DD$ and $(\lambda_k)_{k=1}^\infty\subset(0,1]$
are chosen so that $\phi(0,0)>-\infty$ and $\phi\in\Cal C^\infty(\CC\times\CC_\ast)$, where $\CC_\ast:=\CC\setminus\{0\}$. Note that $G$ is a pseudoconvex Hartogs domain.

Let $c_t:=(0,0,t)\in G$, $t>0$, $z^0:=(b,0,0)\in G$ with $b\neq0$, and let $X^0:=(1,0,0)$.
We will show that
\begin{gather*}
\bbs^{(p)}_G(0,z^0)=0<|b|e^{\phi(0,0)}\leq\bbs^{(p)}_G(c_t,z^0),\\
\bS^{(2p)}_G(0;X^0)=0<e^{\phi(0,0)}\leq\bS^{(2p)}_G(c_t;X^0),\quad 0<t\ll1,
\end{gather*}
which shows that \emph{the functions $\bbs^{(p)}_G(\cdot,z^0)$ and $\bS^{(2p)}_G(\cdot;X^0)$ are not upper semicontinuous at $0$.}

Indeed, the function $G\ni(z_1,z_2,z_3)\overset{v}\tuu(|z_1|e^{\phi(z_2,z_3)})^{1+\eps/p}$ belongs to $\CS^{(p)}_G(c_t)$ for all $\eps>0$ and $t>0$. Hence,
$\bbs^{(p)}_G(c_t,z^0)\geq|b|e^{\phi(0,0)}>0$.
Analogously, the function $G\ni(z_1,z_2,z_3)\overset{v}\tuu|z_1|e^{\phi(z_2,z_3)}$ belongs to $\CS^{(2p)}_G(c_t)$ for all $t>0$.
Hence $\bS^{(2p)}_G(c_t;X^0)\geq\limsup_{\lambda\to0}\frac{v(c_t+\lambda X^0)}{|\lambda|}=e^{\phi(0,t)}\geq e^{\phi(0,0)}>0$.

On the other hand, let $\sqrt[p]{u}\in\CS^{(p)}_G(0)$ (resp.~ $\sqrt[2p]{u}\in\CS^{(2p)}_G(0)$). Since $\CC\times\{a_k\}\times\{0\}\subset G$, we get $u(z_1,a_k,0)=\const(k)$, $z_1\in\CC$, $k\in\NN$. Since $\{0\}\times\CC\times\{0\}\subset G$,
we get $u(0,z_2,0)=\const=u(0)=0$, $z_2\in\CC$. Thus $u(z_1,a_k,0)=0$, $z_1\in\CC$, $k\in\NN$.
Since $u\in\Cal C^p(\{0\})$ (resp.~$u\in\Cal C^{2p}(\{0\})$) we conclude that $u=0$ in $U\times\{0\}$, where $U$ is a neighborhood of $(0,0)$.
Since $\log u\in\PSH(G)$, we get $u(z_1,z_2,0)=0$ for all $(z_1,z_2,0)\in G$. Consequently, $\bbs^{(p)}_G(0,z^0)=0$ (resp.~$\bS^{(2p)}_G(0;X^0)=0$).
\end{Example}

\begin{Example}
In view of Example \ref{Exam1}, one could expect that perhaps the families $(\bbs^{(p)\ast}_G)_G$ and/or $(\bS^{(2p)\ast}_G)_G$ are holomorphically contractible, where $\bbs^{(p)\ast}_G:=(\bbs^{(p)}_G)^\ast$, $\bS^{(2p)\ast}_G:=(\bS^{(2p)}_G)^\ast$, and ${}^\ast$ denotes the upper semicontinuous regularization. We will prove that unfortunately they are not holomorphically contractible.

Keep the notation from Example \ref{Exam1}. Let
$$
D:=\{(z_1,z_2)\in\CC^2: (z_1,z_2,0)\in G\},\quad D\ni(z_1,z_2)\overset{F}\tuu(z_1,z_2,0)\in G.
$$
Then $\bbs^{(p)\ast}_G(0,z^0)\geq\limsup_{t\to0+}\bbs^{(p)}_G(c_t,z^0)\geq|b|e^{\phi(0,0)}>0$ and
$\bS^{(2p)\ast}_G(0;X^0)\geq\limsup_{t\to0+}\bS^{(2p)}_G(c_t;X^0)\geq e^{\phi(0,0)}>0$.

On the other hand, let $w^0\in D\cap(\CC\times\DD)$ and let $\sqrt[p]{u}\in\CS^{(p)}_D(\{w^0\})$ (resp.~ $\sqrt[2p]{u}\in\CS^{(2p)}_D(\{w^0\})$).
Since $\CC\times\{a_k\}\subset D$, we get $u(z_1,a_k)=\const(k)$, $z_1\in\CC$, $k\in\NN$. Since $\{0\}\times\CC\subset D$, we get $u(0,z_2)=\const=u(0,0)$, $z_2\in\CC$. Thus $u(z_1,a_k)=\const$, $z_1\in\CC$, $k\in\NN$.
Since $u\in\Cal C^p(\{w^0\})$ (resp.~$u\in\Cal C^{2p}(\{w^0\})$) we conclude that $u=0$ in $U\times\{0\}$, where $U$ is a neighborhood of $w^0$. Hence, since $\log u\in\PSH(G)$, we get $u(z_1,z_2)=0$ for all $(z_1,z_2)\in D$. Consequently,
$\bbs^{(p)}_D=0$ on $(D\cap(\CC\times\DD))\times D$ (resp.~$\bS^{(2p)}_D=0$ on $(D\cap(\CC\times\DD))\times\CC^2$). In particular, $\bbs^{(p)\ast}_D(0,(b,0))=0$ (resp.~$\bS^{(2p)\ast}_D(0;(1,0))=0$) and therefore
\begin{align*}
\bbs^{(p)\ast}_G(F(0,0),F(b,0))&>0=\bbs^{(p)\ast}_D((0,0),(b,0)),\\
\bS^{(2p)\ast}_G(F(0,0);F'(0,0)(1,0))&>0=\bS^{(2p)\ast}_D((0,0);(1,0)).
\end{align*}
\end{Example}

\begin{Example}\label{Exam2}
For $n\geq2$ and $\alpha=(\alpha_1,\dots,\alpha_n)\in\RR^n\setminus\{0\}$ let
$$
\BD_\alpha:=\{z\in\CC^n(\alpha): |z^\alpha|:=|z_1|^{\alpha_1}\cdots|z_n|^{\alpha_n}<1\},
$$
where $\CC^n(\alpha):=\{(z_1,\dots,z_n)\in\CC^n: \forall_{j\in\{1,\dots,n\}}: (\alpha_j<0 \Longrightarrow z_j\neq0)\}$.
Note that $\BD_\alpha$ is a pseudoconvex Reinhardt domain. For $a=(a_1,\dots,a_n)\in\BD_\alpha$ define
\begin{gather*}
\Xi(a):=\{j\in\{1,\dots,n\}: \alpha_j>0,\;a_j=0\},\\
r(a):=\begin{cases} 1 \casesep \sigma(a)=0 \\ \sum\limits_{j\in\Xi(a)}\alpha_j\casesep \sigma(a)\geq1\end{cases},\quad
\sigma(a):=\#\Xi(a),\quad
\mu(a):=\min\{\alpha_j: j\in\Xi(a)\}.
\end{gather*}
Note that if $\sigma(a)=1$, then $r(a)=\mu(a)$.

The following results are known (cf.~\cite{JarPfl2013}, \S\S\;6.2, 6.3, and \cite{JarPfl2018}, Theorem 1).

\bu If $\alpha_1,\dots,\alpha_n\in\ZZ$ are relatively prime, then
\begin{gather*}
\bm_{\BD_\alpha}(a,z)=\bm_{\DD}(a^\alpha,z^\alpha),\quad
\bg_{\BD_\alpha}(a,z)=\big(\bm_{\DD}(a^\alpha,z^\alpha)\big)^{1/r},\\
\bA_{\BD_\alpha}(a;X)=\Big(\bgamma_\DD(a^\alpha;
\tfrac1{r!}\prod_{j\notin\Xi(a)}a_j^{\alpha_j}\cdot\prod_{j\in\Xi(a)}X_j^{\alpha_j})\Big)^{1/r},\quad r=r(a),\\
\bbs_{\BD_\alpha}(a,z)=\begin{cases}\bm_{\DD}(a^\alpha,z^\alpha)\casesep \sigma(a)=0\\
|z^\alpha|^{1/\mu(a)}\casesep \sigma(a)\geq1\end{cases},\\
\bS_{\BD_\alpha}(a;X)=\begin{cases}\bA_{\BD_\alpha}(a;X)\casesep \sigma(a)\leq1\\
0\casesep \sigma(a)\geq2\end{cases},\quad a,z\in\BD_\alpha,\; X\in\CC^n.\\
\end{gather*}

\bu If $\alpha\notin\RR\cdot\ZZ^n$, then
\begin{gather*}
\bm_{\BD_\alpha}\equiv0,\quad
\bg_{\BD_\alpha}(a,z)=\begin{cases} 0 \casesep \sigma(a)=0\\
|z^\alpha|^{1/r} \casesep \sigma(a)\geq1 \end{cases},\\
\bA_{\BD_\alpha}(a;X)=\begin{cases}0\casesep \sigma(a)=0\\ \big(\prod_{j\notin\Xi(a)}|a_j|^{\alpha_j}\cdot\prod_{j\in\Xi(a)}|X_j|^{\alpha_j}\big)^{1/r}\casesep \sigma(a)\geq1
\end{cases},\quad r=r(a),\\
\bbs_{\BD_\alpha}(a,z)=\begin{cases}0 \casesep \sigma(a)=0\\
|z^\alpha|^{1/\mu(a)}\casesep \sigma(a)\geq1\end{cases},\\
\bS_{\BD_\alpha}(a;X)=\begin{cases}\bA_{\BD_\alpha}(a;X)\casesep \sigma(a)=1\\
0\casesep \sigma(a)\neq1\end{cases},\quad a, z\in\BD_\alpha,\; X\in\CC^n.
\end{gather*}

In particular, if $n=3$ and $\alpha=(1,2,2)$, then
$\bbs_{\BD_\alpha}((0,0,0),z)=|z^\alpha|$ and  $\bbs_{\BD_\alpha}((1/k,0,0),z)=|z^\alpha|^{1/2}$, $k\in\NN$.
Thus, the function $\bbs_{\BD_\alpha}(\cdot, z^0)$ is not upper semicontinuous at $(0,0,0)$ for all $z^0=(z^0_1,z^0_2,z^0_3)\in\BD_\alpha$
with $z^0_1z^0_2z^0_3\neq0$.

Notice that using the above effective formulas one may easily construct many other counterexamples.
\end{Example}

\begin{Example}
Keep the notation from Example \ref{Exam2}. Assume that $\alpha_1,\dots,\alpha_n\in\RR_\ast:=\RR\setminus\{0\}$, $a_1\cdots a_s\neq0$, $a_{s+1}=\dots=a_n=0$,
$s:=n-\sigma(a)$. In particular, $\alpha_{s+1},\dots,\alpha_n>0$.

First observe that if $\sigma(a)\leq1$, then $\bg^{p+\eps}_{\BD_\alpha}(a,\cdot)\in\Cal C^p(\{a\})$ and consequently $\bbs^{(p)}_{\BD_\alpha}(a,\cdot)=\bg_{\BD_\alpha}(a,\cdot)$
(Remark \ref{Rem1}\eqref{Rem1d}). Similarly, if $\sigma(a)\leq1$, then $\bg^{2p}_{\BD_\alpha}(a,\cdot)\in\Cal C^\infty(\{a\})$ and consequently
$\bS^{(2p)}_{\BD_\alpha}(a;\cdot)=\bA_{\BD_\alpha}(a;\cdot)$ (Remark \ref{Rem1}\eqref{Rem1e}). Problems start when $\sigma(a)\geq2$.
We do not know effective formulas for $\bbs^{(p)}_{\BD_\alpha}(a,\cdot)$, $p\neq2$, and $\bS^{(2p)}_{\BD_\alpha}(a;\cdot)$, $p\geq2$. To illustrate problems we discuss some particular cases.

\medskip

\bu Assume that $\sigma(a)\geq1$ and $k_j:=\frac{p\alpha_j}{r(a)}\in\NN$, $j=s+1,\dots,n$. Then
$$
\bg^{2p}_{\BD_\alpha}(a,z)=\prod_{j=1}^s|z_j|^{2p\alpha_j/r(a)}\cdot\prod_{j=s+1}^n|z_j|^{2k_j}
$$
and consequently $\bg^{2p}_{\BD_\alpha}(a,\cdot)\in\Cal C^\infty(\{a\})$ which gives
$\bS^{(2p)}_{\BD_\alpha}(a;\cdot)=\bA_{\BD_\alpha}(a;\cdot)$. For example, if $n=2$, $\alpha=(1,1)$, $a=(0,0)$, then
$\bS^{(4k)}_{\BD_\alpha}(a;X)=|X_1X_2|^{1/2}$, $k\in\NN$.

\medskip

\bu Assume that $\sigma(a)\geq1$ and there exists a $j_0\in\{s+1,\dots,n\}$ such that $\frac{2p\alpha_{j_0}}{r(a)}\notin\NN$.  Then $\bS^{(2p)}_{\BD_\alpha}(a;\cdot)\equiv0$.

Indeed, we may assume that $j_0=n$. Let $r:=r(a)$ and let $k\in\NN_0$ be such that $k<\frac{2p\alpha_n}r<k+1$. In view of Remark \ref{Rem1}\eqref{Rem1b} we have to prove that
$u^{(2p)}(a)\equiv0$ for all $\sqrt[2p]{u}\in\CS^{(2p)}_{\BD_\alpha}(a)$. Fix such a $u$ and suppose that $u^{(2p)}(a)(X^0)\neq0$ for some $X^0\neq0$. We have
\begin{multline*}
\Big(\frac1{(2p)!}|u^{(2p)}(a)(X)|\Big)^{1/2p}\leq\bS^{(2p)}_{\BD_\alpha}(a;X)\\
\leq\bA_{\BD_\alpha}(a;X)=\Big(\prod_{j=1}^s|a_j|^{\alpha_j}\cdot\prod_{j=s+1}^n|X_j|^{\alpha_j}\Big)^{1/r}.
\end{multline*}
Write $u^{(2p)}(a)(X_1^0,\dots,X_{n-1}^0,tX_n^0)=A_dt^d+\dots+A_0$, $t\in\RR$, with $A_d\neq0$. We have $|A_dt^d+\dots+A_0|\leq\const|t|^{2p\alpha_n/r}$, $t\in\RR$. Taking $t\too\infty$ we get $d\leq k$. On the other hand, taking $t\too0$ we get $A_d=0$; a contradiction.

For example let $n=2$, $\alpha=(q,1)$,  $a=(0,0)$, where
$$
0<q\notin\{\tfrac{2p-k}k: k=1,\dots,2p-1\}\cap\{\tfrac{k}{2p-k}: k=1,\dots,2p-1\}.
$$
Then $\bS^{(2p)}_{\BD_\alpha}(a;\cdot)\equiv0$.

\medskip

\bu As a consequence we conclude that for every $s\in\{0,\dots,n-2\}$ there exists a set $C_s$ dense in $\RR_\ast^s\times\RR_{>0}^{n-s}$ ($\RR_{>0}:=(0,+\infty)$) such that for any $\alpha\in C_s$,
 $a\in\BD_\alpha\cap(\CC_\ast^s\times\{0\}^{n-s})$, and $p\in\NN$ we have $\bS^{(2p)}_{\BD_\alpha}(a;\cdot)\equiv0$.

Indeed, we may put
$$
C_s:=(\RR_\ast^s\times\RR_{>0}^{n-s})\setminus\bigcup_{\substack{p,k\in\NN:\; k<2p\\ j\in\{s+1,\dots,n\}}}\{\alpha\in\RR^n: 2p\alpha_j=k(\alpha_{s+1}+\dots+\alpha_n)\}.
$$
\end{Example}

Now we turn to discuss a special case where $G\subset\CC^n$ is a complete $n$-circled domain
(Example \ref{Exam2a}).

\begin{Example}\label{Exam2a} Let $G\subset\CC^n$ be a  complete $n$-circled domain,
i.e.~for any $z=(z_1,\dots,z_n)\in G$ and $\lambda=(\lambda_1,\dots,\lambda_n)\in\overline\DD^n$, the point $\lambda\cdot z:=(\lambda_1z_1,\dots,\lambda_nz_n)$ belongs to
$G$.

(a) Since $\bbs^{(p)}_G(0,\cdot)\leq\bg_G(0,\cdot)$ and the Green function is upper semicontinuous,
the function $\bbs^{(p)}_G(0,\cdot)$ is continuous at $0$.

\medskip

(b) The function $\bbs^{(p)}_G(0,\cdot)$ is upper semicontinuous in the domain $G\setminus\BV_0$, where
$\BV_0:=\{(z_1,\dots,z_n)\in\CC^n: z_1\cdots z_n=0\}$.

Indeed, let $M:=\{a\in G: \bbs^{(p)}_G(0,\cdot) \text{ is not upper semicontinuous at }a\}$.
Since $\bbs^{(p)}_G(0,\cdot)$ is invariant under $n$-rotations (i.e.~under mappings $G\ni z\tuu \lambda\cdot z\in G$,
$\lambda\in\TT^n$, where $\TT:=\partial\DD$), the set $M$ is also invariant under $n$-rotations.  It is known that $M$ is pluripolar,
i.e.~there exists a $v\in\PSH(\CC^n)$, $v\not\equiv-\infty$, such that $M\subset v^{-1}(-\infty)$ (cf.~\cite{Kli1991}, Theorem 4.7.6).
Suppose that $a=(a_1,\dots,a_n)\in M\setminus\BV_0$. Then $v(\lambda\cdot a)=-\infty$ for all $\lambda\in\TT^n$.
Consequently, by the maximum principle for plurisubharmonic functions,  $v(z_1,\dots,z_n)=-\infty$ for all $|z_j|\leq|a_j|$, $j=1,\dots,n$. Hence $v\equiv-\infty$; a
contradiction.

\medskip

(c) Let $a=(0,\dots,0,a_{s+1},\dots,a_n)=:(0,b)\in G\cap\BV_0$, $1\leq s\leq n-1$, $a_{s+1}\cdots a_n\neq0$.
Define $D:=\{\zeta\in\CC^{n-s}: (\underset{s\times}{\underbrace{0,\dots,0}},\zeta)\in G\}$.
Note that $D$ is a complete $(n-s)$-circled domain
with $b\in D$. Let $\Fh_D$ denote the Minkowski functional of $D$ ($\Fh_D(\zeta):=\inf\{1/t: t>0,\;t\zeta\in D\}$, $\zeta\in\CC^{n-s}$). Observe that $\Fh_D$ is continuous (because $D$ is $(n-s)$-circled).

\emph{Assume that $\bbs^{(p)}_D(0,b)=\Fh_D(b)$}. Then the function $\bbs^{(p)}_G(0,\cdot)$ is upper semicontinuous at $a$.

Indeed, let $0<R<1$ and $k>0$ be such that $\Fh_D(b)<R$ and $\|b\|<k$. Note that $\{\zeta\in D: \Fh_D(\zeta)<R,\;\|\zeta\|<k\}\subset\subset
D$. Consequently, there exists an $\eps>0$ such that $U:=\{(z',z'')\in\CC^n: \|z'\|<\eps,\;\Fh_D(z'')<R,\;\|z''\|<k\}\subset G$. Then for $z=(z',z'')\in U$ we have
$$
\bbs^{(p)}_G(0,z)\leq\bbs^{(p)}_U(0,z)\leq\bg_U(0,z)\leq\max\Big\{\frac{\|z'\|}\eps, \frac{\Fh_D(z'')}R, \frac{\|z''\|}k\Big\}.
$$
Hence
$$
\limsup_{z\to a}\bbs^{(p)}_G(0,z)\leq\limsup_{(z',z'')\to(0,b)}
\max\Big\{\frac{\|z'\|}\eps, \frac{\Fh_D(z'')}R, \frac{\|z''\|}k\Big\}=\max\Big\{\frac{\Fh_D(b)}R, \frac{\|b\|}k\Big\}.
$$
Letting $R\too1$ and $k\too+\infty$ we get $\limsup_{z\to a}\bbs^{(p)}_G(0,z)\leq\Fh_D(b)$.

On the other side, since the projection $\CC^s\times\CC^{n-s}\ni(z',z'')\tuu z''\in D$ is well-defined, we get
$\Fh_D(b)=\bbs^{(p)}_D(0,b)\leq\bbs^{(p)}_G(0,a)$.

\medskip

(d) Observe that $\bbs^{(p)}_D(0,b)=\Fh_D(b)$ in the case where $D$ is convex. If $s=n-1$, then $D$ is either a disc or the
whole $\CC$. Thus, if $s=n-1$, then  the function $\bbs^{(p)}_G(0,\cdot)$ is upper semicontinuous at each point $a\in\BV_0$
of the form $a=(0,\dots,0,a_j,0,\dots,0)\in G$.

\medskip

(e) Consequently, if $n=2$, then the function $\bbs^{(p)}_G(0,\cdot)$ is globally upper semicontinuous.

\medskip

(f) If $G$ is bounded, then the function $\bbs^{(p)}_G(0,\cdot)$ is globally upper semicontinuous.

Indeed, we proceed by induction on $n\geq2$. The case $n=2$ is solved in (e). Suppose the result is true for $n-1\geq2$.  Let $a=(a_1,\dots,a_n)\in G\cap\BV_0$ (see (b)). We may assume that $a_{n-1}\neq0$, $a_n=0$. Define $D:=\{z'\in\CC^{n-1}: (z',0)\in G\}$; $D$ is a bounded complete $(n-1)$-circled domain. Thus, by the inductive assumption, $\bbs^{(p)}_D(0,\cdot)$ is upper semicontinuous. Since $G$ is bounded, for every $0<r<1$ with $a\in rG$ there exists an $\eps>0$ such that $(rD)\times\DD(\eps)\subset\subset G$. Suppose that $G\subset\DD^n(R)$ and let $\eta>0$ be such that
$rR|\frac{z_n}{z_{n-1}}|<\eps$ for $z\in U:=\{z=(z',z_n)\in a+\DD^n(\eta): z'\in rD\}$. For $z\in U$ consider the holomorphic mapping $F_z:rD\too G$, $F_z(w):=(w,w_{n-1}\frac{z_n}{z_{n-1}})$. We have $\bbs^{(p)}_G(0,F_z(w))\leq\bbs^{(p)}_{rD}(0,w)=\bbs^{(p)}_D(0,w/r)$. In particular,
$\bbs^{(p)}_G(0,z)=\bbs^{(p)}_G(0,F_z(z'))\leq\bbs^{(p)}_D(0,z'/r)$. Thus $\limsup_{z\to a}\bbs^{(p)}_G(0,z)\leq\limsup_{z\to a}\bbs^{(p)}_D(0,z'/r)=\bbs^{(p)}_D(0,a'/r)$. Letting $r\too1-$ (and using once again the upper semicontinuity of $\bbs^{(p)}_D(0,\cdot)$ we get
$\limsup_{z\to a}\bbs^{(p)}_G(0,z)\leq\bbs^{(p)}_D(0,a')\leq\bbs^{(p)}_G(0,a)$ (cf.~(c)).

\medskip

\emph{Note that if $n\geq3$ and $D$ is unbounded, then it is not know whether the function $\bbs^{(p)}_G(0,\cdot)$ is globally upper semicontinuous}.
\end{Example}

\section{Increasing domains property}

Let $(G_k)_{k=1}^\infty$ be sequence of domains in $\CC^n$ such that $G_k\nearrow G$, i.e.~$G_k\subset G_{k+1}$, $k\in\NN$, and
let $G=\bigcup_{k=1}^\infty G_k$. It is known that if $\CF\in\{\CM, \CK\}$, then $d_{G_k}^\CF\searrow d_G^\CF$
and $\delta_{G_k}^\CF\searrow\delta_G^\CF$ (cf.~\cite{JarPfl2013}, Propositions 2.7.1(a), 4.2.10(a)). We will show that
\emph{this is not true for $\CF=\CS^{(p)}$}.

\begin{Example}\label{Exam3}
Let
$$
\phi_k(\lambda):=\sum_{s=2}^k\frac1{s^2}\log\Big|\lambda-\frac1s\Big|,\quad k\geq2,\quad
\phi(\lambda):=\sum_{s=2}^\infty\frac1{s^2}\log\Big|\lambda-\frac1s\Big|,\quad |\lambda|<\frac12.
$$
Observe that $\phi_k\in\PSH$ and $\phi_k\searrow\phi$. Moreover, $\phi_k\in\Cal C^\infty(\tfrac1k\DD)$. Define
\begin{gather*}
G_k:=\{(z_1,z_2)\in\CC^2: |z_1|<1/2,\;|z_2|e^{\phi_k(z_1)}<1\},\\ G:=\{(z_1,z_2)\in\CC^2: |z_1|<1/2,\;|z_2|e^{\phi(z_1)}<1\}.
\end{gather*}
Note that $G_k$ is a Hartogs domain in $\CC^2$, $k\geq2$, and $G_k\nearrow G$.
For each $k\geq2$ the function
$G_k\ni(z_1,z_2)\tuu (|z_2|e^{\phi_k(z_1)})^{1+\eps/p}$ belongs to $\CS^{(p)}_{G_k}((0,0))$, $\eps>0$. Hence
$\bbs^{(p)}_{G_k}((0,0),(0,z_2))\geq|z_2|e^{\phi_k(0)}\geq|z_2|e^{\phi(0)}$ for $|z_2|<e^{-\phi(0)}$.

Analogously, since the function
$G_k\ni(z_1,z_2)\tuu|z_2|e^{\phi_k(z_1)}$ belongs to $\CS^{(2p)}_{G_k}((0,0))$, we get
$\bS^{(2p)}_{G_k}((0,0);(0,X_2))\geq |X_2|e^{\phi(0)}$ for $X_2\in\CC$ and $k\geq2$.

Now let $\sqrt[p]{u}\in\CS^{(p)}_G((0,0))$. Since $\{1/s\}\times\CC\subset G$, the Liouville type theorem for subharmonic functions gives
$u(1/s,z_2)=\const(s)=:c_s$, $s\geq2$, $z_2\in\CC$. Since $u(0,0)=0$, we conclude that $c_s\too0$. Since $u$ is continuous near
$(0,0)$, we get $u(0,z_2)=\lim_{s\to+\infty}u(1/s,z_2)=\lim_{s\to+\infty}c_s=0$, $|z_2|\ll1$. Hence, since $\log u\in\PSH(G)$, we have
$u(0,z_2)=0$ for all $|z_2|<e^{\phi(0)}$. Consequently, $\bbs^{(p)}_G((0,0),(0,z_2))=0$, $|z_2|<e^{\phi(0)}$,
and $\bS^{(p)}_G((0,0);(0,X_2))=0$, $X_2\in\CC$.
\end{Example}

\section{Relations between $(\bm_G, \bbs_G, \bg_G)$ and $(\bgamma_G, \bS_G, \bA_G$)}

We will discuss the following two problems.
\emph{Find a pseudoconvex domain $G\subset\CC^n$, $a\in G$, and $z^0\in G$ (resp.~$X_0\in\CC^n$) such that}
\begin{gather*}
\bm_G(a,z^0)<\bbs_G(a,z^0)<\bg_G(a, z^0)\\
\text{(resp. } \bgamma_G(a;X_0)<\bS_G(a;X_0)<\bA_G(a;X_0)\text{)}.
\end{gather*}

\begin{Example}\label{Exam4}
If $\alpha_1,\dots,\alpha_n\in\ZZ$ are relatively prime, $\sigma(a)\geq2$, and $\mu(a)\geq2$, then the domain $G=\BD_\alpha$
(cf.~Example \ref{Exam2}) is an example of a pseudoconvex domain (unfortunately, unbounded)
such that $\bm_G(a,\cdot)<\bbs_G(a,\cdot)<\bg_G(a,\cdot)$ on $\BD_\alpha\setminus\BV_0$.
\emph{It is not know whether there exists a bounded pseudoconvex domain with this property.}
\end{Example}

\begin{Example}\label{Exam5}
Let $G\subset\CC^n$ be a balanced domain (i.e.~$\overline\DD\cdot G=G$) and let $\Fh_G(z)$
be the \emph{Minkowski functional of $G$}. It is known that $G=\{z\in\CC^n: \Fh_G(z)<1\}$.
Moreover, $\bg_G(0,\cdot)=\Fh_G$ in $G$ $\Longleftrightarrow$ $\bA_G(0;\cdot)\equiv\Fh_G$ $\Longleftrightarrow$
$G$ is pseudoconvex $\Longleftrightarrow$ $\log\Fh_G\in\PSH(\CC^n)$ (cf.~\cite{JarPfl2013}, Proposition 4.2.10(b)).

Let $\wdht G$ be the convex envelope of $G$. It is known that $\wdht G$ is also balanced and
$\Fh_{\wdht G}=\sup\{q: q:\CC^n\too[0,+\infty) \text{ is a $\CC$-seminorm with } q\leq\Fh_G\}$.
Moreover (cf.~\cite{JarPfl2013}, Proposition 2.3.1(d)), $\bgamma_G(0;\cdot)\equiv\Fh_{\wdht G}$.
Thus, if $G$ is pseudoconvex, then $\bgamma_G(0;\cdot)=\Fh_{\wdht G}\geq\bS_G(0;\cdot)$ and hence
$\bgamma_G(0;\cdot)\equiv\bS_G(0;\cdot)\equiv\Fh_{\wdht G}\leq\Fh_G\equiv\bA_G(0;\cdot)$. Consequently, we get the following result.

If $G$ is a balanced pseudoconvex non-convex domain, then
$$
\Fh_{\wdht G}\equiv\bgamma_G(0;\cdot)\equiv\bS_G(0;\cdot)\myfrac{\leq}{\not\equiv}\bA_G(0;\cdot)\equiv\Fh_G.
$$
\end{Example}
In particular, the result solves the problem formulated in Example 4.2.17 from \cite{JarPfl2013}.

\begin{Example}\label{Exam6}
Keep the notation from Example \ref{Exam1}. Then
$$
\bgamma_G(c_t;X_0)<\bS^{(2p)}_G(c_t;X_0)=\bA_G(c_t;X_0)=e^{\phi(0,t)},\quad p\in\NN,\; 0<t\ll1.
$$

Indeed, the function $G\ni(z_1,z_2,z_3)\overset{v}\tuu|z_1|e^{\phi (z_2,z_3)}$ is of the class $\CS^{(2p)}_G(c_t)$, which gives
$$
\bS^{(2p)}_G(c_t;X_0)\geq\limsup_{\lambda\to0}\frac{v(c_t+\lambda X_0)}{|\lambda|}=e^{\phi(0,t)}>0,\quad t>0.
$$
Observe that the mapping $e^{-\phi(0,t)}\DD\ni\lambda\overset{F}\tuu(\lambda,0,t)\in G$ is well-defined. Hence, using the
holomorphic contractibility, we get
$$
\bA_G(c_t;X_0)=\bA_G(F(0);F'(0)(X_0))\leq\bA_{e^{-\phi(0,t)}\DD}(0;1)=e^{\phi(0,t)}.
$$
Thus $\bS^{(2p)}_G(c_t;X^0)=\bA_G(c_t;X^0)=e^{\phi(0,t)}\geq e^{\phi(0,0)}>0$, $t>0$.

Now, to get the result it suffices to show that $\bgamma_G((0,0,0);X^0)=0$ and then use the continuity of $\bgamma_G(\cdot;X^0)$.
For, let $f\in\CO(G,\DD)$, $f(0,0,0)=0$. Since $\{0\}\times\CC^2\subset G$, the Liouville theorem implies that $f(0,\cdot,\cdot)=\const$.
Since $f(0,0,0)=0$, we get $f(0,\cdot,\cdot)\equiv0$. Since $\CC\times\{a_k\}\times\{0\}\subset G$, we get $f(\cdot,a_k,0)=\const(k)$.
Thus $f(\cdot,a_k,0)\equiv0$. Since the sequence $(a_k)_{k=1}^\infty$ is dense in $\DD$, we conclude that $f=0$ on
$(\CC\times\DD\times\{0\})\cap G$. Thus $f(z_1,0,0)=0$ provided that $|z_1|<e^{-\phi(0,0)}$. Hence $f'(0,0,0)(X^0)=0$ and so
$\bgamma_G((0,0,0);X^0)=0$.
\end{Example}

\section{Derivative}

Recall that for $\CF\in\{\CM, \CK\}$ we have $\delta^\CF_G(a;X)=\limsup_{\lambda\to0}\frac{d^\CF_G(a,a+\lambda X)}{|\lambda|},\quad a\in G,\;X\in\CC^n$ (cf.~\eqref{eq1} and \eqref{eq2}). It is an open problem whether
$$
\bS^{(2p)}_G(a;X)=\limsup_{\lambda\to0}\frac{\bbs^{(2p)}_G(a,a+\lambda X)}{|\lambda|},\quad a\in G,\;X\in\CC^n.
$$
Observe that
\begin{multline*}
\bS^{(2p)}_G(a;X)=\sup\Big\{\limsup_{\lambda\to0}\frac{v(a+\lambda X)}{|\lambda|}: v\in\CS^{(2p)}_G(a)\Big\}\\
\leq\limsup_{\lambda\to0}\frac{\bbs^{(2p)}_G(a,a+\lambda X)}{|\lambda|}\leq\limsup_{\lambda\to0}\frac{\bg_G(a,a+\lambda X)}{|\lambda|}=\bA_G(a;X),
\end{multline*}
so the problem is trivial if $\bS^{(2p)}_G(a;X)=\bA_G(a;X)$.

\bibliographystyle{amsplain}