Quantifying alignment in carbon nanotube yarns and similar two-dimensional anisotropic systems

The uniaxial orientational order in a macromolecular system is usually speci-fied using the Hermans factor which is equivalent to the second moment of the system's orientation distribution function (ODF) expanded in terms of Legendre polynomials. In this work, we show that for aligned materials that are two-dimensional (2D) or have a measurable 2D intensity distribution, such as carbon nanotube (CNT) textiles, the Hermans factor is not appropriate. The ODF must be expanded in terms of Chebyshev polynomials and therefore, its second moment is a better measure of orientation in 2D. We also demonstrate that both orientation parameters (Hermans in three dimensional (3D) and Chebyshev in 2D) depend not only on the respective full-width-at-half-maximum of the peaks in the ODF but also on the shape of the fitted functions. Most importantly, we demonstrate a method to rapidly estimate the Chebyshev orientation parameter from a sample's 2D Fourier power spectrum, using an analysis program written in Python which is available for open access. As validation examples, we use digital photographs of dry spaghetti as well as scanning electron microscopy images of direct-spun carbon nanotube fibers, proving the technique's applicability to a wide variety of fibers and images.


| INTRODUCTION
The physical properties of anisotropic materials such as carbon nanotube (CNT) yarns, 1 carbon fibers or drawn polymer fibers and films are highly correlated with the degree of internal orientation, which can be increased by post-treatments such as mechanical drawing. For example, stretched (and reinforced) electro-spun cellulose fibers show great improvements in tensile strength and thermal stability, and reduced bursting. 2 The use of oriented dye molecules in luminescent solar concentrators increases emission into preferred waveguide modes, thereby increasing their optical quantum efficiency. 3 CNT buckypapers, films and fibers aligned in-situ or by post-processing exhibit considerable enhancements in their electrical and thermal conductivities, carrier mobility, current carrying capacities, tensile strengths and elastic moduli. [4][5][6][7][8][9][10][11][12] The degree of internal alignment is therefore considered an important metric to understand structure-property relations in macroscopic materials. Popular methods to estimate alignment rely on measuring some form of anisotropy in material properties along and perpendicular to the alignment axis such as: birefringence and linear dichroism 13,14 ; fluorescence polarization 15,16 ; and (specifically for CNTs) polarized Raman spectroscopy, 4,14,17 or anisotropy in electrical and thermal conductivities. 4,12 However, anisotropy in material properties is a consequence of structural anisotropy, that is, it is not a direct measurement of macromolecular orientation. A true measure of orientational order in a bulk sample is obtained when the orientation of all sub-elements is averaged with respect to a principal axis of reference (or 'director' following the nomenclature of liquid crystal literature 18 ). The number of sub-elements, which may be molecules, crystals, polymer blocks or macromolecules, is usually large. Hence, their orientation is described by a continuous function known as an orientational distribution function (ODF) which, when averaged, gives quantifiable parameters from the moments of the ODFthese are called orientational order parameters (OPs). Hermans and co-workers pioneered the work on orientation parameters [19][20][21] in the late 1930sa brief history of these advances up to the present day is provided in Appendix I. The Hermans factor, which is equivalent to the average of the second Legendre polynomial, hP 2 i, is now widely accepted as an orientational order parameter suitable for systems with uniaxial symmetry. It is usually calculated from crystallite ODFs determined using wide/small angle X-ray diffraction (WAXS/SAXS). [22][23][24] Due to their high penetration depths in organic materials, X-rays can give information about 3D orientation distributions; however, to fully analyze such distributions, the traditional pole figure method is preferred. 23,25 By contrast, most WAXS/ SAXS plots are not truly 3D since they are obtained by projecting the 3D intensities onto a 2D surface of film/ detector. Nonetheless, the two are mathematically related. 26 Our interest lies in the determination of the degree of uniaxial orientation in aligned CNT films, mats, or textiles, as again, the alignment is key to several physical properties. The X-ray diffraction method has also been widely applied in CNTs, but this technique is relatively slow to apply, particularly, when large number of samples need to be analyzed. An approach involving the analysis of spatial frequencies in scanning electron micrographs (SEM) by Fourier transforms, in development since at least 1998, offers an alternate, easy and fast technique for determining orientation. 14,[27][28][29][30] This technique involves taking circumferential scan of pixel intensities from the Fourier transform of an SEM image to obtain an ODF which has been used to calculate Hermans orientation parameter in CNT films. 30,31 However, bulk CNT mats or textiles are in most cases formed by the assembly of CNT layers via a wind-up process such that CNTs are restricted to individual layers. Spin-coated CNT films and those produced by similar methods are also limited to a single plane on the surface of substrates except for vacuum-filtered films where CNTs might project out-of-the-plane. The application of a 3D orientation parameter like hP 2 i to CNT mats may therefore give a misleading measure of orientation, especially when the ODF is obtained from the Fourier analysis of an SEM image. As an alternative, our group introduced the idea of using a 2D orientation factor based on the Chebyshev polynomial, T 2 , initially for deposited CNT films, 27 and recently for layered CNT films. 32 In the present work, our effort is two-fold. First, we discuss the 2D orientation parameter in full mathematical detail, deriving an ODF that aptly describes orientation of macromolecules in uniaxial planar systems and showing that its second moment is the average of Chebyshev polynomial, hT 2 i. We discuss the necessity for using a function to fit experimental data and how the orientation parameters relate to shape and full-width-at-half-maximum (FWHM) of the fitted functions. Secondly, we demonstrate a program that rapidly estimates Chebyshev orientation parameter from SEM images using the Fourier transform approach. While this article focuses on the specific case of CNT films, we anticipate that this program will be useful to researchers working on various anisotropic films; it is available freely for download from our Github page 33 or directly from authors on request.

