Distinguishing noisy crystalline structures using bond orientational order parameters

The bond orientational order parameters originally introduced by Steinhardt \emph{et. al.} [Phys. Rev. B \textbf{28}, 784 (1983)] are a common tool for local structure characterization in soft matter studies. Recently, Mickel \emph{et. al.} [J. Chem. Phys. \textbf{138}, 044501 (2013)] highlighted problems of the bond orientational order parameters due to the ambiguity of the underlying neighbourhood definition. Here we show the difficulties of distinguish common structures like FCC- and BCC-based structures with the suggested neighbourhood definitions when noise is introduced. We propose a simple improvement to the neighbourhood definition that results in robust and continuous bond orientational order parameters with which we can accurately distinguish crystal structures even when noise is present.


Introduction
A great benefit of studies using soft matter systems is observability on the local scale.Thermodynamic and statistical processes like phase transitions and crystallization can be studied on the level of individual constituents, like solid particles in colloids, dusty plasma and granular media, gas bubbles in foams or droplets in emulsion [1,2,3,4,5,6,7].
Steinhardt and co-workers introduced the bond orientational order parameters q l [8] as a useful local order metric in such studies.The arrangement structure of an individual constituent is quantified by the angles among its 'bonds', i.e. the directions to the neighbouring constituents (see Sec. 3 for details).The q l are rotation invariant and are sensitive to symmetries in the bond angles, thus comparing measured values to values of ideal lattices allows for identifying crystalline arrangements on the level of individual constituents.
These Steinhardt bond orientational order parameters are widely used [9,10,11], but are not uniquely defined due to the ambiguity of neighborhood.Here we continue a discussion by Mickel et al. on different definitions and quantification of neighbourhood [10], and show an optimized definition in the presence of noise.Noise implies that the measured positions of the constituents fluctuate around the equilibrium, mean or real positions.Such noise may arise from a finite resolution of the experimental technique and the resulting uncertainty in measured positions.Intrinsic noise also emerges from thermal fluctuations of positions and from polydispersity of constituent shape and interactions [2,12].The measured arrangement thus will usually exhibit a distribution of local arrangements of the constituents.The local order metric reflects this distribution and extraction of information on underlying crystal structures thus needs additional measures.For example, the analysis was constricted to statistical decomposition of the measured distributions of the q l into the distributions of q l obtained numerically from different noisy crystal structures in order to retrieve information on the apparent fractions of the respective crystal [2,13,14,15].Assignment of crystal structures to q l of individual constituents in studies with noise relies on the distance of the individual q l to the q l of the ideal crystal lattice site.Sensitivity of this approach was improved by using more than one q l and mapping regions in the q l − q m -plane to certain crystal structures or by using information on the neighbours of the neighbours [9,16,17].
However, the bond orientational order parameter has to rely on a neighbourhood definition.This quantification of neighbourhood itself is affected by the noise to a degree which depends on the used definition, as we highlight in Sec. 4. Consequently, measured distributions of q l may not reflect the distribution of local arrangements and structure identification turns ambiguous.Mickel et al. have discussed different neighbourhood definitions and proposed an unambiguous and continuous morphometric neighbourhood definition (see Sec. 3) [10].We suggest a simple optimization of the morphometric neighbourhood definition.This suggested definition minimizes detrimental effects of noise and facilitates discrimination between the common crystal structures like body-centered cubic (BCC), face-centered cubic (FCC) or hexagonally closepacked (HCP) lattices even in the presence of noise.We describe in Sec. 2 how we generate the test data.In Sec. 3 the Steinhardt bond orientational parameters and the modifications introduced by Mickel et al. are discussed.We also motivate our optimization in this section.The effects of noise and the difficulty of discriminating BCC-, FCC-and HCP-based structures are highlighted in Sec. 4. We show that the improved version of the morphometric neighbourhood definition enables the desired structure differentiation among the common crystal structures with noise.

