Bayesian Modelling of Pattern Formation from One Snapshot of Pattern

Natsuhiko Yoshinaga [email protected] WPI - Advanced Institute for Materials Research, Tohoku University, Sendai 980-8577, Japan MathAM-OIL, AIST, Sendai 980-8577, Japan Satoru Tokuda Research Institute for Information Technology, Kyushu University, Kasuga 816-8580, Japan MathAM-OIL, AIST, Sendai 980-8577, Japan

Abstract

Partial differential equations (PDE) have been widely used to reproduce patterns in nature and to give insight into the mechanism underlying pattern formation. Although many PDE models have been proposed, they rely on the pre-request knowledge of physical laws and symmetries, and developing a model to reproduce a given desired pattern remains difficult. We propose a novel method, referred to as Bayesian modelling of PDE (BM-PDE), to estimate the best dynamical PDE for one snapshot of a target pattern under the stationary state without ground truth. We apply BM-PDE to nontrivial patterns, such as quasi-crystals (QCs), a double gyroid and Frank Kasper structures. By using the estimated parameters for the approximant of QCs, we successfully generate, for the first time, three-dimensional dodecagonal QCs from a PDE model. Our method works for noisy patterns and the pattern synthesised without the ground truth parameters, which are required for the application toward experimental data.

The design of structures of materials is one of the most important issues in various fields of physical science, as their structures are related to their physical properties. The structures are often characterised by periodic or quasi-periodic order. These ordered structures, which we call pattern, are ubiquitous in nature ranging from fluid convectionswift:1977 to the microphase separation of block copolymersbates:1990 ; Harrison:2000 and atomic and molecular crystalsElder:2002 ; Provatas:2011 . Surprisingly, the same pattern appears in different systems with completely different length scalescross:1993 . Complex patterns such as quasi-crystal (QC), double gyroid (DG) and Frank Kasper (FK) phases appear not only in metallic alloyFrank:1959 ; Shechtman:1984 but also in soft materials such as block copolymersbates:1990 ; Bates:2019 , biomaterialsScherer:2013 , surfactantsYue:2016 , liquid crystalsZeng:2004 and colloidal assembliesHynninen:2007 .

Understanding a generic mechanism of symmetry selection is an important step to understand structural pattern formation. A continuum approach, using nonlinear partial differential equations (PDEs), is useful for this purposecross:1993 ; Provatas:2011 . For example, it was shown that at least two length scales are necessary for the formation of QCs by using a phenomenological modelLifshitz:1997 . These studies have clarified generic pictures on a specific pattern. Still, when encountering a new pattern, we need to find a governing equation and parameters. They require a sophisticated guess and trial-and-error. Therefore, it is necessary to develop a systematic method to estimate the best model for the target pattern at hand.

Recent developments of imaging techniques give us various structural information, and thus it is desired to understand pattern formation and the design principle of a desired structure Kalinin:2015 . There are two challenging issues to consider the inverse structural design applied to the real-world data. First, the structural data is often stationary because a high-resolution time-dependent image is hard to acquire. Experimentally obtained structure is stable, or at least meta-stable. We want to estimate the best dynamical model to reproduce the target pattern as a (meta-)stable structure from stationary data. Second, there is no ground truth of the target structure; the true model to reproduce the experimental data is not available. Still, finding a surrogate model or phenomenological model is helpful to clarify underlying mechanisms of the structural formationKennedy:2001 . This issue is called model inadequacy and is one of the biggest challenges in model estimationSmith:2013 .

Automatic discovery of a model equation is a recent key topic in data-driven scienceDaniels:2015 ; Brunton:2016 ; Brunton::2019 , and several methods have been proposed for PDEs of time-series data Voss:2004 ; MUeLLER:2004 ; Rudy:2017 ; Zhao:2020 . These approaches are designed to handle time-series data and cannot be applied to estimation for a stationary pattern, because it does not have information about trajectories from their initial conditions. In addition, the previous studies focus on the estimation for the data with ground truth. The data is generated typically from numerical results of a known model with known parameters, and then, those parameters are estimated from the data. Therefore, the two issues, the estimation for stationary data and data without ground truth, remain elusive even in a technical aspect of model discovery.

To overcome these issues, in this study, an inverse problem is formulated to reproduce a given snapshot of a pattern, after which we propose a method, referred to as Bayesian modelling of partial differential equations (BM-PDE), to identify the best dynamical PDE model and its parameters. We apply our method to the problem without ground truth. We prepare the target pattern according to crystallographic symmetry, which is independent of any candidate models, and then perform an estimation of the best model. We demonstrate the BM-PDE for complex patterns such as QC, DG, and FK A15 patterns. We also demonstrate that from the estimated parameters, three-dimensional dodecagonal QCs can be generated. The success shows a potential application of BM-PDE to understand a mechanism of structural formation of novel materials.

Basic Formula

We consider a pattern (or crystalline structures) expressed by the scalar density field $\psi({\bf x})$ . An example of a two-dimensional dodecagonal QC (DDQC) is shown in Fig. 1(a,b). Higher density spots may be considered as a position of particles. We estimate a dynamical model to reproduce a target pattern $\psi^{*}({\bf x})$ as a stable pattern at the steady state $\psi_{s}({\bf x})$ of a nonlinear partial differential equation $\partial_{t}\psi({\bf x})=f_{\mu}[\psi({\bf x})]$ (Fig. 1(c-e)). If the PDE and its parameters $\mu$ are ground truth for the target pattern, $\partial_{t}\psi^{*}({\bf x})=f_{\mu}[\psi^{*}({\bf x})]=0$ is satisfied. Our target pattern is one snapshot and has information only about the stationary state (Fig. 1(d-g)). Its transient pattern from the initial state to the stationary state is not available. We assume the stationary state is stable in the sense that the pattern is generated from a broad range of the initial conditions (Fig. 1(f,g)). This assumption is natural when the target pattern is obtained from an experimental result; the pattern should be reproducible.

At a first look, the estimation for the stationary data is impossible. When $f_{\mu}[\psi^{*}({\bf x})]=0$ is a true model, we may have a series of equally true models such as $f_{\mu}[\psi^{*}({\bf x})]^{2}=0$ , and $\left(f_{\mu}[\psi^{*}({\bf x})]+1\right)f_{\mu}[\psi^{*}({\bf x})]=0$ . Therefore, the estimation is not unique. Nevertheless, as we will see later, the estimation of a model that reproduces the target pattern as a stable structure plays a role as regularisation (see also Supplementary Sec. LABEL:SM.sec.relatedworks seeSM for the comparison with other approaches) . We should note that our problem is not the parameter estimation for $f_{\mu}[\psi^{*}({\bf x})]=0$ in which the stability of the stationary state is not guaranteed.

To see the difficulty of estimating stationary data, it is instructive to consider the state-space model, widely used in data assimilationEvensen:2009 ; Law:2015 . The cost function $E$ consists of measurement (observation) and model errors, and is expressed as

\displaystyle E[\mu,\psi(\bf x)]

\displaystyle=\frac{1}{2}\|\psi^{*}({\bf x})-\psi({\bf x})\|+\frac{1}{2}\|\partial_{t}\psi({\bf x})-f_{\mu}(\psi({\bf x}))\|.

(1)

In many cases, the norm $\|\cdot\|$ is chosen to be the square norm. If the observation contains an error, we have to estimate both the parameters $\mu$ and the state $\psi({\bf x})$ . When the model represents a deterministic system, the state is described by its initial condition $\psi_{0}({\bf x})$ , and accordingly the cost function becomes $E[\mu,\psi_{0}({\bf x})]$ . When the observation does not contain noise, the first term in equation (1) vanishes, and the problem falls into a simple regression (see also Supplementary Sec. LABEL:SM.sec.relatedworks seeSM ). In the conventional data assimilation, both $\mu$ and $\psi_{0}$ are estimated by minimising $E[\mu,\psi_{0}({\bf x})]$ Law:2015 . However, in the problem of pattern formation, the specific initial condition to produce the pattern is not of our interest. Moreover, there are many initial conditions that asymptotically reproduce the same pattern (Fig. 1(f)). Therefore, the estimation under a given stationary target pattern cannot be unique. In the case of time-series data, the estimation is possible because each trajectory is from a different initial condition. In our method, we sample the stationary pattern by solving the model under each parameter and each realisation from the random initial condition. By marginalising the initial conditions, we may obtain a unique estimation.

The basic structure of our estimation is schematically shown in Fig. 1(h). Parameters are estimated from the posterior distribution under the target pattern, whereas the best model is estimated from the log marginal likelihood. For a given model $m_{i}$ and parameters $\mu$ in a PDE, the stationary pattern is uniquely determined under each initial condition, $\psi_{0}({\bf x})$ (Fig. 1(c)). We treat an initial pattern as a latent variable that is marginalised using a random variable for the initial condition. This is because the obtained pattern may be translationally shifted or rotated by changing an initial pattern (see Fig 1(c)). We quantify the similarity between two patterns $\psi_{1}$ and $\psi_{2}$ by the distance between them defined as

\displaystyle E[\psi_{1},\psi_{2}]

\displaystyle=\left|{\bf\Psi}[\psi_{1}]-{\bf\Psi}[\psi_{2}]\right|^{2}.

(2)

Our target pattern is ordered and has many invariants; the pattern must be identified under change by translation and rotation, and also the pattern does not change by the action of the symmetry group that the pattern has (see Fig. 1(a,b)). Here we introduce order parameter ${\bf\Psi}[\psi({\bf x})]=\left(\Psi_{1}[\psi({\bf x})],\ldots\Psi_{l_{0}}[\psi({\bf x})]\right)$ , which maps the pattern onto the feature space and eliminates the redundant information of the ordered pattern $\psi({\bf x})$ due to symmetry (see Fig. 1(c) and Methods). The distance defined by equation (2) identify the patterns up to symmetry transformation thanks to the order parameter.

Bayesian modelling

