Stability of three-dimensional icosahedral quasicrystals in multi-component systems

TL;DR

This paper investigates the stability of 3D icosahedral quasicrystals in multi-component systems using a coupled-mode Swift-Hohenberg model and a projection method, demonstrating their stability and extending the understanding of quasicrystal stabilization mechanisms.

Contribution

It introduces a systematic computational approach to analyze quasicrystal stability in multi-component systems, extending the length-scale interaction mechanism.

Findings

3D icosahedral quasicrystals are stable in the model.

The projection method accurately computes free energies.

Multi-component interactions stabilize quasicrystals.

Abstract

The relative stability of three-dimensional icosahedral quasicrystals in multi-component systems has been investigated based on a coupled-mode Swift-Hohenberg model with two-length-scales. A recently developed projection method, which provides a unified numerical framework to study periodic crystals and quasicrystals, is used to compute free energies to high accuracy. Compared with traditional approaches, the advantage of the projection method has been also discussed detailedly. A rigorous and systematical computation demonstrates that three-dimensional icosahedral quasicrystal, two-dimensional decagonal quasicrystal are stable phases in such a simple multi-component coupled-mode Swift-Hohenberg model. The result extends the multiple length-scales interaction mechanism which can stabilize quasicrystals from single-component to multi-component systems.

Figures

Tables

Table 1. Table 1: The minimal integer L 𝐿 L for desired E D ​ A subscript 𝐸 𝐷 𝐴 E_{DA} . The size of computational domain is [ 0 , L ⋅ 2 ​ π ) 2 superscript 0 ⋅ 𝐿 2 𝜋 2 [0,L\cdot 2\pi)^{2} in CAM.

$E_{DA}$	0.166879	0.091809	0.067405	0.055535	0.037427
$L$	126	204	1288	2084	3372

Table 2. Table 2: A comparison of the convergent free energy calculated by PM and CAM with different numbers of discrete points N 𝑁 N . The model parameters set as c = 80 , τ = 0 , t = 0 , t 0 = − 0.3 , g 0 = − 0.3 , g 1 = 2.2 , g 2 = 2.2 , d 0 = d 1 = d 2 = d 3 = 0 formulae-sequence 𝑐 80 formulae-sequence 𝜏 0 formulae-sequence 𝑡 0 formulae-sequence subscript 𝑡 0 0.3 formulae-sequence subscript 𝑔 0 0.3 formulae-sequence subscript 𝑔 1 2.2 formulae-sequence subscript 𝑔 2 2.2 subscript 𝑑 0 subscript 𝑑 1 subscript 𝑑 2 subscript 𝑑 3 0 c=80,\tau=0,t=0,t_{0}=-0.3,g_{0}=-0.3,g_{1}=2.2,g_{2}=2.2,d_{0}=d_{1}=d_{2}=d_{3}=0 .

$N$	PM	CAM($L=126$)
10	-4.296833632029e-02	-4.287269226697e-02
12	-4.296833659056e-02	-4.287269226627e-02
14	-4.296833659212e-02	-4.287269226626e-02
16	-4.296833659224e-02	-4.287269226626e-02
18	-4.296833659225e-02	-4.287269226626e-02
20	-4.296833659225e-02	-4.287269226626e-02
22	-4.296833659225e-02	-4.287269226626e-02

Equations

$$F[{\phi }_1,{\phi }_2,\cdots ,{\phi }_m]$$

$$=$$

$$\displaystyle\tau_{\sigma_{i_{1}}\sigma_{i_{2}}\cdots\sigma_{i_{n}}}[\sigma_{i_{1}}\varphi_{i_{1}}(\mathbf{r})][\sigma_{i_{2}}\varphi_{i_{2}}(\mathbf{r})]\cdots[\sigma_{i_{n}}\varphi_{i_{n}}(\mathbf{r})]\,d\mathbf{r}\Bigg{\}},$$

$$\displaystyle F[\psi,\phi]=\;\frac{1}{V}\bigg{\{}$$

$$+\int d\mathbf{r}(\tau {\psi }^2+g_0{\psi }^3+{\psi }^4+t{\varphi }^2+t_0{\varphi }^3$$

$$+{\varphi }^4-g_1{\psi }^2\varphi -g_2\psi {\varphi }^2+d_0\psi \varphi$$

$$\displaystyle+d_{1}\psi^{2}\phi^{2}+d_{2}\psi^{3}\phi+d_{3}\psi\phi^{3})\bigg{\}}.$$

$$\frac{\delta F}{\delta\psi(\mathbf{r})}\bigg{|}_{\psi^{*}}=\;0,~{}~{}~{}~{}~{}~{}~{}~{}\frac{\delta F}{\delta\phi(\mathbf{r})}\bigg{|}_{\phi^{*}}=\;0.$$

