# Conversion of stable crystals to metastable crystals in a solution by
periodic change of temperature

Hiroyasu Katsuno<EMAIL_ADDRESS>Institute of Low Temperature
Science, Hokkaido University, Kita-19, Nishi-8, Kita-ku, Sapporo, Hokkaido,
060-0819, Japan  Makio Uwaha  Science Division, Center for General Education,
Aichi Institute of Technology, 1247 Yachigusa, Yakusa-cho, Toyota, Aichi
470-0392, Japan

###### Abstract

Using a Becker-Döring type model including cluster incorporation, we study the
possibility of conversion of stable crystals to metastable crystals in a
solution by a periodic change of temperature. At low temperature, both stable
and metastable crystals grow by coalescence with abundant clusters. At high
temperature, a large amount of small clusters produced by the dissolution of
crystals inhibits the dissolution of crystals, and the imbalance in the amount
of crystals increases. By repeating this process, the periodic temperature
change can convert stable crystals into metastable crystals.

##  I Introduction

Among many possible crystal structures, the most stable structure is realized
and metastable structures disappear in equilibrium. Exceptional cases are
chiral crystals which have two stable structures that are thermodynamically
equivalent. When the most stable structure is chiral, it is inevitable that
both left- and right-handed crystals exist in equilibrium. In 2005, Viedma
demonstrated the conversion of a racemic mixture of chiral sodium chlorate
crystals into homochiral crystals by grinding crystals in a solution[1]. The
chirality conversion by grinding is also possible in conglomerate forming
chiral molecules that racemize in a solution[2, 3]. In contrast with Ostwald
ripening observed in the process of relaxation to equilibrium, it is called
Viedma ripening (VR). Several years later, some researchers showed a similar
phenomenon with periodic temperature change (temperature cycling: TC) of a
solution. The boiling of a solution with powder crystals is a simple way to
realize TC [4, 5]. More quantitative experiments were performed with
controlled temperature change of a solution[6, 7]. One of the important
features of VR and TC is the exponential amplification of an initial small
crystal enantiomeric excess during the conversion.

From the theoretical analysis[8] of chemical reactions that bring complete
homochirality, it is likely that important factors to realize a homochiral
state in VR/TC are also nonlinear autocatalysis and recycling of the product.
Nonlinear autocatalysis increases the asymmetry in the ratio of two
enantiomers and the product is recycled to produce the dominant enantiomer. In
VR/TC experiment, the enhancement of dissolution of crystals by grinding
crystals/increasing temperature corresponds to the recycling process. Then,
the question is what processes in crystallization correspond to the nonlinear
autocatalytic process. Several mechanisms for VR have been proposed such as
the chiral cluster incorporation in crystallization[9], the catalytic surface
reaction[10], the mutual inhibition[11] based on the Frank model[12], and the
secondary nucleation by the shear stress of the flow[13]. The chiral cluster
incorporation is adopted for the explanation of VR in various approaches: the
simple rate equation[9, 14, 15, 16, 17, 18], the time change of the crystal
size distribution such as a Backer-Döring (BD) model[19, 20, 21, 22, 23] and a
population balance (PB) model[24], and Monte Carlo simulation[25, 26]. The
cluster incorporation mechanism developed in BD type models also explains the
chirality conversion by TC for achiral and chiral molecules[22, 27]. For the
crystallization of chiral molecules, simple molecular incorporation without
clusters described by a PB model reproduces the chirality conversion in TC[28,
29]. There is no common view on the mechanism of TC at the moment[30].

While the essential process in TC for chirality conversion is still not very
clear, the chirality conversion by TC is completely different from the
relaxation to equilibrium. The results of TC can be interpreted as the system
chooses one state from two thermodynamically equivalent states. The final
state is determined by the initial state and dynamics of the system. Here, a
simple question arises: is it possible to choose an energetically unfavorable
state from two thermodynamically inequivalent states in TC? A recent
experiment with the combination of VR and TC shows an extraordinary phenomenon
that stable racemic crystals are converted to metastable chiral crystals[31].
In a numerical study, the conversion to the metastable phase crystals by VR is
confirmed[32]. These results show that the conversion of the phase of crystals
(hereafter, we call phase conversion) is possible by the simple methods, and
that even metastable phase crystals win if the growth rate of metastable
crystals becomes larger than that of stable phase crystals.

In this paper, we study the phase conversion by TC in a solution which
contains the stable and the metastable phase crystals. We discuss the
mechanism of the phase conversion in detail.

##  II Model

Figure 1:  Schematics of our model. “Crystals” means clusters whose size is
larger than $i_{\rm s}$ and do not coagulate with other “crystals”.

For the investigation of the phase conversion by TC, we use a generalized BD
model with incorporation of clusters in crystallization. The original BD model
describes the change of the cluster size distribution through monomers[33, 34,
35]. Our generalization of the BD model are i) introduction of the stable and
the metastable crystal phases and ii) cluster incorporation into crystals of
the same phase as schematically shown in Fig. 1. In our model, the stable and
metastable phases are characterized by the two parameters: the equilibrium
number of monomers (solubility)111 In the original BD model, the time change
of the number density is described. Here, we call the number density the
number for simplicity.  and the interfacial energy. Although two sets of
parameters are necessary for the stable and the metastable phases, we assume
the interfacial energies of both phases are the same. This setting ensures
that the metastable phase crystals are always energetically unfavorable for
any given size. The solubilities of the stable phase $n_{1,\rm s}^{\rm eq}$
and the metastable phase $n_{1,\rm m}^{\rm eq}$ are related to the difference
in the chemical potential: $\Delta\mu_{\rm m}-\Delta\mu_{\rm s}=k_{\rm
B}T\ln\left(n_{\rm 1,s}^{\rm eq}/n_{\rm 1,m}^{\rm eq}\right)<0$, where
$\Delta\mu_{\rm\gamma}$ is the difference in chemical potential between the
$\gamma$-phase crystal and the solution, $k_{\rm B}$ is the Boltzmann
constant, and $T$ is temperature. Hereafter, we call $n_{1,\rm s}^{\rm
eq}/n_{1,\rm m}^{\rm eq}(<1)$ the solubility ratio.