| Orientational geometry in 3D and 2D
The orientation of a macromolecule in three dimensions is usually defined using three Euler angles (ϕ, θ, ψ) formed between two Cartesian coordinate systems XYZ and x 0 y 0 z 0the former is taken as the laboratory/sample frame while the latter is fixed in the macromolecule (see Figure 1(a)). Here, ϕ is the azimuthal angle formed between X and the projection of z 0 axis on the XY plane varying from 0 to 2π, θ is the polar angle formed between the Z and z 0 axes varying from 0 to π, and ψ is the 2π rotation of the macromolecule about its z 0 axis. The orientation of all such macromolecules in the sample is then defined by the orientation distribution function (ODF), f (ϕ, θ, ψ). This function gives the probability of finding a macromolecule with an orientation in the range (ϕ, θ, ψ) to (ϕ + dϕ, θ + dθ, ψ + dψ). 34,35 Since the ODF is a 3D probability density, it must be a positive function, that is, f(ϕ, θ, ψ) ≥ 0 and normalized to unity, as shown in Equation (1). The term sinθdϕdθdψ in this equation is the volume element for a sphere of unit radius in the spherical polar coordinate system. 36 Relevant to this work are cases of high symmetry, such as films or drawn fibers, where the macromolecules are not only symmetric about their own z 0 axis but also oriented relative to some preferred direction, usually denoted as the sample Z axis. This is known as uniaxial symmetry. Depending on the geometry, the following subcategories are possible: (a) 3D fiber symmetrywhere the macromolecules are symmetrically distributed in an azimuthal direction, as illustrated in Figure 1(b). The ODF f(ϕ, θ, ψ)then reduces to f(θ) 3D .
(b) 3D planar symmetryrepresented in Figure 1(c); here, the macromolecules mainly lie in the XZ plane, but azimuthal rotation (out-of-plane) of the macromolecules is possible in principle, so f(ϕ, θ, ψ) = f(ϕ, θ) 3D . As a special case, when the layers are all approximately parallel, ϕ $ 0 and so the ODF reduces to f(θ) 3D . However, since the ODF is still described in 3D, it retains the sinθ term in Equation (1). An example of this situation is that of Xray diffraction where although a 3D intensity distribution is projected onto a 2D plane, diffraction occurs in 3D. 26 (c) 2D planar symmetryalso known as uniaxial planar distribution 37 ; here, as shown in Figure 1(d), the macromolecules are confined such that the measurement of the ODF strictly comes from a 2D plane. For example, a 2D Fourier transform of an SEM image is equivalent to diffraction from a plane. Specification of orientation in this case does not require the three Euler angles, just the angle ϕ between the Z and z 0 axes is sufficient; the ODF can be written as f(ϕ) 2D and the reference direction Z is now the molecular plane. Importantly, the 2D area element does not include the sinθ term in Equation (1) which alters how the ODF is defined in 2D. Rewriting the normalization condition for the 3D and 2D cases, we have: For 3D, For 2D,

