Detection of phonon and phason modes in intrinsic colloidal quasicrystals by reconstructing their structure in hyperspace

Quasicrystals may be constructed as projections or sections of higher-dimensional crystals [28]. Conversely, embedding the physical space into hyperspace naturally restores the symmetry between phonons and phasons as hyperspace excitations on an equal footing. In the following we consider for concreteness a two-dimensional decagonal quasicrystal that is derived from the five-dimensional integer lattice. Our technique naturally adapts to quasicrystals of other symmetry and dimensionality.

In order to construct the decagonal quasicrystal, we employ the Penrose basis vectors ${{\vec{p}}_{i}}$ that are also known from the Penrose tiling [36]. Two of these vectors, ${{\left(\,{{\vec{p}}_{1}}\right)}_{j}}=\cos \left(2j\pi /5\right)$ and ${{\left(\,{{\vec{p}}_{2}}\right)}_{j}}=\sin \left(2j\pi /5\right)$ for $j=0\cdots 4$, are chosen to span the physical space. Orthogonal to ${{\vec{p}}_{1}}$ and ${{\vec{p}}_{2}}$, the two-dimensional complementary space is spanned by ${{\left(\,{{\vec{p}}_{3}}\right)}_{j}}=\cos \left(4j\pi /5\right)$ and ${{\left(\,{{\vec{p}}_{4}}\right)}_{j}}=\sin \left(4j\pi /5\right)$. The fifth basis vector, ${{\left(\,{{\vec{p}}_{5}}\right)}_{j}}=1$, corresponds to the body diagonal of the hypercubic lattice. Projections along ${{\vec{p}}_{5}}$ do not lead to new degrees of freedom because ${{\vec{p}}_{5}}$ is commensurate with the hyperlattice. Hence, ${{\vec{p}}_{5}}$ will not play an important role in our analysis. The length of the physical-space projection of a primitive hypercubic lattice vector is ${{a}_{\parallel}}:={{(5/2)}^{-1/2}}$, and corresponds to the long length scale in the quasicrystal. The short length scale is given by ${{a}_{\parallel}}/\tau$ where $\tau =\left(\sqrt{5}+1\right)/2$ is the golden ratio.

In order to construct the two-dimensional quasicrystal, only particles are considered that are within a so-called acceptance domain. For the decagonal case, we choose an acceptance domain of decagonal shape in the ${{\vec{p}}_{3,4}}$-plane with a circumradius of ${{a}_{\parallel}}/3$, and infinite extent in ${{\vec{p}}_{5}}$-direction. The acceptance domain is attached to each lattice point of the five-dimensional integer lattice such that the whole structure is periodic. The resulting decagonal quasicrystal is similar to laser-induced quasicrystals [21–24], because the decagonal acceptance domain does closely resemble the equipotential surfaces of the laser field potential energy (circles in the ${{\vec{p}}_{3,4}}$-plane). The polygonal truncation of the acceptance domain corresponds to enforcing a minimal-distance criterion between two neighbouring particles in physical space. This relationship between projected and laser-induced quasicrystals was discussed by Jagannathan and Duneau for the case of octagonal structures [37].

In simulations, we employ finite periodic approximants to the aperiodic quasicrystal: the irrational Penrose vectors are approximated such that the physical-space box is spanned by multiples of lattice vectors of the hyperlattice, i.e. the modified Penrose vectors satisfy ${{\left({{L}_{x}}\,{{\vec{p}}_{1}}\right)}_{j}}={{m}_{j}}$ and ${{\left({{L}_{y}}{{\vec{p}}_{2}}\right)}_{j}}={{n}_{j}}$, where m_j and n_j are integers for all j. Unless otherwise noted, our examples we will use a system with m₁  =  13 and n₁  =  21 (${{L}_{x}}/{{L}_{y}}\approx 1.902\,51$). The excitation-free approximant in this simulation box consists of N  =  1686 particles. A larger system, with m₁  =  8 and n₁  =  55 (N  =  2728, ${{L}_{x}}/{{L}_{y}}\approx 0.448\,67$), is used in section 3.3.