$$\left\{\;\begin{aligned} &\begin{aligned} \frac{\partial\psi}{\partial\alpha}=\;-\frac{\delta F}{\delta\psi}=\;&-c(\nabla^{2}+1)^{2}\psi-2\tau\psi-3g_{0}\psi^{2}\\ &-4\psi^{3}+2g_{1}\psi\phi+g_{2}\phi^{2}{\color[rgb]{0,0,0}-d_{0}\phi}\\ &{\color[rgb]{0,0,0}-2d_{1}\psi\phi^{2}-3d_{2}\psi^{2}\phi-d_{3}\phi^{3}},\end{aligned}\\ &\begin{aligned} \frac{\partial\phi}{\partial\alpha}=\;-\frac{\delta F}{\delta\phi}=\;&-c(\nabla^{2}+q^{2})^{2}\phi-2t\phi-3t_{0}\phi^{2}\\ &-4\phi^{3}+2g_{2}\psi\phi+g_{1}\psi^{2}{\color[rgb]{0,0,0}-d_{0}\psi}\\ &{\color[rgb]{0,0,0}-2d_{1}\psi^{2}\phi-d_{2}\psi^{3}-3d_{3}\psi\phi^{2}}.\end{aligned}\end{aligned}\right.$$

$$\psi (\mathbf{r})=\sum _{\mathbf{h}\in {\mathbb{Z}}^n}\hat{\psi }(\mathbf{h})e^{i[(\mathcal{P}\cdot \mathbf{Bh}{)}^T\cdot \mathbf{r}]},\mathbf{r}\in {\mathbb{R}}^d,$$

$$\psi (\mathbf{r})=\sum _{\mathbf{k}}\hat{\psi }(\mathbf{k})e^{i{\mathbf{k}}^T\cdot \mathbf{r}},\mathbf{r}\in {\mathbb{R}}^d,$$

$$\varphi (\mathbf{r})=\sum _{\mathbf{k}}\hat{\varphi }(\mathbf{k})e^{i{\mathbf{k}}^T\cdot \mathbf{r}},\mathbf{r}\in {\mathbb{R}}^d,$$

$$\left\{\;\begin{aligned} &\begin{aligned} \frac{\partial\widehat{\psi}(\mathbf{k})}{\partial\alpha}=&-c(-|\mathbf{k}|^{2}+1)^{2}\widehat{\psi}(\mathbf{k})-2\tau\widehat{\psi}(\mathbf{k})-3g_{0}\widehat{\psi^{2}}(\mathbf{k})\\ &-4\widehat{\psi^{3}}(\mathbf{k})+2g_{1}\widehat{(\psi\phi)}(\mathbf{k})+g_{2}\widehat{\phi^{2}}(\mathbf{k}){\color[rgb]{0,0,0}-d_{0}\widehat{\phi}(\mathbf{k})}\\ &{\color[rgb]{0,0,0}-2d_{1}\widehat{(\psi\phi^{2})}(\mathbf{k})-3d_{2}\widehat{(\psi^{2}\phi)}(\mathbf{k})-d_{3}\widehat{\phi^{3}}(\mathbf{k})},\end{aligned}\\ &\begin{aligned} \frac{\partial\widehat{\phi}(\mathbf{k})}{\partial\alpha}=&-c(-|\mathbf{k}|^{2}+q^{2})^{2}\widehat{\phi}(\mathbf{k})-2t\widehat{\phi}(\mathbf{k})-3t_{0}\widehat{\phi^{2}}(\mathbf{k})\\ &-4\widehat{\phi^{3}}(\mathbf{k})+2g_{2}\widehat{(\psi\phi)}(\mathbf{k})+g_{1}\widehat{\psi^{2}}(\mathbf{k}){\color[rgb]{0,0,0}-d_{0}\widehat{\psi}(\mathbf{k})}\\ &{\color[rgb]{0,0,0}-2d_{1}\widehat{(\psi^{2}\phi)}(\mathbf{k})-d_{2}\widehat{\psi^{3}}(\mathbf{k})-3d_{3}\widehat{(\psi\phi^{2})}(\mathbf{k})}.\end{aligned}\end{aligned}\right.$$

$$\left\{\;\begin{aligned} &\widehat{\psi^{2}}(\mathbf{k})=\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=1\\ \mathbf{k}_{1}+\mathbf{k}_{2}=\mathbf{k}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\psi}(\mathbf{k}_{2}),\\ &\widehat{\phi^{2}}(\mathbf{k})=\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}=\mathbf{k}\end{array}$}}\widehat{\phi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2}),\\ &\widehat{\psi^{3}}(\mathbf{k})=\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=1\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=\mathbf{k}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\psi}(\mathbf{k}_{2})\widehat{\psi}(\mathbf{k}_{3}),\\ &\widehat{\phi^{3}}(\mathbf{k})=\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=\mathbf{k}\end{array}$}}\widehat{\phi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3}),\\ &\widehat{(\psi\phi)}(\mathbf{k})=\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=1,|\mathbf{k}_{2}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}=\mathbf{k}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})\\ &{\color[rgb]{0,0,0}\widehat{(\psi\phi^{2})}(\mathbf{k})=\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=1,|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=\mathbf{k}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3})}\\ &{\color[rgb]{0,0,0}\widehat{(\psi^{2}\phi)}(\mathbf{k})=\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=1,|\mathbf{k}_{3}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=\mathbf{k}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\psi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3})}.\end{aligned}\right.$$

$$\left\{\;\begin{aligned} &\begin{aligned} (\widehat{\psi})_{\text{BDF2}}=\;&-c(-|\mathbf{k}|^{2}+1)^{2}\widehat{\psi}_{\alpha+\delta\alpha}-2\tau\widehat{\psi}_{\alpha+\delta\alpha}\\ &-3g_{0}\widehat{\bar{\psi}^{2}}_{\alpha+\delta\alpha}-4\widehat{\bar{\psi}^{3}}_{\alpha+\delta\alpha}+2g_{1}\widehat{(\bar{\psi}\bar{\phi})}_{\alpha+\delta\alpha}\\ &+g_{2}\widehat{\bar{\phi}^{2}}_{\alpha+\delta\alpha}{\color[rgb]{0,0,0}-d_{0}\widehat{\bar{\phi}}_{\alpha+\delta\alpha}-2d_{1}\widehat{(\bar{\psi}\bar{\phi}^{2})}_{\alpha+\delta\alpha}}\\ &{\color[rgb]{0,0,0}-3d_{2}\widehat{(\bar{\psi}^{2}\bar{\phi})}_{\alpha+\delta\alpha}-d_{3}\widehat{\bar{\phi}^{3}}_{\alpha+\delta\alpha}},\end{aligned}\\ &\begin{aligned} (\widehat{\phi})_{\text{BDF2}}=\;&-c(-|\mathbf{k}|^{2}+q^{2})^{2}\widehat{\phi}_{\alpha+\delta\alpha}-2t\widehat{\phi}_{\alpha+\delta\alpha}\\ &-3t_{0}\widehat{\bar{\phi}^{2}}_{\alpha+\delta\alpha}-4\widehat{\bar{\phi}^{3}}_{\alpha+\delta\alpha}+2g_{2}\widehat{(\bar{\psi}\bar{\phi})}_{\alpha+\delta\alpha}\\ &+g_{1}\widehat{\bar{\psi}^{2}}_{\alpha+\delta\alpha}{\color[rgb]{0,0,0}-d_{0}\widehat{\bar{\psi}}_{\alpha+\delta\alpha}-2d_{1}\widehat{(\bar{\psi}^{2}\bar{\phi})}_{\alpha+\delta\alpha}}\\ &{\color[rgb]{0,0,0}-d_{2}\widehat{\bar{\psi}^{3}}_{\alpha+\delta\alpha}-3d_{3}\widehat{(\bar{\psi}\bar{\phi}^{2})}_{\alpha+\delta\alpha}},\end{aligned}\\ &(\bar{\psi}\bar{\phi})_{\alpha+\delta\alpha}=\bar{\psi}_{\alpha+\delta\alpha}\bar{\phi}_{\alpha+\delta\alpha},\\ &(\bar{\psi}^{2}\bar{\phi})_{\alpha+\delta\alpha}=\left(\bar{\psi}_{\alpha+\delta\alpha}\right)^{2}\bar{\phi}_{\alpha+\delta\alpha},\\ &(\bar{\psi}\bar{\phi}^{2})_{\alpha+\delta\alpha}=\bar{\psi}_{\alpha+\delta\alpha}\left(\bar{\phi}_{\alpha+\delta\alpha}\right)^{2},\end{aligned}\right.$$

$$F$$

$$\displaystyle+g_{0}\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=1\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=\mathbf{0}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\psi}(\mathbf{k}_{2})\widehat{\psi}(\mathbf{k}_{3})$$

$$\displaystyle+\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=|\mathbf{k}_{4}|=1\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}+\mathbf{k}_{4}=\mathbf{0}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\psi}(\mathbf{k}_{2})\widehat{\psi}(\mathbf{k}_{3})\widehat{\psi}(\mathbf{k}_{4})$$

$$\displaystyle+t\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}=\mathbf{0}\end{array}$}}\widehat{\phi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})$$

$$\displaystyle+t_{0}\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=\mathbf{0}\end{array}$}}\widehat{\phi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3})$$

$$\displaystyle+\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=|\mathbf{k}_{4}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}+\mathbf{k}_{4}=\mathbf{0}\end{array}$}}\widehat{\phi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3})\widehat{\phi}(\mathbf{k}_{4})$$

$$\displaystyle-g_{1}\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=1,|\mathbf{k}_{3}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=\mathbf{0}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\psi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3})$$

$$\displaystyle-g_{2}\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=1,|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=\mathbf{0}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3})$$

$$\displaystyle{\color[rgb]{0,0,0}+d_{0}\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=1,|\mathbf{k}_{2}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}=\mathbf{0}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})}$$

