# Deposition of quantum dots in a capillary tube

Yong Lin Kong  Department of Mechanical and Aerospace Engineering, Princeton
University, Princeton, NJ 08544  François Boulogne  Department of Mechanical
and Aerospace Engineering, Princeton University, Princeton, NJ 08544
Hyoungsoo Kim  Department of Mechanical and Aerospace Engineering, Princeton
University, Princeton, NJ 08544  Janine Nunes  Department of Mechanical and
Aerospace Engineering, Princeton University, Princeton, NJ 08544  Jie Feng
Department of Mechanical and Aerospace Engineering, Princeton University,
Princeton, NJ 08544  Howard A. Stone  Department of Mechanical and Aerospace
Engineering, Princeton University, Princeton, NJ 08544

###### Abstract

The ability to assemble nanomaterials, such as quantum dots, enables the
creation of functional devices that present unique optical and electronic
properties. For instance, light-emitting diodes with exceptional color purity
can be printed via the evaporative-driven assembly of quantum dots.
Nevertheless, current studies of the colloidal deposition of quantum dots have
been limited to the surfaces of a planar substrate. Here, we investigate the
evaporation-driven assembly of quantum dots inside a confined cylindrical
geometry. Specifically, we observe distinct deposition patterns, such as
banding structures along the length of a capillary tube. Such coating behavior
can be influenced by the evaporation speed as well as the concentration of
quantum dots. Understanding the factors governing the coating process can
provide a means to control the assembly of quantum dots inside a capillary
tube, ultimately enabling the creation of novel photonic devices.

##  1 Introduction

Nanomaterials enable the creation of functional devices that possess unique
properties. For instance, cadmium selenide-zinc sulfide quantum dots (QDs) [1]
have a narrow emission spectra that enable the creation of light-emitting
diodes with exceptional color purity [2]. Further, such nanomaterials can be
dispersed into solvent to form solution-processable inks, which can be
integrated with coating or printing processes to create photonic devices on a
two-dimensional substrate [3] or on a three-dimensional construct [4]. The
integration of electronics with a variety of constructs could enable new
applications, such as the realization of next-generation biomedical devices
[5] or on-body electronics [6]. However, most of the studies focus on the
coating of nanomaterials onto an open surface. Yet, in a three-dimensional
configuration, the total area of a confined surface can be significantly
larger than its exposed surface. More importantly, in many biologically
evolved or synthetically engineered systems, such as the transport of oxygen
and nutrient in a blood vessel or a DNA analysis in a microfluidic chip, many
of the critical processes occur inside a confined region rather than on its
open surface. Hence, the ability to coat the surface within a confined space
with nanomaterials could enable the integration of functional devices for new
applications.

Unfortunately, the confined nature prohibits the coating of such surfaces by
using conventional lithography [7] or many of the novel processing strategies
that were developed for open surfaces such as dip-pen nanolithography [8],
transfer [9], or 3D printing [4]. On the other hand, evaporative self-assembly
is a promising method in which the nanoparticles can be assembled into a
variety of patterns by manipulating parameters such as the evaporation rate
[10], Marangoni effects, [11] etc. For instance, on a planar surface, the self
assembly of micrometer size band structures have been observed in the
literature [12, 13, 14, 15, 16, 17]. Such patterns are believed to originate
from the pinning of the contact line of a droplet during the evaporation
process, which causes the accumulation of the particles near the contact line
in a so-called coffee ring effect [18, 19]. As the suspension evaporates, the
contact line depins before it is pinned again at the next position, and in
this sequential way, stripes or multiple irregular rings are formed [20].

In contrast, experiments that involve the active retraction of a substrate or
meniscus may be described as an active set-up [21], and have been shown to be
able to generate highly regular bands in a planar configuration. The
withdrawal speed of the substrate from a liquid bath of colloidal solution can
be controlled to create a highly regular array of band structures [14]. The
formation of bands are described as a consequence of the breakage of the
meniscus as it is being stretched during the evaporation [14].The substrate
withdrawal speed in such configurations plays a critical role in determining
the coating behavior, which can typically be classified into an evaporative
regime for small capillary number and Landau-Levich regime for larger
capillary number, where the capillary number represents the ratio of viscous
to interfacial stresses. [22, 23, 24, 25]. The readers may refer to a recent
review article on patterned deposition [21] for a comprehensive list of
references regarding the coating of open surfaces.

Nevertheless, significantly fewer studies have been conducted to investigate
the coating of a confined space. The confinement of a closed space imposes a
limitation on the evaporation rate. Previous studies have been performed in
capillary tubes with silica particles [26, 27, 28, 15], polystyrene beads [29,
28], gold nanoparticles [28], tobacco mosaic virus [30] and proteins [31]. Our
study differs in three major aspects. First, the majority of the previous
studies were done with structurally simple micron-size particles that are
spherical in shape. In contrast, typical functional nanoparticles such as QDs
have a more complex structure with a polyhedral shape. For instance, (CdSe)ZnS
QD is a core-shell particle that consists of a small inorganic semiconductor
core and a wider-bandgap inorganic semiconductor shell that is passivated with
organic ligands such as octadecylamine. The organic ligands confer solubility
for the QDs in a non-polar solvent, but they also play an important role that
affects the self-assembly during the deposition process [32]. Secondly, most
of the previous studies in a confined geometry have not examined the real-time
dynamics of the deposition. In our study, such observations are possible as
the intrinsic photoluminescence with narrow width of emission of QDs enables
the real-time visualization of the nano-particle coating process without the
need to tag with fluorescent polymers, which potentially affects the
deposition behavior. Further, QDs have been shown to exhibit fluorescence
lifetimes that are 10 - 100 times longer than organic fluorescent polymers
[33]. This enable reliable visualization and consistent measurement of the
coated QDs pattern based on the photoluminescence intensity. Third, the
ability of visualizing the particles enables the study of nanoparticle
coatings at significantly lower concentrations than previous studies. The
volume fractions that have been studied here are approximately 10 - 100 times
(0.02 - 0.4 vol %) lower than the range that was typically studied in other
reports in the literature. This limit is important for the purpose of
nanoelectronic printing where a thin layer deposition of nanoparticles is
required [2, 4].

##  2 Experimental Section

Figure 1: Experimental schematic and CdSe/ZnS quantum dots (QDs). (A) QDs
suspended in toluene are loaded to initial length $\ell$ via capillary forces
into a glass capillary tube of length $\ell_{g}$. The right end of the tube is
then sealed and the tube is placed on a horizontal surface with an open end
(left side). (B) Cross-section of the tube at the open end. Toluene evaporates
from the open end (diameter $d$) and the meniscus moves towards the closed end
with velocity $u_{m}$. The distance between the meniscus and the tube’s open
end is defined as $x_{m}$. (C) Transmission electron microscope image of the
QDs, which have polygonal shapes in cross section. (D) The QDs have an average
size distribution around 8 nm.

The fabrication of functional devices inside a confined geometry requires the
controlled deposition of the nanoparticles inside the surface. Our study aims
to understand these coating processes that can ultimately be used to create
such functional devices inside a closed construct such as a cylindrical
channel. In our experiments, a glass capillary tube (Wiretrol, Drummond
Scientific Company, Broomall, PA) with diameter $d$ and length, $\ell_{g}$, of
75 mm is first held vertically with a linear stage (PT1, Thorlabs, Newton, NJ)
and immersed into a suspension of QDs. The particles used were cadmium
selenide-zinc sulfide QDs (QSP-600, Ocean Nanotech, San Diego, CA). The QDs
are dispersed in toluene (Sigma-Aldrich, St. Louis, MO) at concentration
$c_{0}$ from 0.05 mg/ml to 1.00 mg/ml ($\approx$ 0.02 vol% to 0.40 vol%). The
suspension fills the tube to a length $\ell_{0}$. The other end is sealed with
hot melt adhesive before the capillary tube is carefully withdrawn from the
suspension.

The capillary tube is then placed on a horizontal surface, with a loaded
suspension of a length $\ell$, as shown in Figure 1A. Toluene evaporates from
the open end and the meniscus moves towards the closed end with velocity,
$u_{m}$, as sketched in Figure 1B. The distance between the meniscus and the
tube’s open end is defined as $x_{m}$. In contrast, any gas trapped at the
closed end is saturated with toluene. Hence, the meniscus at the closed end
will have a negligible movement in comparison with the meniscus at the open
end. This effect isolates the evaporation area to the meniscus at the open end
of the tube. In comparison, previous studies of the colloidal deposition
inside capillary tubes have been done by vertically [29, 34] immersing the
tube into suspensions and allowing both the bath and the liquid in the tube to
evaporate simultaneously, where the movement of the liquid is determined by
the rate of evaporation from both the open end of the capillary tube as well
as the bath reservoir.

Here, the samples are illuminated with ultra-violet light and the images are
acquired with a digital single-lens reflex camera (D5200, Nikon, Tokyo, Japan)
with a macro lens (85 mm f/3.5G Micro Nikkor Lens, Nikon, Tokyo, Japan) that
is attached to a custom-made stage. The experiments are conducted in a glove
box covered with dark cloth to minimize background noise for subsequent image
processing. The temperature on the stage of the samples is found to be stable
and consistent at 26$\pm$1∘C throughout the experiments, which is slightly
above the room temperature due to the heat from the light source.

The QDs are characterized with transmission electron microscopy (CM100, FEI,
Hillsboro, OR). As described above, CdSe/ZnS QDs are core-shell nanoparticles
with polygonal shape as shown in Figure 1C. By measuring the diameter of the
particle obtained by assuming circularity, the TEM image shows the mean size
distribution of 8 nm in Figure 1D. In addition to the real-time images
acquired with the camera, the deposited QD features are characterized with
confocal microscopy (TCS SP5, Leica Microsystems, Wetzlar, Germany). To
examine the inner surface, we first cleaved the tubes into smaller sections of
approximately 3 mm in length with a diamond scriber. The sectioned tube is
then attached to a double sided carbon tape before the capillary tube is
scribed along the length of the tube with a diamond scriber. The scribed
capillary tube is then cracked open with a metal rod, where the outer surface
of the capillary tube pieces would be attached on the double sided carbon
tape. A fluorescent image is taken to identify the region of interest prior to
imaging with Scanning Electron Microscopy (Quanta 200 FEG ESEM, FEI,
Hillsboro, OR), and Atomic Force Microscopy (Dimension NanoMan, Veeco,
Plainview, NY).

##  3 Results and Discussion

Figure 2: Evaporatively-driven motion of the meniscus inside the capillary
tubes. (A) Time series images show the column length $\ell$ as a function of
time for a 0.86 mm diameter capillary tube (interval = 160 minutes); the scale
bar is 1 mm. (B) The column length $\ell$ versus time for different tube
diameters, $d$ = 0.49, 0.60 and 0.86 mm, as toluene evaporates; $c_{0}$ = 0.50
mg/ml. The inset shows the movement of the meniscus follow $t^{1/2}$, where
$\ell_{0}$ is the initial column length. (C) The speed, $u_{m}$, of the
meniscus as it travels along the capillary (x); $c_{0}$ = 0.50 mg/ml. The
speed decreases with $x_{m}$. The evaporation is limited by the diffusion of
the toluene from the meniscus front to the opening of the tube, and can be
described with a theoretical model derived from Fick’s first law without the
need of fitting parameters (solid red line).

When QD suspensions are illuminated with ultra-violet light, the image
processing of the photoluminescence intensity enables the tracking of the
evaporatively-driven motion of the menisci inside the capillary tubes such as
the images shown in Figure 2A. Figure 2B shows measurements of the meniscus
positions versus time $t$ for three different tube diameters. The amount of
toluene lost ($x_{m}=\ell_{0}-$ $\ell$) scales linearly with $t^{1/2}$ as
plotted in the inset, which suggests that the evaporation is limited by the
diffusion of the toluene molecules to the opening of the tube. We hypothesize
that the drying dynamics is controlled by the vapor phase. To examine this
hypothesis, we adapted the derivation of Erbil et. al [35] to describe our
experimental results, where the system is approximated with Fick’s first law.
We assume that the Péclet number, which represents the balance between
convection and diffusion, is small, i.e. Pe = $\frac{x_{m}^{2}}{D_{0}t}\ll 1$.
Hence, for the vapor phase concentration, $c_{v}$, this limit corresponds to a
quasi-stationary state at the region away from the tube opening such that
$\frac{\partial c_{v}}{\partial t}=0$.

$u_{m}=\frac{{\rm d}x_{m}}{{\rm d}t}=-\frac{D_{0}}{\rho}\frac{\partial
c_{v}}{\partial x},$ (1)