Our goal is to find the most probable model $\hat{m}$ described by a PDE and its parameters $\hat{\mu}$ for a given target pattern $\psi^{*}({\bf x})$ . We also want to quantify the uncertainty of the estimation. To achieve this, we use the cost function $E[\psi_{s},\psi^{*}]$ , or called the energy, expressed by the order parameter ${\bf\Psi}$ , and compute the distance from the target pattern, ${\bf\Psi}^{*}={\bf\Psi}[\psi^{*}({\bf x})]$ , to the numerically generated stationary pattern for each model and parameter set, ${\bf\Psi}[\psi_{s}({\bf x})]$ . Our purpose is not to estimate specific initial states $\psi_{0}$ for the target pattern $\psi^{*}$ , but to estimate the best model that could generate patterns similar to $\psi^{*}$ independent of the initial state. Therefore, our best parameter set $\hat{\mu}$ is defined by the mean of the marginal probability distribution under a model $m$

\displaystyle p(\mu\mid{\bf\Psi}^{*},\beta,m)=\int p(\psi_{0})p(\mu\mid{\bf\Psi}^{*},\psi_{0},\beta,m)d\psi_{0}.

(3)

The integral over the initial conditions $\psi_{0}$ implies that the posterior distribution of the parameters is chosen so that the estimated parameters can generate the target pattern from a wide range of the initial conditions. We may avoid the parameters that can generate the target pattern only for a specific initial condition (Fig. 1(f,g)). Following Bayes’ theorem, the posterior distribution under a fixed $\psi_{0}$ is given by

\displaystyle p(\mu\mid{\bf\Psi}^{*},\psi_{0},\beta,m)

\displaystyle=\frac{p({\bf\Psi}^{*}\mid\psi_{0},\mu,\beta,m)p(\mu\mid m)}{p({\bf\Psi}^{*}\mid m,\beta)}.

(4)

The likelihood is represented by $p({\bf\Psi}^{*}\mid\mu,m,\beta)\propto e^{-\beta E[\psi,\psi^{*}]}$ , where the inverse temperature $\beta$ is associated with variance of the observation noise. This likelihood implies that the error in the measurement is given by ${\bf\Psi}^{*}={\bf\Psi}+\bm{\xi}$ with the Gaussian noise $\bm{\xi}$ with zero mean and its variance $\beta^{-1}$ . The prior distributions $p(\psi_{0})$ and $p(\mu\mid m)$ are assumed to the uniform distribution. The normalisation factor $p({\bf\Psi}^{*}\mid\beta,m)$ , or the log marginal likelihood (free energy) $F(\beta,m)\equiv-\log p({\bf\Psi}^{*}\mid\beta,m)$ , is one of the criteria of model selection and hyperparameter estimation mackay1992bayesian ; Bishop:2006 . Both $p({\bf\Psi}^{*}\mid\beta,m)$ and $F(\beta,m)$ play a role as a probability density of the models and hyperparameters (see Methods). Therefore, our best model $m$ and inverse temperature $\beta$ are both determined by maximising $p({\bf\Psi}^{*}\mid\beta,m)$ , or equivalently, minimising $F(\beta,m)$ Tokuda:2017 .

target patterns without ground truth

In this work, we consider two types of target patterns; (i) a numerical solution of the PDE model that we use, and (ii) a pattern synthesised by superposition of plane waves as equation (14). The former has ground truth, whereas the latter does not. The synthesised pattern is independent of our PDE models, and therefore, it is not necessarily a solution for the PDEs. We will demonstrate that our approach still works for the problem without ground truth. This problem is an intermediate step between the estimation of the problem with ground truth and experimental data. We do not know the true PDE for experimental data. Still, we would like to estimate the best PDE model to explain the data to understand the mechanism of the structural pattern formation. Our approach has the advantage in that we do not know a priori the model that reproduces the target pattern. In most of model identification, the problem with ground truth has been used.

To synthesise the target structure by equation (14), positions of the peaks and their amplitude in the Fourier space are required (Fig.2(a)). These quantities can be directly measured in the scattering experiments in two dimensions and three dimensions. The positions and the amplitude in the Fourier space can be reproduced from the structure factor, as shown in Fig.2. We use this information to synthesise the target. We demonstrate the estimation of PDEs for the Frank-Kasper A15 structure taken from the experimental dataJayaraman:2018 (Fig.2(b)). The detailed method to synthesise the target structure is shown in Methods.

The estimation for phase field crystal models

We demonstrate our method for a class of phase-field crystal models. Specifically, we make a family of the models by changing mainly the linear part of equation (5). Our approach is not limited to the choice of this family. In fact, we can replace equation (5) with any other families of models. A possible extension of this approach, including the estimation of nonlinear terms, is discussed in Supplementary Sec.LABEL:sec.family.model seeSM . We consider a model, called $m\in M$ , expressed by a nonlinear PDE of the form

\displaystyle\partial_{t}\psi({\bf x})

\displaystyle=\mathcal{L}_{\mu}^{(m)}\psi({\bf x})+\mathcal{N}\left[\psi({\bf x})\right]

(5)

with a set of parameters $\mu$ . The PDE is decomposed into two parts. The linear term is expressed by the linear operator, $\mathcal{L}_{\mu}^{(m)}$ , acting on $\psi({\bf x})$ . Because we are interested in patterns in bulk, not affected by boundaries, we use periodic boundary conditions.

We make a family of models, $M=\{m_{i}\}_{i=1,2,\ldots,i_{\rm max}}$ . In model $m_{i}$ the linear operator has $i$ peaks in its spectrum (see Fig. 3(f-h)). Our family of models is designed to have a physical interpretation that the system has $i$ length scales for model $m_{i}$ because peaks in the spectrum correspond to the number of length scales. Each length scale is characterised by its wavelength $q_{i}$ and the value of its spectrum at the wavelength $a_{i}$ ( see Fig. 3(f-h)). We also use the mean density $\bar{\psi}$ and system sizes in each direction as parameters. Note that $\bar{\psi}$ is identical to the parameter for the nonlinear $\psi^{2}$ term in equation (5).

two-dimensional QC with ground truth

To give better insight into the BM-PDE, we first focus on an example of a two-dimensional QC with 12-fold symmetry (DDQC) shown in Fig. 1(a). This pattern has been studied using a model with two length scalesLifshitz:1997 . The target pattern in this section is numerically produced with a set of parameters $\mu^{*}$ . The two-length-scale model is used, that is, $m^{*}=m_{2}$ .

For $m=m_{2}$ , the cost function $E[\psi_{s},\psi^{*}]$ decreases during the sampling from the posterior distribution, and the generated stationary patterns from equation (5) converge to QCs which are similar to the target pattern (Fig.3(b), see also Supplementary Fig.LABEL:SM_fig2D seeSM ). The estimated parameters well agree with the parameters that we used to generate the target pattern (see Fig.3(d) and Supplementary Table.LABEL:estimated.parameters.numerical.12fold seeSM ). The estimated length scale is $q_{1}=0.52$ with the second length scale $q_{0}=1$ . The ratio between them agrees to $q_{0}/q_{1}=2\cos(\pi/12)\approx 1.9319$ , which is the known value to generate this patternLifshitz:1997 . The BM-PDE automatically estimates this ratio starting from uniform prior distribution of the wavelength. The estimation also works for other parameters (Supplementary Table.LABEL:estimated.parameters.numerical.12fold and Supplementary Fig.LABEL:SM_fig2D seeSM ).

Using the estimated parameters, we may generate a pattern similar to the target pattern from uniform random initial density (see the inset of Fig.3(b)). To see the quality of the estimation, we measure the steady distribution of the cost function. Figure 3(b) shows that there are two distinct states: one has a higher cost function $E\gtrsim 10^{2}$ and the other has a lower cost function $E\lesssim 10^{2}$ . The latter corresponds to QC, whereas hexagonal patterns mainly dominate the former. The gap between the two states indicates that the QC patterns require a high resolution in the parameter search. The large step in the parameter space cannot find the optimal parameters because their range is narrow in the prior range of the parameters as in the inset of Fig.3(d). Therefore, the conventional gradient method with a fixed step size either cannot find the QC when the step size is large, or impractical when it is small. The BM-PDE use hierarchical sampling, such as replica exchange Monte Carlo (REMC) Geyer:1991 ; Hukushima:1996 . In this method, the knowledge of the step size corresponding to the observation noise $\beta$ is not a prerequisite for the parameter search.

The same algorithm is performed for the models with one length scale $m=m_{1}$ and three length scale $m=m_{3}$ . It is not surprising that the cost function $E$ of the one-length-scale model is much higher than that of $m_{2}$ because it cannot reproduce a QC pattern (Fig. 3(a)). In fact, the model $m=m_{1}$ with the estimated parameters produces hexagonal patterns rather than QCs. The three-length-scale model does reproduce QCs, which are comparable to the target pattern (Fig. 3(c)). However, our Bayesian estimation selects the two-length-scale model. Figure 3(e) demonstrates that the minimum free energy, that is, the marginal likelihood, is lower for $m=m_{2}$ . Therefore, we estimate that this QC pattern is described by the interaction of two modes with different length scale.

various target patterns without ground truth

The BM-PDE is not restricted to the estimation of the parameters that are used to generate the target pattern. Using a two-dimensional DDQC, we demonstrate that the BM-PDE successfully estimates the best model and approximated parameters for the target pattern without ground truth. The DDQC is synthesised by the superposition of 12 plane waves in equation (14) (see Methods). The pattern is similar to numerically produced QC used in the previous section (see Fig. 1(b) and Figure 4(c)), but in the current case, there are no ground truth parameters and a true model. The target pattern can only approximately be accessed by one of the models in equation (5). In contrast with the numerically produced pattern, estimated parameters do not reproduce exactly the same pattern as the target pattern, and therefore the cost function is relatively high (Supplementary Fig. LABEL:SM12fold seeSM ). Nevertheless, both two-length-scale and three-length-scale models reproduce DDQC patterns. The estimated parameters reproduce the inherent ratio of the length scales $q_{0}/q_{1}=1.948\approx 2\cos(\pi/12)$ . The marginal likelihood indicates that the two length scale is favourable ( Figure 4(c)).