$$\displaystyle{\color[rgb]{0,0,0}+d_{1}\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=1,|\mathbf{k}_{3}|=|\mathbf{k}_{4}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}+\mathbf{k}_{4}=\mathbf{0}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\psi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3})\widehat{\phi}(\mathbf{k}_{4})}$$

$$\displaystyle{\color[rgb]{0,0,0}+d_{2}\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=1,|\mathbf{k}_{4}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}+\mathbf{k}_{4}=\mathbf{0}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\psi}(\mathbf{k}_{2})\widehat{\psi}(\mathbf{k}_{3})\widehat{\phi}(\mathbf{k}_{4})}$$

$$\displaystyle{\color[rgb]{0,0,0}+d_{3}\sum_{\mbox{\tiny$\begin{array}[]{c}|\mathbf{k}_{1}|=1,|\mathbf{k}_{2}|=|\mathbf{k}_{3}|=|\mathbf{k}_{4}|=q\\ \mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}+\mathbf{k}_{4}=\mathbf{0}\end{array}$}}\widehat{\psi}(\mathbf{k}_{1})\widehat{\phi}(\mathbf{k}_{2})\widehat{\phi}(\mathbf{k}_{3})\widehat{\phi}(\mathbf{k}_{4})}.$$

$$F_{IQC,L}=$$

$$+120t_0{\hat{\varphi }}^3+3330{\hat{\varphi }}^4-120g_1{\hat{\psi }}^2\hat{\varphi }$$

$$\displaystyle-120g_{2}\widehat{\psi}\widehat{\phi}^{2}{\color[rgb]{0,0,0}+2100d_{1}\widehat{\psi}^{2}\widehat{\phi}^{2}}$$

$$\displaystyle{\color[rgb]{0,0,0}+1440d_{2}\widehat{\psi}^{3}\widehat{\phi}+1440d_{3}\widehat{\psi}\widehat{\phi}^{3}},$$

$$F_{DQC,L}=$$

$$\displaystyle-20g_{1}\widehat{\psi}^{2}\widehat{\phi}-20g_{2}\widehat{\psi}\widehat{\phi}^{2}{\color[rgb]{0,0,0}+140d_{1}\widehat{\psi}^{2}\widehat{\phi}^{2}}$$

$$\displaystyle{\color[rgb]{0,0,0}+60d_{2}\widehat{\psi}^{3}\widehat{\phi}+60d_{3}\widehat{\psi}\widehat{\phi}^{3}},$$

$$F_{L+B[1],L}=$$

$$\displaystyle{\color[rgb]{0,0,0}+8d_{1}\widehat{\psi}^{2}\widehat{\phi}^{2}},$$

$$F_{L+B[2],L}=$$

$$\displaystyle{\color[rgb]{0,0,0}+8d_{1}\widehat{\psi}^{2}\widehat{\phi}^{2}},$$

$$F_{L2,L}=2\tau {\hat{\psi }}^2+6{\hat{\psi }}^4,$$

$$F_{HEX,L}=6{\hat{\psi }}^2+90{\hat{\psi }}^4,$$

$$\mathbf{k}=p_1{\mathbf{e}}_1^{\ast }+p_2{\mathbf{e}}_2^{\ast }+\cdots +p_d{\mathbf{e}}_d^{\ast }.$$

$$\varphi (\mathbf{r})=\sum _{\mathbf{k}}\hat{\varphi }(\mathbf{k})e^{i\cdot (L\mathbf{k})\cdot \mathbf{r}/L},\mathbf{r}\in [0,2\pi L{)}^d.$$

Full text

Keywords

Icosahedral quasicrystals, Two-length-scales, Coupled-mode Swift-Hohenberg model,

Phase diagrams, Projection method

PACS

61.44.Br, 61.44.Fw, 64.60.Bd, 02.70.Hm

DOI

I INTRODUCTION

Quasicrystals Levine1984 are a class of important ordered materials possessing quasiperiodic positional order and long-range orientational order between periodic crystals and amorphous materials. The mathematical description of quasicrystals can track back to the work done independently by Meyer Meyer1972 and Penrose penrose1974role in early 1970. The first quasicrystal, actually a three dimensional (3D) icosahedral quasicrystal (IQC), was discovered by Shechtman in a rapidly-quenched Al-Mn alloy until 1982 Shechtman1984 . Since the first discovery of quasicrystals, the quasiperiodic long-range order has been found in a large number of metallic alloys Steurer2009 ; Steurer2004 ; Tsai2008 , and also in a host of soft-matter systems Zeng2004 ; Zhang2012 . Among these discoveries, icosahedral symmetric quasicrystals are the most frequently found, concretely in more than one hundred different metallic alloys Tsai2003 ; Tsai2008 . However, the thermodynamic stability of quasicrystals in multi-component systems, including IQCs, remains a challenge Henley2006 ; Lifshitz2007 , mainly due to the lack of appropriate theoretical models and high-precision numerical methods.

Theoretical approaches to investigating the stability of an ordered phase, including periodic and quasiperiodic one, often involve minimizing an appropriate free energy functional of the system and comparing the free energy of different candidate structures. A systematic examination of the relative stability of quasicrystals requires the availability of suitable free energy functionals and accurate methods to evaluate the free energy of phases with quasicrystalline orders. Nowadays, various coarse-grained free energy functionals have been proposed to investigate the phase behavior of physical systems, especially the thermodynamic stability of ordered phases including quasicrystals and periodic crystals Alexander1978 ; Bak1985 ; Jaric1985 ; Kalugin1985 ; Gronlund1988 ; Lifshitz1997 ; Subramanian2016 ; Jiang2017 . Some of these works show the possibility of stable IQCs based on Landau-type free energy functionals with one order parameter Subramanian2016 ; Jiang2017 . Inspired by the Alexander-McTague theory Alexander1978 , Mermin and Troian Mermin1985 introduced a second order parameter to obtain stable IQCs. On the other hand, Swift and Hohenberg Swift1977 explicitly added a positive-definite gradient term into the free energy functional to represent the effect of characteristic length-scales. This idea was rapidly extended to more than one characteristic length-scale and widely utilized to explore the quasicrystalline order Mermin1985 ; Mueller1994 ; Lifshitz1997 ; Dotera2007 ; Archer2013 ; Jiang2016 . In particular, Dotera Dotera2007 extended the Mermin-Troian model to study an ABC star copolymer system with an incompressible condition. His work shows the emergence and stability of two-dimensional (2D) decagonal quasicrystals (DQCs). Recently, Jiang et al. Jiang2016 have constructed a two-component coupled-mode Swift-Hohenberg (CMSH) model with two length-scales to explore the thermodynamic stability of the periodic and quasiperiodic structures. According to their work, the 2D decagonal and dodecagonal quasicrystals are both stable when the ratio of two-length-scales is appropriately chosen. In the present work, we will study the relative stability of 3D IQCs in multi-component systems based on the CMSH model.