where $c_{v}$ is the vapor concentration, $D_{0}$ is the vapor phase diffusion
coefficient ($8.11\times 10^{-6}$ m2/s) [35] and $\rho$ is the liquid phase
density of toluene ($8.61\times 10^{2}$ kg/m3). Further, we can also assume
that the toluene vapor concentration is zero at the outlet. We can then
approximate

$\frac{\partial c_{v}}{\partial
x}=\frac{c_{e}-0}{x_{m}-0}=\frac{c_{e}}{x_{m}}=\frac{P_{sat}M}{RTx_{m}},$ (2)

where we have used the ideal gas law with $P_{sat}$, the saturated vapor
pressure at the absolute temperature of 299 K ($4.00\times 10^{3}$ Pa), M the
molecular weight ($9.21\times 10^{-2}$ kg/mol), R the ideal gas constant (8.31
J/mol K), and T the absolute temperature (299 K). Since the vapor pressure of
the toluene is significant, we then incorporate Stefan’s general diffusion
theory [36] as a correction factor, where
$P=P_{T}\ln\left(\frac{P_{T}}{P_{T}-P_{sat}}\right)$ with $P_{T}$ is the total
pressure ($1.01\times 10^{5}$ Pa). Combining equations (1) and (2) yields

$u_{m}=\frac{K}{x_{m}},$ (3)

where $K=\frac{D_{0}MP_{T}}{\rho
RT}\ln\left(\frac{P_{T}}{P_{T}-P_{sat}}\right)$. The integration of equation
(3) yields $\ell(t)=\ell_{0}-\sqrt{2Kt}$ where $\ell(t)=\ell_{0}-x_{m}$. This
behavior, $u_{m}\propto x_{m}^{-1}$, describes the reduction of evaporation
speed with distance, $x_{m}$, as shown in Figure 2C. We note that the
prefactor is independent of the capillary tube diameter, as shown by the solid
red line in Figure 2C. The prediction (the solid red line in Figure 2C), which
includes the prefactor K, matches our experimental data in a region away from
the tube entrance where Pe = $\frac{x_{m}^{2}}{D_{0}t}\approx 10^{-4}\ll 1$.
This supports our hypothesis that the evaporation is limited by the diffusion
of the evaporated toluene molecules to the opening of the tube. This drying
dynamics is different from the confined unidirectional drying with a high
volume fraction of particles [37, 38, 39].

Near the opening of the tube, we observed that the QDs self-assembled into
micrometer wide bands and this feature is regular over a millimeter length
scale, similar to the other reports in the literature as discussed in the
introduction. This pattern is visible without a microscope, and is
particularly clear at higher concentrations of QDs (for instance, 0.50 mg/ml),
as shown in the example in Figure 3A. We have also verified that there is no
observable quenching effect throughout the duration of the experiments (see
Supporting Information). We then used a fluorescence confocal microscope to
further examine the region coated by particles (e.g. Figure 3B,C). The narrow
emission spectrum of QDs provides a clear visualization of the coated region;
we note that the color tunability of QDs based on the change of core-size
would enable the study to be extended to include a mixture QDs of different
sizes.

Figure 3: QDs bands inside a capillary tube and the concentration dependence
of the band width. (A) Image of bands formed from the QDs as toluene
evaporates from the open end on the left. $c_{0}$ = 0.50 mg/ml. (B,C) Confocal
microscopy image showing QD band arrays along the tubes; with initial
concentration, $c_{0}$, of 0.50 mg/ml (0.20 vol%) and 1.00 mg/ml (0.40 vol%).
We defined $w$ as the width of the band. (D) The graph shows the average
width, measured from five different samples, of a band located at $x_{m}$ = 3
mm. Figure 4: QD bands as a function of distance from the tube’s opening. (A)
The width of the band, w, and (B) the distance between two bands measured
along the capillary tube for three different initial concentrations, $c_{0}$.
(C) Confocal image shows the transition from bands to a fully coated region
($c_{0}$ = 0.50 mg/ml, $d$ = 0.90 mm). On the fully coated region, the
patterns perpendicular to the band direction are formed from the cracking of
the deposited QDs.  Figure 5: Growth of a single band. (A) Time series images
show the growth of a single band ($d=0.86\text{ mm},c_{0}=0.50\text{ mg/ml}$).
(B) The graph shows the relationships between the width, $w_{m}$, and
wavelength, $\lambda$ (see inset for the definition).

As a means of controlling the width of the band, we examined the dependence of
the band width with the initial concentration. To isolate the effect of
concentration, we analyzed the width of the bands formed 3 mm away from the
tube opening, where the speed of the meniscus movement is approximately 0.4
$\mathrm{\SIUnitSymbolMicro m}$$\cdot$$\text{s}^{-1}$. We found that the band
width increases approximately linearly with an increase of initial
concentration, as shown in Figure 3D. This observation is consistent with the
increase of the particle flux to the pinned contact line. Larger stripes can
be formed at a specific location before the meniscus retracted as more
particles are being deposited in the same position in the same amount of time.
This result is also consistent with the observations by Watanabe et al. [14]
and Kaplan et al. [40] in a planar geometry, where a dependence of the band
width with the concentration of the particles was reported.

Figure 6: QDs bands inside a capillary tube. (A) Confocal microscopy image
showing QD band arrays along the tubes; concentration $c_{0}$ = 0.05 mg/ml
(0.02 vol%). (B) Secondary electron microscopy image (orange box in A) shows
the morphology and regularity of the band array, with a band width of 5
$\mathrm{\SIUnitSymbolMicro m}$. (C) Atomic force microscopy image of one of
the bands; color represents the thickness of the QD film. Inset shows the
profilometry measurement extracted along the path ($x_{r}$) of the red arrow
on the AFM image. The band thickness increases linearly as the meniscus
travels with $x_{r}$ (4 $\mathrm{\SIUnitSymbolMicro m}$ to 6
$\mathrm{\SIUnitSymbolMicro m}$), before a sharp decrease as the meniscus
slipped to form the next band.

These particulate bands have also been observed in other studies with a
confined cylindrical configuration with polystyrene beads [29], tobacco mosaic
virus [30], proteins [31] and polymers [41]. However, it is important to note
that the origin of such formation also differs based on the experimental
configuration. For instance, in Abkarian et al.’s experimental configuration
[29], the bands are formed when the contact lines are depinned due to the
relative increase of column height as liquid evaporates faster from the bath.
In our experimental configuration, the band formation of the meniscus is only
affected by the evaporation rate at the tube opening and the local
concentration of the QDs as the meniscus is controlled by the evaporation of
the toluene, and only one end of the capillary tube is open. This
configuration enables us to understand the key factors that determines the
features of the band formation.

Next, we examined the dependence of the band structures on the distance from
the tube opening. We found that in general, the width of the band increases as
the distance from the opening increases, as shown in Figure 4A. This
observation could be affected by two factors: first, the contact line speed
decreases as the distance from the tube open end ($x_{m}$) increases, as
described in Figure 2. Secondly, the concentration of the suspension increases
as toluene evaporates. The quantitative relationships between these factors
and the deposited patterns will be described in a future study. We have also
examined the distance between two bands and its relationship with the location
of the tube, $x$. The distance is measured between the two center points of
the bands. It is found that with the exception of the 0.05 mg/ml
concentration, most of the spacing between the two bands remain approximately
constant. As the width grows to be larger than the spacing between the bands,
the bands overlap and form a fully coated film, as shown in Figure 4C.

A closer inspection of the confocal images of the bands shows a subtle
undulation at the right side of the band (for instance, at the bands of Figure
3B). Similar observations can also be seen in the bands reported in the
literature [42]. However, the growth of the fingers of the band are not
described in detail, partly due to the challenge of visualizing the deposition
process in real-time without fluorescence microscopy. In our experiment, the
high photoluminescence efficiency of quantum dots enables the direct
observation of the band growth with a fluorescence confocal microscope (see
Supporting Information, Movie 1). A time series of images shows the growth of
a single band in Figure 5A. Fingers can be seen during the growth and the
pattern moves laterally along the band as it grows. In most cases, the fingers
are smoothed before the meniscus front has slipped. In general, the width
grows approximately linearly with time while the wavelength increases in the
beginning of the growth but saturates at the end, as shown in the graph in
Figure 5B. The plotted values represent the average of the minimum width per
period, $w_{m}$, and average wavelength, $\lambda$, (see inset at Figure 5B
for the definition) that are observed in the time series images during the
growth of a single band as shown in Figure 5A.

To examine the morphology of the QD bands, the capillary tube is cleaved and
examined with a scanning electron microscope (SEM). An array of regular bands
with 5 $\mathrm{\SIUnitSymbolMicro m}$ width is observed (Figure 6B), as well
as some dewetted tails in the gap between the bands. SEM images also reveal
cracks on the bands. The origin of cracks may originate from the contraction
of the QDs through capillary pressure, which is being constrained by the
adhesion of the QD film with the glass substrate [43, 44]. Further, SEM images
also show a smaller crack wavelength at the beginning of the band than the end
of the band, which suggests that the deposit thickness increases in the
direction of the drying [37].

Indeed, examining one of the QD bands with atomic force microscopy verified
the above observation. Profilometry data shows that the deposit thickness
increases linearly as the band grows to a relative thickness of 90 nm before a
sharp decrease of thickness (Figure 6C). This thickness profile suggests that
the band grows gradually when particles accumulate at the pinned contact line,
before the meniscus slips to the next position as the toluene evaporates. This
observation is also consistent with reports in the literature, which used
spherical polystyrene particles [29, 15] and a similar thickness profile of
the bands was measured. On the other hand, the thickness profile with Tobacco
Mosaic Virus rods [30] (300 nm $\times$ 18 nm) is observed to generate a band
structure with monolayer thickness. The differences could be due to the aspect
ratio of the particles, which affects the packing structure when the particles
accumulate at the pinned contact line.

##  4 Conclusion

We presented a method to study the deposition of functional nanoparticles in a
confined geometry. The intrinsic photoluminescence property of QDs enables
clear observations of the coating processes without the interference from
fluorescent polymers. The experimental platform also enables the real-time
measurement of the evaporation speed as well as the visualization of the
growth of bands, even at a low volume fraction of QDs. Decoupling the
evaporation rate and the concentration in the liquid phase requires further
investigations. For instance, the direct measurement of concentration is
possible by calibrating the volume fraction of QDs with its photoluminescence
intensity. These understandings could provide key insights on the coating of
nanoparticles in a confined geometry and ultimately enable the creation of
functional devices inside such constructs.

### acknowledgement

The authors thank Dr. Sangwoo Shin, Jesse Ault, Dr. Jason Wexler and Prof.
Michael C. McAlpine for helpful discussions throughout the project. Y.L.K.
thanks the Department of Mechanical and Aerospace Engineering, Princeton
University for partial support. F.B. acknowledges that the research leading to
these results received funding from the People Programme (Marie Curie Actions)
of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA
grant agreement 623541.