We summarise the results of a variety of patterns in Figure 4. For each target pattern, we can reproduce visually similar patterns, and the most probable number of length scales. In two-dimensional systems, stripe and hexagonal are the two most popular patterns under one length scale. These patterns are obtained from the conventional phase-field crystal model Elder:2004 . The marginal likelihood calculated in the BM-PDE indeed estimates that one length scale is favourable ( Figure 4(a,b)). The difference between stripe and hexagonal patterns appears in the mean density $\bar{\psi}$ . A quadratic nonlinear term is necessary to reproduce hexagonal patterns, and this implies that it appears at $|\bar{\psi}|\gg 0$ Elder:2004 . When $\bar{\psi}\simeq 0$ , stripe patterns appear. The estimated values of the mean density are consistent with the results from the phase diagram reported in literatureElder:2004 ; we obtain the estimated mean density $\hat{\bar{\psi}}\simeq-0.23$ and $\hat{\bar{\psi}}=-0.05$ for the hexagonal and stripe target patterns, respectively (see also Supplementary Fig. LABEL:SMhist.stripe and LABEL:SMhist.hexagonal seeSM ). The BM-PDE automatically estimates appropriate parameters from an artificially synthesised snapshot of the target pattern.

double gyroid and Frank Kasper patterns

We discuss an application of the BM-PDE to nontrivial three-dimensional patterns. The target patterns are double gyroid (DG), shown in Fig. 4(d) and Frank Kasper (FK) A15, shown in Fig. 4(e). The DG structure has two separate domains, each of which has degree-three verticesScherer:2013 . The DG structure was predicted theoreticallyPodneks:1996 ; Shi:1999 and numerically Zhang:2008 , but mainly discussed in similar but a different class of models of block copolymersnonomura:2001 . FK A15 has been studied in metallic alloy and soft materialsBates:2019 . For example, the self-consistent field theory, designed to describe block copolymers, is capable of reproducing this pattern, but to our knowledge, this structure has not been reported within the framework of PFC.

In three dimensions, the order parameter may be defined similarly to that in two dimensions (see Methods). Structural candidates in three dimensions are far richer than two-dimensional patterns. In fact, during the sampling, cylindrical patterns with a hexagonal lattice, BCC, and other patterns appear. These patterns can be stable under a certain region of parameter spaceJaatinen:2010 . The BM-PDE can reproduce both the DG and FK A15 patterns in all models ( $m=m_{1},m_{2},m_{3}$ ). The generated structure is the same as the target pattern, which is evident from the peaks in the Fourier space ( Figure 4(d,e) and see also Supplementary Sec. LABEL:SM.sec.gyroid seeSM ). The real space structures of the two patterns are overlapped by translation. The log marginal likelihood in Fig. 4(d,e) shows that the target patterns of DG and FK A15 are expressed most likely by the two length scale and one length scale, respectively.

The DG structure has two length scales with their ratio $1.15$ (see also Supplementary Sec. LABEL:SM.sec.gyroid seeSM ). In the one-length-scale model, by taking $a_{0}\gg 0$ , several modes with slightly different length scales may become unstable so that the double gyroid pattern is realised. In fact, it was discussed that DG appears not at $a_{0}=\epsilon\approx 0$ with the small number $\epsilon$ , but at $a_{0}\gtrsim\mathcal{O}(1)$ Podneks:1996 ; Shi:1999 . In addition, this pattern appears between the stripe patterns and cylindrical (hexagonal) patterns in the phase diagram, namely $1\gg|\bar{\psi}|\gtrsim 0$ . On the other hand, in the model of block copolymers, the FK A15 phase appears near BCC patterns in the phase diagram Bates:2019 . This suggests that $|\bar{\psi}|\gg 0$ to obtain FK A15. The estimated mean density $\bar{\psi}$ well agrees with this tendency. (Supplementary Figs. LABEL:SMhist.gyroid and LABEL:SMFKA15 seeSM ). The two models $m_{1}$ and $m_{2}$ have a comparable probability for both patterns. The DG pattern chooses $m_{2}$ possibly because $m_{2}$ has a wider range of spectrum amplitude in the phase diagram. The target structure of FK A15 has ellipsoidal domains. The generated structures with the estimated parameters are also ellipsoidal (see Supplementary Fig.LABEL:SMsphericity seeSM for quantitative results on sphericity). Such deformation is also reported for FK A15 structure made of block copolymersLee:2014 .

We should note that the target structure of the FK A15 is taken from the experimental result Jayaraman:2018 . From the X-ray scattering data, the position and amplitude of peaks in the Fourier space are identified (Fig.2(b)). We use the positions of the peaks and their amplitude in the data in Jayaraman:2018 .

Robustness against noise

We hypothesise that the robust estimation of the model and parameters for the target pattern without ground truth originates from robustness against noise. To see this, we add Gaussian white noise to the target pattern in Fig. 1(a) and study its effect on the estimation of the parameters. We focus on the two-length-scale model $m=m_{2}$ and the estimation of wavelength because this parameter has a narrow acceptable range. For the DDQC, the wavenumber is required to be close to $0.51$ , and otherwise, the pattern cannot be generated because the mode coupling between two length scales does not occur. Figure 5(a) shows that even when the amplitude of the noise is about $20\%$ of the amplitude of the density field of the target pattern, the estimation of the parameter successfully works. Beyond the noise amplitude, the error of the estimated wavenumber becomes large, and the fraction of the patterns different from the target pattern increases.

We compare the BM-PDE with the standard regression methods in which parameters $\mu$ are estimated by minimising the cost function $\|\partial_{t}\psi-f(\psi;\mu)\|$ under an appropriate norm (see Methods). This method relies on the state noise added in the dynamical equation, and thus, does not give a good estimate for the measurement noiseVoss:2004 (see also Supplementary Sec.LABEL:SM.sec.relatedworks seeSM ). In fact, Fig. 5(a) shows that the estimated wavenumber deviates from ground truth even under $0.1\%$ noise amplitude.

In contrast with the noiseless target pattern, which has its log marginal likelihood $F[\beta]$ monotonically decreasing as $\beta$ increases in Fig. 3(e), the log marginal likelihood for the target pattern with noise has a minimum at the finite $\beta$ as in Fig. 5(b). The minimum demonstrates the optimal noise corresponding to the noise amplitude in the target patternTokuda:2017 . The noise level $\beta^{-1}$ at the minimum of log marginal likelihood increases as the amplitude of noise in the target increases. At the large noise amplitude $\gtrsim 30\%$ , the minimum $\beta$ reaches a comparable value with the cost function at the gap between two structures in Fig. 3(b): QC and hexagonal patterns. Thus, DDQC can no longer be reproduced. Interestingly, the minimum of the log marginal likelihood is also observed for the target pattern synthesised by the function of equation (14) (Fig. 3(c)). In this case, there is no ground truth parameter that reproduces exactly the same pattern as the target pattern. The optimal $\beta$ at the minimum of $F[\beta]$ describes the deviation of the target pattern from the patterns that the models can reach. Without the Bayesian inference, we cannot successfully estimate the uncertainty, which plays a similar role to noise.

The BM-PDE works for the damaged target structure in which some spatial information is lost. We demonstrate the estimation for the damaged data in Supplementary Sec.LABEL:SM.sec.noise seeSM . The BM-PDE does not rely on the time and spatial derivatives of the target structure. Therefore, the damage in the data can be randomly distributed, namely the information of neighbouring data point is not necessary. This is not the case if the information of spatial derivatives is necessary for the estimation.

Three-dimensional dodecagonal pattern expressed by a PDE

Using the estimated parameters, we may investigate additional physical insights of the model. To demonstrate it, we consider DDQC in three dimensions (Fig. 6(b,c)). Although the icosahedral QCs have been found in the PDE model using two-length-scale PFCSubramanian:2016 , other three-dimensional structures are largely unexplored Damasceno:2017 . The axisymmetric dodecagonal structure has twelve-fold symmetry along the axis ( $c$ -axis) of axisymmetry (Fig. 6(b,c)).

To generate this structure, larger system size is required, but it is computationally demanding. Moreover, as we will discuss later, we found that the kinetics of the formation of the DDQC in three dimensions is fundamentally different from other structures. Because the FK A15 is known as the simplest approximant of the DDQCIacovella:2011 , it is fair to assume that the DDQC may appear near the estimated parameters of the Frank-Kasper A15. Therefore, we focus on the estimated parameters, and solve the estimated models with larger systems size in a longer time scale.

The structure discussed in the previous sections, including the DDQC in two dimensions, may appear from a random initial condition by quenching. The random initial structure reorganises to the desired pattern without noise. In addition, the system reaches its stationary state before $t=10^{2}-10^{3}$ . Note that the time scale $t=1$ corresponds to the relaxation time of the structure when $a_{i}=1$ . The formation of the DDQC in three dimensions is entirely different from the two-dimensional structure. It requires not only a larger system size but also an annealing process. Without the state noise, the initial random structure forms micelle-like isotropic structures (Fig.6(a)), which is frozen. This is particularly the case for the system size $N>64$ . To overcome the trap at the metastable structure, we add the state noise into the PDE model (see Methods), and anneal the structure by decreasing the noise amplitude. We found that the annealing process is very slow; the DDQC appear around $t\gtrsim 5\times 10^{4}$ .

Figure 6(b) shows the generated structure for the system size $N=256$ . The generated structure is periodic along the $c$ -axis. We denote the other two axes as the $a$ -axis and the $b$ -axis, and in the Fourier space, these axes are denoted by $k_{a}$ , $k_{b}$ , and $k_{c}$ . As in Fig. 6(c), the structure has twelve-fold symmetry in the plane perpendicular to the $c$ -axis. The twelve-fold symmetry arises from the layered structure in which two layers with the hexagonal symmetry rotate $\pi/6$ with each other (Fig. 6(b)). The centre of the dodecagonal ring is located between the two layers. There are 62 main peaks in the Fourier space. Along the symmetry axis, there are five 12-fold rings and two spots along the $c$ -axis. From the peaks in the Fourier space shown in Fig. 6(c), it is evident that there are two length scales $q_{0}$ and $q_{1}$ whose ratio is $q_{1}/q_{0}=2/\sqrt{5}$ . This value is, in fact, close to the estimated ratio of parameters $\hat{q}_{1}/\hat{q}_{0}=0.889$ . The layered structure with the dodecagonal symmetry has been reported in particle based simulationsDzugutov:1993 ; Damasceno:2017 . We have tried both one-length-scale and two-length-scale models with the estimated parameters for the Frank-Kasper A15 structure. All the results shown in Fig. 6 are by the two-length-scale model. The one-length-scale model cannot reproduce the quasi-crystal. This is natural because more than two length scales are necessary to form QCs.