Besides a proper free energy functional to describe multi-component systems, examining the thermodynamic stability of quasicrystals requires accurate and efficient methods to evaluate the free energy of various ordered phases. Due to the lack of translational symmetry, the computation of quasicrystals is harder to carry out within a finite domain as done for periodic crystals. In the literature, utilizing a large periodic structure to approximate a quasicrystal is a commonly used method in the study of the quasicrystalline order Lifshitz1997 ; Subramanian2016 ; Quandt1999 ; Skibinsky1999 ; Dotera2006 . The method actually obtains a crystalline approximant, therefore, it is named the crystalline approximant method (CAM). From the viewpoint of numerical computation, CAM is based on the approximation of irrational numbers by integers or rational numbers, corresponding to well-known Diophantine approximation (DA) problem in the number theory Meyer1972 ; Jiang2014 . Because of the existence of DA, CAM has to be implemented in a very large computational region which means an unacceptable computational amount if small DA error is required. Furthermore, it has been verified that the gap between the free energy of quasicrystals and their corresponding approximants cannot be vanished in any finite computational region Meyer1972 . In order to avoid the approximation error, an alternative approach is proposed to calculate quasicrystals based on the fact that quasiperiodic lattices can be generated by a cut-and-project method from higher-dimensional periodic lattices Meyer1972 . This method provides a basic framework to study quasicrystals, originally proposed by Meyer in studying the relationship between harmonic analysis and algebraic numbers Meyer1972 . More recently, Jiang and Zhang proposed a projection method (PM) to obtain the density profile of quasicrystals and evaluate their energy density on high accuracy. The PM shows that the Fourier spectrum of a quasiperiodic structure can be embedded into a higher-dimensional crystallographic point packing set of corresponding periodic structure. Consequently, quasicrystals can be computed precisely in a higher-dimensional space and then be recovered by projecting the higher-dimensional reciprocal lattice vectors back to the original Fourier space through a projection matrix. This method can avoid the DA error effectively. As a particular case, it can be further used to investigate periodic crystals by setting the projection matrix as an identity matrix. From this perspective, the PM can calculate free energy of quasicrystals and periodic crystals with the same accuracy. Therefore, the PM provides a unified computational framework to study the relative stability of quasicrystals and periodic crystals.

In the present work, we will mainly explore the formation and the thermodynamic stability of quasicrystals and periodic crystals in multi-component systems based on the CMSH model. A careful comparison of the free energy of abundant possible ordered phases, including 3D IQCs, 2D DQCs and periodic crystals, leads to confirm the relative stability and construct phase diagrams. Theoretical results predict that IQCs can emerge and occupy a thermodynamic stable region in the phase diagram within the multi-component CMSH model.

II THEORETICAL FRAMEWORK

Since the discovery of quasicrystals, a large number of theoretical studies have been carried out to investigate their symmetry, structure characterization physical properties and thermodynamic stability. Current theoretical results demonstrate that the formation and stability of the quasicrystalline phases may be characterized by two or more length-scales interaction potential Mermin1985 ; Mueller1994 ; Lifshitz1997 ; Dotera2007 ; Archer2013 ; Jiang2016 ; Subramanian2016 . Furthermore, since the formation of quasicrystals in multi-component systems could be hardly described by only one order parameter, more order parameters should be considered in a coarse-grained free energy functional. Combining with the above factors, the CMSH model Jiang2016 was proposed to exploit the formation and the thermodynamic stability of quasicrystals. In order to examine the relative stability of different ordered phases, the free energy of candidate structures, corresponding to local minima of the free energy functional, should be calculated precisely. With an appropriate representation of quasicrystals, the PM provides a unified numerical framework to study periodic and quasiperiodic structures with the same precision Jiang2014 . The rest of this section will give a brief introduction of the CMSH model and the PM for quasicrystals and periodic crystals.

II.1 Coupled-mode Swift-Hohenberg model

The original CMSH model only considers two-component systems with two-length-scales Jiang2016 . In this paper, we will give a general framework of multi-component systems with multi-length-scales. To be more concrete, the free energy functional of the CMSH model for a $m$-component system can be written as

[TABLE]

where $\varphi_{j}$ is the $j$-th order parameter. $V$ represents the system volume. $\sigma_{i_{j}}\in\{0,1\}$, $j=1,2,\cdots,n$. $q_{j}>0$ denotes the $j$-th characteristic length scale. $c>0$ is an energy penalty factor to ensure $m$ characteristic length-scales. The cross-terms in Eq.(1) demonstrate the $n$-interactions among $m$ components. $\tau_{\sigma_{i_{1}}\sigma_{i_{2}}\cdots\sigma_{i_{n}}}$ is interaction intensity related to physical conditions, such as temperature, pressure, and physical or chemical properties of materials. Comparing with the common Landau theory of phase transition, the crucial feature of CMSH model is the occurrence of multi-length-scales $q_{j}$, $j=1,2,\cdots,m$. It should be noted that the number of characteristic multi-length-scales could be unequal to $m$.

In current work, we consider a two-component system with two-length-scales. At the same time, the highest-degree of interactions is fourth. More specifically, the general CMSH model Eq.(1) is simplified by the following free energy functional,

[TABLE]

Order parameters $\psi(\mathbf{r})$ and $\phi(\mathbf{r})$ correspond to the density profile of two-component systems such as metallic alloys or soft matters. $\tau$, $g_{0}$, $t$, $t_{0}$, $g_{1}$, $g_{2}$, $d_{0}$, $d_{1}$, $d_{2}$ and $d_{3}$ are all interaction parameters. $q_{j}$, $j=1,2$ represents the characteristic length-scale. According to the theories of Alexander and MacTague Alexander1978 , and Lifshitz and Petrich Lifshitz1997 , the cubic terms, associated with triangle interactions, in the above free energy functional play a crucial role in stabilizing periodic and quasiperiodic structures since they can reduce the value of free energy. The ratio $q=q_{2}/q_{1}$ is related to the symmetry of ordered structures. It should be pointed out that one of the length-scales can be always taken as a unit wavelength (as specified by $q_{1}=1$) and the second length-scale is specified by $q_{2}=q$. The difference of the two-component CMSH model from a single-component incommensurate phase-field-crystal model, such as Lifshitz-Petrich model Lifshitz1997 , is the two-length-scales potential occurring in different order parameters. Moreover, in consideration of an incompressibility condition Dotera2007 ; Jiang2016 , the two-component CMSH model can be utilized to describe three-component systems, such as Al-Cu-Fe alloys Tsai1987 ; Ishimasa1988 , Zn-Mg-Sc alloys Kaneko2001 , Al-Rh-Si alloys Koshikawa2003 , Au-Al-Yb alloys Ishimasa2011 , $etc$. It should be remarked that the CMSH model is a phenomenological model. If the model parameters can match a given experimental system, then the CMSH model could be used to explain the corresponding physical system.

Theoretically, the ordered patterns including periodic and quasiperiodic structures correspond to local minima of the free energy functional with respect to order parameters $\psi$ and $\phi$ in the above energy functional. Accordingly, the order parameters $\psi^{*}$ and $\phi^{*}$ located in the equilibrium state are the minima of the free energy density functional, which means

[TABLE]

In order to find the equilibrium state, the Allen-Cahn dynamic equation is utilized to minimize the free energy functional and yields

[TABLE]

The variable $\alpha$ does not represent time but a parameter controlling the iteration steps.

II.2 Projection Method