For the generalization ii), the acceptable incorporation process depends on
the cluster size $i$. Small clusters, whose size is smaller than $i_{\rm s}$,
can be incorporated into clusters/crystals (we call small clusters simply
clusters). Large clusters, whose size is larger than $i_{\rm s}$, can not be
incorporated into other crystals (we call large clusters simply crystals) .

In the numerical calculation, we investigate the time change of the number
$n_{i,\gamma}$, where $\gamma=\textrm{m}$ and $\gamma=\textrm{s}$ represent
the metastable phase and the stable phase. The change of $n_{i,\gamma}$ by
incorporation of a monomer or other clusters is given by
$\sigma_{i,j}n_{i,\gamma}n_{j,\gamma}$, where the incorporation rate
$\sigma_{i,j}$ is proportional to the collision cross-section of the size $i$
and $j$, i.e., $\sigma_{i,j}=ai^{2/3}j^{2/3}$ with a coefficient $a$. The
change by dissolution is $\lambda_{i,j}^{\gamma}n_{i,\gamma}$, where the
dissolution rate $\lambda_{i,j}^{\gamma}$ of $j$-mer from $i$-mer is
determined from the detailed balance condition:

$\displaystyle\lambda_{i,j}^{\gamma}n_{i,\gamma}^{\rm
eq}=\sigma_{i-j,j}n_{i-j,\gamma}^{\rm eq}n_{j,\gamma}^{\rm eq},$ (1)

where $n_{i,\gamma}^{\rm eq}$ is the equilibrium Boltzmann distribution and we
have assumed that $\sigma_{i,j}$ is independent of $\gamma$ for simplicity.
The Boltzmann distribution is given by the solubility $n_{1,\gamma}^{\rm eq}$
and the free energy of an $i$-mer

$\displaystyle n_{i,\gamma}^{\rm eq}=n_{1,\gamma}^{\rm
eq}\exp[-\bar{\alpha}(i^{2/3}-1)],$ (2)

where $\bar{\alpha}$ corresponds to the interfacial energy divided by $k_{\rm
B}T$. Change of the number of monomers by the incorporation/dissolution
process is written as

$\displaystyle\frac{\partial n_{1}}{\partial
t}=-2\sigma_{1,1}(n_{1})^{2}+\sum_{\gamma={\rm s},{\rm
m}}\left[-\sum_{j=2}^{i_{\rm
max}-1}\sigma_{1,j}n_{1}n_{j,\gamma}+2\lambda_{2,1}^{\gamma}n_{2,\gamma}+\sum_{j=2}^{i_{\max}-1}\lambda_{j+1,j}^{\gamma}n_{j+1,\gamma}\right],$
(3)

The number of monomers is denoted by $n_{1}(=n_{1,\rm m}=n_{1,\rm s})$ because
monomers are common. For simplicity, it is assumed that a dimer formed by the
coalescence of monomers belongs to the stable or the metastable phase with the
equal probability. Change of the number of clusters with $2\leq i\leq i_{\rm
s}$ is written as

$\displaystyle\frac{\partial n_{i,\gamma}}{\partial
t}=\sum_{j=1}^{[(i+1)/2]}\left(\sigma_{i-j,j}^{\prime}n_{i-j,\gamma}n_{j,\gamma}-\lambda_{i,j}^{\gamma}n_{i,\gamma}\right)-\sum_{j=1}^{i_{\rm
max}-i}\left[\left(\sigma_{i,j}+\delta_{i,j}\sigma_{i,i}\right)n_{i,\gamma}n_{j,\gamma}-\lambda_{i+j,j}^{\gamma}n_{i+j,\gamma}\right],$
(4)

where $\sigma_{i,j}^{\prime}$ is defined as