Using the obtained stable DDQC, we also study their fluctuations under the fixed amplitude of the noise. Figure 6(d) demonstrates the reconstruction of a dodecagonal ring. This process is due to the change of the position of the centre along the $c$ -axis to the ring (see the insets of Fig. 6(d)). This reconstruction has also been observed in a particle-based modelDzugutov:1995 . The fluctuation is characterised by the intermediate scattering function, $F({\bf k},t)=\left<\psi_{\bf k}(0)\psi^{\dagger}_{\bf k}(t)\right>$ for the wave vectors corresponding to the peaks in the Fourier space (Fig.6(e)). Here, $\dagger$ indicated complex conjugate. The relaxation time of $F({\bf k},t)$ indicates the diffusive time scale of the structure at ${\bf k}$ (Fig. 6(f)). The peaks other than the main peaks of the DDQC shown in Fig. 6(c) decay quickly in a short time scale of $t\sim 10^{3}$ . Among the main peaks, the larger peaks (cyan in Fig. 6(e)) are stable against noise even in the longer time scale of $t\gtrsim 10^{4}$ , whereas the smaller peaks (magenta in Fig. 6(e)) show diffusive behaviour. This slower diffusion at $t\gtrsim 10^{4}$ corresponds to phason flips. In fact, from the inverse Fourier transformation, we found the peaks with the intermediate amplitude (magenta in Fig. 6(e)) correspond to the centres of the dodecagonal rings. We stress that these analyses of fluctuations are achieved by the estimation of the dynamical PDE from the stationary target pattern in the previous section.

Discussions and Conclusion

In summary, we propose the inverse problem of equation discovery for a stationary pattern without ground truth. We successfully estimate the best PDE for complex patterns. Our approach has several key features compared with the previous inference methods; it is designed for the data of a stationary pattern without ground truth, In addition, the method has generality to apply a wide class of PDE models and target patterns. Here, we summarise these features.

The estimation for stationary data: The previous approaches to estimate the PDE from data are estimation of dynamical equations from dynamical non-stationary data. In this case, the data has information about trajectories of dynamics. In this work, we may estimate dynamical equations from stationary data, which is one snapshot of the stationary pattern. Such estimation looks impossible for two reasons: One is lack of information of transient dynamics. The second reason is that in pattern formation, the dynamical equation is nonlinear, and therefore, $f(\psi)=0$ has multiple stationary solutions. Only one solution from the multiple solutions is not enough to reconstruct $f(\psi)$ uniquely. We take the inverse problem for the stationary pattern as an estimation of a dynamical PDE from marginalised initial conditions (see Supplementary Sec.LABEL:SM.sec.relatedworks.inverse.problem seeSM for comparison to other methods). The stable structure should be generated from a wide range of initial conditions. In the BM-PDE, the estimation of the model reproducing a target pattern as a stable solution by marginalising the initial conditions as in equation (3) (see also Fig. LABEL:fig_otherfamily(j) for demonstration). This is in contrast with the method to estimate time-series data, for example, using data assimilation, in which a single initial condition is estimated from dataEvensen:2009 .

The order parameters play an important role in the BM-PDE. The order parameters can extract symmetries of the pattern and identify the two patterns that are generated with the same parameters but the different initial conditions. This identification is a necessary step to marginalise the initial patterns.

The estimation without ground truth: The BM-PDE not only works for a numerically produced pattern, which has ground truth of parameters, but also for an synthesised pattern by superposition of plane waves of equation (14). In the latter case, our family of models does not include ground truth. To demonstrate the ability to estimate the model for the structure without ground truth, we have used experimental data (Fig.2). Even in this case, we have successfully estimated the best PDE model and its parameters. This is because the BM-PDE quantifies the uncertainty by estimating the optimal noise amplitude, and is robust against noise.

Our method may be used to find a structure that we do not know how to reproduce by a PDE model. We have found the three-dimensional DDQC from the estimated parameters for its approximant. To our knowledge, this is the first three-dimensional DDQC found in the PDE model.

The generality of BM-PDE: In the BM-PDE, the choice of models and the cost function are independent. Therefore, we may replace the models with any other PDE models. To demonstrate this feature, we have shown the estimation of parameters not only in the linear operator but also in nonlinear terms (Supplementary Sec. LABEL:sec.family.Hermite and LABEL:sec.family.nonlinear seeSM ). The choice of the target pattern is also flexible. By comparing two patterns in the space of the order parameter, we can estimate the model that can reproduce a similar pattern in the sense that it has the same symmetry to the target pattern. Thanks to this feature, the BM-PDE works for the noisy target.

Methods

The main framework of BM-PDE consists of three parts: a family of PDE models, characterisation of a structure with order parameter, and statistical inference. Here, we summarise their basic steps.

PDE models

We consider a pattern described by a scalar density field $\psi({\bf x})$ in a box with a periodic boundary condition, and ${\bf x}\in[-L_{x}/2,L_{x}/2]\times[-L_{y}/2,L_{y}/2]$ in two dimensions and ${\bf x}\in[-L_{x}/2,L_{x}/2]\times[-L_{y}/2,L_{y}/2]\times[-L_{z}/2,L_{z}/2]$ in three dimensions. The density field $\psi({\bf x})$ is obtained by an unknown model of a partial differential equation. In this study, we focus on a family of the phase-field crystal (PFC) models given by

\displaystyle\partial_{t}\psi

\displaystyle=\mathcal{L}\psi+\Delta\psi^{3}

(6)

with the linear operator denoted by $\mathcal{L}$ and the nonlinear term is given in the second term. The family is constructed by modifying the linear operator $\mathcal{L}$ so that the system has one or more length scales. The simplest PFC model is a conserved version of the Swift-Hohenberg (SH) equation Provatas:2011 . The equation is given by

\displaystyle\partial_{t}\psi

\displaystyle=\Delta\left[-\epsilon\psi+\left(q_{0}+\Delta\right)^{2}\psi+\psi^{3}\right]

(7)

where the total mass is conserved as

\displaystyle\bar{\psi}

\displaystyle=\frac{1}{V}\int\psi({\bf x})d{\bf x}.

(8)

Here, $V$ is the total volume (area in two dimensions) in the system, and $q_{0}$ is a characteristic wavenumber corresponding to the length scale $2\pi/q_{0}$ . The parameter $\epsilon$ controls whether the uniform state $\psi({\bf x})=\bar{\psi}$ is stable $\epsilon\lesssim 0$ or unstable $\epsilon\gtrsim 0$ . The precise value of the threshold is dependent on other parameters and a type of patterns. The PFC equation reproduces a stripe (also called lamellar or smectic) and hexagonal patterns in two dimensionsPismen:2006 ; Elder:2004 , and a lamellar, hexagonal cylinder, BCC, and hexagonal closed packing patterns Zhang:2008 ; Jaatinen:2010 . A finite mean density $\bar{\psi}$ plays a role as the quadratic nonlinear term in equation (7). This can be seen by subtracting the mean density as $\psi\rightarrow\psi+\bar{\psi}$ in equation (7) as,

\displaystyle\partial_{t}\psi

\displaystyle=\Delta\left[\left(-\epsilon+3\bar{\psi}^{2}\right)\psi+\left(q_{0}+\Delta\right)^{2}\psi+3\bar{\psi}\psi^{2}+\psi^{3}\right].

(9)

The PFC equation has a single characteristic length at which a real part of the eigenvalue is positive (or, at least, negative but close to zero). The linear spectrum is shown in Fig. LABEL:SMspectrum. To extend equation (7) for the arbitrary number of length scales, we use $m$ characteristic wavenumbers, $q_{0},q_{1},\cdots,q_{m-1}$ , and the values of the spectral, $a_{0},a_{1},\cdots,a_{m-1}$ , at the wavenumber $k=q_{i}$ . A family of our models is conveniently described by the equations for the Fourier transform of the density field, that is, $\tilde{\psi}({\bf k})=\mathcal{F}[\psi({\bf x})]$ . Our models corresponding to equation (6) are given by

\displaystyle\partial_{t}\tilde{\psi}({\bf k})

\displaystyle=\mathcal{L}_{k}\tilde{\psi}({\bf k})+\mathcal{F}\left[\Delta\psi({\bf x})^{3}\right],

(10)

and its linear operator is expressed in the Fourier space as

\displaystyle\mathcal{L}_{k}

\displaystyle=a_{0}S_{0}(k)+a_{1}S_{1}(k)+\cdots+{\color[rgb]{0,0,0}a_{m-1}S_{m-1}(k)}+k^{2}\left({\color[rgb]{0,0,0}q_{0}}^{2}-k^{2}\right)^{2}\left(q_{1}^{2}-k^{2}\right)^{2}\cdots\left({\color[rgb]{0,0,0}q_{m-1}^{2}}-k^{2}\right)^{2}.

(11)

The function $S_{i}(k)$ for $i\in[0,m-1]$ is chosen so that the coefficient $a_{i}$ corresponds to the peak as a function of $k$ of $\mathcal{L}_{k}$ at $k=q_{i}$ (see Fig. 3(f-h)). Since we may freely choose a unit of length scale, we fix to be $q_{0}=1$ when $i\geq 2$ in $m=m_{i}$ . This implies that we impose the length scale $2\pi/q_{i}$ . The concrete form of $S_{i}(k)$ is shown in Supplementary Information Sec.LABEL:SM.sec.PDE.models seeSM . The parameter $s_{0}$ describes the sharpness of the peaks. To make the spectrum sharp enough, we chose the parameter to be $s_{0}=100$ for the one-length and two-length models, and $s_{0}=2000$ for the three-length model. Both $a_{i}$ and $q_{i}$ are chosen as parameters to be optimised. We have other parameters such as the system size $L$ in each direction and the mean density of a pattern $\bar{\psi}$ . In this study, we use the periodic boundary condition. A set of parameters is thus $\mu=\{L_{\alpha},\bar{\psi},a_{0},q_{0},\ldots,a_{l-1},q_{l-1}\}$ with $\alpha\in[1,\ldots,d]$ in space dimension $d$ .