Due to the spatial periodicity, the computation of periodic crystals can be carried out within a unit cell with periodic boundary conditions. On the contrary, this method is not suitable for quasicrystals because they are quasiperiodic in at least one direction. For the quasicrystalline order, an efficient method is *the projection method *(PM) based on the fact that a $d$-dimensional quasicrystal is a combination of a class of exponentials \big{\{}e^{i\mathbf{k}^{T}\mathbf{r}}\big{\}}, $\mathbf{r}\in\mathbb{R}^{d}$, $\mathbf{k}$ belongs to a countable set Jiang2014 . From the higher-dimensional description Steurer2009 , any $d$-dimensional quasicrystal can be embedded into an $n$-dimensional periodic structure ($n\geqslant d$). Specifically, the $n$-dimensional reciprocal vectors can be spanned by a set of primitive vectors $\mathbf{b}_{i}$, which are the primitive reciprocal vectors in the $n$-dimensional reciprocal space, with integer coefficients. It means any reciprocal vector of an $n$-dimensional periodic structure can be expressed as $\mathbf{H}=\mathbf{Bh}$, where $\mathbf{B}=(\mathbf{b}_{1},\dots,\mathbf{b}_{n})\in\mathbb{R}^{n\times n}$ is the primitive reciprocal lattice and $\mathbf{h}\in\mathbb{Z}^{n}$. Then the wave vector $\mathbf{k}$ in $d$-dimensional physical space is obtained by the projection, $\mathbf{k}=\mathcal{P}\cdot\mathbf{H}$, where $\mathcal{P}$ is a projection matrix of $d\times n$-order. The dimensionality $n$ is determined by the rotational symmetry of the quasicrystal. In particular, it can be obtained by an additive Euler function hiller1985crystallographic . In the view of numerical computation, the projection matrix can be directly confirmed through the representation of the primitive basis vectors Jiang2014 .

In detail, using the $n$-dimensional periodic structure and the projection matrix, any $d$-dimensional quasiperiodic function $\psi(\mathbf{r})$ can be expanded as

[TABLE]

where the Fourier coefficient $\widehat{\psi}(\mathbf{h})$ can be easily obtained by using the $n$-dimensional $L^{2}$-inner product, $\widehat{\psi}(\mathbf{h})=\left<\tilde{\psi}(\tilde{\mathbf{r}}),e^{i[(\mathcal{P}\cdot\mathbf{Bh})^{T}\cdot\tilde{\mathbf{r}}]}\right>$, with $\{\widehat{\psi}(\mathbf{h})\}_{\mathbf{h}\in\mathbb{Z}^{n}}\in\ell^{2}(\mathbb{Z}^{n})$ and $\tilde{\mathbf{r}}=\sum_{i=1}^{n}c_{i}\mathbf{a}_{i},c\in[0,1)$. $\mathbf{a}_{i}\in\mathbb{R}^{n}$, $i=1,2,\dots,n$, is the reciprocal primitive vector which satisfies the dual relationship, $\mathbf{a}_{i}\cdot\mathbf{b}_{j}=2\pi\delta_{ij}$. Furthermore, the function $\tilde{\psi}(\tilde{\mathbf{r}})$ is the inverse Fourier transform of the Fourier coefficient $\widehat{\psi}(\mathbf{h})$. From the expansion Eq.(5), the $d$-dimensional quasiperiodic structure can be also treated as a hyperplane of an $n$-dimensional periodic structure whose orientation is determined by the projection matrix $\mathcal{P}$. In order to describe the position of the quasilattice in $d$-dimensional Fourier space, the notation $\mathbf{k}$ is utilized to replace $\mathcal{P}\cdot\mathbf{Bh}$ in Eq.(5). With this notation, the projection method has the following form,

[TABLE]

where $\mathbf{k}=\sum_{i=1}^{n}h_{i}(\mathcal{P}\mathbf{b}_{i})\in\mathbb{R}^{d}$. Despite being similar to the common Fourier series, it should be emphasized that the distribution of $\mathbf{k}$ is not a periodic lattice. Similarly, the order parameter $\phi(\mathbf{r})$ in multi-component systems can be expanded as follows,

[TABLE]

where $\mathbf{k}$ has the same definition as Eq.(6).

For a given structure of interest, the reciprocal lattice vectors are determined by its symmetry, and the optimal coefficients are obtained by minimizing the free energy functional. Based on these expansions, we can compute quasicrystals in the high-dimensional space and then project it onto the lower-dimensional physical space. Since the reciprocal lattice in the $n$-dimensional space is periodic, the computation can be performed on a uniform mesh grid. As a particular case, a $d$-dimensional periodic structure can be described by the PM when the projection matrix is set as a $d$-order identity matrix. In this sense, the PM becomes the common Fourier-spectral method. Therefore, this perspective provides a unified computational framework of periodic and quasiperiodic structures with the same accuracy. More significantly, the PM can evaluate the free energy of the corresponding ordered structure to high accuracy.

In practice, inserting the generalized Fourier expansions Eqs.(6) and (7) into the dynamic equations (4), we can obtain the following iterative equations expressed by Fourier coefficients,

[TABLE]

In this expression, the quadratic, cubic and cross terms are given by,

[TABLE]

From these expressions, it is clear that the nonlinear (quadratic, cubic and cross) terms in Eq.(8) are $n$-dimensional convolutions in the reciprocal space. A direct evaluation of these nonlinear terms is extremely expensive. Instead, these terms are simple multiplication in the $n$-dimensional real space. The pseudospectral method takes advantage of this observation by evaluating the gradient terms in the Fourier space and the nonlinear terms in the real space by performing the efficient Fast Fourier Transformation (FFT) algorithm. Thus it provides an efficient technique to find the solutions of dynamic equations.

In order to solve the time-dependent equations (8) numerically, it is necessary to apply a time discretization scheme. In this work, we propose a second-order precision method by combining the second-order Adam-Bashforth with Lagrange extrapolation approach (BDF2-LE) to discretize the dynamic equations. In particular, the scheme can be written as

[TABLE]

where $(\cdot)_{\text{BDF2}}=\dfrac{3(\cdot)_{\alpha+\delta\alpha}-4(\cdot)_{\alpha}+(\cdot)_{\alpha-\delta\alpha}}{2\delta\alpha}$, and $\bar{(\cdot)}_{\alpha+\delta\alpha}$ can be calculated by the Lagrange extrapolation approach $\bar{(\cdot)}_{\alpha+\delta\alpha}=2(\cdot)_{\alpha}-(\cdot)_{\alpha-\delta\alpha}$.

II.3 Two Modes Approximation Method

From the observation of the CMSH model (2), it is clear that the penalty factor $c$ affects the emergence of nonzero Fourier vectors. In order to systematically investigate the ideal and actual phase behavior, we consider the hard constraint ($c\to+\infty$) and the soft constraint (finite $c$). Under the hard constraint, $\widehat{\psi}$ (resp. $\widehat{\phi}$) should be strictly restricted on the circle with radius $1$ (resp. $q$), otherwise, the value of the free energy functional (2) will be infinite. Therefore, it is only required to analyze the entropy part. This is the key idea of the two modes approximation method (TMAM). In particular, under the hard constraint, the expansion terms of Eqs.(6) and (7) are finite for any given symmetric structure. Inserting Eqs.(6) and (7) into the CMSH model (2), we obtain the following expression with respect to Fourier coefficients $\widehat{\psi}(\mathbf{k})$ and $\widehat{\phi}(\mathbf{k})$,

[TABLE]

It should be noted that the summations are all finite in the above expression. Through the TMAM, the free energy functional is turned into a function with a finite number of variables. The energy expressions of all involved ordered structures can be found in Sec. III.3.

III RESULTS AND DISCUSSION

In this section, all experiments were performed on a desktop computer with a 3.20 GHz CPU (i5-6500, 4 processors). All codes were written in MATLAB without parallelization.

III.1 The advantages of PM