## References

  * [1] A P Alivisatos.  Semiconductor Clusters, Nanocrystals, and Quantum Dots.  Science., 271(5251):933–937, 1996.
  * [2] Y Shirasaki, G J Supran, M G Bawendi, and V Bulovic.  Emergence of Colloidal Quantum-dot Light-emitting Technologies.  Nat. Photonics, 7(1):13–23, 2013.
  * [3] Vanessa Wood, Matthew J. Panzer, Jianglong Chen, Michael S. Bradley, Jonathan E. Halpert, Moungi G. Bawendi, and Vladimir Bulovic.  Inkjet-Printed Quantum Dot-Polymer Composites for Full-Color AC-Driven Displays.  Adv. Mater., 21(21):2151–2155, June 2009.
  * [4] Yong Lin Kong, Ian A Tamargo, Hyoungsoo Kim, Blake Johnson, Maneesh K Gupta, Tae-Wook Koh, Huai-An Chin, Daniel A Steingart, Barry P Rand, and Michael C McAlpine.  3D Printed Quantum Dot Light-Emitting Diodes.  Nano Lett., 14(12):7017–7023, October 2014.
  * [5] R H Kim, D H Kim, J L Xiao, B H Kim, S I Park, B Panilaitis, R Ghaffari, J M Yao, M Li, Z J Liu, V Malyarchuk, D G Kim, A P Le, R G Nuzzo, D L Kaplan, F G Omenetto, Y G Huang, Z Kang, and J A Rogers.  Waterproof AlInGaP Optoelectronics on Stretchable Substrates with Applications in Biomedicine and Robotics.  Nat. Mater., 9(11):929–937, 2010.
  * [6] D H Kim, N Lu, R Ma, Y S Kim, R H Kim, S Wang, J Wu, S M Won, H Tao, A Islam, K J Yu, T I Kim, R Chowdhury, M Ying, L Xu, M Li, H J Chung, H Keum, M McCormick, P Liu, Y W Zhang, F G Omenetto, Y Huang, T Coleman, and J A Rogers.  Epidermal Electronics.  Science., 333(6044):838–843, 2011.
  * [7] S Coe, W K Woo, M Bawendi, and V Bulovic.  Electroluminescence from Single Monolayers of Nanocrystals in Molecular Organic Devices.  Nature, 420(6917):800–803, 2002.
  * [8] Richard D Piner, Jin Zhu, Feng Xu, Seunghun Hong, and Chad A Mirkin.  ”Dip-Pen” Nanolithography.  Science., 283(5402):661–663, 1999.
  * [9] J H Ahn, H S Kim, K J Lee, S Jeon, S J Kang, Y G Sun, R G Nuzzo, and J A Rogers.  Heterogeneous Three-dimensional Electronics by Use of Printed Semiconductor Nanomaterials.  Science., 314(5806):1754–1757, 2006.
  * [10] Terry P Bigioni, Xiao-Min Lin, Toan T Nguyen, Eric I Corwin, Thomas a Witten, and Heinrich M Jaeger.  Kinetically Driven Self Assembly of Highly Ordered Nanoparticle Monolayers.  Nat. Mater., 5(4):265–270, April 2006.
  * [11] Hua Hu and Ronald G Larson.  Marangoni Effect Reverses Coffee-ring Depositions.  J. Phys. Chem. B, 110(14):7090–4, April 2006.
  * [12] Eiki Adachi, Antony S. Dimitrov, and Kuniaki Nagayama.  Stripe Patterns Formed on a Glass Surface during Droplet Evaporation.  Langmuir, 11(4):1057–1060, 1995.
  * [13] Jun Xu, Jianfeng Xia, and Zhiqun Lin.  Evaporation-Induced Self-Assembly of Nanoparticles from a Sphere-on-Flat Geometry.  Angew. Chemie Int. Ed., 46(11):1860–1863, 2007.
  * [14] Satoshi Watanabe, Koji Inukai, Shunsuke Mizuta, and Minoru T Miyahara.  Mechanism for Stripe Pattern Formation on Hydrophilic Surfaces by Using Convective Self-assembly.  Langmuir, 25(13):7287–95, July 2009.
  * [15] Hugues Bodiguel, Frédéric Doumenc, and Béatrice Guerrier.  Stick-slip Patterning at Low Capillary Numbers for An Evaporating Colloidal Suspension.  Langmuir, 26(13):10758–10763, 2010.
  * [16] Hyun Suk Kim, Cheol Hee Lee, P. K. Sudeep, Todd Emrick, and Alfred J. Crosby.  Nanoparticle Stripes, Grids, and Ribbons Produced by Flow Coating.  Adv. Mater., 22(41):4600–4604, 2010.
  * [17] Diego Noguera-Marin, Carmen L Moraila-Martinez, Miguel A Cabrerizo-Vilchez, and Miguel A Rodriguez-Valverde.  Transition from stripe-like patterns to a particulate film using driven evaporating menisci.  Langmuir, 30(25):7609–7614, 2014.
  * [18] Robert D. Deegan, Olgica Bakajin, Todd F. Dupont, Greb Huber, Sidney R. Nagel, and Thomas A. Witten.  Capillary Flow as The Cause of Ring Stains from Dried Liquid Drops.  Nature, 389(6653):827–829, October 1997.
  * [19] Robert Deegan, Olgica Bakajin, Todd Dupont, Greg Huber, Sidney Nagel, and Thomas Witten.  Contact Line Deposits in an Evaporating Drop.  Phys. Rev. E, 62(1):756–765, July 2000.
  * [20] Leonid Shmuylovich, Amy Q. Shen, and Howard A. Stone.  Surface Morphology of Drying Latex Films: Multiple Ring Formation.  Langmuir, 18:3441–3445, 2002.
  * [21] Uwe Thiele.  Patterned Deposition at Moving Contact Lines.  Adv. Colloid Interface Sci., 206:399–413, April 2014.
  * [22] Maël Le Berre, Yong Chen, and Damien Baigl.  From Convective Assembly to Landau-Levich Deposition of Multilayered Phospholipid Films of Controlled Thickness.  Langmuir, 25(5):2554–2557, 2009.
  * [23] F Doumenc and B Guerrier.  Drying of a Solution in a Meniscus: A Model Coupling the Liquid and the Gas Phases.  Langmuir, 26(17):13959–13967, 2010.
  * [24] Marco Faustini, Benjamin Louis, Pierre A Albouy, Monika Kuemmel, and David Grosso.  Preparation of Sol-Gel Films by Dip-Coating in Extreme Conditions.  J. Phys. Chem. C, 114(17):7637–7645, 2010.
  * [25] Guillaume Berteloot, Adrian Daerr, François Lequeux, and Laurent Limat.  Dip coating with colloids and evaporation.  Chem. Eng. Process., 68:69 – 73, 2013.
  * [26] H. Wang, X. Li, H. Nakamura, M. Miyazaki, and H. Maeda.  Continuous Particle Self-Arrangement in a Long Microcapillary.  Adv. Mater., 14(22):1662–1666, November 2002.
  * [27] Hiroyuki Nakamura, Xianying Li, Hongzhi Wang, Masato Uehara, Masaya Miyazaki, Hazime Shimizu, and Hideaki Maeda.  A simple method of self assembled nano-particles deposition on the micro-capillary inner walls and the reactor application for photo-catalytic and enzyme reactions.  Chem. Eng. J., 101(1):261–268, 2004.
  * [28] Chih-Jung Kuo, King Yiu Lam, Srikar Vengallatore, and Yoav Peles.  Colloidal self-assembly on internal surfaces of partially sealed microchannels.  J. Micromech. Microeng., 16(8):1667, 2006.
  * [29] Manouk Abkarian, Janine Nunes, and Howard A Stone.  Colloidal Crystallization and Banding in a Cylindrical Geometry.  J. Am. Chem. Soc., 126(19):5978–9, May 2004.
  * [30] Yuan Lin, Elizabeth Balizan, L Andrew Lee, Zhongwei Niu, and Qian Wang.  Self-assembly of Rodlike Bio-nanoparticles in Capillary Tubes.  Angew. Chem. Int. Ed. Engl., 49(5):868–72, January 2010.
  * [31] Yuan Lin, Zhaohui Su, Elizabeth Balizan, Zhongwei Niu, and Qian Wang.  Controlled Assembly of Protein in Glass Capillary.  Langmuir, 26(15):12803–9, August 2010.
  * [32] Henk N.W. Lekkerkerker and Remco Tuinier.  Colloids and the Depletion Interaction.  Springer Netherlands, 1 edition, 2011.
  * [33] Ute Resch-Genger, Markus Grabolle, Sara Cavaliere-Jaricot, Roland Nitschke, and Thomas Nann.  Quantum dots versus organic dyes as fluorescent labels.  Nat. Methods, 5(9):763–775, 2008.
  * [34] Wenhua Guo, Ming Wang, Wei Xia, and Lihua Dai.  Evaporation-induced Self-assembly of Capillary Cylindrical Colloidal Crystal in a Face-centered Cubic Structure with Controllable Thickness.  J. Mater. Res., 27(13):1663–1671, May 2012.
  * [35] H. Yildirim Erbil and Yonca Avci.  Simultaneous Determination of Toluene Diffusion Coefficient in Air from Thin Tube Evaporation and Sessile Drop Evaporation on a Solid Surface.  Langmuir, 18(13):5113–5119, June 2002.
  * [36] T. R. Marrero and E. A. Mason.  Gaseous Diffusion Coefficients.  J. Phys. Chem. Ref. Data, 1(1):3, 1972.
  * [37] C. Allain and L. Limat.  Regular Patterns of Cracks Formed by Directional Drying of a Collodial Suspension.  Phys. Rev. Lett., 74(15):2981–2984, 1995.
  * [38] KJ Wallenstein and WB Russel.  The theory of delamination during drying of confined colloidal suspensions.  J. Phys. Condens. Matter, 23(19):194104, 2011.
  * [39] Pierre Lidon and Jean-Baptiste Salmon.  Dynamics of unidirectional drying of colloidal dispersions.  Soft Matter, 10(23):4151–4161, 2014.
  * [40] C. Nadir Kaplan, Ning Wu, Shreyas Mandre, Joanna Aizenberg, and L. Mahadevan.  Dynamics of evaporative colloidal patterning.  Physics of Fluids, 27(9), 2015.
  * [41] Yingjuan Sun, Guihua Xiao, Yuan Lin, Zhaohui Su, and Qian Wang.  Self-assembly of Large-scale P3HT Patterns by Confined Evaporation in the Capillary Tube.  RSC Adv., 5(26):20491–20497, 2015.
  * [42] Zhiqun Lin.  Controlled Evaporative Assembly of Polymers from Confined Solutions.  J. Polym. Sci. Part B Polym. Phys., 48(24):2552–2557, December 2010\.
  * [43] A Groisman and E Kaplan.  An experimental study of cracking induced by desiccation.  Europhys. Lett., 25(6):415, 1994.
  * [44] Raymond C. Chiu and Michael J. Cima.  Drying of granular ceramic films: Ii, drying stress and saturation uniformity.  J. Am. Chem. Soc., 76(11):2769–2777, 1993.


# Monotonic Gaussian process for physics-constrained machine learning with
materials science applications

Anh Tran Corresponding author: Anh Tran<EMAIL_ADDRESS> Scientific
Machine Learning
Sandia National Laboratories
Albuquerque, NM 87185
Email<EMAIL_ADDRESS>
Kathryn Maupin  Optimization and Uncertainty Quantification
Sandia National Laboratories
Albuquerque, NM 87185
Email<EMAIL_ADDRESS>
Theron Rodgers
Computational Materials & Data Science
Sandia National Laboratories
Albuquerque  NM 87185
Email<EMAIL_ADDRESS>

###### Abstract

Physics-constrained machine learning is emerging as an important topic in the
field of machine learning for physics. One of the most significant advantages
of incorporating physics constraints into machine learning methods is that the
resulting model requires significantly less data to train. By incorporating
physical rules into the machine learning formulation itself, the predictions
are expected to be physically plausible. Gaussian process (GP) is perhaps one
of the most common methods in machine learning for small datasets. In this
paper, we investigate the possibility of constraining a GP formulation with
monotonicity on three different material datasets, where one experimental and
two computational datasets are used. The monotonic GP is compared against the
regular GP, where significant reduction in the posterior variance is observed.
The monotonic GP is strictly monotonic in the interpolation regime, but in the
extrapolation regime, the monotonic effect starts fading away as one goes
beyond the training dataset. Imposing monotonicity on the GP comes at a small
accuracy cost, compared to the regular GP. The monotonic GP is perhaps most
useful in applications where data is scarce and noisy, and monotonicity is
supported by strong physical evidence.

##  1 Introduction

Physics-constrained machine learning is an important toolbox in the machine
learning era, with many possible applications, including those in biology,
materials science, astrophysics, and finance. Breakthroughs in deep learning
have been recognized as some of the most revolutionary contributions to
science, but its Achilles heel is the size of the dataset required to train
huge deep learning architectures.

However, this requirement cannot always be satisfied, especially in fields
involving experimentation, where data can be resource-intensive, in terms of
both time and labor. In this context, physics-constrained machine learning
arises as an alternative solution, with equivalent or better accuracy, while
requiring significantly less data to train. By incorporating some fundamental
physics rules into the machine learning framework, the behavior of the machine
learning predictions could be more physically plausible. One common constraint
is the monotonicity, where the predictions are expected to behave
monotonically with respect to one or several input variables.

Materials science is a field where theory has overwhelmed in the last several
decades and centuries. The process-structure-property relationship is perhaps
one of the most well-studied topics, in terms of theoretical analysis.
Existing well before the computer era, material scientists tended to agree
that mathematics and physics are the right tools to bridge the gap between
inputs and quantities of interest. With the emergence of computers,
supercomputers, and subsequently deep learning, the main tool has changed
significantly, starting with the Materials Genome Initiative [1]. However,
brute-force applications of deep learning has suffered from the lack of data
in materials science, where physics-constrained machine learning seems to be a
more promising candidate. Furthermore, materials science theory suggests that
there are many simple physical rules that one can take advantage of when
constructing a machine learning framework, for example, the Hall-Petch
relationship [2]. It should be noted that many of these rules or simple
equations are involve monotonicity. If these rules are successfully exploited,
then a physics-constrained machine learning could be achieved with much less
data, while retaining the accuracy performance.