In order to make a pattern, higher-order spatial derivatives are necessary. Polynomial expansion in terms of the wavenumber $k$ instead of equation (11) may be available to make several length scales (see Fig. 1(e)). Nevertheless, it suffers from the large value of the coefficient of each order in the polynomials because we may not use prior distribution to confine the parameter space (see also Supplementary Sec. LABEL:SM.sec.polynomial.expansion seeSM ).

Numerical simulations of PDE models

Numerical simulations of the PDEs are performed using the pseudo-spectral method in which the linear terms are computed in the Fourier space, and the nonlinear terms are computed in real space. Since our PDEs contain higher-order derivatives, we use the operator-splitting methodCox:2002 ; Aranson:2015 . Both real and Fourier spaces are discretised into $N^{d}$ meshes in $d$ -dimensional space. Instead of changing the system size $L_{i}$ in each dimension $i\in[1,d]$ under the periodic domain, we change the mesh size $dx_{i}$ so that the system size becomes $L_{i}=(N-1)dx_{i}$ .

The number of mesh points is fixed to be $N=128$ in two dimensions and $N=32$ in three dimensions. The larger $N$ is better in terms of accuracy, but computational time is scaled roughly as $\mathcal{O}(N^{d})$ . In REMC, we need to simulate it in $N_{\rm rep}$ replicas, and therefore, the choice of $N$ is made by the balance between accuracy and realistic computational time. In addition, the larger system size suffers from longer relaxation time and a higher probability that topological defects such that dislocations and disclinations appear. We set the total number of steps to be $10^{5}$ with a time step $dt=0.01$ . We confirmed this is enough to obtain the stationary patterns studied in this work, but may be changed depending on the pattern of interest. Note that statistical inference in this study is entirely independent of the algorithm of numerical simulations solving a PDE. An efficient algorithm would improve the performance of estimation, and one may replace the numerical scheme suit for one’s purpose.

Formation of three-dimensional DDQC

To generate the DDQC in three dimensions, we have to add noise in equation (5),

\displaystyle\partial_{t}\psi({\bf x},t)

\displaystyle=\mathcal{L}_{\hat{\mu}}^{(m)}\psi({\bf x},t)+\mathcal{N}[\psi({\bf x},t)]+\gamma(t)\xi({\bf x},t).

(12)

Note that the noise is state noise, and nothing to with the observation noise. In the annealing process, the amplitude of the noise $\gamma(t)$ is decreased in time. The annealing schedule is chosen as $\gamma(t)=0.198/\log(t+\tau_{a})$ with $\tau_{a}=10^{4}$ . This choice ensures $\gamma(t)\sim 1/\log t$ for large $t$ Geman:1984 . The Gaussian white noise with the zero mean $\left<\xi({\bf x},t)\right>=0$ and its variance

\displaystyle\left<\xi({\bf x},t)\xi({\bf x}^{\prime},t^{\prime})\right>

\displaystyle=-2\Delta\delta({\bf x}-{\bf x}^{\prime})\delta(t-t^{\prime})

(13)

is used. The Laplacian in equation (13) ensures the conservation of the density and ensures to reach the equilibrium state. We may also use other statistical properties of the noise as long as the density conservation is ensured. For example, we have tested $\left<\xi({\bf x},t)\xi({\bf x}^{\prime},t^{\prime})\right>=2\delta({\bf x}-{\bf x}^{\prime})\delta(t-t^{\prime})$ with $\int\xi({\bf x}^{\prime},t^{\prime})d{\bf x}=0$ . We have confirmed in all cases the DDQC appears when the amplitude of the noise $\gamma$ is decreased.

Fluctuation of the DDQC is studied by fixed $\gamma$ . We set it to be $\gamma=0.171$ . Note that the DDQC is destroyed for the noise with the amplitude of $\gamma^{2}=0.198$ .

Target patterns and characterisation of patterns

Target Pattern

A target pattern $\psi^{*}({\bf x})$ is prepared in two ways. One is a numerical solution of equation (5) for given parameters under a given model. The resulting pattern is numerically transformed into the Fourier space, and then the order parameter $\boldmath{\Psi}^{*}_{l}$ is calculated.

The second way is completely independent of the models. The target pattern is expressed as a density field $\psi^{*}({\bf x})$ that is a superposition of the cosine function. The simplest case is a stripe pattern, which is described only by one wave in one direction in a two-dimensional space. Similarly, a hexagonal pattern is expressed by two-dimensional waves in three directions. The target pattern is expressed as

\displaystyle\psi^{*}({\bf x})

\displaystyle=\sum_{i}b_{i}\cos({\bf q}^{*}_{i}\cdot{\bf x})

(14)

where the wave vectors ${\bf q}^{*}_{i}$ are chosen at the position appropriate to express symmetries of the target pattern. The amplitude of each mode $b_{i}$ is also chosen properly. We numerically make the Fourier transform of equation (14) to obtain $\hat{\psi}^{*}({\bf k})=\mathcal{F}[\psi^{*}(\bf x)]$ and calculate the Fourier spectrum $|\hat{\psi}^{*}({\bf k})|$ from which we obtain the order parameter. The Fourier transform of equation (14) defined in the infinite domain is expressed by the superposition of the delta function at the position of ${\bf q}^{*}_{i}$ . Nevertheless, the numerical Fourier transform in the bounded domain results in peaks smeared around ${\bf q}^{*}_{i}$ . To remove the artefact, we set $|\psi^{*}({\bf k})|=0$ except for the region $|\psi^{*}({\bf k})|>\alpha{\rm max}|\psi^{*}({\bf k})|$ . Here, the value of $\alpha$ is chosen so that the peaks of the minimal height are left. We choose $\alpha=0.6$ for the two-dimensional target patterns, whereas $\alpha=0.01$ for the three-dimensional target patterns.

The dodecagonal QC pattern in two dimensions is synthesised by $\psi=\sum_{i=1}^{12}\cos({\bf q}^{*}_{i}\cdot{\bf x})$ in which the wave vectors ${\bf q}^{*}_{i}$ are chosen at the position of the vertices of the hexagon with a radius $|{\bf q}^{*}_{1}|=2\pi/\sqrt{2+\sqrt{3}}$ and the hexagon with a radius $|{\bf q}^{*}_{2}|=2\pi$ rotated by $\pi/12$ . The DG pattern is expressed by 24 wave vectors of ${\bf q}^{*}=(\pm 2,\pm 1,\pm 1)$ and 12 wave vectors of ${\bf q}^{*}=(\pm 2,\pm 2,0)$ with their permutation along the $x,y,z$ directionsSchnering:1991 ; Yamada:2004 . The amplitude of the latter wave vector is $\sqrt{8/6}\simeq 1.15$ times longer than the former waves. The FKA15 pattern is expressed by 24 wave vectors ${\bf q}^{*}=(\pm 2,\pm 1,0)$ , 24 wave vectors ${\bf q}^{*}=(\pm 2,\pm 1,\pm 1)$ , 6 wave vectors of ${\bf q}^{*}=(\pm 2,0,0)$ with their permutation along the $x,y,z$ directionsImperorClerc:2012 .

order parameters

To assure translational invariance, we use a spectrum of the Fourier transform of the pattern and expand it by the basis functions, each of which expresses certain point group symmetries. In two dimensions, the basis functions reflect $n$ -fold rotational symmetry, as shown in Fig. 1(b), whereas in three dimensions, spherical harmonics expansion is used. Projection of the Fourier spectrum of the pattern onto the basis function is given by $A_{l,\pm}$ in two dimensions, and $A_{lm}$ in three dimensions.

The order parameter is a rotational invariant form of the quantity $A_{lm}$ with $l\in[0,l_{0}]$ and $m\in\{\pm 1\}$ in two dimensions and $m\in[-l,l]$ in three dimensions. In two dimensions, $A_{l,\pm 1}$ is described by

\displaystyle A_{l,\pm 1}\left[\psi\right]

\displaystyle=\int|\psi({\bf k})|\begin{pmatrix}\cos l\theta_{k}\\ \sin l\theta_{k}\end{pmatrix}d{\bf k}^{2}

(15)

where $+1$ $(-1)$ corresponds to $\cos l\theta_{k}$ $(\sin l\theta_{k})$ , respectively, and $\theta_{k}$ is a polar angle in the Fourier space. We use the Fourier transform of the pattern as

	$\displaystyle\tilde{\psi}({\bf k})$	$\displaystyle=\int\psi({\bf x})e^{i{\bf k}\cdot{\bf x}}d{\bf x}$		(16)
	$\displaystyle\psi({\bf x})$	$\displaystyle=\int_{\bf k}\tilde{\psi}({\bf k})e^{-i{\bf k}\cdot{\bf x}}$		(17)

where the volume in the Fourier space is $\int_{\bf k}=\frac{1}{(2\pi)^{d}}d^{d}{\bf k}$ . We denote $A_{l,\pm}$ in the vector form as ${\bf A}_{l}=(A_{l,+},A_{l,-})$ . The maximum mode is denoted by $l_{0}$ . In three dimensions, $A_{lm}$ is given by

\displaystyle A_{lm}\left[\psi\right]

\displaystyle=\int|\psi({\bf k})|Y_{l}^{m}(\theta_{k},\varphi_{k})d{\bf k}^{3}

(18)

with spherical harmonics $Y_{l}^{m}(\theta_{k},\varphi_{k})$ in the spherical coordinates of the Fourier space $(k,\theta_{k},\varphi_{k})$ . Note that $m$ in the superscript of $Y_{l}^{m}(\theta_{k},\varphi_{k})$ and subscript of $A_{lm}$ should not be confused by $m$ describing a model in $M$ . The zeroth mode $l=0$ corresponds to the mean amplitude of $\tilde{\psi}({\bf k})$ , which is independently considered by $\bar{\psi}$ . We, therefore, use the sum for $l\in[1,l_{0}]$ in the cost function. The maximum mode is denoted by $l_{0}$ . We use the convention of spherical harmonics