High accurate numerical methods are crucial to theoretically investigating the relative stability of ordered structures. The CAM is a widely used approach to examine quasiperiodic structures, however, it inevitably accompanies the DA error, $E_{DA}$ Meyer1972 ; Jiang2014 , which dominates the computational precision. In order to avoid $E_{DA}$, the PM is a reliable approach Jiang2014 . In this subsection, we will demonstrate that the PM has advantages over the CAM. The main idea of CAM is using a large periodic structure to approximate a quasicrystal. For a desired computational precision, CAM always requires a large computational region to reduce $E_{DA}$. More details about CAM can be found in Appendix A. Recently, based on the observation that a quasiperiodic phase can be embedded in a higher-dimensional periodic structure, an efficient method, i.e., PM, has been formulated in the Fourier space. Comparing with CAM, PM not only requires less computational cost, but also computes real quasicrystalline structures, as well as their free energy to high accuracy.

In the following, we will take the 2D DQC as an example to demonstrate the correctness and efficiency of PM in detail, the morphology of 2D DQC can be found in FIG. 5 (a). The 2D DQC can be embedded into a $4$-dimensional periodic structure. Therefore the calculation of PM is implemented in 4D space while that of CAM is carried out in 2D space. In order to compare the computational cost of both approaches, we set $N$ discrete points in the length $2\pi$ of each direction. In practice, the computational region of PM is $[0,2\pi)^{4}$ in 4D space, but that of CAM should be $[0,L\cdot 2\pi)^{2}$, where the value of $L\in\mathbb{Z}$ is dependent on the desired precision of $E_{DA}$. Therefore the number of discrete points of PM and CAM is $N^{4}$ and $(LN)^{2}$, respectively. As FIG. 10 (in Appendix A) shows, $E_{DA}$ does not decay monotonically as $L$ increases. The minimal integer $L$ of desired $E_{DA}$ for 2D DQC is listed in TAB. 1.

It should be pointed out that the PM computes real DQCs, while the CAM merely calculates corresponding crystalline approximants.

TAB. 2 gives the free energy of DQCs calculated by PM and by CAM ($L=126$, $E_{DA}=0.166879$) with different discretizations when $c=80$, $\tau=0$, $t=0$, $t_{0}=-0.3$, $g_{0}=-0.3$, $g_{1}=2.2$, $g_{2}=2.2$ and $d_{0}=d_{1}=d_{2}=d_{3}=0$. Applying CAM to compute the crystalline approximants of DQCs, the minimum computational region is $[0,126\times 2\pi)^{2}$. It was found that two approaches obtain different convergent energy values as $N$ increases. While $N\geq 18$, the free energy computed by both methods has $13$ significant digits.

Subsequently, we will demonstrate the correctness of PM in evaluating free energy. Subjected to the hard constraint $c\rightarrow\infty$, the nonzero Fourier modes are restricted on circles with radii $1$ and $q$. Under the symmetric assumption, we can analytically obtain the exact free energy function of any desired pattern, such as the DQC, by the TMAM chaikin1995principles ; Lifshitz1997 ; Jiang2015 ; Jiang2017 . Theoretically, for an ordered structure, the value of free energy obtained by a correct numerical method should converge to that under the hard constraint case as $c$ increases. To observe this computational phenomenon, we apply PM and CAM with $L=126$ ($E_{DA}=0.166879$), $L=204$ ($E_{DA}=0.091809$) to compute the DQC with same parameters and increasing penalty factor $c$ and observe the energy values. $16$ discrete points, corresponding to $16$ Fourier basis functions, in the length $2\pi$ of each direction are used in both approaches. The green line denotes the free energy $F_{c\to\infty}$ under hard constraint whose energy expression is given by Eq.(13). It is obvious that the free energy is heavily dependent on $E_{DA}$. When the $E_{DA}$ decreases from $0.166879$ to $0.091809$, the intersection point where the free energy computed by the CAM nearly equals to $F_{c\to\infty}$ is from about $c=500$ to $c=1800$. Therefore in order to obtain high accurate free energy by CAM, it is required to decrease the $E_{DA}$ with increasing computational area. However, the free energy obtained by CAM is still divergent from $F_{c\to\infty}$ with the growth of $c$ as long as $E_{DA}$ does not vanish. In other words, CAM may not obtain the correct free energy when $c$ is large enough. The essential reason is that the CAM cannot capture the real spectrum points as FIG. 11 shows. Correspondingly the energy of interaction potential part diverges quickly as $c$ becomes larger, as shown in FIG. 1 (b). On the contrary, the free energy and the interaction potential part obtained by the PM converge to the hard constraint case with the increase of $c$. Therefore, numerical results demonstrate that the PM can evaluate the free energy of quasicrystals correctly.

Furthermore, we compare the computational amount of two numerical methods in the above simulations. FIG. 2 presents the CPU time of PM and CAM ($L=126$, $L=204$) with respect to different $N$. It is intelligible that the blue line ($L=204$) is always above the red one ($L=126$) due to the formidable growth of computational domain $[0,L\cdot 2\pi)^{2}$. Meanwhile, $N$ is from $10$ to $22$, PM (black line) costs much less computational amount than CAM regardless of $L=126$ or 204. The reason is attributed to the implementation of PM in a high-dimensional unit cell $[0,2\pi)^{4}$. The degree of freedom of PM is $N^{4}$, while that of CAM is $\left(L\cdot N\right)^{2}$. The computational cost of PM is much smaller than that of CAM, since $L$ grows faster than $N$.

III.2 Candidate ordered phases

Based on the PM, minimizing the free energy functional (2) of the CMSH model allows us to investigate the equilibrium phase behavior of quasicrystals and periodic crystals. In what follows, we will focus on the occurrence and stability of 3D IQCs and related ordered structures, thus we will set the ratio between two characteristic length scales as $q=q_{2}/q_{1}=2\cos(\pi/5)$ Subramanian2016 ; Jiang2017 . Due to the choice of the ratio $q$, the 2D DQC is also considered as a possible equilibrium structure. In the PM, 3D IQCs are required to be embedded into $6$-dimensional periodic phases, and DQCs can be projected from the $4$-dimensional periodic structures. The projection matrices have been given in our previous work Jiang2014 ; Jiang2015 ; Jiang2016 ; Jiang2017 . In addition, a large number of related 2D periodic structures are also contained as the candidate structures in our study. In the following simulations, the $n$-dimensional Fourier space is discretized by $16$ basis functions along each direction. The total number of variables $16^{n}$ is enough to determine the relative stability of candidate phases according to the analysis in Sec. III.1.

In order to quickly obtain desired ordered structures, the initial configurations are important in numerical simulations. For the two-length-scales of the CMSH model, the initial nonzero Fourier vectors are chosen as FIG. 3 shows. The central black dot represents the origin of the Fourier space and the blue (resp. red) dots surrounding it are endpoints of wave vectors in $\mathbf{k}_{\psi}$ (resp. $\mathbf{k}_{\phi}$). Among these initial configurations, the first one is the IQC. The blue (resp. red) points represent $30$ diffraction points with icosahedral symmetry located on the spherical surface of radius $1$ (resp. $q$). Edges on spherical surfaces with radii $1$ and $q$ are indicated by the magenta and luminous green lines, respectively. For the other 2D patterns, the blue points stand for diffraction points located on the circle with radius $1$, while the red ones located on the circle with radius $q$.