The Gaussian process (GP) stands out among other machine learning approaches
and has been one of the most popular methods for small datasets across
multiple disciplines. In particular, it is well known that GPs are often used
as underlying surrogate models for Bayesian optimization, which is an
efficient active machine learning method and very popular among materials
scientists. For example, Tallman et al. [3, 4] used a GP as a surrogate model
to bridge microstructure crystallography and homogenized materials properties,
where crystal plasticity finite element models are used to construct the
database, and Yabansu et al. [5] applied GPs to capture the evolution of
microstructure statistics. Tran and Wildey [6, 7] also used a GP as a
surrogate model for structure-property relationship and solved a stochastic
inverse problem using the acceptance-rejection sampling method. More examples
can be found in [8, 9].

One of the main advantages of using a GP is its flexibility and, therefore,
its ability to adopt extensions. Fernández-Godino et al. [10] proposed a novel
approach to constrain symmetry on GPs. Linear operator constraints can be
enforced through a novel covariance function, as demonstrated by Jidling et
al. [11]. Agrell [12] presented an approach to constrain GPs such that a set
of linear transformations of the process are bounded. Furthermore, Lange-
Hegermann [13] constrained multi-output GP priors to the solution set of a
system of linear differential equations subjected to boundary conditions. A
comprehensive survey study of constrained GPs, including, but not limited to,
bound constraints, non-negativity, monotonicity, linear partial differential
operator constraints, and boundary conditions of a partial differential
equation can be found in Swiler, et al. [14].

The present work focuses on bound-type constraints. Riihimäki and Vehtari [15]
proposed enforcing monotonicity to constrain GPs using binary classification,
where the derivative is classified as either monotonically increasing or
decreasing with respect to a set of specific input variables. Golchi et al.
[16] extended this approach in a deterministic computer experiments setting
and sampled from the exact joint posterior distribution rather than an
approximation. More recently, Ustyuzhaninov et al. [17] presented a
monotonicity constraint in a Bayesian non-parametric setting based on
numerical solutions of stochastic differential equations, and Pensoneault et
al. [18] presented a novel approach to construct a non-negative GP. Tan [19]
proposed a B-spline model that linearly constrains the model coefficients to
achieve monotonicity, which converges to the projection of the GP in
$\displaystyle L_{2}$. Chen et al. [20] proposed an optimal recovery method
that is equivalent to maximum a posterior estimation for a GP constrained by a
partial differential equation.

In this paper, we adopt the monotonic GP formulation from Riihimäki and
Vehtari [15] and apply it to three different datasets, both in interpolation
and extrapolation regimes, to survey the advantages and disadvantages between
the regular GP and the monotonic GP. We focus on experimental and
computational materials science applications where monotonicity is physically
expected. Materials science, specifically experimental materials science, is
well-known for data scarcity due to its resource-intensive experiments. An
overview of Gaussian Processes and the incorporation of monotonic constraints
therein are discussed in Section 2. Three examples of how classical and
monotonic GPs compare in performance are given in Sections 4, 5, and 6. An
overall discussion and conclusion are given in Section 7.

##  2 Overview of Gaussian Processes

###  2.1 Nomenclature

The symbols used in this paper are as follows:

  * •

$\displaystyle\mathbf{X}=\\{\mathbf{x}^{(i)}\\}_{i=1}^{N}$: training dataset
of size $\displaystyle N$, $\displaystyle\mathbf{X}\in\mathbb{R}^{N\times D}$

  * •

$\displaystyle\mathbf{x}=[\mathbf{x}_{1},\dots,\mathbf{x}_{D}]$: training
input of $\displaystyle D$ dimensionality

  * •

$\displaystyle\mathbf{X}_{m}$: derivative inducing points,
$\displaystyle\mathbf{X}_{m}\in\mathbb{R}^{m\times D}$

  * •

$\displaystyle\mathbf{y}$: noisy observation
$\displaystyle\mathbf{y}(\mathbf{x})=\mathbf{f}(\mathbf{x})+\varepsilon$,
$\displaystyle\varepsilon\sim\mathcal{N}(0,\sigma^{2})$

  * •

$\displaystyle\mathbf{f}$: noiseless output
$\displaystyle\mathbf{f}(\mathbf{x})$

  * •

$\displaystyle\mathbf{f}^{\prime}$: first derivative of
$\displaystyle\mathbf{f}$

  * •

$\displaystyle x^{*}$: testing input

  * •

$\displaystyle f^{*}$: testing output

  * •

$\displaystyle N$: number of data points

  * •

$\displaystyle M$: number of inducing data points for derivatives

  * •

$\displaystyle D$: input dimensionality

  * •

$\displaystyle i$, $\displaystyle j$: dummy data point index

  * •

$\displaystyle d$, $\displaystyle g$, $\displaystyle h$: dummy dimensionality
index

###  2.2 Classical Gaussian processes

A thorough mathematical explanation and derivation of Gaussian processes can
be found in the canonical text by Rasmussen and Williams [21], but we provide
a brief overview here.

Let us assume that we can model the true process, $\displaystyle\mathbf{y}$,
with a zero-mean GP,

$\mathbf{f}|\mathbf{\mathbf{X}}\sim\mathcal{N}(\mathbf{0},\mathbf{K}_{\mathbf{f},\mathbf{f}}),$
(1)

where the entries in the covariance matrix
$\displaystyle\mathbf{K}_{\mathbf{f},\mathbf{f}}$ are defined by a covariance
function. Although there are many covariance functions to choose from, in this
paper, we focus on the squared exponential covariance function

${\text{Cov}\left[f^{(i)},f^{(j)}\right]}=\mathbf{K}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})=\eta^{2}\exp\left[-\frac{1}{2}\sum_{d=1}^{D}\rho_{d}^{-2}(x^{(i)}_{d}-x^{(j)}_{d})^{2}\right],$
(2)

where $\displaystyle\eta$ and
$\displaystyle\mathbf{\rho}=\\{\rho_{1},\dots,\rho_{D}\\}$ are hyper-
parameters, representing the signal variance and length-scale (also known as
correlation length) parameters, respectively. The observations
$\displaystyle\mathbf{y}$ are then given by

$\mathbf{y}|\mathbf{f}\sim\mathcal{N}(\mathbf{f},\sigma^{2}\mathbf{I}),$ (3)

where $\displaystyle\sigma^{2}$ is the intrinsic noise variance.

At a desired input prediction location, $\displaystyle x^{*}$, the posterior
prediction distribution, also called the testing distribution, is a Gaussian
distribution, where the posterior mean and posterior variance are,
respectively,

${\mathbb{E}\left[f^{*}|x^{*},\mathbf{y},\mathbf{X},\mathbf{\theta}\right]}=\mathbf{K}_{\mathbf{*},\mathbf{f}}[\mathbf{K}_{\mathbf{f},\mathbf{f}}+\sigma^{2}\mathbf{I}]^{-1}\mathbf{y}$
(4)

and

${\mathbb{V}\left[f^{*}|x^{*},\mathbf{y},\mathbf{X},\mathbf{\theta}\right]}=\mathbf{K}_{\mathbf{*},\mathbf{*}}-\mathbf{K}_{*,\mathbf{f}}[\mathbf{K}_{\mathbf{f},\mathbf{f}}+\sigma^{2}\mathbf{I}]^{-1}\mathbf{K}_{\mathbf{f},\mathbf{*}}.$
(5)

Here, the hyper-parameters are collected into
$\displaystyle\mathbf{\theta}=\\{\eta,\mathbf{\rho},\sigma\\}$.

The GP is trained by solving for the hyper-parameters
$\displaystyle\mathbf{\theta}$ that maximize the log-marginal likelihood

$\begin{array}[]{lll}\log
p(\mathbf{y}|\mathbf{X},\mathbf{\theta})&=&-\frac{1}{2}\mathbf{y}^{\top}[\mathbf{K}_{\mathbf{f},\mathbf{f}}+\sigma^{2}\mathbf{I}]^{-1}\mathbf{y}\\\
&&-\frac{1}{2}\log{|\mathbf{K}_{\mathbf{f},\mathbf{f}}+\sigma^{2}\mathbf{I}|}-\frac{N}{2}\log{(2\pi)}\end{array}$
(6)

In Equation 6, the first term, so-called the “data fit” term, quantifies how
well the model fits the data described in the Mahalanobis distance. The second
term, called the “complexity” term, quantifies the model complexity where
smoother covariance matrix with smaller determinant is encouraged [21]. The
last term indicates that the likelihood of data tends to decrease with larger
data sets [22].

It is noteworthy that the derivative operator is a linear operator; therefore,

${\mathbb{E}\left[\frac{\partial f^{(i)}}{\partial
x_{d}^{(i)}}\right]}=\frac{\partial{\mathbb{E}\left[f^{(i)}\right]}}{\partial
x_{d}^{(i)}},$ (7) ${\text{Cov}\left[\frac{\partial f^{(i)}}{\partial
x_{d}^{(i)}},f^{(j)}\right]}=\frac{\partial}{\partial
x_{d}^{(i)}}{\text{Cov}\left[f^{(i)},f^{(j)}\right]},$ (8)

and

${\text{Cov}\left[\frac{\partial f^{(i)}}{\partial x_{d}^{(i)}},\frac{\partial
f^{(j)}}{\partial x_{g}^{(j)}}\right]}=\frac{\partial^{2}}{\partial
x_{d}^{(i)}\partial x_{g}^{(j)}}{\text{Cov}\left[f^{(i)},f^{(j)}\right]}.$ (9)

For the squared exponential kernel described in Equation 2, Equations 8 and 9
become

$\begin{array}[]{lll}{\text{Cov}\left[\frac{\partial f^{(i)}}{\partial
x_{g}^{(i)}},f^{(j)}\right]}&=&-\eta^{2}\exp\left(-\frac{1}{2}\sum_{d=1}^{D}\rho_{d}^{-2}(x_{d}^{(i)}-x_{d}^{(j)})^{2}\right)\\\
&&\rho_{g}^{-2}\left((x_{g}^{(i)}-x_{g}^{(j)})\right),\end{array}$ (10)

and

$\begin{array}[]{lll}{\text{Cov}\left[\frac{\partial f^{(i)}}{\partial
x_{d}^{(i)}},\frac{\partial f^{(j)}}{\partial
x_{h}^{(j)}}\right]}&=&\eta^{2}\exp\left(-\frac{1}{2}\sum_{d=1}^{D}\rho_{d}^{-2}(x_{d}^{(i)}-x_{d}^{(j)})^{2}\right)\\\
&&\rho_{g}^{-2}\left(\delta_{gh}-\rho_{h}^{-2}(x_{g}^{(i)}-x_{g}^{(j)})(x_{h}^{(i)}-x_{h}^{(j)})\right),\end{array}$
(11)

respectively, where $\displaystyle\delta_{gh}=1$ if $\displaystyle g=h$ and 0
otherwise.

For an arbitrary testing point $\displaystyle x^{*}$, the derivatives with
respect to the dimension of the posterior mean and posterior variance are,
respectively,

${\mathbb{E}\left[\frac{\partial f^{*}}{\partial
x^{*}_{d}}\right]}=\frac{\partial\mathbf{K}_{*,\mathbf{f}}}{\partial
x^{*}_{d}}[\mathbf{K}_{\mathbf{f},\mathbf{f}}+\sigma^{2}\mathbf{I}]^{-1}\mathbf{y},$
(12)

and

$\begin{array}[]{lll}{\mathbb{V}\left[\frac{\partial f^{*}}{\partial
x^{*}_{d}}\right]}=\frac{\partial^{2}\mathbf{K}_{*,*}}{\partial
x^{*}_{d}\partial x^{*}_{d}}-\frac{\partial\mathbf{K}_{*,\mathbf{f}}}{\partial
x^{*}_{d}}[\mathbf{K}_{\mathbf{f},\mathbf{f}}+\sigma^{2}\mathbf{I}]^{-1}\frac{\partial\mathbf{K}_{\mathbf{f},*}}{\partial
x^{*}_{d}}.\end{array}$ (13)

###  2.3 Monotonic Gaussian processes

In this work, we adopt the monotonic GP formulation from Riihimäki and Vehtari
[15], which has been implemented in the GPstuff toolbox [23, 24]. To constrain
the classical GP to be monotonic, Riihimäki and Vehtari proposed a block
covariance matrix structure similar to that of many multi-fidelity approaches,
e.g. [25, 26, 27], and approximated the posterior using expectation
propagation [28], which is arguably more efficient than Laplace’s method. For
the sake of completeness, we assume a zero-mean GP and briefly summarize the
formulation. Readers are referred to the original work of Riihimäki and
Vehtari [15] for further details. Roughly speaking, the crux of the approach
is to augment the covariance matrix $\displaystyle\mathbf{K}$ using a block
structure, as is done in multi-fidelity and gradient-enhanced GP. Now, the
difference is that the augmented block encodes the binary classification of
whether the GP is supposed to be monotonically increasing or decreasing at
some locations. Naturally, the formulation of GP lends itself into binary
classification problem with the linear logistic regression or the probit
regression, taken from the cumulative density function of a standard normal
distribution $\displaystyle\Phi(z)=\int_{-\infty}^{t}\mathcal{N}(t|0,1)dt$.