The hyperlattice can be subject to spatially non-uniform displacements. If the displacements associated with an excitation occur in the direction of physical space, such modes are called phonons, and phasons otherwise. In periodic approximants, all modes can be decomposed into discrete modes that are the harmonics of the simulation box, such as displayed in figure 1. In the following, we also refer to modes by their Miller indices $\left({{h}_{x}},{{h}_{y}}\right)$, by their wavevector $\vec{k}=\left(2\pi {{h}_{x}}/{{L}_{x}},2\pi {{h}_{y}}/{{L}_{y}}\right)$ or by the corresponding wavelength $\lambda =2\pi /|\vec{k}|$.

Zoom In Zoom Out Reset image size

Each mode carries a polarisation consisting of two phononic (u_1,2 along ${{\vec{p}}_{1,2}}$) and two phasonic (w_1,2 along ${{\vec{p}}_{3,4}}$) components. The modes are subsumed into a four-dimensional displacement field $\left({{u}_{1}},{{u}_{2}},{{w}_{1}},{{w}_{2}}\right)$. For each hyperlattice point, the acceptance domain is shifted by the local value of this field, figuratively stretching (phononic) and buckling (phasonic components) the physical plane. While phonon modes displace particles in the physical space directions, phason modes shift the cut plane in a complementary-space direction. A phasonic excitation (see bold dashed line in figure 2) alters the pattern of hits and misses between the cut plane and the acceptance domain. The physical expression is the insertion or deletion of particles. In the case of small phasonic gradients, a deletion is typically associated with an insertion at another nearby position (coloured squares in figure 2). In the physical-space view, this gives the impression of a particle jumping. Such phasonic flips occur along well-defined displacement vectors which can be derived from the hyperspace geometry. The orientational discreteness of those physical-space motion vectors reflects in the small number of ten possible jump directions (i.e. colours) in figure 1. In this rather extreme example, about one in four particles undergoes a flip as a result of the phasonic excitation.

Zoom In Zoom Out Reset image size

In this section we discuss the analysis of phonons and phasons in a quasicrystal by reconstructing its embedding into hyperspace. In addition to the analysis of phonons, which can be straightforwardly carried out in physical space, this will give indications about the phasonic distortions underlying the observed flip patterns.

Basically, the reconstruction is achieved by reversing the projection procedure described in the previous section. For any two-dimensional particle configuration derived from a decagonal quasicrystal with known origin, one can construct an embedding into hyperspace as explained in the following. First, each particle i of the test pattern is embedded into five-dimensional space at the location ${{\vec{R}}_{i}}={{x}_{i}}\,{{\vec{p}}_{1}}+{{y}_{i}}\,{{\vec{p}}_{2}}$ parametrised by its physical-space coordinates $\left({{x}_{i}},{{y}_{i}}\right)$.

The particle i is then assigned to the hyperlattice site ${{\vec{L}}_{i}}$ closest to its position. In order to eliminate the trivial ${{\vec{p}}_{5}}$ component from the embedding, we project the five-dimensional integer lattice onto a subspace orthogonal to ${{\vec{p}}_{5}}$. The hyperspace points that are assigned to a lattice point ${{\vec{L}}_{i}}$ correspond to a four-dimensional Voronoi cell around that point. The residual parts of coordinates ${{\vec{R}}_{i}}-{{\vec{L}}_{i}}$ carry, separated into the physical and complementary directions, the information about phononic ($\vec{u}$) and phasonic ($\vec{w}$) displacements, respectively. In this way, we are able to reconstruct four-dimensional occupations from two-dimensional information, as long as deviations from an ideal embedded crystal are small compared to the extent of the acceptance domain, or rather the extension of the four-dimensional Voronoi cell.