Applying effective numerical methods and appropriate initial values outlined above to the dynamic equation (4), we can obtain abundant equilibrium structures, as shown in FIG. 4 and FIG. 5. Among these phases, the diffraction patterns and the physical morphologies of 3D IQCs are given in FIG. 4, whose initial nonzero Fourier vectors are illustrated with FIG. 3 (a). These nonzero Fourier vectors with diffraction intensity larger than $1\times 10^{-4}$ are given in FIG. 4 (a). Except for a host of small black spots, the pattern is entirely consistent with the initial configuration of IQC in the Fourier space. Projecting the 3D diffraction pattern of IQCs onto appropriate planes, we can found two-, three-, and five-fold symmetric axes in the spectra as shown in FIG. 4 (a1)-(a3). These diffraction patterns are consistent with the experimental results of metallic alloys, such as Al-Cu-Fe, Zn-Mg-Sc, Al-Rh-Si alloys Tsai1987 ; Ishimasa1988 ; Kaneko2001 ; Koshikawa2003 ; Maezawa2004 ; Honma2007 ; Ishimasa2011 . Correspondingly, FIG. 4 (b) displays its morphology by using the blue and orange colors to stand for components $\psi$ and $\phi$, respectively. The individual morphologies of $\psi$ and $\phi$ are given in FIG. 4 (b1)-(b2).

The morphologies of other candidate phases are shown in FIG. 5 in which (a) gives 2D DQC, (b) and (c) present lamellar phase with alternating beads but different density distribution (L+B[1], L+B[2]), (d) and (e) are hexagonal (HEX) and lamellar (LAM) structures, respectively. The initial non-zero Fourier vectors of these ordered patterns are given in FIG. 3 (b)-(f) in sequence.

III.3 Stability of IQCs and phase diagrams

In the subsection, we will investigate the relative stability of these candidate ordered structures in two-component CMSH model under the hard constraint (the limiting case $c\to+\infty$) and the soft constraint (finite $c$). Subjected to the hard constraint, the Fourier vectors of $\psi$ (resp. $\phi$) should be strictly restricted on the circle with radius $1$ (resp. $q$). Otherwise, the interaction energy becomes infinity. As described in Sec. II.3, for these candidate patterns whose basic wave vectors $\mathbf{k}$ can be found in FIG. 3, their energy expressions are given in Eqs. (12)-(17) under the hard constraint.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\widehat{\psi},\widehat{\phi}\in\mathbb{R}$ stand for the Fourier coefficients of the order parameters $\psi$ and $\phi$, respectively.

From these analytical expressions, it is apparent that the parameters $t_{0}$ and $g_{0}$ only exist in the expression of IQCs. It means that $t_{0}$ and $g_{0}$ can affect the relative stability of IQCs explicitly. In order to analyze the influence, we first select different values of $t_{0}$, $g_{0}$, fix $q=2\cos(\pi/5)$, $g_{1}=2.2$, $g_{2}=2.2$, $d_{0}=d_{1}=d_{2}=d_{3}=0$, and leave $\tau$ and $t$ as free parameters to observe the phase behavior of these candidate patterns. FIG. 6 (a)-(e) present phase diagrams of $[t_{0},g_{0}]=[0.7,0.7]$, $[t_{0},g_{0}]=[0.0,0.0]$, $[t_{0},g_{0}]=[-0.2,-0.2]$, $[t_{0},g_{0}]=[-0.3,-0.3]$, and $[t_{0},g_{0}]=[0.0,-0.2]$. The stable region of the different candidate phase is denoted by the corresponding color. MIX denotes the coexist region of L+B[1] and LAM in which the free energy of two phases is indistinguishable. In FIG. 6 (a), the stable area of IQCs in the phase diagram is almost invisible. When $t_{0}=0$ and $g_{0}=0$, IQCs have an obvious stable region as shown in FIG. 6 (b). It means that the three-body interactions, corresponding to the cubic terms, can affect the stability of IQCs. We further lower the values of $t_{0}$ and $g_{0}$, the stable region of IQCs grows rapidly, as FIG. 6 (c)-(d) show. It conducts that negative numbers of $t_{0}$ and $g_{0}$ are more beneficial to form and stabilize IQCs than the positive case. The phase diagram with the unequal model parameters $t_{0}$ and $g_{0}$ is also given in FIG. 6 (e). Comparing with FIG. 6 (b), we find that the black region extends to the part above the diagonal $\tau=t$ as $g_{0}$ descends. On the contrary, FIG. 6 (e) has less area of the black part under the diagonal than the figure (c) due to the growth of $t_{0}$. As the model parameter $\tau$ or $t$ becomes negative, the absolute value of the cubic term becomes greater and leads to stronger three-body interactions. Since IQCs have the largest number of the basic Fourier wave vectors, three-body interactions make more contributions to IQCs than the other candidate phases. Then we investigate the influence of the last four terms through setting $t_{0}=-0.3,g_{0}=-0.3,q=2\cos(\pi/5),g_{1}=2.2,g_{2}=2.2,d_{0}=d_{1}=d_{2}=d_{3}=1$, and still letting $\tau$ and $t$ be free parameters to construct the phase diagram, as shown in FIG. 6 (f). Comparing with the subfigure (d), the stable regions of DQCs and IQCs dwindle sharply and LAM is invisible in the phase diagram. It shows that the four-body interactions can affect the stability of quasicrystals as well as periodic crystals.

We have explored the phase behavior of the two-component system described by the CMSH model under the hard constraint with different parameters. Subsequently, we turn to study the relative stability of the system under the soft constraint (finite $c$) using the highly accurate PM. Due to the expansive computational cost, we select model parameters, $c=80$, $q=2\cos(\pi/5)$, $t_{0}=-0.3$, $g_{0}=-0.3$, $g_{1}=2.2$, $g_{2}=2.2$ and $d_{0}=d_{1}=d_{2}=d_{3}=0$, as an example. The parameters are the same as the hard constraint case of FIG. 6 (d) except for the constraint factor $c$. The corresponding phase diagram ($c=80$) is given in FIG. 7. It is obvious that the phase diagram of $c=80$ is very different from the hard constraint case of FIG. 6 (d). The discrepancies can be ascribed to more complex behaviors of interaction terms in multi-component systems. As shown in FIG. 7, DQC and L+B[1] occupy relatively large areas which are almost equal to the stable region of IQCs, but visibly smaller than the hard constraint case. Under the soft constraint, the relative stability of IQCs has been significantly intensified and its stable region almost doubled in size. The reason is mainly attributed to the appearance of high-order non-zero Fourier modes which can form more cubic interactions to stabilize IQCs. Accordingly, the regions of DQC and L+B[1] have to shrink in the phase diagram. The phenomenon is different from the single component system in which the 2D DQC becomes metastable under soft constraint Jiang2017 . At the same time, the stable area of the L+B[2] phase remains unchanged. The HEX phase appears in the phase diagram and replaces the MIX region in FIG. 6 (d).

The biggest distinction between the soft and the hard constraints is the emergence of high-order Fourier modes when $c$ is finite which do not locate on the circles with radii $1$ and $q$. To further investigate the phase behavior, it is useful to analyze the contribution to the free energy from different Fourier modes. We split the Fourier modes into two parts: fundamental modes and higher-harmonics. The fundamental modes are located on the circles with radii $1$ and $q$, as shown in FIG. 3, and the higher-harmonics are the other Fourier modes. The fundamental part of energy is defined by the contribution of the fundamental Fourier modes to the free energy, while the higher-harmonic energy is the remainder when subtracting the fundamental part from the total energy.