$\displaystyle\sigma_{i,j}^{\prime}=\left\\{\begin{array}[]{cc}\sigma_{1,1}/2&\textrm{if}\
\ i=j=1\\\ \sigma_{i,j}&\textrm{otherwise},\end{array}\right.$ (7)

$\delta_{i,j}$ is the Kronecker delta, and $i_{\rm max}$ is the maximum size
of crystals preset in our numerical calculation. Change of the number of
crystals with $i>i_{\rm s}$ is

$\displaystyle\frac{\partial n_{i,\gamma}}{\partial
t}=\sum_{j=1}^{\textrm{min}\\{i_{\rm
s},[(i+1)/2]\\}}\left(\sigma_{i-j,j}n_{i-j,\gamma}n_{j,\gamma}-\lambda_{i,j}^{\gamma}n_{i,\gamma}\right)-\sum_{j=1}^{i_{\rm
s},i+j\leq i_{\rm
max}}\left(\sigma_{i,j}n_{i,\gamma}n_{j,\gamma}-\lambda_{i+j,j}^{\gamma}n_{i+j,\gamma}\right).$
(8)

We periodically change temperature $T$ of the system. The temperature profile
is simple: $T$ is constant and low in the first half of the period $P$ and
high in the second half. The solubility $n_{1,\gamma}^{\rm eq}$ at low
temperature is smaller than that at high temperature. The effective
interfacial energy $\bar{\alpha}$ at low temperature is larger than
$\bar{\alpha}$ at high temperature. A weak temperature dependence of the
constant $a$ in the incorporation rate is ignored.

The parameters in our numerical calculation are as follows: The constant in
the incorporation rate is $a=1$. In the low temperature, the solubility of the
metastable crystals is $n_{\rm 1,m}^{\rm eq}=10^{-2}$ and the effective
interfacial energy is $\bar{\alpha}=5$. In the high temperature, $n_{\rm
1,m}^{\rm eq}=1.5\times 10^{-2}$ and $\bar{\alpha}=0.5$. The period of the
temperature cycle is $P=10$. The total mass is conserved and normalized:
$n_{1}+\sum_{i>1,\gamma}in_{i,\gamma}=1$. The maximum crystal size is $i_{\rm
max}=1000$, and the maximum cluster size is $i_{\rm s}=10$. The initial
condition is that all monomers are in crystals and they have the same size
$i=500$.

##  III Numerical results & discussion

The state of the system is indicated by the relative crystal mass difference
between the stable and the metastable crystals

$\displaystyle\phi=\frac{M^{\rm m}-M^{\rm s}}{M^{\rm m}+M^{\rm s}},$ (9)

where the mass of $\gamma$-phase crystals is defined by

$\displaystyle M^{\rm\gamma}=\sum_{i>i_{\rm s}}^{i_{\rm max}}in_{i,\gamma}.$
(10)

We call $\phi$ the excess parameter in the present paper.

We demonstrate the conversion of the stable crystals to the metastable
crystals by TC when the difference in solubility is small. The time change of
the excess parameter $\phi$ with the solubility ratio $n_{\rm 1,\rm s}^{\rm
eq}/n_{\rm 1,m}^{\rm eq}=0.99$ is shown in Fig. 2. When the initial condition
is $\phi(0)=0.05$ (open triangle in Fig. 2), all stable crystals are converted
to metastable crystals. When the initial condition is $\phi(0)=0.01$ (filled
inverted triangle in Fig. 2), all metastable crystals are converted to the
stable crystals.

In the previous work[32], it was shown that the phase conversion of the stable
crystals to the metastable crystals is possible by VR. Using similar parameter
values of the solubilities and the interfacial energies (see Fig. 4 in [32]),
the conversion by TC is faster than that by VR. In an experiment of the
chirality conversion, the conversion time of TC is about $20$ times shorter
than that of VR[37]. Our result is consistent with the experimental result
although direct correspondence between the theoretical and the experimental
parameter values cannot be made.

Figure 2:  Time change of the excess parameter $\phi$. Open triangle:
$\phi(0)=0.05$, filled inverted triangle: $\phi(0)=0.01$. The dotted line
shows $\phi=0$.

In our result shown in Fig. 2, the initial excess parameter $\phi(0)$ is
amplified exponentially in time. The behavior of the excess parameter is
similar to that of the crystal enantiomeric excess of chiral crystals during
VR and TC. However, the metastable phase crystals can win only when $\phi(0)$
is above a critical value.

Fig. 3 shows the crystal mass distributions $in_{i,\gamma}$ at
$t=500,1000,1500$. The corresponding time change of the excess parameter
$\phi$ is shown with open triangles in Fig. 2. Red (light gray) and blue (dark
gray) areas represent mass distributions of the metastable crystals and the
stable crystals. The initial crystal distribution is $in_{i,\rm
m}=0.525\delta_{i,500}$ and $in_{i,\rm s}=0.475\delta_{i,500}$ with the
initial excess parameter $\phi(0)=0.05$. In the initial relaxation ($t\leq
50$), the crystals dissolve rapidly due to undersaturation. Then, the
distribution takes a monotonically decreasing form at the end of high
temperature period as shown in Fig. 3(a). Thereafter, the mass difference
increases as the whole distribution gradually spreads [Fig. 3(b)]. Finally,
the metastable crystals dominate the system, and the stable crystals disappear
completely [Fig. 3(c)].

Figure 3:  The crystal mass distribution at (a) $t=500$ ($\phi=0.18$), (b)
$t=1000$ ($\phi=0.65$), and (c) $t=1500$($\phi=1$). The initial distribution
is $in_{i,\gamma}=\delta_{i,500}(1\pm\phi)/2$ with $\phi=0.05$. Red (light
gray) and blue (dark gray) areas represent the masses of metastable crystals
and stable crystals.  Figure 4:  Time change of (a) the excess parameter
$\phi$, (b) the incorporation $\Sigma_{100}^{\gamma}$ and the dissolution
$\Lambda_{100}^{\gamma}$ of a crystal of the size $i=100$ [see Eq. (11)] in
the 51st period. The low temperature period is $500\leq t<505$, and the high
temperature period is $505\leq t<510$. The dotted line shows the magnitude of
the excess parameter at $t=500$. Open diamonds and filled squares represent
$\Sigma_{100}^{\gamma}$ for the metastable crystals and the stable crystals,
respectively. Open inverted triangles and filled triangles represent
$\Lambda_{100}^{\gamma}$ for a metastable crystal and a stable crystal,
respectively.

To investigate the process in one temperature cycle, the time change of the
excess parameter in the 51st period is plotted in Fig. 4(a). The low
temperature period is $500\leq t<505$, and the high temperature period is
$505\leq t<510$. The dotted line shows the magnitude of the excess parameter
$\phi$ at the beginning of the 51st period, and, at the end of the cycle,
$\phi$ slightly increases compared to the beginning. Both type of crystals
grow thanks to the supersaturation of monomers/clusters which have been
supplied in the high temperature period. It is necessary to calculate the
supersaturation in order to determine which crystal is more favorable for
growth. However, the conventional definition of supersaturation using the
number of monomers is inadequate because of the presence of cluster
incorporation. Instead, we use the growth rate of a given size crystal as the
degree of saturation. The growth rate $v_{i}^{\gamma}$ of a size $i$-mer of
the phase $\gamma$ is defined by the rate of incorporation/dissolution of
monomers and clusters:

$\displaystyle v_{i}^{\gamma}$ $\displaystyle=\sum_{j=1}^{i_{\rm
s}}j\left(\sigma_{i,j}n_{j,\gamma}-\lambda_{i,j}^{\gamma}\right)\equiv\Sigma_{i}^{\gamma}-\Lambda_{i}^{\gamma},$
(11)

where $\Sigma^{\gamma}_{i}$ and $\Lambda^{\gamma}_{i}$ represent the sum of
the incorporation terms $j\sigma_{i,j}n_{j}^{\gamma}$ and the sum of the
dissolution terms $j\lambda_{i,j}^{\gamma}$, respectively. The conventional
relation between the supersaturation and the growth rate is obtained when
$i_{\rm s}=1$.

Figure 4 (b) shows $\Sigma_{100}^{\gamma}$ (squares and diamonds) and
$\Lambda_{100}^{\gamma}$ (triangles and inverted triangles) in the same 51st
period. Open and filled symbols represent the data for the metastable and the
stable phase crystals, respectively. As the growth rate is the difference
between the incorporation $\Sigma_{100}^{\gamma}$ and the dissolution
$\Lambda_{100}^{\gamma}$, the stable and the metastable crystals grow at low
temperature ($500\leq t<505$) and dissolve at high temperature ($505\leq
t<510$). In the low temperature period, there is no significant difference in
the growth rates. In the high temperature period, crystals dissolve and the
excess parameter increases. Fig. 4(b) shows that the metastable crystals are
more difficult to dissolve than the stable crystals because the incorporation
contribution of the metastable crystal $\Sigma^{\rm m}_{100}$ is more than
that of the stable crystal $\Sigma^{\rm s}_{100}$ while the dissolution
contributions $\Lambda_{100}^{\gamma}$ are not so different. This is because
the term $\Lambda^{\gamma}_{100}$ related to dissolution is independent of the
cluster size distribution since we have assumed that the solubility of the
metastable phase is not much different from that of the stable phase. The term
$\Sigma^{\gamma}_{100}$ related to growth depends on the number of monomers
and clusters. At the beginning of the high temperature period, the amount of
the metastable clusters produced by the dissolution is larger than that of the
stable clusters because of the abundant metastable crystals. As a result, the
excess parameter $\phi$ increases at high temperature excessively compensating
the loss at low temperature.

Figure 5:  Flow diagram from various initial conditions ($M^{\rm m}+M^{\rm
s}=1$) with $\phi(0)=0.2,0.28,0.3,0.4$ for the solubility ratio $n_{1,\rm
s}^{\rm eq}/n_{1,\rm m}^{\rm eq}=0.9$. The two stable fixed points are located
at $(0.69,0)$ (open red circle) and $(0,0.73)$ (filled blue square). The
unstable fixed point and the critical value of the parameter are located at
$(0.39,0.17)$ (magenta cross) and $(0.645,0.355)$ (filled black triangle),
respectively.  Figure 6:  The phase conversion diagram shows the final phase
for the given solubility ratio and the initial excess parameter $\phi(0)$. In
the red area and the blue area, the final phase is metastable and stable,
respectively. Black triangles, open circles, blue squares and yellow crosses
represent $\phi_{\ast}$, $\phi^{\rm m}_{\rm s}$, $\phi_{\rm s}^{\rm s}$ and
$\phi_{\rm u}$, respectively. Lines are guides for the eyes.

We search the parameter range where the energetically unfavorable metastable
phase can win the stable phase in TC. Figure 5 shows the flow diagram in a
mass space at the beginning of low temperature period for the solubility ratio
$0.9$. In our initial condition, all molecules are in crystals: $M^{\rm
m}+M^{\rm s}=1$. When the TC starts, both types of crystals dissolve and the
masses of crystals decrease. After a rapid initial relaxation, passing by
either side of the unstable fixed point, the system reaches one of the stable
fixed points. From the two flow lines of $\phi(0)=0.28$ and $0.30$, the
critical value of the parameter $\phi_{\ast}=0.29$ with respect to the initial
condition is found. If $\phi(0)>\phi_{\ast}$, all stable crystals disappear
and only metastable crystals remain. The final state is characterized by the
excess parameter $\phi(\infty)\equiv\phi^{\rm m}_{\rm s}=1$. If
$\phi(0)<\phi_{\ast}$, the metastable crystals disappear and the stable
crystals remain, that is, $\phi(\infty)\equiv\phi_{\rm s}^{\rm s}=-1$. The
crystal mass at the stable fixed point of the stable phase is larger than that
of the metastable phase owing to the solubility. The two flow lines with
$\phi(0)=0.28$ and $0.30$, are parallel during the initial relaxation, but
their directions become opposite at the unstable fixed point $(0.39,0.17)$.
The unstable fixed point in the mass space is characterized by the excess
parameter $\phi_{\rm u}=0.39$.

From the flow diagram for various solubility ratios, the phase conversion
diagram is obtained as shown in Fig. 6. Black triangles show the critical
value $\phi_{\ast}$ in the initial condition. The critical value $\phi_{\ast}$
increases linearly on decreasing the solubility ratio. However, the critical
value $\phi_{\ast}$ vanishes around the solubility ratio $n_{\rm 1,s}^{\rm
eq}/n_{\rm 1,m}^{\rm eq}\simeq 0.74$ because the stable fixed point $\phi^{\rm
m}_{\rm s}$ and the unstable fixed point $\phi_{\rm u}$ disappear. The stable
crystals prevail in the final state starting from any initial condition if the
solubility ratio is less than the value: $n_{\rm 1,s}^{\rm eq}/n_{\rm
1,m}^{\rm eq}\leq 0.74$, because the stable fixed point related to the stable
phase crystal $\phi_{\rm s}^{\rm s}$ is always present. This behavior of
$\phi$ is the subcritical pitchfork bifurcation.

##  IV Summary

We investigated the conversion from the stable phase crystals to the
metastable phase crystals by simple temperature cycling with the use of the
generalized BD model. During TC, crystals grow at low temperature and dissolve
at high temperature. When crystals dissolve, the majority metastable crystals
yield a larger amount of clusters, which prevent the dissolution of the
majority crystals through the large incorporation term $\Sigma_{i}^{\gamma}$.
At the end of one cycle, the excess parameter is amplified. From the flow
diagram in the mass space for various solubility ratios, the phase conversion
diagram is constructed. If the solubility ratio $n_{\rm 1,s}^{\rm eq}/n_{\rm
1,m}^{\rm eq}$ decreases, the initial relative amount of the metastable
crystals required for the metastable phase to become dominant is large. The
critical initial excess parameter $\phi_{\ast}$ increases linearly and makes a
jump to $\phi_{\ast}=1$ because the stable fixed point related to the
metastable phase crystals vanishes. The behavior of $\phi$ corresponds to the
subcritical pitch fork bifurcation.

To the best of our knowledge, there is no experimental result that directly
matches our scheme. We believe that the experimental verification would be
possible if the difference in solubility between stable and metastable
crystals is small. Although the VR+TC experiment of aspartic acid crystals has
succeeded in converting stable racemic crystals into metastable chiral
crystals[31], the experimental system is different from our model due to the
molecular compositions of the stable and the metastable phase crystals.
Theoretical study of the phase conversion of aspartic acid crystals is now
underway.

###### Acknowledgements.

This work is supported by JSPS KAKENHI Grant Number JP18K03500 and JP21K03379.

## References

  * Viedma [2005] C. Viedma, Phys. Rev. Lett. 94, 065504 (2005).
  * Noorduin  _et al._ [2008] W. L. Noorduin, T. Izumi, A. Millemaggi, M. Leeman, H. Meekes, W. J. Van Enckevort, R. M. Kellogg, B. Kaptein, E. Vlieg, and D. G. Blackmond, J. Am. Chem. Soc. 130, 1158 (2008).
  * Sögütoglu  _et al._ [2015] L.-C. Sögütoglu, R. R. Steendam, H. Meekes, E. Vlieg, and F. P. Rutjes, Chem. Soc. Rev. 44, 6723 (2015).
  * El-Hachemi  _et al._ [2009] Z. El-Hachemi, J. Crusats, J. M. Ribo, and S. Veintemillas-Verdaguer, Cryst. Growth. Des. 9, 4802 (2009).
  * Viedma and Cintas [2011] C. Viedma and P. Cintas, Chem. Commun. 47, 12786 (2011).
  * Suwannasang  _et al._ [2013] K. Suwannasang, A. Flood, C. Rougeot, and G. Coquerel, Cryst. Growth. Des. 13, 3498 (2013).
  * Suwannasang  _et al._ [2016] K. Suwannasang, A. E. Flood, and G. Coquerel, Cryst. Growth Des. 16, 6461 (2016).
  * Saito and Hyuga [2004] Y. Saito and H. Hyuga, J. Phys. Soc. Jpn. 73, 33 (2004).
  * Uwaha [2004] M. Uwaha, J. Phys. Soc. Jpn. 73, 2601 (2004).
  * Saito and Hyuga [2008] Y. Saito and H. Hyuga, J. Phys. Soc. Jpn. 77, 113001 (2008).
  * Saito and Hyuga [2010] Y. Saito and H. Hyuga, J. Phys. Soc. Jpn. 79, 083002 (2010).
  * Frank [1953] F. C. Frank, Biochem. Biophys. Acta 11, 459 (1953).
  * Cartwright  _et al._ [2007] J. H. E. Cartwright, O. Piro, and I. Tuval, Phys. Rev. Lett. 98, 165501 (2007).
  * Uwaha [2008] M. Uwaha, J. Phys. Soc. Jpn. 77, 083802 (2008).
  * McBride and Tully [2008] J. M. McBride and J. C. Tully, Nature 452, 161 (2008).
  * Steendam  _et al._ [2015] R. R. Steendam, J. Dickhout, W. J. Van Enckevort, H. Meekes, J. Raap, F. P. Rutjes, and E. Vlieg, Cryst. Growth Des. 15, 1975 (2015).
  * Spix  _et al._ [2016] L. Spix, A. H. Engwerda, H. Meekes, W. J. van Enckevort, and E. Vlieg, Cryst. Growth Des. 16, 4752 (2016).
  * Katsuno and Uwaha [2017] H. Katsuno and M. Uwaha, Phys. Rev. E 95, 062804 (2017).
  * Uwaha and Katsuno [2009] M. Uwaha and H. Katsuno, J. Phys. Soc. Jpn. 78, 023601 (2009).
  * Wattis [2011] J. A. Wattis, Orig. Life Evol. Biosph. 41 (2011).
  * Iggland and Mazzotti [2011] M. Iggland and M. Mazzotti, Cryst. Growth Des. 11, 4611 (2011).
  * Blanco  _et al._ [2013] C. Blanco, J. Crusats, Z. El-Hachemi, A. Moyano, S. Veintemillas-Verdaguer, D. Hochberg, and J. M. Ribó, ChemPhysChem 14, 3982 (2013).
  * Blanco  _et al._ [2015] C. Blanco, J. M. Ribó, and D. Hochberg, Physical Review E 91, 022801 (2015).
  * Iggland  _et al._ [2014] M. Iggland, R. Müller, and M. Mazzotti, Cryst. Growth Des. 14, 2488 (2014).
  * Katsuno and Uwaha [2009] H. Katsuno and M. Uwaha, J. Cryst. Growth 311, 4265 (2009).
  * Ricci  _et al._ [2013] F. Ricci, F. H. Stillinger, and P. G. Debenedetti, J. Chem. Phys. 139, 12629 (2013).
  * Katsuno and Uwaha [2016] H. Katsuno and M. Uwaha, Phys. Rev. E 93, 013002 (2016).
  * Bodák  _et al._ [2018] B. Bodák, G. M. Maggioni, and M. Mazzotti, Cryst. Growth Des. 18, 7122 (2018).
  * Bodák  _et al._ [2019] B. Bodák, G. M. Maggioni, and M. Mazzotti, Cryst. Growth Des. 19, 6552 (2019).
  * Uwaha and Katsuno [2022] M. Uwaha and H. Katsuno, J. Cryst. Growth 598, 126873 (2022).
  * Viedma and Ortiz [2021] C. Viedma and J. E. Ortiz, Isr. J. Chem. 61, 758 (2021).
  * Katsuno and Uwaha [2021] H. Katsuno and M. Uwaha, J. Phys. Soc. Jpn. 90, 044001 (2021).
  * Becker and Döring [1935] R. Becker and W. Döring, Ann. Phys.(Leipzig) 24, 719 (1935).
  * Markov [2016] I. V. Markov,  _Crystal growth for beginners: fundamentals of nucleation, crystal growth and epitaxy_ (World scientific, 2016).
  * Uwaha [2010] M. Uwaha, in  _AIP Conference Proceedings_ , Vol. 1270 (American Institute of Physics, 2010) pp. 155–172.
  * Note [1] In the original BD model, the time change of the number density is described. Here, we call the number density the number for simplicity.
  * Belletti  _et al._ [2022] G. Belletti, J. Schuurman, H. Stinesen, H. Meekes, F. P. Rutjes, and E. Vlieg, Cryst. Growth Des. 22, 1874 (2022).


title = {Algebraic surfaces}, series = {Universitext}, note = {Translated from
the 1981 Romanian original by Vladimir Ma\c{s}ek and revised by the author},
publisher = {Springer-Verlag, New York}, date = {2001}, pages = {xii+258},
isbn = {0-387-98668-5}}

  * [5]