In effect, we assign each physical particle to a point in hyperspace. The phononic offset is fully recovered, and all the residuals ${{\vec{u}}_{i}}$ indeed collapse to the origin if no phonons are present. The phasonic residuals $\vec{w}$ fall into the acceptance domain in complementary space, and fill it densely in the infinite-system limit. Displays of the residual phasonic coordinates along ${{\vec{p}}_{3}}$ and ${{\vec{p}}_{4}}$ are given as insets of figure 1 and the systems in section 3 with the central decagonal domain containing all particles which have not been flipped (black dots).

In the hyperspace picture, a phasonic flip corresponds to an assignment of a particle to an adjacent lattice site, so that its complementary-space coordinates lie outside the acceptance domain of the original lattice site (coloured squares in inset of figure 1). Conversely, the physical directions of the embedding contain the phononic part of particle motions between an idealised template and the analysed structure. This decoupling between phononic and phasonic movements effectively separates the two modes of hyperspace displacement. We exploit this to identify phasonic flips.

We now proceed to describe our analysis method for the system as a whole, to extract amplitudes of collective phononic and phasonic excitations. We apply the downhill simplex optimisation algorithm [38] to determine the amplitudes of a set of modes which have been chosen in advance. Our choice are plane waves with amplitudes and polarisations given by $\vec{u}$ and $\vec{w}$ for the phononic and phasonic part of the distortion space, respectively. In order to avoid mixing amplitude-type and phase-type variables, we use as orthogonal basis the harmonic functions $\cos \left(\vec{k}\centerdot \vec{r}\right)$, $\sin \left(\vec{k}\centerdot \vec{r}\right)$.

The reconstruction of a phononic displacement field by embedding it into the hyperspace is lossless. Therefore, we may extract the field directly by applying the downhill simplex optimisation to an energy functional that is proportional to the sum of phononic misfits (excluding jumped particles).

In contrast, the reconstruction of phason modes uses a more sophisticated energy functional: given the current parameter set, a trial system is constructed by the method described in section 2.2. Its particle placements are evaluated with penalties for false predictions. With the goal of an approximative recovery of thermal systems in mind, we design the energy functional to be a continuous function of the distortion amplitudes. Rather than a hard binary hit-or-miss criterion, we choose the error function, as a function of distance from the edge of the acceptance domain, to penalise badly-placed particles. The downhill simplex algorithm now minimises the difference in distribution of particles between the probed system and trial configurations. During the optimisation run, the sigmoid width of the error function is reduced, starting from a value of $0.03{{a}_{\parallel}}$, reduced step-wise by an overall factor of two.

The result of these two schemes is a quantitative interpretation of how the real-space distortions and phasonic jumps can be understood by deformations of the physical plane within its embedding into hyperspace.

We study thermally activated phononic and phasonic excitations in colloidal quasicrystals by employing Brownian dynamics simulations. In the simulations, we integrate the overdamped Langevin equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \gamma \frac{\text{d}}{\text{d}t}\vec{r}={{{\vec{F}}}_{\text{int}}}+{{{\vec{F}}}_{T}} \end{array}\end{eqnarray}$$

with friction coefficient γ and position $\vec{r}$ of the particle in two dimensions. The random thermal force ${{\vec{F}}_{T}}$ complies with $\langle {{\vec{F}}_{T}}\rangle =0$ and $\langle {{F}_{T,k}}(t)\centerdot {{F}_{T,l}}\left({{t}^{\prime}}\right)\rangle =2\gamma {{k}_{\text{B}}}T{{\delta}_{kl\delta}}\left(t-{{t}^{\prime}}\right)$ where k_B is the Boltzmann constant and T the temperature of the system. The internal force ${{\vec{F}}_{\text{int}}}$ results from a pair potential V(r) such that ${{\vec{F}}_{\text{int}}},{{}_{i}}\,\,={\sum}_{i\ne j}{{\vec{F}}_{ij}}$ with ${{\vec{F}}_{ij}}=-\vec{\nabla}V\left({{r}_{ij}}\right)$, ${{r}_{ij}}=\,\,|{{\vec{r}}_{i}}-{{\vec{r}}_{j}}|$ denotes the distance between two particles i and j.