| ODFs in 3D and 2D
We now proceed to discuss the form of the ODFs in 3D and 2D. Van Gurp showed that a suitable ODF can be constructed either by a trial-and-error method or by expressing the 3D ODF as a series of even Legendre polynomials P n (cosθ). 34 The coefficients or moments of this expansion, a n , can be shown to be averages of the Legendre polynomials themselves, that is, a n ¼ 2nþ1 2 P n cosθ ð Þ h i. The expression for the 3D ODF can then be written as shown in Equation (4) (see sections S1, S4, S5 in SI for definitions and detailed derivation).
To understand the above expression, we consider an expansion of up to the first two even moments as shown in Equation (5) and simplify it further by substituting P ¼ i . Compared to the function cos 2 θ, the ODF in (5) satisfies the conditions for being a 3D ODF, that is, it can be normalized to 1 (see section S6 in SI) and is positive for P in the range, À0.5 ≤ P ≤ 1.
However, Equation (5) does not satisfy the normalization condition for being a 2D ODF as given by Equation (3) (see section S6 in SI for a proof)so it is not a valid 2D ODF. Below, we describe how a suitable 2D ODF can be derived following van Gurp's series expansion method. 34 First, using trial and error, we guess an ODF of the form described in Equation (6) that is normalizable in 2D.
Equation (6), unlike (5), does not contain the second Legendre polynomial, (1/2)(3cos 2 ϕ À 1), rather, a term (2cos 2 ϕ À 1) which resembles the second Chebyshev polynomial, T 2 (cosϕ). Therefore, we express the 2D ODF in a series expansion of even Chebyshev polynomials of the first kind as shown in Equation (7). The coefficients of this expansion can be obtained using the orthogonality conditions of Chebyshev polynomials and using the definition of 2D average of a function (see sections S2-S4 in SI for definitions and detailed derivations). Analogous to the 3D case, the coefficients turn out to be the average of Chebyshev polynomials as shown in Equation (8) which leads to final form of the 2D ODF in Equation (9).
The simple ODF in (6) can be shown to be a special case of (9) when the latter is expanded up to first two even terms, with T 0 = 1, T 2 = 2cos 2 ϕ À 1 and T = 2hT 2 i. It is to be noted, however, that the ODFs in Equation (5) and (6) represent broad distributions; higher order coefficients, hT n i or hP n i, with n going up to 100 may be required to describe and differentiate highly aligned samples. 34,38

| Orientational order parameters in 3D and 2D
Since all orientational information is contained within the moments of the ODF, in practice, knowing the true ODF is not essentialan indirect measurement in the form of some intensity distribution, I(θ) or I(ϕ), is sufficient. As discussed before, intensity distributions for CNT films can be obtained via SAXS (3D distribution) or via our preferred method, Fourier transform of an SEM image (2D distribution). Once the intensity distribution is known, orientation order parameters (OPs) which are the moments of the ODF, hP n (cosθ)i or hT n (cosϕ)i can be calculated. The most widely used parameter is hP 2 (cosθ)i (or hP 2 i for short) which is known as the Hermans orientation parameter (HOP).
For the 2D case, we recommend hT 2 (cosϕ)i (or hT 2 i) which we call the planar or Chebyshev orientation parameter (COP). These two OPs can be calculated using Equations (10) to (13).
The denominators in Equations (11) and (13) are normalization factors introduced so that the measured intensity distributions satisfy the normalization condition for ODFs. We reiterate that the HOP and COP are moments of the 3D and 2D ODFs, respectively. So, using a 3D orientation distribution (e.g., from X-ray) to calculate the COP or using a 2D orientation distribution (e.g., from SEM) to calculate the HOP will lead to under-or overestimation of the orientational orderthe two operations are therefore strictly incorrect.
The HOP has values in the range À0.5 ≤ hP 2 i ≤ 1 for macromolecules aligned perpendicular to and along the reference direction, respectively. By comparison, the COP has values in the range À1 ≤ hT 2 i ≤ + 1 for macromolecules arranged perpendicular to and along the reference direction, respectively. The Hermans and Chebyshev parameters are zero when the macromolecules are randomly orientated with respect to the reference axis (corresponding to the horizontal axis in an SEM image), and also if all the molecules are perfectly aligned at 54.7 or 45 , respectively, to the reference axis. However, this ambiguity can be resolved by defining hP 2 i or hT 2 i to be the highest eigenvalue of the orientation tensor, that is, the orientation is defined with respect to a new axis which is taken along the direction with the highest intensity in the Fourier space. With this change of reference axis, the orientation parameter values now represent alignment of the macromolecules. For example, when all macromolecules are oriented in the 45 direction, and the 45 direction is defined as the reference direction, then hT 2 i is unity, which is what we expect when all macromolecules are oriented along the same direction. We, therefore, recognize "alignment" of macromolecules in a sample as the degree of orientation and not the direction of orientation. In the case of continuously drawn CNT films, textiles or fibers, the predominant orientation originates from the direction of drawing which is the fiber/ textile axis. The need for a rotated reference axis would only originate from misaligned mounting in the SEM.