BirkarCaucherZhangDe-QiEffectivity of iitaka fibrations and pluricanonical
systems of polarized pairsPubl. Math. Inst. Hautes Études
Sci.1232016283–331ISSN 0073-8301@article{BZ16, author = {Birkar, Caucher},
author = {Zhang, De-Qi}, title = {Effectivity of Iitaka fibrations and
pluricanonical systems of polarized pairs}, journal = {Publ. Math. Inst.
Hautes \'{E}tudes Sci.}, volume = {123}, date = {2016}, pages = {283–331},
issn = {0073-8301}}

  * [7]
CaoJunyanAlbanese maps of projective manifolds with nef anticanonical
bundlesEnglish, with English and French summariesAnn. Sci. Éc. Norm. Supér.
(4)52201951137–1154ISSN 0012-9593@article{Cao19, author = {Cao, Junyan}, title
= {Albanese maps of projective manifolds with nef anticanonical bundles},
language = {English, with English and French summaries}, journal = {Ann. Sci.
\'{E}c. Norm. Sup\'{e}r. (4)}, volume = {52}, date = {2019}, number = {5},
pages = {1137–1154}, issn = {0012-9593}}

  * [9]
CampanaFrédéricCaoJunyanMatsumuraShin-ichiProjective klt pairs with nef anti-
canonical divisorAlgebr. Geom.820214430–464ISSN 2313-1691@article{CCM21,
author = {Campana, Fr\'{e}d\'{e}ric}, author = {Cao, Junyan}, author =
{Matsumura, Shin-ichi}, title = {Projective klt pairs with nef anti-canonical
divisor}, journal = {Algebr. Geom.}, volume = {8}, date = {2021}, number =
{4}, pages = {430–464}, issn = {2313-1691}}

  * [11]
CasciniPaoloTanakaHiromuXuChenyangOn base point freeness in positive
characteristicEnglish, with English and French summariesAnn. Sci. Éc. Norm.
Supér. (4)48201551239–1272ISSN 0012-9593@article{CTX15, author = {Cascini,
Paolo}, author = {Tanaka, Hiromu}, author = {Xu, Chenyang}, title = {On base
point freeness in positive characteristic}, language = {English, with English
and French summaries}, journal = {Ann. Sci. \'{E}c. Norm. Sup\'{e}r. (4)},
volume = {48}, date = {2015}, number = {5}, pages = {1239–1272}, issn =
{0012-9593}}

  * [13]
ChenYifeiZhangLeiThe subadditivity of the kodaira dimension for fibrations of
relative dimension one in positive characteristicsMath. Res.
Lett.2220153675–696ISSN 1073-2780@article{CZ15, author = {Chen, Yifei}, author
= {Zhang, Lei}, title = {The subadditivity of the Kodaira dimension for
fibrations of relative dimension one in positive characteristics}, journal =
{Math. Res. Lett.}, volume = {22}, date = {2015}, number = {3}, pages =
{675–696}, issn = {1073-2780}}

  * [15]
EjiriShoWhen is the albanese morphism an algebraic fiber space in positive
characteristic?Manuscripta Math.16020191-2239–264ISSN
0025-2611@article{Ejiri-2019, author = {Ejiri, Sho}, title = {When is the
Albanese morphism an algebraic fiber space in positive characteristic?},
journal = {Manuscripta Math.}, volume = {160}, date = {2019}, number = {1-2},
pages = {239–264}, issn = {0025-2611}}

  * [17]
EkedahlTorstenFoliations and inseparable morphismstitle={Algebraic geometry,
Bowdoin, 1985}, address={Brunswick, Maine}, date={1985}, series={Proc. Sympos.
Pure Math.}, volume={46, Part 2}, publisher={Amer. Math. Soc., Providence,
RI}, ISBN 0-8218-1480-X1987139–149@article{Ekedahl-1987, author = {Ekedahl,
Torsten}, title = {Foliations and inseparable morphisms}, conference =
{title={Algebraic geometry, Bowdoin, 1985}, address={Brunswick, Maine},
date={1985}, }, book = {series={Proc. Sympos. Pure Math.}, volume={46, Part
2}, publisher={Amer. Math. Soc., Providence, RI}, }, isbn = {0-8218-1480-X},
date = {1987}, pages = {139–149}}

  * [19]
EjiriShoPatakfalviZsoltThe demailly–peternell–schneider conjecture is true in
positive characteristicarXiv preprint arXiv:2305.021572023@article{EP23,
author = {Ejiri, Sho}, author = {Patakfalvi, Zsolt}, title = {The
Demailly–Peternell–Schneider conjecture is true in positive characteristic},
journal = {arXiv preprint arXiv:2305.02157}, year = {2023}}

  * [21]
FlorisEnricaLazićVladimirA travel guide to the canonical bundle
formulatitle={Birational geometry and moduli spaces}, series={Springer INdAM
Ser.}, volume={39}, publisher={Springer, Cham}, 202037–55@article{FL20, author
= {Floris, Enrica}, author = {Lazi\'{c}, Vladimir}, title = {A travel guide to
the canonical bundle formula}, conference = {title={Birational geometry and
moduli spaces}, }, book = {series={Springer INdAM Ser.}, volume={39},
publisher={Springer, Cham}, }, date = {2020}, pages = {37–55}}

  * [23]
FujinoOsamuMoriShigefumiA canonical bundle formulaJ. Differential
Geom.5620001167–188ISSN 0022-040X@article{FM00, author = {Fujino, Osamu},
author = {Mori, Shigefumi}, title = {A canonical bundle formula}, journal =
{J. Differential Geom.}, volume = {56}, date = {2000}, number = {1}, pages =
{167–188}, issn = {0022-040X}}

  * [25]
FultonWilliamIntersection theoryErgebnisse der Mathematik und ihrer
Grenzgebiete. 3. Folge. A Series of Modern Surveys in Mathematics [Results in
Mathematics and Related Areas. 3rd Series. A Series of Modern Surveys in
Mathematics]22Springer-Verlag, Berlin1998xiv+470ISBN 3-540-62046-XISBN
0-387-98549-2@book{FultonIT, author = {Fulton, William}, title = {Intersection
theory}, series = {Ergebnisse der Mathematik und ihrer Grenzgebiete. 3. Folge.
A Series of Modern Surveys in Mathematics [Results in Mathematics and Related
Areas. 3rd Series. A Series of Modern Surveys in Mathematics]}, volume = {2},
edition = {2}, publisher = {Springer-Verlag, Berlin}, date = {1998}, pages =
{xiv+470}, isbn = {3-540-62046-X}, isbn = {0-387-98549-2}}

  * [27]
GrothendieckA.Éléments de géométrie algébrique. iv. Étude locale des schémas
et des morphismes de schémas. iiiInst. Hautes Études Sci. Publ.
Math.281966255ISSN 0073-8301@article{GorIV-III, author = {Grothendieck, A.},
title = {\'{E}l\'{e}ments de g\'{e}om\'{e}trie alg\'{e}brique. IV. \'{E}tude
locale des sch\'{e}mas et des morphismes de sch\'{e}mas. III}, journal =
{Inst. Hautes \'{E}tudes Sci. Publ. Math.}, number = {28}, date = {1966},
pages = {255}, issn = {0073-8301}}

  * [29]
HaconChristopher D.PatakfalviZsoltZhangLeiBirational characterization of
abelian varieties and ordinary abelian varieties in characteristic $p>0$Duke
Math. J.168201991723–1736ISSN 0012-7094@article{HPZ19, author = {Hacon,
Christopher D.}, author = {Patakfalvi, Zsolt}, author = {Zhang, Lei}, title =
{Birational characterization of Abelian varieties and ordinary Abelian
varieties in characteristic $p>0$}, journal = {Duke Math. J.}, volume = {168},
date = {2019}, number = {9}, pages = {1723–1736}, issn = {0012-7094}}

  * [31]
IgusaJun-ichiA fundamental inequality in the theory of picard varietiesProc.
Nat. Acad. Sci. U.S.A.411955317–320ISSN 0027-8424@article{Igusa55-ineq, author
= {Igusa, Jun-ichi}, title = {A fundamental inequality in the theory of Picard
varieties}, journal = {Proc. Nat. Acad. Sci. U.S.A.}, volume = {41}, date =
{1955}, pages = {317–320}, issn = {0027-8424}}

  * [33]
IitakaShigeruAlgebraic geometryNorth-Holland Mathematical Library24An
introduction to birational geometry of algebraic varietiesSpringer-Verlag, New
York-Berlin1982x+357ISBN 0-387-90546-4@book{Iit82, author = {Iitaka, Shigeru},
title = {Algebraic geometry}, series = {North-Holland Mathematical Library},
volume = {24}, note = {An introduction to birational geometry of algebraic
varieties}, publisher = {Springer-Verlag, New York-Berlin}, date = {1982},
pages = {x+357}, isbn = {0-387-90546-4}}

  * [35]
JiLenaWaldronJoeStructure of geometrically non-reduced varietiesTrans. Amer.
Math. Soc.3742021128333–8363ISSN 0002-9947@article{JW21, author = {Ji, Lena},
author = {Waldron, Joe}, title = {Structure of geometrically non-reduced
varieties}, journal = {Trans. Amer. Math. Soc.}, volume = {374}, date =
{2021}, number = {12}, pages = {8333–8363}, issn = {0002-9947}}

  * [37]
KawamataYujiroSubadjunction of log canonical divisors. iiAmer. J.
Math.12019985893–899ISSN 0002-9327@article{Kaw98, author = {Kawamata, Yujiro},
title = {Subadjunction of log canonical divisors. II}, journal = {Amer. J.
Math.}, volume = {120}, date = {1998}, number = {5}, pages = {893–899}, issn =
{0002-9327}}

  * [39]
KeelSeánBasepoint freeness for nef and big line bundles in positive
characteristicAnn. of Math. (2)14919991253–286ISSN 0003-486X@article{Keel99,
author = {Keel, Se\'{a}n}, title = {Basepoint freeness for nef and big line
bundles in positive characteristic}, journal = {Ann. of Math. (2)}, volume =
{149}, date = {1999}, number = {1}, pages = {253–286}, issn = {0003-486X}}

  * [41]
LazarsfeldRobertPositivity in algebraic geometry. iErgebnisse der Mathematik
und ihrer Grenzgebiete. 3. Folge. A Series of Modern Surveys in Mathematics
[Results in Mathematics and Related Areas. 3rd Series. A Series of Modern
Surveys in Mathematics]48Classical setting: line bundles and linear
seriesSpringer-Verlag, Berlin2004xviii+387ISBN 3-540-22533-1@book{Lazarsfeld,
author = {Lazarsfeld, Robert}, title = {Positivity in algebraic geometry. I},
series = {Ergebnisse der Mathematik und ihrer Grenzgebiete. 3. Folge. A Series
of Modern Surveys in Mathematics [Results in Mathematics and Related Areas.
3rd Series. A Series of Modern Surveys in Mathematics]}, volume = {48}, note =
{Classical setting: line bundles and linear series}, publisher = {Springer-
Verlag, Berlin}, date = {2004}, pages = {xviii+387}, isbn = {3-540-22533-1}}

  * [43]
LiuQingAlgebraic geometry and arithmetic curvesOxford Graduate Texts in
Mathematics6Translated from the French by Reinie Erné; Oxford Science
PublicationsOxford University Press, Oxford2002xvi+576ISBN
0-19-850284-2@book{LiuAGAC, author = {Liu, Qing}, title = {Algebraic geometry
and arithmetic curves}, series = {Oxford Graduate Texts in Mathematics},
volume = {6}, note = {Translated from the French by Reinie Ern\'{e}; Oxford
Science Publications}, publisher = {Oxford University Press, Oxford}, date =
{2002}, pages = {xvi+576}, isbn = {0-19-850284-2}}

  * [45]
PatakfalviZsoltWaldronJoeSingularities of general fibers and the lmmpAmer. J.
Math.14420222505–540ISSN 0002-9327@article{PW22, author = {Patakfalvi, Zsolt},
author = {Waldron, Joe}, title = {Singularities of general fibers and the
LMMP}, journal = {Amer. J. Math.}, volume = {144}, date = {2022}, number =
{2}, pages = {505–540}, issn = {0002-9327}}

  * [47]
QueenClifford S.Non-conservative function fields of genus one. iArch. Math.
(Basel)221971612–623ISSN 0003-889X@article{Queen-Non-conservative-I, author =
{Queen, Clifford S.}, title = {Non-conservative function fields of genus one.
I}, journal = {Arch. Math. (Basel)}, volume = {22}, date = {1971}, pages =
{612–623}, issn = {0003-889X}}

  * [49]
SchröerStefanOn fibrations whose geometric fibers are nonreducedNagoya Math.
J.200201035–57ISSN 0027-7630@article{Sch10, author = {Schr\"{o}er, Stefan},
title = {On fibrations whose geometric fibers are nonreduced}, journal =
{Nagoya Math. J.}, volume = {200}, date = {2010}, pages = {35–57}, issn =
{0027-7630}}

  * [51]
TanakaHiromuInvariants of algebraic varieties over imperfect fieldsTohoku
Math. J. (2)7320214471–538ISSN 0040-8735@article{Tanaka21I, author = {Tanaka,
Hiromu}, title = {Invariants of algebraic varieties over imperfect fields},
journal = {Tohoku Math. J. (2)}, volume = {73}, date = {2021}, number = {4},
pages = {471–538}, issn = {0040-8735}}

  * [53]
WangYuanOn the characterization of abelian varieties for log pairs in zero and
positive characteristicEnglish, with English and French summariesAnn. Inst.
Fourier (Grenoble)72202262515–2539ISSN 0373-0956@article{Wang-2022, author =
{Wang, Yuan}, title = {On the characterization of abelian varieties for log
pairs in zero and positive characteristic}, language = {English, with English
and French summaries}, journal = {Ann. Inst. Fourier (Grenoble)}, volume =
{72}, date = {2022}, number = {6}, pages = {2515–2539}, issn = {0373-0956}}

  * [55]
WitaszekJakubOn the canonical bundle formula and log abundance in positive
characteristicMath. Ann.38120213-41309–1344ISSN 0025-5831@article{Wit21,
author = {Witaszek, Jakub}, title = {On the canonical bundle formula and log
abundance in positive characteristic}, journal = {Math. Ann.}, volume = {381},
date = {2021}, number = {3-4}, pages = {1309–1344}, issn = {0025-5831}}

  * [57]
ZhangLeiAbundance for non-uniruled 3-folds with non-trivial albanese maps in
positive characteristicsJ. Lond. Math. Soc. (2)9920192332–348ISSN
0024-6107@article{Zha19, author = {Zhang, Lei}, title = {Abundance for non-
uniruled 3-folds with non-trivial Albanese maps in positive characteristics},
journal = {J. Lond. Math. Soc. (2)}, volume = {99}, date = {2019}, number =
{2}, pages = {332–348}, issn = {0024-6107}}

  * [59]