The particles in our system interact through the Lennard-Jones–Gauss potential [33, 34]

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{V}_{\text{LJG}}}(r)}{{{\epsilon}_{0}}}={{\left(\frac{{{r}_{\text{LJ}}}}{r}\right)}^{12}}-2{{\left(\frac{{{r}_{\text{LJ}}}}{r}\right)}^{6}}-{{\epsilon}_{\text{G}}}\exp \left(-\frac{{{\left(r-{{r}_{\text{G}}}\right)}^{2}}}{2r_{\text{LJ}}^{2}{{\sigma}^{2}}}\right) \end{array}\end{eqnarray}$$

with the energy unit ${{\epsilon}_{0}}$. For the parameters ${{r}_{\text{G}}}=1.52$, ${{\epsilon}_{G}}=1.8$ and ${{\sigma}^{2}}=0.02$, V_LJG describes a double-well potential with two incommensurate length scales (see figure 3), favouring local decagonal symmetry and supporting quasicrystalline structures [33].

Zoom In Zoom Out Reset image size

Our simulation quantities are dimensionless: energies are given relative to ${{\epsilon}_{0}}$, and temperatures T in terms of ${{\epsilon}_{0}}/{{k}_{\text{B}}}$. The time unit is given by the Brownian time ${{\tau}_{\text{B}}}=r_{\text{LJ}}^{2}\gamma /{{\epsilon}_{0}}$. We choose $\text{d}t={{10}^{-4}}{{\tau}_{\text{B}}}$ for a single Brownian dynamics time step. Furthermore, we apply a cut-off radius ${{r}_{\text{c}}}=2.5{{r}_{\text{LJ}}}$ and subtract the cutoff force $F\left({{r}_{\text{c}}}\right)=-7.83\centerdot {{10}^{-3}}$ to ensure continuity. We determine the unit distance ${{r}_{\text{LJ}}}=0.6554{{a}_{\parallel}}$ as the length scale of the pair potential such that—in accordance with the barostat used in references [33, 34]—a pressure-free system at T  =  0 is achieved. Note there is a substantial detuning of the Lennard-Jones minimum relative to the pair distances of the actual system (figure 3), in contrast to the near-perfect match of the Gauss minimum.

We obtain initial positions of the particles from the projection method using a periodic approximant. The systems are heated up and held at a constant temperature for $t=10{{\tau}_{\text{B}}}$, so that thermal equilibrium is reached. We remain well below the melting temperature ${{T}_{\text{M}}}=0.56\pm 0.02$ from [34]. Prior to analysis, systems are instantaneously quenched: we continue the Brownian dynamics simulation at T  =  0 for a duration of $t=5{{\tau}_{\text{B}}}$ to allow the particles to find the equilibrium positions in their local environment, i.e. to freeze out phononic thermal noise.

We first determine the number of phasonic flips which are introduced by a phasonic excitation with respect to the reference state. The phasonic excitation is constructed by choosing a set of mutually orthogonal harmonics with respect to the box size. To each individual mode $\vec{k}$, a random amplitude vector $\vec{w}\left(\vec{k}\right)$ and a random phase are assigned. We quantify the total magnitude of the excitation W by the sum of the absolute values of the amplitudes of the individual modes, i.e. $W={\sum}_{{\vec{k}}}|\vec{w}\left(\vec{k}\right)|$.

Figure 4 (left) displays the phasonic flip fraction, i.e. the number of flipped particles divided by the total number of particles, as a function of the total amplitude W of the phasonic excitation. We find that the flip fraction increases linearly with the total amplitude and therefore with the amplitude of phasonic strain. Under the premise that all particle jumps in a given configuration are caused by phasonic excitations, this relationship could be used to estimate the total phasonic amplitude from the number of flips. Surprisingly, the number of wavevectors in the mode set seems to be insignificant, i.e. it is irrelevant if the total amplitude is contained in a single harmonic or in arbitrarily many modes.

Zoom In Zoom Out Reset image size