Monotonicity is imposed at $\displaystyle M$ inducing locations
$\displaystyle\mathbf{X}_{m}\in\mathbb{R}^{M\times D}$. At the location
$\displaystyle\mathbf{x}^{(i)}\in\mathbf{X}_{m}$, the derivative of
$\displaystyle\mathbf{f}$ is non-negative with respect to the input dimension
$\displaystyle d_{i}$. The key idea is to assume a probit likelihood at the
location $\displaystyle\mathbf{x}^{(i)}$ as

$p\left(m_{d_{i}}^{(i)}\Bigg{|}\frac{\partial f^{(i)}}{\partial
x_{d_{i}}^{(i)}}\right)=\Phi\left(\frac{\partial f^{(i)}}{\partial
x_{d_{i}}^{(i)}}\frac{1}{\nu}\right),$ (14)

where

$\Phi(z)=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{z}{\sqrt{2}}\right)\right]=\int_{-\infty}^{z}\mathcal{N}(t|0,1)dt$
(15)

is the cumulative distribution function of the standard normal distribution.
Conceptually, imposing the probit likelihood is intrinsically similar to
classification using GP (cf. Section 3.6 [21] and Kuss et al. [29] for binary
classification {-1,+1} using GP). Put simply, Equation 14 models the
probability of “penalized” monotonicity with the parameter $\displaystyle\nu$,
in the same way one would in binary classification [29], where
$\displaystyle\nu=10^{-6}$ is recommended.

Assume that at $\displaystyle\mathbf{X}_{m}$, the function
$\displaystyle\mathbf{f}$ is known to be monotonic. The joint prior for
$\displaystyle\mathbf{f}$ and its derivatives
$\displaystyle\mathbf{f}^{\prime}$ is given by

$p(\mathbf{f},\mathbf{f}^{\prime}|\mathbf{X},\mathbf{X}_{m})=\mathcal{N}\left(\mathbf{f}_{\text{joint}}|\mathbf{0},\mathbf{K}_{\text{joint}}\right),$
(16)

where

$\mathbf{f}_{\text{joint}}=\begin{bmatrix}\mathbf{f}\\\
\mathbf{f}^{\prime}\end{bmatrix},\quad\mathbf{K}_{\text{joint}}=\begin{bmatrix}\mathbf{K}_{\mathbf{f},\mathbf{f}}&\mathbf{K}_{\mathbf{f},\mathbf{f}^{\prime}}\\\
\mathbf{K}_{\mathbf{f}^{\prime},\mathbf{f}}&\mathbf{K}_{\mathbf{f}^{\prime},\mathbf{f}^{\prime}}\end{bmatrix}.$
(17)

Using the Bayes rule, the joint posterior is

$p(\mathbf{f},\mathbf{f}^{\prime}|\mathbf{y},\mathbf{m})=\frac{1}{Z}p(\mathbf{f},\mathbf{f}^{\prime}|\mathbf{X},\mathbf{X}_{m})p(\mathbf{y}|\mathbf{f})p(\mathbf{m}|\mathbf{f}^{\prime}),$
(18)

where

$p(\mathbf{m}|\mathbf{f}^{\prime})=\prod_{i=1}^{M}\Phi\left(\frac{\partial
f^{(i)}}{\partial x_{d_{i}}^{(i)}}\frac{1}{\nu}\right).$ (19)

Since the posterior is analytically intractable, local likelihood
approximations are given by the expectation propagation (EP) algorithm,
allowing the approximation of the posterior distribution in Equation 18

$\begin{array}[]{lll}p(\mathbf{f},\mathbf{f}^{\prime}|\mathbf{y},\mathbf{m})&\approx&q(\mathbf{f},\mathbf{f}^{\prime}|\mathbf{y},\mathbf{m})\\\
&=&\frac{1}{Z_{\text{EP}}}p(\mathbf{f},\mathbf{f}^{\prime}|\mathbf{X},\mathbf{X}_{m})p(\mathbf{y}|\mathbf{f})\prod_{i=1}^{M}t_{i}(\tilde{Z}_{i},\tilde{\mu}_{i},\tilde{\sigma}^{2}_{i}),\end{array}$
(20)

where $\displaystyle
t_{i}(\tilde{Z}_{i},\tilde{\mu}_{i},\tilde{\sigma}^{2}_{i})=\tilde{Z}_{i}\mathcal{N}(f^{\prime}_{i}|\tilde{\mu}_{i},\tilde{\sigma}^{2}_{i})$
are local likelihood approximations with site parameters
$\displaystyle\tilde{Z}_{i},\tilde{\mu}_{i},\tilde{\sigma}$ from the EP
algorithm.

The approximate posterior is analytically tractable as a product of Gaussian
distributions and can be simplified to

$q(\mathbf{f},\mathbf{f}^{\prime}|\mathbf{y},\mathbf{m})=\mathcal{N}(\mathbf{f}_{\text{joint}}|\mathbf{\mu},\Sigma),$
(21)

where

$\mathbf{\mu}=\Sigma\tilde{\Sigma}_{\text{joint}}^{-1}\tilde{\mathbf{\mu}}_{\text{joint}},\quad\Sigma=[\mathbf{K}_{\text{joint}}^{-1}+\tilde{\Sigma}_{\text{joint}}^{-1}]^{-1},$
(22)

and

$\tilde{\mathbf{\mu}}_{\text{joint}}=\begin{bmatrix}\mathbf{y}\\\
\tilde{\mathbf{\mu}}\end{bmatrix},\quad\tilde{\Sigma}_{\text{joint}}=\begin{bmatrix}\sigma^{2}\mathbf{I}&0\\\
0&\tilde{\Sigma}\end{bmatrix},$ (23)

where $\displaystyle\tilde{\mu}$ is the vector of site means
$\displaystyle\tilde{\mu}_{i}$, and
$\displaystyle\tilde{\Sigma}=\text{Diag}[\tilde{\sigma}_{i}^{2}]_{i=1}^{M}$.
The approximation for the logarithm of the marginal likelihood is computed as

$\begin{array}[]{lll}\log
Z_{\text{EP}}&=&-\frac{1}{2}\log|\mathbf{K}_{\text{joint}}+\tilde{\Sigma}_{\text{joint}}|\\\
&&-\frac{1}{2}\tilde{\mathbf{\mu}}_{\text{joint}}^{\top}[\mathbf{K}_{\text{joint}}+\tilde{\Sigma}_{\text{joint}}]^{-1}\tilde{\mathbf{\mu}}_{\text{joint}}\\\
&&+\sum_{i=1}^{M}\frac{(\mu_{-i}-\tilde{\mu}_{i})^{2}}{2(\sigma^{2}_{-i}+\tilde{\sigma}_{i}^{2})}\\\
&&+\sum_{i=1}^{M}\log\Phi\left(\frac{\mu_{-i}}{\nu\sqrt{1+\sigma^{2}_{-i}/\nu^{2}}}\right)\\\
&&+\frac{1}{2}\sum_{i=1}^{M}\log(\sigma^{2}_{-i}+\tilde{\sigma}^{2}_{i}),\end{array}$
(24)

where $\displaystyle\mu_{-i}$ and $\displaystyle\sigma^{2}_{-i}$ are the
parameters of the cavity distribution in EP. By introducing $\displaystyle M$
monotonic inducing point, the cost complexity for optimizing the log-marginal
likelihood increases from $\displaystyle\mathcal{O}\left(N^{3}\right)$ to
$\displaystyle\mathcal{O}\left((N+M)^{3}\right)$. The posterior mean and
posterior variance of the testing distribution are, respectively, given by

${\mathbb{E}\left[f^{*}|x^{*},\mathbf{y},\mathbf{X},\mathbf{m},\mathbf{X}_{m}\right]}=\mathbf{K}_{*,\text{joint}}[\mathbf{K}_{\text{joint}}+\tilde{\Sigma}_{\text{joint}}]^{-1}\tilde{\mu}_{\text{joint}}$
(25)

and

${\mathbb{V}\left[f^{*}|x^{*},\mathbf{y},\mathbf{X},\mathbf{m},\mathbf{X}_{m}\right]}=\mathbf{K}_{*,*}-\mathbf{K}_{*,\text{joint}}[\mathbf{K}_{\text{joint}}+\tilde{\Sigma}_{\text{joint}}]^{-1}\mathbf{K}_{*,\text{joint}}.$
(26)

##  3 Numerical examples

###  3.1 Example 1: Noisy logistic regression

Figure 1: Comparison between the regular and monotonic GPs in the
homoscedastic noisy logistic regression example. Training and testing data
points are shown as black dots and magenta squares, respectively. By enforcing
the monotonicity, the monotonic GP has a better approximation compared to the
regular GP.

In this example, we consider the 1d monotonic function

$y=\frac{1}{1+e^{-x}}+\varepsilon,$ (27)

where $\displaystyle\varepsilon\sim\mathcal{N}(0,0.1^{2})$ on [-3,3] with 10
samples, as shown in Figure 1. Figure 1 shows the comparison between the
regular (dashed red) and monotonic GP (solid blue), where the training dataset
is plotted as black dots, and testing dataset is plotted as magenta squares.
The noise $\displaystyle\varepsilon$ is substantial enough to corrupt the
monotonicity of the underlying function. Despite the fact that the noise is
homoscedastic, i.e. the variance of the noise is constant throughout the
bounded domain, the regular GP has not been able to capture the function
correctly. On the other hand, the monotonic GP enforces the monotonicity while
ignoring the noise, which results in a better approximation of the true
function, which is also plotted as magenta squares as testing dataset.

###  3.2 Example 2: Heteroscedastic Hall-Petch relationship

Figure 2: Comparison between the regular and monotonic GPs in the
heteroscedastic settings. The noise $\displaystyle\varepsilon(d)$ increases as
$\displaystyle d$ increases. It is observed that the regular GP is strongly
influence.

Following Counts et al. [30], we adopt the Hall-Petch relationship in
polycrystalline materials from Fernandes and Vieira [31] with a
heteroscedastic noise proportional to $\displaystyle d$,

$\sigma_{Y}=16.47+0.0000288\frac{1}{(10^{-6}\cdot d)^{1.3}}+\varepsilon(d).$
(28)

where $\displaystyle\varepsilon(d)\sim\mathcal{N}(0,2.2\cdot 10^{10}\cdot
d^{3})$, to model the finite-size effect of the representative volume element
[6]. $\displaystyle d$ is considered on the [15$\displaystyle\mu$ m,
350$\displaystyle\mu$ m]. For a fixed-size representative volume element, when
$\displaystyle d$ is large, the number of grain reduces, and therefore, the
noise of $\displaystyle\sigma_{Y}$ increases. The same argument applies when
$\displaystyle d$ is small. Now, the variance of the Monte Carlo estimator
roughly decays at the rate of $\displaystyle N^{-1}$ where $\displaystyle N$
is the number of grains or samples; in this hypothetical example, to examine
the behavior of the monotonic GP, we simply model the volume of the grain as
proportional to $\displaystyle d^{3}$ and therefore, the noise term with zero-
mean and $\displaystyle\mathcal{O}(d^{3})$ variance is considered. Obviously,
the variance of the noise term increases as $\displaystyle d$ increases, which
is consistent to what we have observed in the literature [6].

Figure 2 shows the comparison between the regular and monotonic GP with 20
training data points, plotted as black dots. The regular GP overfits the
training dataset, which results in an oscillatory behavior toward the large
$\displaystyle d$ region, which is unphysical. The monotonic GP, on the other
hand, correctly identifies the trend and monotonically decreases as
$\displaystyle d$ increases, which results in a better approximation compared
to the regular GP.

##  4 Case study: Fatigue life prediction under multiaxial loading

###  4.1 Dataset

Training | Testing
---|---
$\displaystyle\sigma_{a}$ [MPa] | $\displaystyle N_{\text{exp}}$ |  $\displaystyle\sigma_{a}$ [MPa] | $\displaystyle N_{\text{exp}}$
674 | 2908 | 427 | 77730
558 | 8115 | 400 | 113900
556 | 10035 | 411 | 117275
507 | 17012 | 403 | 144264
483 | 19955 | 390 | 192920
505 | 20595 | 391 | 198992
498 | 23780 | 379 | 243816
490 | 25913 | 366 | 376815
484 | 28045 | 369 | 396987
474 | 51430 | 345 | 406800
469 | 52000 | 342 | 1252208
475 | 66200 | 335 | 1444998
|  | 335 | 1528487
Table 1: Training and testing fatigue datasets for S355N steel, adopted from
Karolczuk and Słoński [32].