Generation of structures and tessellation
We numerically mimic the effect of the various noise mechanisms mentioned in Sec. 1 in the positional data.We start by creating ideal HCP, FCC and BCC lattices with ≈ 4500 lattice sites by placing a site at integer multiples of the primitive vectors for the corresponding crystal structure.We then apply normally distributed noise to every Cartesian component of each lattice site to simulate different degrees of noise.The normal distributions from which the x-, y-and z-displacements are drawn have a mean of 0 and standard deviations from σ = 0.01d up to σ = 0.1d in 10 logarithmically spaced steps, where d is the lattice constant.We note that normally distributed noise applied to each component results in a Maxwell-Boltzmann distribution of the absolute radial displacements.The root mean square displacement is related to the componentwise standard deviation via  1: q 4 and q 6 values for noise-free hexagonal closepacked (HCP), face-centered cubic (FCC) and bodycentered cubic (BCC) lattices.The q l -values of the BCC lattice sites depend sensitively on the definition of neighbourhood, in contrast to the values for FCC and HCP lattice sites.q l values are given for neighbourhood definitions using a fixed threshold of 1.1d and using facet-weighting α of order 0, 1, and 2.
An common criterion for predicting melting of a crystal is the Lindemann criterion [18], which states that thermal motion in a crystal leads to melting once the root mean square displacement exceeds a certain fraction δ l of the lattice constant.Typical values of this critical fraction are δ l ≈ 0.22 [19], which is still higher than the highest fraction we tested of 171.Therefore, crystals can realistically exhibit deviations as high as the ones we investigate here.
Our protocol leads to 10 homogeneous structures with increasing levels of noise.Figure 1 shows the radial distribution functions of the BCC-and FCC-based structures with noise with σ = 0.01d and σ = 0.036d.The structures exhibit long-ranged positional correlations and distinguishable local environments even with noise, only the correlation peaks become increasingly broadened by the random displacements.
We use the Voro++-code for a Voronoi tesselation of the test structures [20].The code provides vertices, edges and area of the facets of the Voronoi cells of the sites.We include all sites for the calculation of the bond orientational order but we exclude in our results the boundary sites, i.e. those that share a facet with with the lattice boundary, and any site sharing facet with a boundary site.This exclusion of two layers of sites is necessary when working with BCC-based structures, since the Voronoi facets of the sites are determined both by the nearest and the next-nearest neighbours (see Fig. 2. The data shown is thus for about 3400 sites each.
These test structures are consequently evaluated using the Steinhardt bond orientational order parameter using different neighbourhood definitions.

The order metric and the neighbourhood definition
The Steinhardt bond orientational parameters q l are rotation invariant metrics which are sensitive to periodicity in the orientations between neighbouring sites [8].The q l Fig. 2: The bond orientational order of the ideal noisefree lattice sites for the four different neighbourhood definitions, given as position in the q 4 -q 6 plane (compare Tab. 1).The position of the BCC lattice site shifts closer or further from the close-packed lattice sites, depending on whether a fixed threshold ("t") or a facet-weighting α of order 0, 1, and 2 is used in the neighbourhood definition.The inset shows the outline of the Voronoi cell of an FCC (left) and BCC (right) lattice site.are calculated for any site i in the medium according to: where Z i is the number of neighbours of the test site i, Y m l are the spherical harmonics functions of degree l and order m, x j is the vector connecting the test site with its jth neighbour and θ and φ are the angles of x j in spherical coordinates.
An obvious method to determine the neighbours is to count all site positions within a cutoff distance r c as neighbours [1,2,14,15].This approach is easy to implement, but introduces the parameter r c .This parameter is either set arbitrarily [1,14] or is defined by the first minimum of the radial distribution function g(r) [2,15].A problem with this definition of r c is apparent from Fig. 1.The first minimum of g(r) may readily be at distances smaller or larger than the shell of the next nearest neighbours in the BCC structure, depending on the level of noise.The usage of a fixed r c to determine neighbourhood also implies a discontinuous quantification, as infinitesimal shifts in position can produce or annihilate a bond to a neighbor.
A parameter-free determination of neighbourhood is based on the Delaunay-triangulation.The Delaunay-triangulation creates the dual graph of the Voronoi-tesselation, thus every site that shares a facet of the Voronoi cell of the test site is counted as a neighbour [21].As Mickel et al. pointed out, this method of triangulation is parameterfree, but still not continuous [10].Infinitesimal changes in the site positions can lead to forming or disappearance of facets and thus of bonds to neighbours, resulting in a discontinuous change in the value of q l .They instead suggested a morphometric neighbourhood, i.e. a neighbourhood relation that is weighted with the area of the Voronoi facet associated with the bond: Here A f is the area of the facet f associated with neighbour j, A is the total surface area of the Voronoi cell, and the sum is over the ensemble of facets of the Voronoi cell F(i) of site i.In this definition the q l are unambiguous, parameter-free and continuous.Infinitesimal positional changes will create or annihilate facets with small areas and consequently the new or vanished bonds will only negligibly contribute to q l .The original formula for the q l (Eq.1), with using a Delauney-triangulation to determine neighbourhood, and eq. 2, using the morphometric neighbourhood definition with the additional facet-weighting, can be written into the general formula The exponent α of the facet area equals 0 for the Delaunay triangulation, which is just normalized by the number of bonds, and is 1 for the morphometric facet-weighting.The normalization A(α) of each bond reads The impact of the neighbourhood definition on the value of q l can be nicely illustrated using this parametrization.We use the q 4 and q 6 bond orientational order parameters.These parameters are most sensitive to the cubic and hexagonal symmetry, respectively, as we investigate here and are commonly used in studies on crystallization [1,2,14,15].The q l of a lattice site in the ideal crystalline lattices are summarized in Tab. 1 and Fig. 2. We compare the q 4 -and q 6 -values obtained using a distance threshold r c of 1.1 times the lattice constant d, and using the Voronoi-tesselation-based methods with an area-weighting exponent α of 0 an 1.The q l -values of the FCC and HCP lattice sites do not depend on the method of defining and weighting neighbourhood since they share only equal facets with the 12 nearest neighbours.The BCC lattice sites, in contrast, change their q l -values by up to an order of magnitude with the different neighbourhood definitions, and the q 4 is more affected than the q 6 .The strong dependence of the q l of the BCC lattice sites on the neighbourhood definition can be rationalized with their Voronoi cell (see inset in Fig. 2).The Voronoi cells of the BCC lattice share facets with the 8 nearest, but also smaller facets with the 6 nextnearest neighbours.The 6 next-nearest neighbours are not counted as bonds when using a threshold, are counted equally when using Delauney-(α=0)-weighting and are counted to less extent with morphometric α=1-weighting.
We suggest that even better results can be obtained by using an α=2 facet-weighting of the bonds.This is motivated primarily by the behaviour of q 4 when considering the BCC lattice.When a threshold is used such that only the 8 nearest neighbours are considered, a high fourfold symmetry compared to HCP and FCC is found (see Tab. 1).If instead the next nearest neighbours are also taken into account and weighted equally, such as for α=0, the fourfold symmetry is nearly completely destroyed.Heuristically, the more the close neighbourhood is emphasized, the easier it is to distinguish BCC by its fourfold symmetry from FCC and HCP.An α=2 facetweighting emphasizes the contribution from the larger facets assigned to the nearest neighbours compared to the facets of the next-nearest neighbours, thus is closer to using a threshold than α=1 facet-weighting.Still the q l stay a continuous measure as with the α=1-weighting.The values of q 4 and q 6 of the ideal lattice sites are given in Fig. 2 and Tab. 1.The difference among the q 4 and the q 6 of the ideal FCC and the ideal BCC lattice is larger using the α=2-weighting than using the α=1-weighting, which facilitates structure identification even in noisy crystals.

Application to noisy crystals
The dependence of the q l -values on the weighting of the neighbourhood may be considered a nuisance in treatment of ideal structures.The numeric value of the q l shift, but still stay unique for the different crystal structures.This changes when noise is present in the structures.
Noise has severe impact on the q l -values of the sites.We exemplify the effect of the noise when using the three established neighbourhood metrics in Fig. 3. None of the three neighbourhood metrics can be considered satisfactorily.A splitting of the q 4 -q 6 -values of the individual lattice sites of the noisy BCC-structure into isolated clouds can be observed in Fig. 3, a), where the threshold-definition of neighbourhood is used.The individual distributions of q 4 and q 6 on the axes of the graph show that the clouds of the BCC-structure already start to overlap with the cloud of the FCC-structure at this moderate noise of 3.6%.The q 4 -q 6 -values of the noisy FCC-structure splits up into individual clouds in a similar fashion when using a Delauney triangulation to determine neighbourhood (Fig. 3, b)).These clouds readily overlap with the values of the ideal HCP lattice and the cloud of the noisy BCC-structure for this moderate noise.Thus the q 4 -and the q 6 -values in a) and b) give the impression that the sampled structures consist of individual sub-populations with distinct local structures, in contrast to the homogeneous creation protocol and the continuously distributed displacements used to mimic the noise.
The reason for this behaviour is the discretized neighbourhood-discrimination.With noise, a neighbour may be just below or above the threshold distance to be counted as a neighbour and the q l -value jumps between distinct values (Fig. 3, a)).Similarly, in Fig. 3, b), minute displacements of the lattice sites leads to creation or annihilation of facets, and the q l -value again jumps between distinct values.The BCC structure is particularly sensitive to the first mechanism because of the close distance between the nearest neighbours (r = 1d) and the second nearest neighbours (r = 2/ √ 3d ≈ 1.15d).Any threshold chosen between those values will inevitably lead to misidentified neighbours once the first and second peak of the radial distribution function overlap (see Fig. 1,a)).The FCC structure is more sensitive to the second mechanism because its Voronoi cell contains six vertices which split and form new facets upon infinitesimal displacements as described by Troadec et al. [22].
The α=1 facet-weighting of the neighbourhood is continuous, contrary to the first two definitions (Fig. 3 c)).Emergence of a new facet with displacement of a lattice site creates a new neighbour with zero weight, and vanishing facets continuously contribute less to the q l -value.However, even with this neighbourhood-weighting correct characterization of noisy crystal structures is difficult.The noisy BCC and FCC structures are constructed to be different, as can be seen from the radial distribution function (Fig. 1).But this structural difference cannot be extracted from the q 4 -and q 6 -values, as the distributions of both parameters overlap.
The result of a q 4 -q 6 -analysis using the α=2 facetweighted neighbourhood definition is displayed in Fig. 4. The analysed sites of FCC and BCC structures with noise generate two continuous clouds of points in the q 4 -q 6plane.This behaviour correctly displays the underlying homogeneous structures.In addition, the clouds of the BCC-and the FCC-structure are now clearly distinguishable, in contrast to the other neighbourhood definitions.Such a q 4 -q 6 -analysis consequently would allow for the analysis of crystallinity just by analyzing the distance of the noisy q l values from the values obtained for the ideal lattices.
The optimization of sensitivity of the q 4 -q 6 -analysis using an α=2-weighting can be seen for the whole tested range of noise in Figure 5.The mean value of all q 4 and q 6 values obtained from the structures at a given noise level are displayed.The distance among the q 4 -q 6 -value of the FCC-based structure and the BCC-based structure at each noise level changes continuously with noise and becomes minimal with the largest noise values tested.The distance among q 4 -q 6 -positions of the FCC-and the BCC-based structure, which indicates the distinguishability of the two structures, is larger for the α=2-weighting (BCC2) compared to α=1-weighting (BCC1).Interestingly, the noise mainly affects and lowers the q 6 -value, while the neighbourhood definition mainly affected the q 4value of the BCC lattice sites (see Fig. 2).The analysis of structures is not limited to q 4 and q 6 .These are commonly used, as they are especially sensitive to cubic or hexagonal symmetry, respectively.Mickel et al.
suggested that the q 2 bond orientational order parameter could be of particular interest in studies on crystallinity.The q 2 parameter vanishes by any kind of mirror symmetry of the Voronoi cell [10].It is thus not as useful as the q 4 and the q 6 in distinguishing between different crystalline structures.However, a vanishing q 2 serves as a good indicator of the presence of crystallinity.
We illustrate this behaviour in Fig. 5, where the change of average q 2 with noise is displayed.q 2 grows with increasing noise for all the neighbourhood definitions and vanishes for vanishing noise.The growth of q 2 with noise is not continuous and monotonous in general, but with the α=2-weighting.Care has to be taken especially with applying the α=0-weighting neighbourhood definition, i.e.Delaunay triangulation, to the FCC structure.The q 2 parameter does not approach 0 continuously as this metric exhibits a discontinuous jump as soon as infinitesimal noise is present (see discussion above and in Troadec et al. [22]).