Figure 4 (right) displays the flip fraction found after annealing by using Brownian dynamics simulations of duration $t=10{{\tau}_{\text{B}}}$ as a function of the inverse temperature. The flip fraction exhibits Arrhenius behaviour $\propto \exp \left(- \Delta E/{{k}_{\text{B}}}T\right)$ with $\Delta E\approx 1.43{{\epsilon}_{0}}$ and is compatible with a thermal activation of phasonic waves.

In this section we test the reliability of recovering a single, known phasonic excitation from a flip pattern that was constructed using the method described in section 2.2. An example of such a system is given in figure 1. We now apply the optimisation algorithm of section 2.3, restricted to the known phasonic wavevector. As can be seen in figure 2, the recovery result (blue bold line) agrees well with the original phason wave (black dashed line). The top panel of figure 5 displays where particles were flipped (black squares) due to the synthesised phason. In addition, we show where the best-fit reconstructed phason wave predicts the deletion (blue circle) and insertion (red) of particles. In $94.0\, \%$ of the created and $94.7\, \%$ of the deleted particle positions, the fit agrees with the synthetic test system, hence the majority of flips has been correctly interpreted. False predictions usually occur if particles in hyperspace are close to the border of the acceptance domain (indicated by a lighter colour in the top pane of figure 5). Part of these discrepancies are due to the deliberate smoothening of the energy functional that we use for optimisation.

Zoom In Zoom Out Reset image size

The quality of recovery of synthetic phasons depends on the wavelength of the phason. As figure 6 demonstrates, the recovery improves with the wavelength of the phasonic wave. The spatial resolution of a phason mode is limited by the typical distance between particles, making short-wavelength phasons less well-defined. Phasons with wavelengths $\lambda >10{{a}_{\parallel}}$ are discovered with good accuracy.

Zoom In Zoom Out Reset image size

We now explore what kind of excitations occur due to thermal fluctuations and how we can analyse them. We use excitation-free quasicrystals as starting configurations of Brownian dynamics simulations at constant temperature T. Noticeable numbers of phasonic flips start to occur above T  =  0. 2 (compare with right-hand panel of figure 4). The vast majority of these flips persist when the system is quenched to T  =  0 for analysis. These phasonic flips support, even at T  =  0, a static phononic displacement field. Below T  =  0.35, we find that the phonon field is typically dominated by a transversal (1, 0) mode aligned with the long dimension of the simulation box. Figure 7 illustrates this effect: after annealing for $t=10{{\tau}_{\text{B}}}$ at T  =  0.300 the final configuration differs by 82 phasonic flips from the initial configuration. The flips enable a phononic strain field that does not level out at T  =  0.

Zoom In Zoom Out Reset image size

We now use our optimisation procedure to extract the amplitudes of the six longest-wavelength modes ($\lambda >13{{a}_{\parallel}}$). For the system shown in figure 7, we find that the majority of strain lies within the transversal (1, 0) mode with amplitude ${{u}_{2}}(1,0)=0.053{{a}_{\parallel}}$, while the other modes vanish within fit accuracy ($|u\left(\vec{k}\right)|<{{10}^{-2}}\times {{u}_{2}}(1,0)$). The best-fit results are consistent between successive runs of the optimisation with different random initial parameters. A system with a larger number of particles is shown in figure 8. In this case, a different approximant is used (see section 2.1), and the quasicrystal is rotated by 90 degrees with respect to the long dimension of the box. Again, a transversal phonon mode dominates the displacement field, but in this larger system, it is accompanied by additional long-wavelength harmonics.

Zoom In Zoom Out Reset image size

Next, we want to interpret flips as manifestation of a continuous phasonic displacement field. Again, we consider the system of figure 7. In the optimisation procedure, we use the same modes as for the phonon fit, but with phasonic polarisations. In the bottom pane of figure 9, we display the best-fit phasonic displacement field within these modes. No mode clearly dominates the phason field. The total amplitude amounts to $W=0.059(14)\,{{a}_{\parallel}}$, with the largest contribution coming from a (2, 0) mode with $|\vec{w}(2,0)|\,=0.019(5)\,{{a}_{\parallel}}$ and elliptical polarisation. When we interpret the density of flips in the annealed system in the context of section 3.1, we expect the summed phason amplitude at about $0.045{{a}_{\parallel}}$, which is in reasonable agreement with the value of W determined in this section by the optimisation method.