\displaystyle Y_{l}^{m}(\theta,\varphi)

\displaystyle=\sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}}P_{l}^{m}(\cos\theta)e^{im\varphi}

(19)

where $P_{l}^{m}(\cos\theta)$ is associated Legendre polynomial with integers $l$ and $m\in[-l,l]$ . Any continuous function on a unit sphere may be expanded.

We define the order parameter ${\bf\Psi}[\psi({\bf x})]=\{\Psi_{l}[\psi({\bf x})]\}_{l=1}^{l_{0}}$ by a rotationally invariant form of the coefficients, $A_{l,\pm}$ or $A_{lm}$ .

\displaystyle\Psi_{l}=\|{\bf A}_{l}\|

\displaystyle\equiv\sqrt{A_{l,+1}^{2}+A_{l,-1}^{2}}

(20)

in two dimensions, and

\displaystyle\Psi_{l}=\|A_{lm}\|

\displaystyle\equiv\sqrt{\frac{4\pi}{2l+1}}\sqrt{\sum_{m=-l}^{l}(-1)^{m}A_{l,m}A_{l,-m}}

(21)

in three dimensions. Here, the prefactor is included because $(Y_{l}^{0})^{2}+\sum_{m=1}^{l}Y_{l}^{m}(\theta,\varphi)Y_{l}^{m*}(\theta,\varphi)=(2l+1)/(4\pi)$ , and the sum of $|A_{lm}|^{2}$ scales $2l+1$ .

We numerically evaluate $A_{l}$ for patterns $\psi(\bf x)$ and $A^{*}_{l}$ for a target pattern $\psi^{*}({\bf x})$ . Since both real space and Fourier space density fields are expressed by values at a finite number of mesh points, the range of a mode $l$ is truncated at the maximum mode $l_{0}$ . The larger mode extracts a finer structure in the Fourier spectrum, and the structure finer than the mesh size is invalid. We thus take $l_{0}=N$ . Note that for odd $l$ , $A_{lm}\simeq 0$ and therefore the dimension of $\bm{\Psi}$ is $l_{0}/2$ .

Statistical inference

Bayesian formulation

We may extend our model naturally toward the Bayesian formulation, which enables us not only to choose the optimal PDE, i.e. parameters and the number of the characteristic length scales, but also to evaluate their uncertainty. To do this, we assume the order parameter ${\bf\Psi}[\psi_{s}]$ of the stationary pattern $\psi_{s}$ is observed as that of the target pattern $\psi^{*}$ with the additive noise $\xi_{l}$ :

\displaystyle\Psi_{l}^{*}

\displaystyle=\Psi_{l}(\psi_{s};\psi_{0},\mu)+\xi_{l},

(22)

where $\xi_{l}$ is the random variable for each mode $l$ distributed according to zero-mean Gaussian distribution with variance $\beta^{-1}\geq 0$ . The noise $\xi_{l}$ and the corresponding inverse temperature $\beta$ play a role of the uncertainty of the measurement. Here, $\psi_{s}$ and $\psi_{0}$ are the stationary and initial states of $\psi$ in equation (6), respectively, and $\mu$ is the set of parameters. The dependence of $\Psi_{l}[\psi_{s}]$ on $\psi_{0}$ and $\mu$ is explicitly represented by $\Psi_{l}(\psi_{s};\psi_{0},\mu)$ . The assumption equation (22) is equivalently represented as the conditional probability density

\displaystyle p(\Psi_{l}^{*}\mid\psi_{0},\mu,\beta)

\displaystyle=\sqrt{\frac{\beta}{2\pi}}\exp\left\{-\frac{\beta}{2}[\Psi_{l}^{*}-\Psi_{l}(\psi_{s};\psi_{0},\mu)]^{2}\right\}.

(23)

We consider the parameter estimation of $\mu$ by marginalising $\psi_{0}$ . Hereafter, the discretization $\psi_{0}=\{\psi_{0}^{(j)}\}_{j=1}^{N^{d}}$ and the reparametrization $\psi_{0}^{(j)}\rightarrow\bar{\psi}+\psi_{0}^{(j)}$ are also considered, where $\bar{\psi}$ and $\psi_{0}^{(j)}$ are the mean density and the (relative) density at the mesh point $j$ , respectively. By Bayes’ theorem, the conditional joint probability density of $\psi_{0}$ and $\mu$ under given $\{\Psi_{l}^{*}\}_{l=1}^{l_{0}}$ , $\beta$ and the model class $m$ is represented as

	$\displaystyle p(\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\psi_{0},\beta,m)$	$\displaystyle=\frac{p(\mu\mid m)}{p(\{\Psi_{l}\}_{l=1}^{l_{0}}\mid\beta,m)}\prod_{l=1}^{l_{0}}p(\Psi_{l}^{*}\mid\psi_{0},\mu,\beta),$
		$\displaystyle\propto\exp\left[-\frac{\beta}{2}E(\psi^{*},\psi_{s};\psi_{0},\mu)\right]$		(24)

where $m$ denotes the number of the characteristic length scales such as equations (LABEL:SHK1), (LABEL:SHK2) or (LABEL:SHK3). Here, $p(\mu\mid m)$ and $p(\psi_{0})$ are the prior distribution defined as the uniform distribution, and the marginal likelihood $p(\{\Psi_{l}\}_{l=1}^{l_{0}}\mid\beta,m)$ is given by

\displaystyle p(\{\Psi_{l}\}_{l=1}^{l_{0}}\mid\beta,m)

\displaystyle=\left(\frac{\beta}{2\pi}\right)^{\frac{l_{0}}{2}}\int\exp\left[-\frac{\beta}{2}E(\psi^{*},\psi_{s};\psi_{0},\mu)\right]p(\psi_{0})p(\mu\mid m)d\psi_{0}d\mu.

(25)

The dependence of $E[\psi^{*},\psi_{s}]$ on $\psi_{0}$ and $\mu$ is explicitly represented by $E(\psi^{*},\psi_{s};\mu,\psi_{0})$ . Note that we assume the causality

\displaystyle p(\psi_{0},\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\psi_{0},\beta,m)=p(\psi_{0})p(\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\psi_{0},\beta,m),

(26)

which ignores (i) the dependence of $\psi_{0}$ on $\{\Psi_{l}^{*}\}_{l=1}^{l_{0}}$ and $\beta$ , and (ii) the correlation between $\psi_{0}$ and $\mu$ . This assumption reflects our ansatz that $\psi_{0}$ is not uniquely determined only by $\psi^{*}$ (or $\{\Psi_{l}^{*}\}_{l=1}^{l_{0}}$ ). Here $\psi_{0}$ is treated as a latent variable. By marginalising out $\psi_{0}$ , the posterior distribution of $\mu$ is given by

\displaystyle p(\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\beta,m)

\displaystyle=\int p(\psi_{0},\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\psi_{0},\beta,m)d\psi_{0}.

(27)

The posterior mean estimator $\hat{\mu}$ , i.e. the mean of $p(\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\beta,m)$ , is adopted as our best parameter set. The standard deviation of $p(\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\beta,m)$ plays a role of the error in $\hat{\mu}$ .

We consider both the hyperparameter estimation of $\beta$ and model selection of $m$ efron1973stein ; akaike1998likelihood ; mackay1992bayesian ; Bishop:2006 ; Tokuda:2017 . By Bayes’ theorem, the joint probability density of $\beta$ and $m$ under given $\{\Psi_{l}^{*}\}_{l=1}^{l_{0}}$ is represented as

\displaystyle p(\beta,m\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}})

\displaystyle=\frac{p(\{\Psi_{l}^{*}\}_{l=1}^{l_{0}}\mid\beta,m)p(\beta)p(m)}{p(\{\Psi_{l}^{*}\}_{l=1}^{l_{0}})},

(28)

where $p(\{\Psi_{l}^{*}\}_{l=1}^{l_{0}})$ is the normalisation constant. Here, $p(\beta)$ and $p(m)$ are the prior distributions defined as the uniform distribution. The maximum a posteriori estimator, or equivalently the empirical Bayes estimator in this setup, is adopted as the pair of our optimal model and temperature

	$\displaystyle(\hat{\beta},\hat{m})$	$\displaystyle=\mathop{\rm arg~{}max}\limits_{\beta,m}p(\beta,m\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}})$		(29)
		$\displaystyle=\mathop{\rm arg~{}max}\limits_{\beta,m}p(\{\Psi_{l}^{*}\}_{l=1}^{l_{0}}\mid\beta,m).$		(30)

For convenience, the Bayes free energy $F(\beta,m)=-\log p(\{\Psi_{l}\}_{l=1}^{l_{0}}\mid\beta,m)$ is defined. Using the Bayes free energy, we may see the optimal model and temperature $(\hat{\beta},\hat{m})$ minimise $F(\beta,m)$ . If $\partial F/\partial\beta=0$ is satisfied at $\beta=\hat{\beta}$ , then we obtain the self-consistent equation

\displaystyle\hat{\beta}

\displaystyle=\frac{1}{\langle E(\psi^{*},\psi_{s};\psi_{0},\mu)\rangle_{\hat{\beta}}},

(31)

where

\displaystyle\langle\cdots\rangle_{\beta}=\int(\cdots)p(\psi_{0},\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\beta,m)d\psi_{0}d\mu.

(32)

By marginalising out $\beta$ , we can also evaluate the uncertainty of $m$ as the probability

\displaystyle p(m\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}})=\int p(\beta,m\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}})p(\beta)d\beta.

(33)

Note that $p(m_{1}\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}})$ , …, $p(m_{i_{\rm max}}\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}})$ demonstrate the probability of each model $m_{1}$ , …, $m_{i_{\rm max}}$ , respectively, based on the observations $\{\Psi_{l}^{*}\}_{l=1}^{l_{0}}$ .

Setup of a prior distribution

We assume no prior information about parameters and latent variables except their range. The prior density of each variable is defined by the continuous uniform distribution whose support equals the domain of each variable. The prior density of $\psi_{0}$ is defined by