We adopt a subset from Karolczuk and Słoński [32, 33], where the training and
testing datasets are listed in Table 1. The training dataset is taken to be 12
data points with high $\displaystyle\sigma_{a}$, whereas 13 data points with
low $\displaystyle\sigma$ are used as the testing dataset. We are interested
in investigating the extrapolation capability of GPs beyond their training
regime as well as their robustness against noise, which is fairly common in
experimental materials science, particularly in fatigue, where stochasticity
plays a non-trivial role. Much of the stochasticity can be traced back to the
physics of fracture and fatigue, where void nucleation, void growth, and void
coalescence all play important roles in the origin, development, and failure
of a material. Due to this underlying stochasticity, sometimes the data shows
a substantially noisy behavior despite a clear physical monotonic trend. That
is, the higher the stress amplitude is, the lower the fatigue life is. This
behavior is typically modeled by the SN equation as

$\log N_{\text{exp}}=A-B\log\sigma_{a},$ (29)

where $\displaystyle A$ and $\displaystyle B$ are material-dependent
coefficients, $\displaystyle\sigma_{a}$ is the stress amplitude, and
$\displaystyle N_{\text{f}}$ is the experimental fatigue life. Experimental
materials science is resource-intensive [34]; thus, practically speaking,
fatigue experimental data is scarce and may not be strictly monotonic, as
shown in Karolczuk and Słoński [32]. In the present work, the input of the GPs
is $\displaystyle\sigma_{a}$, and the output of the GPs is $\displaystyle\log
N_{\text{exp}}$.

###  4.2 Results

Figure 3: Accuracy comparison of testing dataset predictions between the
regular and monotonic GPs in S355N fatigue dataset, where the monotonic GP
outperforms the regular GP. Figure 4: Comparison between the regular and
monotonic GPs in S355N fatigue dataset. $\displaystyle\mu\pm 1\sigma$
confidence intervals are highlighted. It is observed that the monotonicity
starts to disappear as the extrapolation range goes further to
$\displaystyle\sigma_{a}\leq 375$MPa.

We implement and compare two GP variants: the regular GP and a monotonically-
constrained GP. Both are trained using a squared exponential kernel and the
hyper-parameters are optimized using the L-BFGS method [35] with a maximum of
6000 iterations. Figure 3 compares the accuracy between the two GPs against
the testing dataset. The root-mean-square errors for the regular GP and the
monotonic GP are 183.3026 and 2.0441, respectively. While the regular GP tends
to underestimate the experimental data, the monotonic GP provides more
accurate predictions. As shown in Figure 4, outside the training domain, the
regular GP is unable to capture the general trend and behaves wildly, even
when interpolating the data points. On the contrary, the monotonic GP is able
to capture the monotonic trend and extrapolate to a significant distance
outside the training interval. It is important to note that the monotonic GP
is only strictly monotonic within the training domain and may fail to maintain
monotonicity far beyond the training regime. Indeed, in Figure 4, we see that
the monotonic GP fails to maintain monotonicity when
$\displaystyle\sigma_{a}\lesssim 350$ MPa. Figure 4 also shows a non-monotonic
behavior for the regular GP.

It is also observed that the posterior variance of the monotonic GP is
significantly less than the posterior variance of the classic GP, using the
same training dataset. This holds true for both interpolation and
extrapolation cases. The main reason is that the monotonic GP uses an extra
“pseudo” training dataset by imposing $\displaystyle M$ inducing location in
the input domain, as described in Section 2.3, and a training dataset with
more data points reduces the posterior variance.

##  5 Case study: Potts Kinetic Monte Carlo for grain growth

###  5.1 Model description

The details of the Potts kinetic Monte Carlo (KMC) simulation for grain growth
and its implementation in SPPARKS is described in [36, 37], and is summarized
here for the sake of completeness. In the grain growth simulation, the Potts
model [38] is used to simulate curvature-driven grain growth. Grain
microstructures are represented by an integer value stored at each pixel.
During a timestep (referred to here as a Monte Carlo or MC step), pixels in
the simulation are visited and probabilistically change their grain membership
to a neighboring grain based on the Metropolis algorithm. The probability
$\displaystyle P$ of successful change in grain site orientation is calculated
as

$P=\begin{cases}\exp\left(\frac{-\Delta E}{k_{B}T_{s}}\right),&\text{ if
}\Delta E>0,\\\ 1,&\text{ if }\Delta E\leq 0,\end{cases}$ (30)

where $\displaystyle E$ is the total grain boundary energy calculated by
summing all the neighbor interaction energies, $\displaystyle\Delta E$ can be
regarded as the activation energy, $\displaystyle k_{B}$ is the Boltzmann
constant, and $\displaystyle T_{s}$ is the simulation temperature. In the
basic Potts model, the interaction energy between two pixels belonging to the
same grain is zero, and $\displaystyle E$ is incremented by one for each
dissimilar neighbor. From Equation 30, changes that decrease system energy are
preferred, and the total system energy is monotonically decreased through
grain coarsening. It is worthy to note that the $\displaystyle T_{s}$
simulation temperature is not the real system temperature: $\displaystyle
k_{B}T_{s}$ is an energy that defines the thermal fluctuation, i.e. noise,
presented in the kMC simulation [36]. Increasing the value of $\displaystyle
k_{B}T_{s}$ results in microstructures with increasing grain boundary
roughness.

(a) SPPARKS at 50 MC step.

(b) SPPARKS at 100 MC step.

(c) SPPARKS at 200 MC step.

(d) SPPARKS at 350 MC step.

Figure 5: Grain growth simulation via kinetic Monte Carlo (SPPARKS).

###  5.2 Results

Figure 6: Comparison between regular GP and monotonic GP in extrapolating QoI
in the SPPARKS grain growth simulation for $\displaystyle k_{B}T_{s}=0.65$.
Shaded regions correspond to 90% confidence intervals $\displaystyle\mu\pm
1.645\sigma$. Figure 7: Comparison of testing dataset predictions between the
regular and monotonic GPs in kinetic Monte Carlo SPPARKS grain growth
simulation. The output is the average grain size, whereas the input is Monte
Carlo step $\displaystyle t$. An error bar of $\displaystyle\mu\pm
1.645\sigma$, which corresponds to a 90% confidence interval is shown. The
dashed line $\displaystyle y=x$ is presented as a reference.

Using SPPARKS as the forward model, we build a dataset with $\displaystyle
k_{B}T_{s}\in\\{0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85\\}$ at multiple
snapshots in time $\displaystyle t\in\\{$ 0.000, 56.879, 106.899, 152.482,
206.331, 251.886, 306.530, 353.403, 406.482, 450.629, 502.591, 554.926,
603.345, 653.470, 706.620, 756.473, 800.112, 855.231, 901.232, 950.060,
1005.970, 1054.070, 1102.890, 1155.320, 1200.360, 1251.300, 1302.020,
1353.410, 1400.000$\displaystyle\\}$ Monte Carlo steps. In this case study,
the average grain size, measured in pixel2, is used as the quantity of
interest. Similar to Tran et al. [39], we apply a filter of 50 pixel2 to
eliminate small grains before computing the average grain size. The training
and testing datasets are divided as follows: if $\displaystyle k_{B}T_{s}\leq
0.75$ or if $\displaystyle t<1000$, the data point belongs to the training
dataset; otherwise, it belongs to the testing dataset. Under this condition, a
training dataset of 140 data points and a testing dataset of 92 data points
are obtained. The inputs of the GP are $\displaystyle(k_{B}T_{s},t)$, and the
output of the GP is the average grain size.

It is noted that the data is naturally monotonic, as grains grow in time,
which is shown in Figure 5. At the boundary for the training and testing
datasets, at $\displaystyle t=950.060$ Monte Carlo step, for $\displaystyle
k_{B}T_{s}=0.65$, the average grain size is 1203.6605 pixel2, whereas at
$\displaystyle t=1005.970$ Monte Carlo step, for $\displaystyle
k_{B}T_{s}=0.65$, the average grain size is 1279.0099 pixel2. This is clearly
shown in Figure 7, where the testing predictions are mostly more than 1250
pixel2.

Figure 6 compares the regular GP and the monotonic GP in terms of
extrapolating beyond the training regime for $\displaystyle k_{B}T_{s}=0.65$,
showing a consistent behavior between the monotonic GP and the regular GP. The
posterior means between the regular GP and the classical GP are almost exactly
the same. However, the posterior variance is significantly reduced in the
monotonic GP, which could be explained by $\displaystyle M$ extra inducing
points for monotonicity classification, as described in Section 2.3.

Figure 7 compares the regular GP and the monotonic GP in terms of accuracy
with the error bar of $\displaystyle\mu\pm 1.645\sigma$, showing that the
regular GP ($\displaystyle R^{2}=0.9730$) is slightly better than the
monotonic GP ($\displaystyle R^{2}=0.9661$). It can be explained that imposing
monotonicity to the GP formulation comes at a small cost for accuracy.

##  6 Case study: Strain-rate-dependent stress-strain with crystal plasticity
finite element

###  6.1 Model description

We adopt and extend the dataset from one of our previous studies [6],
generated from Fe-22Mn-0.6C TWIP steel with the dislocation-density-based
constitutive model. The material parameters are as described in Steinmetz et
al. [40] and summarized in Section 6.2.3 and Tables 8 and 9 in [41]. The
constitutive model was validated experimentally by Wong et al [42], and for
the sake of completeness, we briefly summarize the constitutive model here.
The TWIP/TRIP steel constitutive model is parameterized in terms of
dislocation density, $\displaystyle\varrho$, the dipole dislocation density,
$\displaystyle\varrho_{\text{di}}$, the twin volume fraction, $\displaystyle
f_{\text{tw}}$, and the $\displaystyle\varepsilon$-martensite volume fraction,
$\displaystyle f_{\text{tr}}$. A model for the plastic velocity gradient with
contribution of mechanical twinning and phase transformation was developed in
Kalidindi [43] and is given by

$\begin{array}[]{lll}\mathbf{L}_{\text{p}}&=&(1-f^{\text{tot}}_{\text{tw}}-f^{\text{tot}}_{\text{tr}})\big{(}\sum_{\alpha=1}^{N_{\text{s}}}\dot{\gamma}^{\alpha}\mathbf{s}_{\text{s}}^{\alpha}\otimes\mathbf{n}_{\text{s}}^{\alpha}\\\
&&+\sum_{\beta=1}^{N_{\text{tw}}}\dot{\gamma}\mathbf{s}_{\text{tw}}\otimes\mathbf{n}_{\text{tw}}^{\beta}+\sum_{\chi=1}^{N_{\text{tr}}}\dot{\gamma}^{\chi}\mathbf{s}_{\text{tr}}\otimes\mathbf{n}_{\text{tr}}^{\beta}\big{)},\end{array}$
(31)

where $\displaystyle\chi=1,\dots,N_{\text{tr}}$ is the
$\displaystyle\varepsilon$-martensite with volume fraction $\displaystyle
f_{\text{tr}}$ on $\displaystyle N_{\text{tr}}$ transformation systems,
$\displaystyle\mathbf{s}_{\text{s}}^{\alpha}$ and
$\displaystyle\mathbf{n}_{\text{s}}^{\alpha}$ are unit vectors along the shear
direction and shear plane normal of $\displaystyle N_{\text{s}}$ slip systems
$\displaystyle\alpha$, $\displaystyle\mathbf{s}_{\text{tw}}^{\alpha}$ and
$\displaystyle\mathbf{n}_{\text{tw}}^{\alpha}$ are those of $\displaystyle
N_{\text{tw}}$ twinning systems $\displaystyle\beta$, and
$\displaystyle\mathbf{s}_{\text{tw}}^{\alpha}$ and
$\displaystyle\mathbf{n}_{\text{tr}}^{\alpha}$ are those of $\displaystyle
N_{\text{tr}}$ transformation systems $\displaystyle\chi$. The Orowan equation
models the shear rate on the slip system $\displaystyle\alpha$ as

$\dot{\gamma}^{\alpha}=\rho_{e}b_{\text{s}}\nu_{0}\exp\left[-\frac{Q}{k_{B}T}\left\\{1-\left(\frac{\tau_{\text{eff}}^{\alpha}}{\tau_{\text{sol}}}\right)^{p}\right\\}^{q}\right],$
(32)