We take the line $t=0$ in the phase diagram of FIG. 7 as an example to analyze the contribution to the free energy from the fundamental and higher-harmonic parts. FIG. 8 gives the free energy curves of these candidate phases. To observe the tendency of the free energy better, we use the free energy of the IQC phase as the baseline. As FIG. 8 shows, the IQC is the most stable phase among these candidate ordered structures when $-1.70\leqslant\tau\leqslant 0.5$. While the L+B[1] phase is favored in the small range $-2\leqslant\tau\leqslant-1.71$. It is consistent with the phase diagram of FIG. 7.

We split the free energy ($t=0$) into the fundamental part and the higher-harmonic part, as FIG. 9 presents. The corresponding parts of IQCs are still utilized as baselines. Firstly, FIG. 9 (a) gives the energy curves of fundamental modes of these candidate patterns. The IQC phase has the largest number, $60$, of nonzero Fourier modes located on the spherical surfaces with radii $1$ and $q$, however, its fundamental energy is not the lowest among these phases in our considered region $\tau\leqslant 0.5$. The reason is attributed to the fact that the quartic terms of IQC, which increase the value of free energy, are much larger than that of other candidate patterns as Eqs.(12)-(17) express. Among these 2D ordered structures, the DQC is favored by the fundamental modes since it has the larger number, $20$, of nonzero Fourier modes located on the circles with radii $1$ and $q$, thus forming more triangles in the Fourier space and obtaining a lower fundamental energy, as shown in FIG. 9 (a). Meanwhile, there is a big difference between the L+B[1] and L+B[2] phases in the fundamental energy although they both have the same number of Fourier modes with nonzero coefficients. From Eqs.(14)-(15), we find that the analytical expression of L+B[2] can be obtained from the expression of L+B[1] by exchanging $\psi$ and $\phi$ if we treat $\tau$ and $t$ equally without discrimination. It means that the components $\psi$ and $\phi$ play a completely contrary role in two different analytical expressions. More precisely, the L+B[1] phase has more four-body interactions with respect to the component $\phi$ in the fundamental part than the L+B[2] structure does. Based on the fact, we conclude that the component $\phi$ takes a more dominant position than $\psi$ to stabilize these candidate ordered structures with the given model parameters.

On the other hand, since the energy penalty factor $c$ is finite, the higher-harmonics cannot be ignored which have a profound impact on the free energy of the multi-component system, as FIG. 9 (b) shows. The IQC structure has the highest higher-harmonics contribution to the free energy, since its Fourier modes with nonzero coefficients form a large number of triad interactions that dramatically reduce the free energy. In contrast, the higher-harmonics contributions of 2D ordered structures generate fewer numbers of cubic interactions than that of the IQC does. Therefore, the IQC has the lowest higher-harmonic energy. Due to the competition of the fundamental modes and higher-harmonics, we obtain the stable regions of L+B[1] and IQC as shown in FIG. 7.

IV SUMMARY

In summary, we have investigated the thermodynamic stability of 3D IQCs in multi-component systems using a CMSH model with multi-length-scales. In the model system, the characteristic length scales contained in the interaction potential functions are acting on different order parameters through the positive-definite gradient terms. For two-component systems, when the ratio of the two-length-scales is set as the golden ratio, we predicted that the 3D IQCs and 2D DQCs are both stable for hard and soft constraints. Under the hard constraint ($c\rightarrow\infty$), the Fourier modes are strictly restricted on the circles of radii $1$ and $q$. Using the two modes approximation method, we systematically analyze the rule of emergence and stability of quasicrystals and a series of periodic crystals. Under the soft constraint (finite $c$), more higher-order non-zero Fourier modes arise besides those on the circles $|\mathbf{k}|=1$ and $|\mathbf{k}|=q$ which have a profound impact on the phase behavior. Thanks to the high-precision PM approach, it enables us to quantitatively analyze the phenomenon. Using the PM, high accurate free energy and more precise phase boundary have been obtained under the soft constraint. Meanwhile, the advantages of PM over the traditional CAM have been also discussed in detail. In particular, the numerical results about 3D IQCs are phenomenologically consistent with lots of experimental observations in metallic alloys Tsai1987 ; Ishimasa1988 ; Kaneko2001 ; Koshikawa2003 ; Maezawa2004 ; Honma2007 ; Ishimasa2011 and that demonstrates the reliability of the PM. These results provide a good understanding of the stability of 3D IQCs in multi-component systems. The numerical approaches and insights will be helpful for further studying related quasiperiodic physical systems.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (11771368), Hunan Science Foundation of China (2018JJ2376), Youth Project Hunan Provincial Education Department of China (Grant No. 16B257), and Project of Scientific Research Fund of Hunan Provincial Science and Technology Department (2018WK4006).

Appendix A CAM

The principal idea of CAM is using the periodic structures to approximate the quasiperiodic structures. For a $d$-dimensional quasicrystal, its reciprocal lattice vector $\mathbf{k}$ can be expressed by $d$ linearly independent reciprocal vectors, $\mathbf{e}^{*}_{1},\mathbf{e}^{*}_{2},\cdots,\mathbf{e}^{*}_{d}$,

[TABLE]

It is important to note that $\mathbf{k}$ cannot be represented by linear combinations of $\mathbf{e}^{*}_{i}$ with integer-valued coefficients since there always exists, at least, an irrational number $p_{i}\in\mathbb{R}$. However, the quasiperiodic function $\phi(\mathbf{r})$ can be approximately expanded as

[TABLE]

If there exists a rational number $L$ such that $Lp_{i}\in\mathbb{Z}$ or $Lp_{i}$ can be sufficiently close to a series of integers, then

[TABLE]

where $[y]$ represents the nearest integer of the real number $y$. Replacing $[L\mathbf{k}]$ by $L\mathbf{k}$, numerical methods designed for periodic structures can be used to treat quasicrystals. It should be noted that the rational number $L$ depends on these irrational numbers $p_{i}$ due to rotational symmetry and the desired precision of the approximation. Without loss of generality, we can always choose $(1,0,\cdots,0)$ as one of the primitive reciprocal vectors which leads $L$ to be an integer.

How to determine $L$ is the so-called Diophantine Approximation (DA) topic in the number theory which deals with the approximation of real numbers by rational numbers or integers. In FIG. 10, we give DA error, $E_{DA}$, as a function of the integer $L$ with the decagonal quasicrystalline order. According to our experiments, the DQCs can be computed by CAM at least $L=126$. We mark the first integer $126$ and the other integers which have lower $E_{DA}$ in FIG. 10. These marked integers show that $L$ increases very quickly as the desired precision improves slightly (also see TAB. 1).

Then we give the approximation and actual positions of these Fourier spectrum points to analyze the influence of $E_{DA}$, as shown in FIG. 11. The central black point is the origin of coordinates and its surrounding black points represent the actual positions. The red, green, blue, magenta and light blue dots show the approximated positions with $L=126$, $204$, $1288$, $2084$ and $3372$, respectively. It is obvious that the green points, $L=204$, are significantly closer to the real positions than the red ones which correspond to $L=126$. Although the positions of the points with the remaining colors can hardly be distinguished from the authentic ones in FIG. 11 (a), a gap still exists if we zoom in this figure, as shown in FIG. 11 (b). As a matter of fact, the $E_{DA}$ always exists in CAM unless $L\to+\infty$.

Reference