Discussion and Conclusion
The Steinhardt bond orientational order parameters q l are a resourceful way to quantify geometric arrangement on a local scale.The comparability among studies is presently limited by the ambiguity how neighbourhood is defined.We have studied different definitions and weightings of neighbourhood in analyzing noisy crystal structures and suggested an optimization.None of the commonly used neighbourhood definitions would allow for a discrimination of a noisy BCC and a noisy FCC structure, as was highlighted in Fig. 3 c) BCC Fig. 5: Evolution of the average bond orientational order parameter attained from BCC-and FCC-based structures as a function of applied noise.In a) the evolution of the average q 4 and q 6 can be followed.The colour indicates the neighborhood definition (purple: α=1, pink: α=2) and the symbols indicate the noisy crystalline structure.From the values of the ideal, noise-free lattices (large symbols, compare Tab. 1) the q 6 -values continuously decrease, while the q 4 of the two crystal structures become increasingly similar.The α=2-weighting emphasizes the differences in symmetry compared to the α=1-weighting, and the values obtained for the two noisy crystal structures differ more.The dotted line indicates the σ = 0.036-noise level discussed in the figures above.The evolution of the q 2 bond orientational order parameter with noise is displayed for the different neighbourhood definitions in the FCC-based structure, b), and the BCC-based structure, c).q 2 vanishes for both ideal, noise-free crystal lattices, but linearity and monotony of the growth of q 2 with noise depend on the neighbourhood definition.
neighbourhood already allows for an unambiguous definition of neighbourhood and continuous values of the q l .However, this definition weights next-nearest neighbours in non-densest structures too heavily, consequently the symmetry of non-densest structures appears closer to the symmetry of densest structures.
This drawback of the α=1 facet-weighting is removed to large extent by using an α=2 facet-weighting.The larger facets associated with the nearest neighbours are emphasized.By this non-densest structures can be distinguished from densest structures, an improvement not limited to BCC structures.In consequence, noisy crystalline structures can be more easily identified by mapping regions in a q 4 -q 6 diagram to the ideal values (see Fig. 4).
Application of α > 2 weighting is also be possible.The neighbourhood metric stays unambiguous, continuous and emphasizes the nearest neighbours.However, we did not find a measurable improvement regarding the ability to distinguish noisy FCC and BCC structures for higher α values.
The q 2 bond orientational order parameter is useful for distinguishing amorphous and crystalline structures as it vanishes for a crystalline lattice.We show that the growth of q 2 with applied noise is continuous and monotonous for the suggested α=2 facet-weighting.

Fig. 1 :
Fig. 1: Radial distribution functions g(r) of a) BCC-based and b) FCC-bsed structures with noise of σ = 0.01d and σ = 0.036d as they are used in the analysis (see Sec. 2 for details of generating the structures).The distance r is scaled by the lattice constant d.The dashed lines indicate the peak positions of the noise-free ideal lattice.The structures with noise still exhibit long-ranged positional correlations and distinguishable local environments.

Fig. 3 :Fig. 4 :
Fig.3: Bond orientational order parameters of a BCC-based structure (orange) and a FCC-based structure (green) with a noise of σ = 0.036d.In a) the neighbours are determined by a threshold value, in b) a Delauney tringulation is used (α=0-weighting), in c) the Delauney neighbours are additionally facet area-weighted (α=1-weighting).The graphs on the x-and y-axis give the distributions of the q 4 -and the q 6 -parameter individually on a linear scale.The large symbols indicate the q l values of the ideal lattices (see Tab. 1).
. An arbitrary threshold and the Delauney triangulation (α=0) introduce a discretized neighbourhood, and consequently the investigated homogeneous structures appear to bear different distinguishable local structures.The morphometric or α=1 facet-weighting of