where $\displaystyle b_{\text{s}}$ is the length of the slip Burgers vector,
$\displaystyle\nu_{0}$ is a reference velocity, $\displaystyle Q_{\text{s}}$
is the activation energy for slip, $\displaystyle k_{B}$ is the Boltzmann
constant, $\displaystyle T$ is the temperature,
$\displaystyle\tau_{\text{eff}}$ is the effective resolved shear stress,
$\displaystyle\tau_{\text{sol}}$ is the solid solution strength, and
$\displaystyle 0<\rho_{\text{s}}\leq 1$ and $\displaystyle 1\leq
q_{\text{s}}\leq 2$ are fitting parameters controlling the glide resistance
profile. Blum and Eisenlohr [44] model the evolution of dislocation densities,
particularly the generation of unipolar dislocation density and formation of
dislocation dipoles, respectively, as

$\dot{\varrho}=\frac{|\dot{\gamma}|}{b_{\text{s}}}\Gamma_{\text{s}}-\frac{2\hat{d}{b_{\text{s}}}}{\varrho}|\dot{\gamma}|,$
(33)
$\dot{\varrho}_{\text{di}}=\frac{2(\hat{d}-\widecheck{d})}{b_{\text{s}}}\varrho|\dot{\gamma}|-\frac{2\widecheck{d}}{b_{\text{s}}}\varrho_{\text{di}}|\dot{\gamma}|-\varrho_{\text{di}}\frac{4\nu_{\text{cl}}}{\hat{d}-\widecheck{d}},$
(34)

where the glide distance below which two dislocations form a stable dipole is

$\hat{d}=\frac{3Gb_{\text{s}}}{16\pi|\tau|},$ (35)

$\displaystyle\widecheck{d}=D_{a}b_{\text{s}}$ is the distance below which two
dislocations annihilate, and the dislocation climb velocity is

$\nu_{\text{cl}}=\frac{GD_{0}V_{\text{cl}}}{\pi(1-\nu)k_{\text{B}}T}\frac{1}{\hat{d}+\widecheck{d}}\exp\left(-\frac{Q_{\text{cl}}}{k_{\text{B}}T}\right).$
(36)

Here, $\displaystyle D_{0}$ is the pre-factor of self-diffusion coefficient,
$\displaystyle V_{\text{cl}}$ is the activation volume for climb, and
$\displaystyle Q_{\text{cl}}$ is the activation energy for climb. Strain
hardening is described in terms of a dislocation mean free path, where the
mean free path is denoted by $\displaystyle\Gamma$. The mean free path for
slip is modeled as

$\frac{1}{\Gamma_{\text{s}}}=\frac{1}{D}+\frac{1}{\lambda_{\text{s}}}+\frac{1}{\lambda_{\text{tw}}}+\frac{1}{\lambda_{\text{tr}}}$
(37)

where $\displaystyle D$ is the average grain size and

$\frac{1}{\lambda_{\text{s}}^{\alpha}}=\frac{1}{i_{\text{s}}}\left(\sum_{\alpha^{\prime}=1}^{N_{\text{s}}}p^{\alpha\alpha^{\prime}}(\varrho^{\alpha^{\prime}}+\varrho_{\text{di}}^{\alpha^{\prime}})\right)^{1/2},$
(38)
$\frac{1}{\lambda_{\text{tw}}^{\alpha}}=\sum_{\beta=1}^{N_{\text{tw}}}h^{\alpha\beta}\frac{f^{\beta}_{\text{tw}}}{t_{\text{tw}}(1-f^{\text{tot}}_{\text{tw}})},\quad\frac{1}{\lambda_{\text{tr}}^{\alpha}}=\sum_{\chi=1}^{N_{\text{tr}}}h^{\alpha\chi}\frac{f^{\chi}_{\text{tr}}}{t_{\text{tr}}(1-f^{\text{tot}}_{\text{tr}})}.$
(39)

Here, $\displaystyle i_{\text{s}}$ is a fitting parameter, $\displaystyle
t_{\text{tw}}$ is the average twin thickness, and $\displaystyle
t_{\text{tr}}$ is the average $\displaystyle\varepsilon$-martensite thickness.
The mean free path for twinning and for transformation are computed,
respectively, as

$\frac{1}{\Gamma^{\beta}_{\text{tw}}}=\frac{1}{i_{\text{tw}}}\left(\frac{1}{D}+\sum_{\beta^{\prime}=1}^{N_{\text{tw}}}h^{\beta\beta^{\prime}}f_{\text{tw}}^{\beta^{\prime}}\frac{1}{t_{\text{tw}}(1-f^{\text{tot}}_{\text{tw}})}\right),$
(40)
$\frac{1}{\Gamma^{\chi}_{\text{tr}}}=\frac{1}{i_{\text{tr}}}\left(\frac{1}{D}+\sum_{\chi^{\prime}=1}^{N_{\text{tr}}}h^{\chi\chi^{\prime}}f_{\text{tr}}^{\chi^{\prime}}\frac{1}{t_{\text{tr}}(1-f^{\text{tot}}_{\text{tr}})}\right),$
(41)

with $\displaystyle i_{\text{w}}$ and $\displaystyle i_{\text{tr}}$ being
fitting parameters. The nucleation rates for twins and
$\displaystyle\varepsilon$-martensite are
$\displaystyle\dot{N}=\dot{N}_{0}P_{\text{ncs}}P$, where the probability
$\displaystyle P$ that a nucleus bows out to form a twin or
$\displaystyle\varepsilon$-martensite is

$p_{\text{tw}}=\exp\left[-\left(\frac{\hat{\tau}_{\text{tw}}}{\tau}\right)^{p_{\text{tw}}}\right],\quad
p_{\text{tr}}=\exp\left[-\left(\frac{\hat{\tau}_{\text{tr}}}{\tau}\right)^{p_{\text{tr}}}\right],$
(42)

$\displaystyle p_{\text{tw}}$ and $\displaystyle p_{\text{tr}}$ are fitting
parameters. For more details, readers are referred to Roters et al. [41] (cf.
Section 6.2.3) and Wong et al [42].

(a) The microstructure statistical volume element with equiaxed grains
investigated in this case study.

(b) A snapshot of von Mises stress for $\displaystyle D=10^{-6}m$,
$\displaystyle\dot{\varepsilon}=10^{-4}$ s-1 before homogenization.

Figure 8: A crystal plasticity finite element example by DAMASK.

Figure 8 shows the equiaxed microstructure used in this case study, which has
been generated using DREAM.3D [45] under random crystallographic textures.
DAMASK [46, 41] is used as a crystal plasticity finite element, where PETSc
[47, 48] is used as the underlying spectral solver.

###  6.2 Results

Figure 9: Homogenized stress-strain response with different
$\displaystyle(D,\dot{\varepsilon},\varepsilon)$ for a fixed microstructure
shown in Figure 8. $\displaystyle\sigma$ is monotonically increasing when
either $\displaystyle D$ decreases, or $\displaystyle\dot{\varepsilon}$
increases, or $\displaystyle\varepsilon$ increases (within the specified
$\displaystyle\varepsilon$ domain). Figure 10: Comparison of testing dataset
predictions between the regular and monotonic GPs in crystal plasticity finite
element DAMASK. The output is the homogenized stress $\displaystyle\sigma$,
whereas the inputs are average grain size, strain-rate, and strain,
respectively, $\displaystyle(D,\dot{\varepsilon},\varepsilon)$.

In this study, we focus on the effect of strain-rate on the stress-strain
curve with various average grain sizes. The effect of strain-rate has been
experimentally validated in Benzing et al [49], and some numerical studies
using crystal plasticity finite element model have been done as well, such as
those found in Singh et al [50]. The inputs of the GP are the average grain
size, the strain-rate, the strain, i.e.
$\displaystyle(D,\dot{\varepsilon},\varepsilon)$, and the output of the GP is
the stress $\displaystyle\sigma$. We vary independently the average grain size
$\displaystyle D\in\\{1.0\cdot 10^{-6},2.5\cdot 10^{-6},5.0\cdot
10^{-6},1.0\cdot 10^{-5},2.5\cdot 10^{-5}\\}$ m, as well as the strain-rate
$\displaystyle\dot{\varepsilon}\in\\{10^{-5},10^{-4},10^{-3},10^{-2},10^{-1},10^{0},10^{1},10^{2},10^{3},10^{4}\\}$
s-1. Numerically unstable results are removed in the post-process. Some
examples of homogenized stress-strain responses are shown in Figure 9 for
visualization purpose.

The training and testing datasets are divided as follows: if the average grain
size $\displaystyle D$ is less than $\displaystyle 5\cdot 10^{-6}$m or if the
strain-rate $\displaystyle\dot{\varepsilon}=10^{1}$ s-1, the data point
belongs to the training dataset; otherwise, it belongs to the testing dataset.
Under this splitting condition, a training dataset of 340 data points and a
testing dataset of 120 data points are obtained.

Figure 10 shows the accuracy comparison between the regular GP and the
classical GP. The $\displaystyle R^{2}$ coefficient for the regular GP is
0.9895, whereas the $\displaystyle R^{2}$ coefficient for the monotonic GP is
slightly lower, 0.9860. Both GPs do very well, but the monotonic GP is
slightly worse. The regular GP predictions are also more consistent in the
large stress regime $\displaystyle\sigma$.

##  7 Discussion and Conclusion

In this paper, we compare the monotonic GP against the regular GP for various
materials datasets, namely fatigue, grain growth, and homogenized stress-
strain behaviors. The fatigue dataset is the only experimental dataset with
noise, whereas the other two datasets, i.e. grain growth and homogenized
stress-strain, are computationally generated by SPPARKS and DAMASK,
respectively. We test out their predictive capabilities in both interpolation
and extrapolation, but are mainly concerned with extrapolation.

For the fatigue example featured in Section 4, the experimental data does not
adhere strictly to theoretical derivations. This is where the monotonic GP
outperforms the regular GP significantly. In the other two examples, which
featured datasets generated by computational models, both GPs perform very
closely with each other, even though the regular GP performs slightly better
than the monotonic GP. This can be explained by imposing the monotonicity on
the regular GP comes a small cost for accuracy and numerical performance. It
is also noted that the monotonic behavior in the monotonic GP is only
guaranteed within the training domain. Beyond the training regime, the
monotonicity effect imposed by introducing $\displaystyle M$ inducing points
fades away: the further the monotonic GP goes beyond the training regime, the
more likely it is the GP will violate the monotonicity constraints. Last but
not least, the posterior variances of the monotonic GP are significantly
smaller than the posterior variances of the regular GP, while the posterior
means are roughly the same.

Constructing an accurate and physically faithful surrogate model is important
for uncertainty quantification [51]. Computationally efficient surrogate
models are commonly used, for example, in sampling [3, 4] or to solve a
stochastic inverse problem in structure-property relationship [6, 7].
Monotonicity is common in materials science, such as the famous Hall-Petch
relationship [2], i.e.
$\displaystyle\sigma_{y}=\sigma_{0}+k\frac{1}{\sqrt{D}}$. Exploiting such
physical constraints and insights may improve the efficiency of machine
learning in the future.

The objective of this work is to examine the performance of the monotonic GP
compared to the regular GP in different settings, with an emphasis on
materials science. Our conclusion is that the monotonic GP is best suited for
sparse and noisy settings, which is fairly common in experimental materials
science, where the data acquisition process is expensive (and intensive in
terms of both time and labor perspectives) and the data is often corrupted by
the noise induced by the finite-size effect of microstructures.

It is worthy to note that the monotonic GP indeed performs slightly worse in
the second and third case studies. In these case studies, both the training
and testing datasets are already monotonic. With the use of linear basis
function for the global trend, the regular GP correctly picks up the monotonic
behavior, and therefore, is able to make accurate prediction. While the
monotonic GP shares the same settings with the regular GP, there are some
approximations involved in constructing the joint covariance matrix. In
ideally monotonic cases such as the second and third examples, the covariance
matrix should asymptotically reduce to the diagonal block matrix as
$\displaystyle\nu\to 0$