| Dependence of OP on shape and width of the distribution
The intensity distribution as measured over the range of 0 to 360 is a continuous periodic distribution with peaks that are 180 apart. However, it is usual practice to analyze/fit only one of the peaks of the distribution which: (1) is equivalent to treating the boundaries as discontinuous, and (2) leads to a poor estimation of the baseline. As we will see below, the background not only includes noise but also contributions arising from the overlap of neighboring peaks. In this work, therefore, we fit the complete intensity distribution using a cumulative function that includes all individual peaks and a proper linear background (see also section 2.7 for discussion on baseline correction).
While there is a consensus that narrower peaks imply a higher degree of alignment, the dependence of OP on the exact shape of the intensity distribution function (IDF) used to fit the measured intensity profiles is the subject of ongoing discussion. 23,26 Most reports use a Gaussian distribution (GD) to fit intensity profiles, but Lorentzian distributions (LD) and combinations of the two (Pseudo-Voigt) have also been used. 23,26,31,32 Recently, a generalized normal distribution (GND) has been shown to best fit SAXS intensity distributions from multi-walled CNT (MWCNT) forests. 26,39 This function is characterized by an additional fitting parameter called the shape or sharpness factor, β, which takes values of 2 for Gaussian (broad) and 1 for Laplace distributions (sharp), respectively. "Sticky" or cohesive granular materials like MWCNT arrays have β in the range 1.37 to 1.65 depending (linearly) on the film's volumetric density. 39 The functional forms of the Gaussian (GD), Lorentzian (LD), and generalized normal distributions (GND) are presented in Equation (14) to (16).
Here μ is the peak center, γ is the half-width, σ is SD ) and α is a factor that is proportional to the half-width (α ¼ γ= ffiffiffiffiffiffi ffi ln2 β p ); Γ is the gamma function. 26 In contrast to Vainio et al., 26 but in agreement with Gbordzoe et al., 23 we found that OP values depend on the shape of the distribution. We generated model continuous IDFs using GD, LD and GND distributions; two sets of these distributions with full-width-at-half-maximum (FWHM) values of 10 and 100 are shown in Figure 2 (a)-(d). The distributions in Figure 2(a),(c) represent orientation along the reference axis while those in Figure 2 (b),(d) represent alignment perpendicular to the reference axis. The individual peaks within the distributions can overlap, depending on their width and shape of their tails, resulting in nonzero backgrounds (see the trough regions in Figure 2(c),(d)). Therefore, fitting the entire distribution is essential to get accurate estimation of orientational order. Note that the OPs are still determined by integrating the distribution from 0 to 180 , or equivalently 180 to 360 , as shown in Equation (11) and (13); For experimental data, averaging the two halves of the ODF prior to integration improves the signal-to-noise ratio.
OP values, calculated as described above, depend on the FWHM of the peaks in the IDF. The calculated OP values for model IDFs of varying widths are plotted in Figure 2(e). The two branches with positive and negative OP values are obtained when the orientation angle is fixed parallel, and perpendicular to, the reference axis, respectively. Also, note that a FWHM greater than 180 has no physical significance -OP values for such widths are shown simply for completeness. We see that for GD, LD and GND functions, the OP values decrease with increasing peak width; though the actual values depend on the functional shape, with Gaussian distributions giving the highest values and Lorentzian the lowest. The GND function (with β = 1.5) produces values that are intermediate to those of the Gaussian and Lorentzians, closely matching Gaussian values at high orientations and those of Lorentzian at low orientations. In addition, for GNDs, OP values also depend on their shape factors as shown in Figure 2(f) (see also Figure S2(a) in SI for a similar trend in HOP). This dependence of OP on a function's shape suggests that comparison of OP values of different samples holds meaning only if the same fitting function is used and, particularly for GNDs, if the same shape factor is used.
Further, the nonlinearity of OPs with the function's FWHM suggests that the sensitivity of OPs to changes in orientation is lost at high alignments (widths less than 25 to 30 ), particularly when Gaussian and GND line shapes are used to estimate OP. We, therefore, recommend using Lorentzian distributions for estimating hT 2 i. A Gaussian distribution is still the best option for estimating hP 2 i, since extremely narrow line widths are required for Lorentzian distributions to produce a hP 2 i of 1, which is not realistic situation in practice.
When studying isotropic samples, or those with two prominent directions, additional peaks can appear in the intensity profiles, as shown in Figure 3(a). COP values then depend on the relative intensities of the secondary (A 2 ) and primary peaks (A 1 ). As shown in Figure 3(b), COP values decrease with increasing intensities of the secondary peaks; when the intensities are comparable, that is, A 2 /A 1 = 1, even sharp distributions can have an orientation value close to zero. This effect of secondary orientation peaks is stronger in case of hP 2 i, but does not alter higher order OPs such as hT 4 i or hP 4 i (see Figure S2-S4 in SI), as by