\displaystyle p(\psi_{0})=\prod_{j=1}^{N^{d}}\varphi(\psi_{0}^{(j)}),

(34)

where $\varphi(\psi_{0}^{(j)})$ is the continuous uniform distribution whose support is $\psi_{0}^{(j)}\in[-0.1,0.1]$ . For equation (23) with $m$ length scales, the set of parameters is defined by

\displaystyle\mu

\displaystyle=\{dx,dy,\bar{\psi},a_{0},q_{1},a_{1},q_{2},a_{2},\cdots,q_{m-1},a_{m-1}\},

(35)

where ${\rm dim}(\mu)=2m+2$ in two dimensions. In three dimensions, the mesh size along the $z$ -axis $dz$ is added in the parameters. The prior density of $\mu$ is also defined by

\displaystyle p(\mu\mid m)

\displaystyle=\varphi(dx)\varphi(dy)\varphi(\bar{\psi})\varphi(a_{0})\prod_{i=1}^{m-1}\varphi(a_{i})\varphi(q_{i}),

(36)

where $\varphi(dx)$ , $\varphi(dy)$ , $\varphi(\bar{\psi})$ , $\varphi(a_{i})$ , and $\varphi(q_{i})$ are the continuous uniform distributions, whose supports are respectively $dx\in[1-(1/q^{*}N),1+(1/q^{*}N)]$ , $dy\in[1-(1/q^{*}N),1+(1/q^{*}N),]$ , $\bar{\psi}\in[-1,0]$ , $a_{i}\in[-0.2,0.2]$ , and $q_{i}\in[0,1]$ . Here, $2\pi/q^{*}$ is the wavelength that is used to synthesise the target pattern.

We also assume no prior information for model and hyperparameter; the prior distribution of each variable is defined by the discrete uniform distribution. The prior distribution $p(\beta)$ is also defined by the discrete uniform distribution with $\beta\in\{\beta_{\alpha}\}_{\alpha=0}^{N_{\rm rep}-1}$ , where $\beta_{0}=0$ and

\displaystyle\beta_{\alpha}

\displaystyle=10^{\log_{10}\beta_{\rm min}+\frac{\alpha-1}{N_{\rm rep}-1}\log_{10}(\beta_{\rm max}/\beta_{\rm min})}

(37)

for $\alpha\in\{1,2,\cdots,N_{\rm rep}-1\}$ . Here, we set as $N_{\rm rep}=40$ , $\beta_{\rm min}=10^{-3}$ and $\beta_{\rm max}=10^{2}$ . Equation (37) means that discretization of $\beta$ is finer at the large $\beta$ (lower variance of noise). The prior distribution $p(m)$ is defined by the discrete uniform distribution at $m\in\{m_{i}\}_{i=1,2\ldots i_{\rm max}}$ . Each grid point $(m_{i},\beta_{\alpha})$ can be regarded as the candidate of model selection equally possible in prior.

sampling from a posterior distribution with Monte Carlo method

The realisation of $p(\mu\mid{\bf\Psi}^{*},\psi_{0},m,\beta)$ is carried out by Monte Carlo (MC) sampling in the parameter space. For each point in the parameter space, we compute a stationary pattern $\psi_{s}$ in the model of equation (5) under the randomly chosen initial condition $\psi_{0}$ . Then, the parameters are changed stochastically according to the Metropolis criterion defined by the cost function (energy) $E[\psi,\psi^{*}]$ and inverse temperature $\beta$ . We use the REMC method Geyer:1991 ; Hukushima:1996 with different inverse temperatures $\beta$ in parallel. The REMC is an efficient method for the estimation of the optimal variance $\hat{\beta}^{-1}$ of the observation noise because the method enables us to sample parameters simultaneously under various $\beta$ . The method also makes an efficient sampling to avoid a local trap in the parameter space. Bridge sampling meng1996simulating ; gelman1998simulating is used to calculate $F(m,\beta)$ for each $m$ . The error bars of $F(m,\beta)$ are calculated by the bootstrap resampling efron1992bootstrap .

The joint distribution $p(\psi_{0},\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\beta,m)$ is realised by the Gibbs sampling based on the relation of equation (26), i.e. the alternately iterative sampling from $p(\psi_{0})$ and $p(\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\psi_{0},\beta,m)$ . The sampling from $p(\psi_{0})$ simply follows equation (34). The sampling from $p(\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\psi_{0},\beta,m)$ follows the procedure below. First, we solve the model of equation (6) under a given initial state $\psi_{0}$ and parameters $\mu$ of equation (35) for a model $m$ . Then, the similarity of an obtained pattern $\psi_{s}$ as a stationary state and a target pattern $\psi^{*}$ is evaluated by the cost function $E(\psi^{*},\psi_{s};\psi_{0},\mu)$ , describing the distance shown in equation (2) in the space of the order parameter $\Psi_{l}$ (see equation (20) or equation (21)). Changing $\mu$ , we may iterate numerical simulations and evaluation of the similarity between them. Following the Metropolis criterion at an inverse temperature $\beta$ , we compare a current cost function with a cost function in a previous step, and decide a current set of parameters is accepted or not.

By using the replica-exchange Monte Carlo (REMC) method, we sample $\psi_{0}$ and $\mu$ from $p(\psi_{0},\mu\mid\{\Psi_{l}^{*}\}_{l=1}^{l_{0}},\beta_{\alpha},m)$ for $N_{\rm rep}$ replicas in parallel Geyer:1991 ; Hukushima:1996 . At higher temperature (smaller $\beta$ ), the motion of one MC step in the parameter space is large, whereas at the lower temperature (larger $\beta$ ), each motion is small so that it intensively samples parameters near the minimum of the cost function. For every two steps, the parameter sets of neighbouring $\beta$ were exchanged following the Metropolis criterion. This process enables us to sample parameters weighed with likelihood effectively Geyer:1991 ; Hukushima:1996 .

The initial parameter set is sampled from the prior distribution of equation (36). The lowest cost function is typically achieved by $1000-2000$ MC steps. In one MC step, the Gibbs sampling is used to perform motion in the parameter space in all directions one by one. After finding the lowest cost function, we restart REMC from the initial parameter set of the lowest energy state to sample its steady state. This is because the relaxation under smaller $\beta$ is much faster than the larger $\beta$ . After $1000$ MC steps, we cut the initial burn-in steps and compute the statistical quantities after $200$ MC steps. Bridge sampling was used to calculate $F(\beta,m)$ for each $m$ meng1996simulating ; gelman1998simulating . The error bars of $F(\beta,m)$ were calculated by the bootstrap resampling efron1992bootstrap .

Regression method for noisy data

In the BM-PDE, the estimation from data with noise is performed by adding zero-mean Gaussian noise $\xi$ in the target pattern, namely $\psi^{*}({\bf x})\rightarrow\psi^{*}({\bf x})+\xi$ . In addition to this estimation using BM-PDE, another parameter estimation method is tested. We performed parameter estimation using the regression method in which the following cost function was used

\displaystyle E

\displaystyle=\frac{1}{2}\int\left[\psi^{*}-f_{\mu}^{dt}(\psi^{*})\right]^{2}d{\bf x}.

(38)

Here, the target pattern $\psi^{*}$ is numerically forward in $dt$ by the model of equation (5) with a parameter set $\{\mu\}$ as $f_{\mu}^{dt}(\psi(t))=\psi(t+dt)$ . If the target pattern is the stationary solution of the model, namely, if the parameters are ground truth to obtain the target pattern, the cost function must be zero. This approach is philosophically the same as the regression method in previous studies in which the cost function is a difference between the left-hand side (time derivative) and the right-hand side (force to change $\psi$ ), that is, $E=(1/2)\|\partial_{t}\psi-f_{\mu}(\psi)\|$ , under an appropriate norm $\|\cdot\|$ Brunton:2016 ; Schaeffer:2017 . The norm is typically chosen as the $L^{2}$ norm

\displaystyle E

\displaystyle=\frac{1}{2}\int\left[\partial_{t}\psi-{\color[rgb]{0,0,0}f_{\mu}(\psi)}\right]^{2}d{\bf x}.

(39)

In the current system, our model is no longer linear in the parameters, and therefore, we cannot use linear regression (including conventional sparse regression). In order to carry out nonlinear regression, we used the REMC method to minimise the cost function, equation (38). The method is similar to our main algorithm to sample parameters and to estimate the optimal noise by temperature $\beta$ . Following Bayes’ theorem, we estimate the best parameters by the sampled values and their error by the standard deviation.

Code availability

All codes used in this work are written in MATLAB and freely available from the corresponding author upon a request.

Data availability

The data that support the findings of this study are availablefrom the corresponding author on reasonable request.

Acknowledgements.

The authors are grateful to Edgar Knobloch, Yasumasa Nishiura, An-Chuang Shi, and Philippe Marcq for helpful discussions. The authors acknowledges the support by JSPS KAKENHI grant numbers 17K05605, 20H05259, and 20K03874 to N.Y. and 20K19889 to S.T. Numerical simulations in this work were carried out in part by AI Bridging Cloud Infrastructure (ABCI) at National Institute of Advanced Industrial Science and Technology (AIST).

Author contribution

N.Y. and S.T. conceived the research. N.Y. carried out simulations. N.Y. and S.T. analysed the results. N.Y. and S.T. wrote the manuscript.

Materials & Correspondence

Correspondence and requests for materials should be addressed to N.Y.

Competing Interests statement

The authors declare no competing financial interests.

References