Zoom In Zoom Out Reset image size

The phason fit is less reliable than the phonon fit. For different initial parameters, the phason fit yields parameter sets which vary within a standard deviation of about $30 \%$. Furthermore, the top pane of figure 9 demonstrates that the optimisation performs worse when predicting the phasonic flips in the Brownian dynamics simulations than for synthetic phasons. Overall, $32 \%$ of the created and $34 \%$ of the deleted particles agree between the annealed configuration and the recovered phasonic field. Note that wrong predictions are concentrated on particles close to the boundary of the acceptance domain (light coloured circles) while flips of particles well within the acceptance domain are usually described correctly (intense coloured circles).

The fraction of correctly characterised flips in principle could be increased by including phason modes with shorter wavelength into the optimisation process. On the other hand, with an increasing number of parameters, the deviations of the final parameters for different initialisations would further increase. Inherently, long-wavelength harmonics can only account for flip events that occur coherently though the whole system. Hence the optimisation systematically ignores strongly localised flip patterns, e.g. driven by thermal fluctuations. Nevertheless, with our choice of mode selection, we are able to extract the smooth part of the phasonic field, at the cost of ignoring localised effects.

In a continuum picture, phononic and phasonic strain fields are expected to be coupled [3]. In principle, one could conclude from this that a continuous phasonic strain field (evidenced by frozen-in phasonic flips) might give rise to the observed static phononic distortions. Our analysis does not reveal any strong correlation between the observed phonon modes and the dominant phasonic modes. However, as shown above, a large fraction of observed phasonic flips cannot be accounted for without shorter-wavelength phasons. In this limit, a continuum description of the quasicrystal might not be appropriate, such that the phonon/phason coupling is not observed on the atomic scale.

Phasonic strain usually is not relaxed at zero temperature: in order to flip a particle one has to cross an energy barrier. However, thermal fluctuations can in principle enable the system to overcome energy barriers and level out phasonic strain. Therefore, we study the relaxation of synthetic phasonic strain fields during Brownian dynamics annealing.

We start with a synthetic single-mode phason and track the relaxation of its amplitude over Brownian dynamics simulation time. Its amplitude may relax, and its phase and polarisation may drift with respect to the original phason. One test system exhibits a (1, 0) mode corresponding to a wavelength of $40.0{{a}_{\parallel}}$ and a phasonic amplitude of $0.19{{a}_{\parallel}}$. A second system has a (2, 3) mode, i.e. a shorter wavelength of $6.6{{a}_{\parallel}}$ and an amplitude of $0.105{{a}_{\parallel}}$. The relaxation of both systems is studied at two temperatures, T  =  0.2 and T  =  0.3. We take periodic snapshots and fit a phasonic wave to each. The amplitude of the recovered phasons is shown in figure 10 as a function of time. As shown by the colour code, drifts of the phase and changes of the polarisation stay small during our simulation.

Zoom In Zoom Out Reset image size

For the temperatures considered, the magnitude of the (1, 0) mode shows no indications of relaxation in the presence of thermal noise. Some phase or polarisation misfit is determined from the beginning of the simulation, but it stays constant during the simulation. With its shorter wavelength, the (2, 3) mode behaves differently: initially, a fast, relaxation-like process reduces the amplitude. Then, the amplitude stabilises at a saturation level which is lower for higher temperatures. Our results may be rationalised considering the fact that, at the same amplitudes, the phasonic strain becomes more localised for shorter wavelengths, which might be conductive to relaxation. The saturation indicates that depending on the temperature there is a certain amount of phasonic stress that can prevail in the system. In other words, the phasonic stress associated with the saturation amplitudes is not large enough to initiate correlated phasonic flips that would be necessary to further relax the system.