| Estimating OP using "FibreCOP"
We have developed a program for rapidly determining orientation parameters in aligned CNT textiles and similar anisotropic systems. The code, written in Python 3.7, is available from our Github repository. 33 This program determines orientation parameters from SEM images via the following sequence of steps: (1) image analysis which includes cropping, optional de-noising, binarization and thresholding, generating a power spectrum via Fourier transformation, applying a circular mask, and calculating a radial intensity profile to get the ODF; (2) filtering noisy data and peak-fitting ODFs using various functional shapes; (3) calculation of FWHM, COP and HOP (although the latter is not recommended for 2D ODFs), as well as the fourth moments. An option to analyze only intensity distribution data (in the absence of an image) is also provided.
To test the program, first, the model IDFs generated in Figure 2(a),(c) were fed to the program as data (not as an image) and the orientation values obtained from the fitted distributions were compared with those calculated directly from the model. COP (FWHM) values for the fitted Lorentzian IDFs were obtained as 0.848 (10 ) and 0.196 (100 ), respectively, which closely matched those calculated from the model, with rounding errors seen only in the fourth significant digit. Secondly, we generated model images from photographs of spaghetti oriented in different directions (see Figure 4 (row (a)). Results of subsequent analyses steps such as binarization, Fourier transformation, masking and peak fitting are depicted in Figure 4, rows (b)-(e); further details can be found in the Methods section. The peaks in the radial distributions for aligned spaghetti (see the first three graphs in Figure 4(e)) have sharp tips with a wide base which require a narrow Lorentzian and a broad Gaussian, respectively, to produce good visual fits as well as lower AIC (Akaike Information Criterion) values. Note that while comparing different nonlinear models, a fit with lowest AIC value is the best fit (see Methods for detailed explanation). 39,40 The distributions, in total, require four or five peaks for an accurate fit. The overall COP values (see Table 1) for the three aligned samples are in the range 0.74 to 0.8; although these values are high, the broad Gaussian tails lower the values from $0.91 expected for narrow Lorentzians with widths less than 3 . The negative signs for the horizontal (vertical in Fourier space) and diagonal samples are due to the choice of reference axes which in digital Fourier image analysis is along the horizontal direction (3 o'clock position). To obtain only positive OP values, the images can be rotated prior to analysis (see Methods section for further details). In comparison to the aligned samples, the randomly oriented fibers have very low COP value close to zero, as expected for a truly isotropic sample. We now consider a final example of the crisscross configuration, as shown in last column of Figure 4, which too has a hT 2 i or COP value close to zero. The ODF for crisscross configuration is similar to those shown in Figure 3(a), that is, it has secondary orientation peaks and so, the overall OP reduces to zero. To distinguish the crisscross and isotropic cases we must consider the biaxial orientation term, hT 4 i which, as listed in Table 1, is far higher for the crisscross sample. As a final test, we analyzed an SEM image of a CNT film produced by the direct-spinning method (see Methods). The SEM and its Fourier power spectrum are shown in Figure 5(a); the Fourier spectrum is diffuse, suggesting that the sample is not highly aligned. Radial intensity profiles with different fitted IDFs are shown in Figure 5(b) and the calculated orientation values are tabulated in Table 2. Three peaks of equal FWHM, centered at 0 , 180 and 360 are sufficient to fit the intensity profile; the orientation values, however, depend on the type The GND produces the best fit because of the additional free parameter, β, as it captures the sharpness of the tip, width of the peak and curvature of the tail, adequately.
Since the widths and β of the GND are constrained to be the same, the total number of fit parameters is 10, just one more than the other fits; the AIC for this fit is the lowest at 10,700. However, as β can vary between samples (see Figure 2(f)), OPs do not hold comparative value. T A B L E 2 Orientation parameters for the example CNT image in Figure 5(a) represented here as spot a. SEM images for spots b -D, are presented in Figure S5 in SI We, therefore, recommend using Lorentzian distributions as fits for CNTs, in general. SEM images of CNTs are usually captured at brightness (B) and contrast (C) settings visually pleasing to the operator. To test if these settings have any bearing on OP values, we captured images at different B and C percentages, as shown in Figure 5(c), varying between very poor contrast to extremely bright images. COP values (using LD fits) for these images are presented in Figure 5(d). As we can see, the variation in COP values within each image and at different settings is low with an overall value of 0.44 ± 0.01. In comparison, higher variation in COP values is seen at different SEM spots on the same sample (see Figure S5 in SI for images, data is tabulated in Table 2). For example, at one random position, d, in the sample, COP value was very low at 0.39. As a best practice, it is therefore advisable to analyze different regions on the same sample to arrive at an average orientation value.

| Application in CNT films
A classic example of orientation analysis is to study alignment changes in polymers as a function of draw ratio (a measure of stretching). Here, we take the example of CNT fibers grown by direct-spinning method and collected as films at two different spinning rates, 30 m min 1 and 60 m min 1 ; SEM images of these samples are shown in Figure 6(a),(b) and the respective Lorentzian fits are shown in respective insets. COP for the sample spun at 30 m min 1 is 0.44 which increases to 0.54 at 60 m min 1 , with corresponding changes in FWHM from 50.5 to 37.3 , indicating a significant increase in CNT alignment at faster spinning rates. We again note that while the GND gives a better fitting, the free fitting parameter, β, varies between the images. In this example, β takes values 1.4 and 1 for the slower and faster spinning rate samples and the COP values with GND fits turn out to be 0.68 for both images. Therefore, the utility of GNDs in determining OP is limited but is rather useful in understanding confinement or coupling of CNTs, as suggested by Vainio et al. 26

| Further discussion
The FWHM of an ODF as a measure of orientation: It is well accepted that narrower distributions have higher OP values, and this is also seen from Figure 2(e). Several CNT reports already rely on simply using FWHM as it is purportedly more stable and does not require complicated peak fitting operations. 23 However, (1) FWHM does not capture information regarding the peak sharpness or tail features of an ODF, (2) for noisy data, peak fitting is necessary to estimate FWHM, (3) peaks in a continuous ODF may have varying FWHM unless constrained to be equal, and (4) when fitting multiple functions to a peak (Figure 4(e)) or when secondary orientation peaks are present (see Figure 3), the correlation between OP and FWHM is lost. Therefore, FWHM is not always a good measure of macromolecular orientation, although in certain limited cases, it can give an estimate of the anisotropy in samples. Choice of OP range: Figure 2(e) shows that hP 2 i is not symmetric about zero while hT 2 i is, that is, when the alignment direction is along the perpendicular to the reference axis, the range for hP 2 i is compressed between 0 to À0.5, rather than 0 to 1. To get a bigger range of values, it is best to use the positive scale for hP 2 i and as a good practice, for hT 2 i as well. The asymmetric range of hP 2 i also poses a problem when using the Fourier transform approach to obtain IDFs. For a system containing anisotropic objects oriented about a single axis, the Fourier transform will produce an intensity distribution with maxima at 90 to this axis, as seen from Figure 4(c). Therefore, either the image should be rotated prior to analysis or the OP transformed after analysis. Mathematically, under a rotation by π/2, hP 2 i does not transform as h i which is valid only at the extremes. The trigonometrically correct transformations for the various order parameters are: These relations are neither well-known nor intuitive, so we have included a proof in section S8 of the SI. Note that in the rotated transformation, the range of hP n 0 i also changes: the unaligned state is represented by 0.5 and not 0. Therefore, it is not advisable to mathematically transform HOP values from the negative to the positive scale, rather it is best to rotate the image itself before analysis. This rotation can be achieved via postprocessing of the image or capturing the images with the principal axis about the perpendicular direction; our program includes a rotation option to handle this issue. Baseline correction and peak smoothing: As mentioned before, we found that orientation values depend on the proper choice of background. In Figure 2(c),(d), the minima of LD and GND functions are nonzero which is due to the overlap of adjacent peaks. The minima must not be arbitrarily reduced to zero, rather, a linear fit with a nonzero slope must be included in the fitting model and ultimately subtracted from the cumulative fit before calculating OPs. Adding a linear background to the fit model takes the number of fit parameters required for fitting continuous distributions to 11 (or 14 for GNDs): two for the linear fit and three (or four) for LD and GD (or GND) used. By constraining the widths (and the shape factor for GNDs) to be equal, the number of free parameters can be reduced to 9 (or 10 for GNDs). Peak smoothing operations such as Savitzky-Golay 41 are alternatives to fitting which can be applied to the usual noisy distribution data to obtain the form of the IDF. However, smoothing is akin to using a function with an additional fitting parameter which changes between samples. Therefore, values of OP from smoothed ODFs do not hold comparative value and are best treated as upper or lower bounds.

| CONCLUSIONS
We have demonstrated a method to facilitate the estimation of alignment in anisotropic materials, such as carbon nanotube yarns and similar aligned systems, using an analysis program which we have made available online. This method relies on analyzing photographs or scanning electron micrographs of the samples to arrive at an ODF which is used to calculate order parameters. We have shown, using a series expansion method, that the appropriate orientation order parameter for anisotropic samples that are 2D or have a measurable 2D ODF is the Chebyshev order parameter, hT 2 i. When defined with respect to the highest intensity direction in the Fourier space as the reference axis, hT 2 i can be used to obtain a measure of structural alignment. We found that the order parameter values depend on both the width and shape of the function used to fit the continuous ODF, hence the order parameters must always be quoted along with the fitting function used. We have demonstrated the working of our analysis program on both macroscopic (spaghetti) and microscopic (CNT) fibers, therefore, we anticipate that this method will find widespread application in various anisotropic systems.

| METHODS
SEM images were obtained on a FEI Nova SEM at 5 keV. The photographs of spaghetti were obtained using a mobile phone camera with its default settings which could include software-based picture corrections. Carbon nanotube fibers were produced via the direct spinning floating catalyst chemical vapor deposition method; details of procedures can be found in our other works. 32,42 A program, "FibreCOP", to calculate orientation parameter from SEM images/photographs of aligned fibers has been developed. The code is written in Python 3.7 and uses standard libraries including OpenCV, numpy, matplotlib, scipy, pandas, lmfit and tkinter; it has a simple graphical user interface making it intuitive for users with little Python knowledge and is available for download and contributions from our Github page. 33 The main analyses steps of the program are as follows: With the original image size in pixels as reference, the image is first cropped into a square whose side is equal to 256n pixel dimensions, n being the highest integer possible (for example, a digital image with a size of 1536 Â 1026 is cropped to 1024 Â 1024 pixels). The remaining area is scanned by shifting the crop area in a finite number of steps, if needed; each square area is then analyzed and used to calculate an average orientation parameter. For further analyses, the images are converted to a binary format using OTSU or Gaussian thresholding methods (Figure 4, row (b)). 43 2D discrete Fourier transformation of the binary image produces a power spectrum ( Figure 4, row (c)) where the direction of orientation is always perpendicular to alignment direction in the image space. A circular mask is then applied (Figure 4, row (d)) before calculating the radial intensity distribution which is a plot of summed pixel intensities along a radius vs. the angle the radius makes with a reference direction (clockwise starting from a 3 o'clock position). To calculate alignment the program shifts the highest intensity peak to either the 3 or 12 o'clock position. For directions less than 45 , the peak is shifted to the 3 o'clock position (positive COP), but for higher angles, the distribution is shifted to align with 12 o'clock position giving negative values. As a good practice, we recommend using the positive OP range, and therefore, an option to rotate the image prior to analysis is provided in the program. The step size for radial scan can be controlled by the user; we recommend a value less than 1 / step (preferably 0.25 /step)this is a choice that will be guided by need for low computation time against data resolution required. The power spectrum from SEM images is usually noisy with a spike at 45 arising from pixel edges-a median filter is applied to remove this noise prior to peak fitting. We use the lmfit package for peak fitting, so fit statistics such as χ 2 , AIC and BIC (Bayesian information criterion) are outputted by the program for comparing models. AIC provides a method to compare nonlinear models by penalizing the model's goodness of fit for increased number of fit parameters. It is defined as AIC ¼ À2lnL þ 2K, where K is the number of free parameters andL is the maximum likelihood estimator for the model. 40 For linear models, it is also formulated as AIC = NlnSS E + 2K where SS E is residual sum of squares. The relative likelihood of different models can be compared by Δ(AIC) = exp((AIC min À AIC i )/2) where AIC min smallest of the AIC obtained. The program outputs both Hermans and Chebyshev orientation parameters, as well as FWHM and fourth momentshowever, we reiterate the point discussed in main text that calculating HOP from power spectra of SEM images is not a valid approach.
APPENDIX I Historical development of orientation parameters: The development of orientational order parameters began with the work of Hermans and co-workers, in the late 1930s and early 40s. [19][20][21] Studying birefringence in amorphous polymers and X-ray diffraction in crystalline polymers, they showed that (1) a measure of orientation in uniaxial fibers could be obtained from the ratio (f H ) between the differences in polarizability along and perpendicular to the reference fiber axis and the polymer chain axes, respectively, and (2) this ratio is related to the average angle, θ, formed between the fiber axis and a crystallographic direction in the crystallites as f H = 1 À (3/2)hsin 2 θi; the angle itself is determined by averaging the diffracted X-ray intensities. 21,44 Several researchers refined this idea of uniaxial fiber orientation through 1940s and 60sa detail of these developments is presented by White and Spruiell, 37,45 and van Gurp. 34 Notably, Muller in 1941, 34 suggested that the X-ray intensity distribution can be expressed as a function of spherical harmonics, in particular, the Legendre polynomials P n (cosθ). He identified the Hermans factor as being the second Legendre polynomial, P 2 (cosθ) = (1/2)(3hcos 2 θi À 1)this variant is now the widely used form of f H . 34,45 Later researchers were concerned with determining the hcos 2 θi part of f H : for example, Wilchinsky showed how the term hcos 2 θi, which is the 3D average over a coordinate sphere, could be determined from X-ray measurements even in the absence of diffracting planes normal (but available in other directions) to the axis. 46,47 In a series of works in 1960s, Roe and Krigbaum further explored Muller's idea, showing the equivalence between a crystallite ODF and the plane normal distribution obtained from X-ray measurements on the crystallite. 35,48 Roe concluded that if an accurate and complete description of the ODF is not required, then it is sufficient to determine the first few moments of the ODF. 49 The exact number of moments needed to describe the ODF depend on the type of measurement made, as shown by McBrierty and Ward: P 2 is sufficient for optical birefringence, infrared absorption and ultrasonic measurements; while for magnetic resonance, fluorescence polarization and elasticity measurements, moments up to P 4 may be required. X-ray diffraction ODFs give greater number of moments, thus providing a complete description of the orientation, 50,51 particularly if the anisotropy is large. For example, in a study on isotactic polypropylene, higher moments up to P 100 were required when the draw ratio exceeded eight (corresponding to a phase transformation from lamellar to fibrillar). 38 Researchers in 1960-80s were concerned with biaxial orientation in polymer films, a case that might happen in blown films or those stretched biaxially. Stein defined an orientation factor of the form f P = 2hcos 2 θi À 1 for each of the three axes in a coordinate sphere. 52 This was adapted by Samuel, 53 and later, White and Spruiell, 37,45 to represent the orientation of macromolecules restricted to a single plane. The latter authors derived the uniaxial planar/2D orientation factor, f P , using a polarizability tensor-based approach (in the current work, we use a series-expansion method) and further showed that it is a special case of biaxial orientation factors.
The use of Hermans orientation factor and small-angle X-ray scattering profiles for estimating orientation has now become routine. Several groups have also used this method to determine orientation in CNT films. 22,23,26 In comparison, the use of the planar orientation parameter remained limited to the field of liquid crystals where it is known as the 2D scalar nematic order parameter, S . The analysis of scanning electron microscope images using digital computation of the Fast Fourier transform (FFT) triggered the use of 2D orientation parameter in CNT literature. The earliest mention of T 2 as a measure of orientation for CNTs seems to be in a work from our group in 1998. 27 An SEM image of CNTs dispersed on a substrate was Fourier transformed to obtain a power spectrum, which was used to calculate hT 2 i. A work on liquid crystalline solutions of CNTs in 2008 estimated a parameter S (equivalent to hT 2 i) from the dichroic ratio. Unrelated to these works, the method of Fourier transformation of an SEM image found other applications. Ayres et al., in 2008, 28 applied the FFT method to study orientation in polymeric scaffolds made of electrospun fibers. They relied on the normalized height (to a baseline of zero) of ODFs to estimate anisotropy (taller and narrower peaks represent higher alignment). On similar lines and around the same time, a Fourier transform misalignment analysis method was developed by Kratman et al. 29 to measure misalignment in carbon fiber composites. In 2016, our group used the FFT method in conjunction with Hermans factor to estimate alignment in CNT epoxy composites. 30 Our concern at that point with this approach was the validity of a 3D orientation parameter in case of 2D mats. Subsequently, in 2017, our group suggested the use of Chebyshev polynomial hT 2 i for layered CNT films, although without providing a detailed mathematical background. 32 Brandley et al. improved this method further and transformed it into a Fourier mapping technique to analyze orientation in large areas of thin CNT veils. However, they limited their method to determining hP 2 i. 31 The present work, via a series expansion method derives the expression for hT 2 i and explains why it is the appropriate parameter to be used together with the FFT image analysis method. We hope that our program simplifies the procedure so that the technique is adapted widely, both within the CNT community and further afield.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.