(1) Swift, J. & Hohenberg, P. C. Hydrodynamic fluctuations at the convective instability. Phys. Rev. A 15, 319–328 (1977).
(2) Bates, F. S. & Fredrickson, G. H. Block copolymer thermodynamics: Theory and experiment. Annu. Rev. Phys. Chem. 41, 525–557 (1990).
(3) Harrison, C. et al. Mechanisms of ordering in striped patterns. Science 290, 1558–1560 (2000).
(4) Elder, K. R., Katakowski, M., Haataja, M. & Grant, M. Modeling elasticity in crystal growth. Phys. Rev. Lett. 88, 245701 (2002).
(5) Provatas, N. & Elder, K. Phase-field methods in materials science and engineering (John Wiley & Sons, 2011).
(6) Cross, M. & Hohenberg, P. Pattern formation outside of equilibrium. Rev. Mod. Phys. 65, 851–1112 (1993).
(7) Frank, F. C. & Kasper, J. S. Complex alloy structures regarded as sphere packings. II. Analysis and classification of representative structures. Acta Crystallographica 12, 483–499 (1959).
(8) Shechtman, D., Blech, I., Gratias, D. & Cahn, J. W. Metallic phase with long-range orientational order and no translational symmetry. Phys. Rev. Lett. 53, 1951–1953 (1984).
(9) Bates, M. W. et al. Stability of the a15 phase in diblock copolymer melts. Proc. Nat. Acad. Sci. USA 116, 13194–13199 (2019).
(10) Scherer, M. R. J. Gyroid and Gyroid-Like Surfaces, 7–19 (Springer International Publishing, Heidelberg, 2013).
(11) Yue, K. et al. Geometry induced sequence of nanoscale frank–kasper and quasicrystal mesophases in giant surfactants. Proc. Nat. Acad. Sci. USA 113, 14195–14200 (2016).
(12) Zeng, X. et al. Supramolecular dendritic liquid quasicrystals. Nature 428, 157–160 (2004).
(13) Hynninen, A.-P., Thijssen, J. H. J., Vermolen, E. C. M., Dijkstra, M. & van Blaaderen, A. Self-assembly route for photonic crystals with a bandgap in the visible region. Nat. Mater. 6, 202–205 (2007).
(14) Lifshitz, R. & Petrich, D. M. Theoretical model for faraday waves with multiple-frequency forcing. Phys. Rev. Lett. 79, 1261–1264 (1997).
(15) Kalinin, S. V., Sumpter, B. G. & Archibald, R. K. Big-deep-smart data in imaging for guiding materials design. Nat. Mater. 14, 973 (2015).
(16) Kennedy, M. C. & O’Hagan, A. Bayesian calibration of computer models. J. Royal Stat. Soc. B 63, 425–464 (2001).
(17) Smith, R. C. Uncertainty quantification: theory, implementation, and applications, vol. 12 (SIAM, 2013).
(18) Daniels, B. C. & Nemenman, I. Automated adaptive inference of phenomenological dynamical models. Nat. Commun. 6, 8133– (2015).
(19) Brunton, S. L., Proctor, J. L. & Kutz, J. N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Nat. Acad. Sci. USA 113, 3932–3937 (2016).
(20) Brunton, S. L. & Kutz, J. N. Data-driven Science and Engineering: Machine Learning, Dynamical Systems, and Control (Cambridge University Press, 2019).
(21) Voss, H. U., Timmer, J. & Kurths, J. Nonlinear dynamical system identification from uncertain and indirect measurements. Int. J. Bifurcat. Chaos 14, 1905–1933 (2004).
(22) Müller, T. G. & Timmer, J. Parameter identification techniques for partial differential equations. Int. J. Bifurcat. Chaos 14, 2053–2060 (2004).
(23) Rudy, S. H., Brunton, S. L., Proctor, J. L. & Kutz, J. N. Data-driven discovery of partial differential equations. Sci. Adv. 3, e1602614 (2017).
(24) Zhao, H., Storey, B. D., Braatz, R. D. & Bazant, M. Z. Learning the physics of pattern formation from images. Phys. Rev. Lett. 124, 060201 (2020).
(25) Supplementary Information is available on ’http://www.wpi-aimr.tohoku.ac.jp/~yoshinaga/publication/Supplement20210922.pdf’.
(26) Evensen, G. Data assimilation: the ensemble Kalman filter (Springer, 2009).
(27) Law, K., Stuart, A. & Zygalakis, K. Data assimilation (Springer, 2015).
(28) MacKay, D. J. Bayesian interpolation. Neural comput. 4, 415–447 (1992).
(29) Bishop, C. M. Pattern recognition and machine learning I, vol. 1 (springer New York, 2006).
(30) Tokuda, S., Nagata, K. & Okada, M. Simultaneous estimation of noise variance and number of peaks in bayesian spectral deconvolution. J. Phys. Soc. Japan 86, 024001 (2017).
(31) Jayaraman, A. & Mahanthappa, M. K. Counterion-dependent access to low-symmetry lyotropic sphere packings of ionic surfactant micelles. Langmuir 34, 2290–2301 (2018).
(32) Geyer, C. J. Markov chain monte carlo maximumlikelihood. Computing science and statistics: Proceedings of the 23rd Symposium on the Interface, American Statistical Association 156 (1991).
(33) Hukushima, K. & Nemoto, K. Exchange monte carlo method and application to spin glass simulations. J. Phys. Soc. Japan 65, 1604–1608 (1996).
(34) Elder, K. R. & Grant, M. Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals. Phys. Rev. E 70, 051605 (2004).
(35) Podneks, V. E. & Hamley, I. W. Landau- brazovskii theory for the ia3̄d structure. J. Exp. Theor. Phys. Lett. 64, 617–624 (1996).
(36) Shi, A.-C. Nature of anisotropic fluctuation modes in ordered systems. J. Phys. Cond. Mat. 11, 10183 (1999).
(37) Zhang, P. & Zhang, X. An efficient numerical method of landau-brazovskii model. J. Comp. Phys. 227, 5859 – 5870 (2008).
(38) Nonomura, M. & Ohta, T. Kinetics of morphological transitions between mesophases. J. Phys. Cond. Mat. 13, 9089–9112 (2001).
(39) Jaatinen, A. & Ala-Nissila, T. Extended phase diagram of the three-dimensional phase field crystal model. J. Phys. Cond. Mat. 22, 205402 (2010).
(40) Lee, S., Leighton, C. & Bates, F. S. Sphericity and symmetry breaking in the formation of frank–kasper phases from one component materials. Proc. Nat. Acad. Sci. 111, 17723–17731 (2014).
(41) Subramanian, P., Archer, A. J., Knobloch, E. & Rucklidge, A. M. Three-dimensional icosahedral phase field quasicrystal. Phys. Rev. Lett. 117, 075501 (2016).
(42) Damasceno, P. F., Glotzer, S. C. & Engel, M. Non-close-packed three-dimensional quasicrystals. Journal of Physics: Condensed Matter 29, 234005 (2017).
(43) Iacovella, C. R., Keys, A. S. & Glotzer, S. C. Self-assembly of soft-matter quasicrystals and their approximants. Proc. Nat. Acad. Sci. 108, 20935–20940 (2011).
(44) Dzugutov, M. Formation of a dodecagonal quasicrystalline phase in a simple monatomic liquid. Phys. Rev. Lett. 70, 2924–2927 (1993).
(45) Dzugutov, M. Phason dynamics and atomic transport in an equilibrium dodecagonal quasi-crystal. Europhys. Lett. 31, 95–100 (1995).
(46) Pismen, L. M. Patterns and interfaces in dissipative dynamics (Springer, 2006).
(47) Cox, S. & Matthews, P. Exponential time differencing for stiff systems. J. Comp. Phys. 176, 430–455 (2002).
(48) Aranson, I. S. (ed.) Physical Models of Cell Motility (Springer, Cham, 2015).
(49) Geman, S. & Geman, D. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE PAMI 6, 721–741 (1984).
(50) von Schnering, H. G. & Nesper, R. Nodal surfaces of fourier series: Fundamental invariants of structured matter. Z. Phys. B 83, 407–412 (1991).
(51) Yamada, K., Nonomura, M. & Ohta, T. Kinetics of morphological transitions in microphase-separated diblock copolymers. Macromolecules 37, 5762–5777 (2004).
(52) Impéror-Clerc, M. Three-dimensional periodic complex structures in soft matter: investigation using scattering methods. Interface focus 2, 589–601 (2012).
(53) Efron, B. & Morris, C. Stein’s estimation rule and its competitors-an empirical bayes approach. J. Am. Stat. Assoc. 68, 117–130 (1973).
(54) Akaike, H. Likelihood and the bayes procedure. In Selected papers of Hirotugu Akaike, 309–332 (Springer, 1998).
(55) Meng, X.-L. & Wong, W. H. Simulating ratios of normalizing constants via a simple identity: a theoretical exploration. Statistica Sinica 831–860 (1996).
(56) Gelman, A. & Meng, X.-L. Simulating normalizing constants: From importance sampling to bridge sampling to path sampling. Statistical science 163–185 (1998).
(57) Efron, B. Bootstrap methods: another look at the jackknife. In Breakthroughs in statistics, 569–593 (Springer, 1992).
(58) Schaeffer, H. Learning partial differential equations via data discovery and sparse optimization. Proc. Royal Soc. A 473, 20160446 (2017).

Figures and Tables

Refer to caption — Figure 1: Schematic illustration of Bayesian modelling of partial differential equations (BM-PDE). (a,b) The example of the target pattern of a dodecagonal QC produced from a numerical result (a), and its Fourier transform (b). The colour bar indicates $[-1.5,3.5]$ in (a) and $[0,3000]$ in (b). (c) The space of patterns $\psi({\bf x})$ and order parameters $\bm{\Psi}[\psi({\bf x})]$ . The PDE is solved with the initial condition $\psi_{0}$ taken from random variables. For each trajectory, the translational position and orientation of the pattern $\psi_{s}$ would change even under exactly the same parameters and the same model. The order parameter identifies the two patterns by extracting symmetries of the pattern. The distance between two patterns is quantified by $E$ . (d-g) Stationary observation and initial conditions. A pattern may be stable, meta-stable, or unstable. The unstable pattern cannot be generated from an initial condition. The better model (d,f) has broader initial conditions that generate stable patterns, than the worse model (e,g). (h) For each model and each set of parameters, there is a stationary pattern $\psi_{s}$ . The cost function $E[\psi^{*},\psi_{s}]$ (energy) is calculated from the order parameters of the target and generated patterns. From the posterior distribution, the best parameters and their errors are estimated. The distribution of the cost function gives the log marginal likelihood of the model from which the selection of models can be made.