$\mathbf{K}_{\text{joint}}=\begin{bmatrix}\mathbf{K}_{\mathbf{f},\mathbf{f}}&\mathbf{K}_{\mathbf{f},\mathbf{f}^{\prime}}\\\
\mathbf{K}_{\mathbf{f}^{\prime},\mathbf{f}}&\mathbf{K}_{\mathbf{f}^{\prime},\mathbf{f}^{\prime}}\end{bmatrix}\xrightarrow[\nu\to
0]{}\begin{bmatrix}\mathbf{I}_{\mathbf{f},\mathbf{f}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{\mathbf{f}^{\prime},\mathbf{f}^{\prime}}\end{bmatrix},$
(43)

which points to how the step function is approximated. In practice, a small
numerical degradation is expected as $\displaystyle\nu$ cannot be smaller than
machine epsilon. In the monotonic GP approach adopted in this paper, the step
function is approximated by a probit function described in Equation 14. It
might be possible to increase the performance of the monotonic GP by even
further reducing the $\displaystyle\nu$ parameter. This may lead to an
improvement in performance for the monotonic GP with monotonic datasets, but
might also lead to other numerical conditioning problems with respect to the
inverse of the joint covariance matrix
$\displaystyle\mathbf{K}_{\text{joint}}$. We conclude that for datasets that
are already monotonic, applying a monotonic GP formulation may be unnecessary
and may lead to a slight degradation in performance. For datasets that are
sparse and noisy, imposing monotonicity when it is appropriate may result in a
substantial improvement in performance.

## Acknowledgment

This work was supported by the Laboratory Directed Research and Development
program at Sandia National Laboratories, a multimission laboratory managed and
operated by National Technology and Engineering Solutions of Sandia LLC, a
wholly owned subsidiary of Honeywell International Inc. for the U.S.
Department of Energy’s National Nuclear Security Administration under contract
DE-NA0003525.

## References

  * [1] National Science and Technology Council (US), 2011.  Materials Genome Initiative for global competitiveness.  Executive Office of the President, National Science and Technology Council.
  * [2] Cordero, Z. C., Knight, B. E., and Schuh, C. A., 2016.  “Six decades of the Hall–Petch effect–a survey of grain-size strengthening studies on pure metals”.  International Materials Reviews, 61(8), pp. 495–512.
  * [3] Tallman, A. E., Stopka, K. S., Swiler, L. P., Wang, Y., Kalidindi, S. R., and McDowell, D. L., 2019.  “Gaussian-process-driven adaptive sampling for reduced-order modeling of texture effects in polycrystalline alpha-Ti”.  JOM, 71(8), pp. 2646–2656.
  * [4] Tallman, A. E., Swiler, L. P., Wang, Y., and McDowell, D. L., 2020.  “Uncertainty propagation in reduced order models based on crystal plasticity”.  Computer Methods in Applied Mechanics and Engineering, 365, p. 113009.
  * [5] Yabansu, Y. C., Iskakov, A., Kapustina, A., Rajagopalan, S., and Kalidindi, S. R., 2019.  “Application of Gaussian process regression models for capturing the evolution of microstructure statistics in aging of nickel-based superalloys”.  Acta Materialia, 178, pp. 45–58.
  * [6] Tran, A., and Wildey, T., 2020.  “Solving stochastic inverse problems for property-structure linkages using data-consistent inversion and machine learning”.  JOM, 73, pp. 72–89.
  * [7] Tran, A., and Wildey, T., 2021.  “Solving stochastic inverse problems for property-structure linkages using data-consistent inversion”.  In TMS 2021 150th Annual Meeting & Exhibition, Springer, pp. 1–8.
  * [8] Tran, A., Tranchida, J., Wildey, T., and Thompson, A. P., 2020.  “Multi-fidelity machine-learning with uncertainty quantification and Bayesian optimization for materials design: Application to ternary random alloys”.  The Journal of Chemical Physics, 153, p. 074705.
  * [9] Khatamsaz, D., Molkeri, A., Couperthwaite, R., James, J., Arróyave, R., Allaire, D., and Srivastava, A., 2021.  “Efficiently exploiting process-structure-property relationships in material design by multi-information source fusion”.  Acta Materialia, 206, p. 116619.
  * [10] Fernández-Godino, M. G., Balachandar, S., and Haftka, R., 2018.  “On the use of symmetries in building surrogate models”.  Journal of Mechanical Design.
  * [11] Jidling, C., Wahlström, N., Wills, A., and Schön, T. B., 2017.  “Linearly constrained Gaussian processes”.  arXiv preprint arXiv:1703.00787.
  * [12] Agrell, C., 2019.  “Gaussian processes with linear operator inequality constraints”.  arXiv preprint arXiv:1901.03134.
  * [13] Lange-Hegermann, M., 2021.  “Linearly constrained Gaussian processes with boundary conditions”.  In International Conference on Artificial Intelligence and Statistics, PMLR, pp. 1090–1098.
  * [14] Swiler, L. P., Gulian, M., Frankel, A. L., Safta, C., and Jakeman, J. D., 2020.  “A survey of constrained gaussian process regression: Approaches and implementation challenges”.  Journal of Machine Learning for Modeling and Computing, 1(2).
  * [15] Riihimäki, J., and Vehtari, A., 2010.  “Gaussian processes with monotonicity information”.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, pp. 645–652.
  * [16] Golchi, S., Bingham, D. R., Chipman, H., and Campbell, D. A., 2015.  “Monotone emulation of computer experiments”.  SIAM/ASA Journal on Uncertainty Quantification, 3(1), pp. 370–392.
  * [17] Ustyuzhaninov, I., Kazlauskaite, I., Ek, C. H., and Campbell, N., 2020.  “Monotonic Gaussian process flows”.  In International Conference on Artificial Intelligence and Statistics, PMLR, pp. 3057–3067.
  * [18] Pensoneault, A., Yang, X., and Zhu, X., 2020.  “Nonnegativity-enforced Gaussian process regression”.  Theoretical and Applied Mechanics Letters, 10(3), pp. 182–187.
  * [19] Tan, M. H. Y., 2017.  “Monotonic metamodels for deterministic computer experiments”.  Technometrics, 59(1), pp. 1–10.
  * [20] Chen, Y., Hosseini, B., Owhadi, H., and Stuart, A. M., 2021.  “Solving and learning nonlinear PDEs with Gaussian processes”.  Journal of Computational Physics, 447, p. 110668.
  * [21] Rasmussen, C. E., 2006.  Gaussian processes in machine learning.  MIT Press.
  * [22] Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., and de Freitas, N., 2016.  “Taking the human out of the loop: A review of Bayesian optimization”.  Proceedings of the IEEE, 104(1), pp. 148–175.
  * [23] Vanhatalo, J., Riihimäki, J., Hartikainen, J., Jylänki, P., Tolvanen, V., and Vehtari, A., 2012.  “Bayesian modeling with Gaussian processes using the GPstuff toolbox”.  arXiv preprint arXiv:1206.5754.
  * [24] Vanhatalo, J., Riihimäki, J., Hartikainen, J., Jylänki, P., Tolvanen, V., and Vehtari, A., 2013.  “GPstuff: Bayesian modeling with Gaussian processes”.  Journal of Machine Learning Research, 14(Apr), pp. 1175–1179.
  * [25] Tran, A., Wildey, T., and McCann, S., 2020.  “sMF-BO-2CoGP: A sequential multi-fidelity constrained Bayesian optimization for design applications”.  Journal of Computing and Information Science in Engineering, 20(3), pp. 1–15.
  * [26] Yang, X., Zhu, X., and Li, J., 2020.  “When bifidelity meets CoKriging: An efficient physics-informed multifidelity method”.  SIAM Journal on Scientific Computing, 42(1), pp. A220–A249.
  * [27] Xiao, M., Zhang, G., Breitkopf, P., Villon, P., and Zhang, W., 2018.  “Extended co-kriging interpolation method based on multi-fidelity data”.  Applied Mathematics and Computation, 323, pp. 120–131.
  * [28] Minka, T. P., 2001.  “Expectation propagation for approximate Bayesian inference”.  In Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, UAI’01, AUAI Press, pp. 362–369.
  * [29] Kuss, M., Rasmussen, C. E., and Herbrich, R., 2005.  “Assessing approximate inference for binary gaussian process classification.”.  Journal of Machine Learning Research, 6(10).
  * [30] Counts, W. A., Braginsky, M. V., Battaile, C. C., and Holm, E. A., 2008.  “Predicting the Hall–Petch effect in fcc metals using non-local crystal plasticity”.  International journal of plasticity, 24(7), pp. 1243–1263.
  * [31] Fernandes, J., and Vieira, M., 2000.  “Further development of the hybrid model for polycrystal deformation”.  Acta materialia, 48(8), pp. 1919–1930.
  * [32] Karolczuk, A., and Słoński, M., 2022.  “Application of the Gaussian process for fatigue life prediction under multiaxial loading”.  Mechanical Systems and Signal Processing, 167, p. 108599.
  * [33] Karolczuk, A., and Kluger, K., 2020.  “Application of life-dependent material parameters to lifetime calculation under multiaxial constant-and variable-amplitude loading”.  International Journal of Fatigue, 136, p. 105625.
  * [34] Arróyave, R., and McDowell, D. L., 2019.  “Systems approaches to materials design: Past, present, and future”.  Annual Review of Materials Research, 49(1), pp. 103–126.
  * [35] Zhu, C., Byrd, R. H., Lu, P., and Nocedal, J., 1997.  “Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization”.  ACM Transactions on mathematical software (TOMS), 23(4), pp. 550–560.
  * [36] Garcia, A. L., Tikare, V., and Holm, E. A., 2008.  “Three-dimensional simulation of grain growth in a thermal gradient with non-uniform grain boundary mobility”.  Scripta Materialia, 59(6), pp. 661–664.
  * [37] Plimpton, S., Battaile, C., Chandross, M., Holm, L., Thompson, A., Tikare, V., Wagner, G., Webb, E., Zhou, X., Cardona, C. G., et al., 2009.  “Crossing the mesoscale no-man’s land via parallel kinetic Monte Carlo”.  Sandia Report SAND2009-6226.
  * [38] Anderson, M., Grest, G., and Srolovitz, D., 1989.  “Computer simulation of normal grain growth in three dimensions”.  Philosophical Magazine B, 59(3), pp. 293–329.
  * [39] Tran, A., Mitchell, J. A., Swiler, L. P., and Wildey, T., 2020.  “An active-learning high-throughput microstructure calibration framework for process-structure linkage in materials informatics”.  Acta Materialia, 194, pp. 80–92.
  * [40] Steinmetz, D. R., Jäpel, T., Wietbrock, B., Eisenlohr, P., Gutierrez-Urrutia, I., Saeed-Akbari, A., Hickel, T., Roters, F., and Raabe, D., 2013.  “Revealing the strain-hardening behavior of twinning-induced plasticity steels: Theory, simulations, experiments”.  Acta Materialia, 61(2), pp. 494–510.
  * [41] Roters, F., Diehl, M., Shanthraj, P., Eisenlohr, P., Reuber, C., Wong, S. L., Maiti, T., Ebrahimi, A., Hochrainer, T., Fabritius, H.-O., et al., 2019.  “DAMASK–The Düsseldorf Advanced Material Simulation Kit for modeling multi-physics crystal plasticity, thermal, and damage phenomena from the single crystal up to the component scale”.  Computational Materials Science, 158, pp. 420–478.
  * [42] Wong, S. L., Madivala, M., Prahl, U., Roters, F., and Raabe, D., 2016.  “A crystal plasticity model for twinning-and transformation-induced plasticity”.  Acta Materialia, 118, pp. 140–151.
  * [43] Kalidindi, S. R., 1998.  “Incorporation of deformation twinning in crystal plasticity models”.  Journal of the Mechanics and Physics of Solids, 46(2), pp. 267–290.
  * [44] Blum, W., and Eisenlohr, P., 2009.  “Dislocation mechanics of creep”.  Materials Science and Engineering: A, 510, pp. 7–13.
  * [45] Groeber, M. A., and Jackson, M. A., 2014.  “DREAM. 3D: a digital representation environment for the analysis of microstructure in 3D”.  Integrating materials and manufacturing innovation, 3(1), p. 5.
  * [46] Diehl, M., Groeber, M., Haase, C., Molodov, D. A., Roters, F., and Raabe, D., 2017\.  “Identifying structure–property relationships through DREAM.3D representative volume elements and DAMASK crystal plasticity simulations: An integrated computational materials engineering approach”.  JOM, 69(5), pp. 848–855.
  * [47] Abhyankar, S., Brown, J., Constantinescu, E. M., Ghosh, D., Smith, B. F., and Zhang, H., 2018.  “PETSc/TS: A modern scalable ODE/DAE solver library”.  arXiv preprint arXiv:1806.01437.
  * [48] Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W., et al., 2019.  “PETSc users manual”.
  * [49] Benzing, J. T., Poling, W. A., Pierce, D. T., Bentley, J., Findley, K. O., Raabe, D., and Wittig, J. E., 2018.  “Effects of strain rate on mechanical properties and deformation behavior of an austenitic Fe-25Mn-3Al-3Si TWIP-TRIP steel”.  Materials Science and Engineering: A, 711, pp. 78–92.
  * [50] Singh, L., Vohra, S., and Sharma, M., 2021.  “Investigation of strain rate behavior of aluminium and AA2024 using crystal plasticity”.  Materials Today: Proceedings.
  * [51] Acar, P., 2021.  “Recent progress of uncertainty quantification in small-scale materials science”.  Progress in Materials Science, 117, p. 100723.