Evolution of Complex 3D Motions in Spicules

Ubiquitous transverse oscillations observed in spicular waveguides, identified as the kink wave-mode had previously been reported along with periodic structural distortions of the flux tubes, observed as cross-sectional width and associated photometric variations. Previous studies identified these perturbations as the observed signatures of concurrent kink and sausage wave-modes. High-resolution Hα imaging-spectroscopy data from the CRisp Imaging SpectroPolarimeter at the Swedish Solar Telescope are used to analyze the off-limb spicular structures. For the first time, the evolution of the resultant transverse displacement of the flux-tube structure, estimated from the perpendicular velocity components, is analyzed along with longitudinal, cross-sectional width, photometric, and azimuthal shear/torsion variations. The pulse-like nonlinear kink wave-mode shows strong coupling with these observables, with a period-doubling, -tripling aspect, supported by mutual phase relations concentrated around 0° and ± 180 ° . The three-dimensional ensemble of the observed dynamical components revealed complexities pertinent to the accurate identification and interpretation of, e.g., linear/nonlinear, coupled/uncoupled magnetohydrodynamical wave-modes in spicules.


Introduction
Motions in the solar convective zone can excite oscillations in thin magnetic flux concentrations that emerge from inter-granular lanes in the photosphere. This oscillatory behavior may become apparent in the dynamics of chromospheric magnetic flux-tube (MFT) structures observed by imaging-spectroscopy, which are often inferred as magnetohydrodynamic (MHD) wave-mode(s) (Zaqarashvili & Erdélyi 2009). Accurate interpretation helps to quantify the associated wave energy flux that is transferred across the Interface Region (chromosphere-Transition Region-lower corona) and to identify potential dissipation mechanisms. Mathematically, the MFTs are usually modeled as cylindrical waveguides. These structures are highlighted in the limb (ondisk) observations as intensity enhancements (depression), suggestive of higher plasma density when compared to the background. At the solar limb, these thin magnetic structures are often observed as short-lived jets, inclined at some angle to the normal. Depending on the location (on-disk/limb/active-/quiet Sun region) and the observed wavelength, these features are often classified as spicules, mottles, fibrils, straws, Rapid Blue/Red shifted Excursions (RBEs/RREs), and very recently as Slender Ca II Fibrils (SCFs; Gafeira et al. 2017b).
MHD wave theory permits the possibility of an infinite number of orthogonal wave-mode(s) in these structures, categorized in linear domain by azimuthal (m) and axial (k z ) wavenumbers in the thin-tube approximation (tube radius =observed length). However, due to observational limitations, so far only the three low-order fundamental modes (kink, sausage, and Alfvén) were reported in chromospheric features (see review by Verth & Jess 2016). However, it must be noted that the solar chromosphere is a highly dynamic, gravitationally stratified, nonlinear, inhomogeneous environment where plasma-beta (β) varies from large, across unity, to very small. Interaction of MHD wave-modes embedded in localized MFT structures with such inhomogeneities can give rise to physical mechanisms that involve resonant absorption (Hollweg 1988), phase-mixing (Heyvaerts & Priest 1983), mode-coupling (Fazel 2016), and/or conversion.
The bulk motions of spicules have been extensively examined and interpreted in terms of discrete MHD wave-mode(s). Transverse displacements of spicules as a consequence of the m=1 kink mode were reported in spicules (Kukhianidze et al. 2006;De Pontieu et al. 2007b;Ebadi & Ghiassi 2014;Tavabi et al. 2015), in their on-disk counterparts or RBE/RREs (Rouppe van der Voort et al. 2009), in mottles (Kuridze et al. 2012(Kuridze et al. , 2013, and in fibrils (Pietarila et al. 2011). Longitudinal (field-aligned) motions associated with mass flows were reported in mottles (Loughhead 1974), spicules (Pereira et al. 2012), and on-disk RBE/RREs (Sekse et al. 2013). A mix of propagating (upward/downward) and standing transverse waves along the spicule structures was also reported by Okamoto & De Pontieu (2011). Sekse et al. (2013) studied on-disk spicule counterparts and noted transverse displacements of the waveguides along with rotational motions as prime components of spicule dynamics along with plasma flows. Rotational motions were also observed and simply interpreted as the m=0 torsional Alfvén wave by De . Cross-sectional and intensity oscillations in on-disk SCFs were recently reported by Gafeira et al. (2017a) and interpreted as sausage mode waves.
The presence of concurrent wave-modes in localized MFTs was also claimed in a few studies. Transverse oscillations along with cross-sectional width variations interpreted as kink and sausage MHD wave-modes in spicule structures were reported by Jess et al. (2012). Morton et al. (2012), also found ubiquitous signatures of transverse, cross-sectional width, and intensity oscillations in on-disk fibril structures. They postulated that these independently excited wave-modes at MFT foot-point can undergo mode-coupling and can also exchange energies through a mode-conversion mechanism at the b = 1 layer in the solar transition region. Morton et al. (2012) claimed that simultaneous nonlinear kink and sausage waves were present in chromospheric MFT structures. Sharma et al. (2017) showed that the rotational motions in spicules appeared either axisymmetric or non-axisymmetric depending on the line-ofsight with respect to the bulk transverse motion.
In the linear regime, the cross-sectional width and photometric variations are signatures for both m=0 sausage and  m 2 fluting modes in MFT structures. Beyond the small amplitude regime, the nonlinear kink mode can also distort the cross-sectional width/shape and induce stronger density/ intensity variations, resulting in accelerated wave damping due to coupling to higher order fluting modes (Ziegler & Ulmschneider 1997a, 1997bRuderman et al. 2010). Recent numerical simulations (Magyar 2016) have shown for coronal loop structures that nonlinear kink oscillations are coupled to the m=2 fluting mode, causing cross-sectional width and intensity fluctuations with half the period of the kink mode. All of these factors can result in ambiguity in the identification of MHD wave-modes in the Sun's atmosphere.
This paper explores the complex and coupled evolution of plasma three-dimensional (3D) motions in spicules by examining the phase relationships between the resultant tube displacement (x r ), cross-sectional width ( ¢ W ), intensity (I), and azimuthal shear/torsion (x t ) components in spatial and temporal domains.

Observations
The Hα data and the observed spicules used in this study are the same as those reported by Sharma et al. (2017). Data were taken by using the CRisp Imaging SpectroPolarimeter (CRISP) at the Swedish 1m Solar Telescope (Scharmer et al. 2003(Scharmer et al. , 2008 on La Palma. The Hα imaging-spectroscopic data of 07:15-07:48 UT, 2012 June 21 is of Active Region NOAA AR11504, which consisted of two sunspots at the limb position (heliocentric coordinates with respect to the disk center, hereby denoted by: Θ=893, Φ=−250). The AR was scanned using 31 equally spaced line positions with 86 mÅ, steps from −1.376 to +1.29 Å, relative to the line center, along with the additional 4 positions in the far blue wing from −1.376 to −2.064 Å. This data was then further processed using the Multi-Object Multi Frame Blind Deconvolution (MOMFBD; van Noort et al. 2005) image restoration algorithm. Also, standard procedures available in the image pipeline for CRISP data (de la Cruz Rodríguez et al. 2015) including differential stretching and removal of dark-and flat-fielding were implemented. Any line-scan positions were, however, not normalized to the local continuum. The final science-grade data of ∼30 minutes duration had a pixel size of 0 059 (∼43 km), angular-resolution of  0. 13 (∼95 km), and cadence of 7.7 s. According to the Nyquist criterion, this pixel (cadence) resolution allowed us to detect MHD waves in the spicular waveguide with wavelength longer than 86 km (15.4 s), though, the angular-resolution restricts any observed wave behavior below ∼100 km to be unrealistic. The leastsuperimposed spicule structures outlined as high-intensity features were studied with a case study (SP1), presented here for detailed analysis and shown in Figure 1. The respective lifetimes, height, length, and other dynamical parameters are all listed in Sharma et al. (2017).

Longitudinal (Field-aligned) Motions
Longitudinal motions associated with mass upflows and downflows were estimated by tracking the visible apex of the spicule structure in intensity images. It is hereby assumed that the apparent changes in apex location at each time frame is solely due to the plasma motion and not due to any thermal or seeing effects. The apex of the inclined spicule structure is visually enhanced with respect to the background emission by the use of an unsharp-mask (USM) procedure for correct estimation of the height. The temporal evolution of the spicule apex (marked as yellow dot) is shown in Figure 1.

Transverse Displacement
Transverse motions, as observed in the plane-of-sky (POS) and the line-of-sight (LOS) domains by imaging-spectroscopy, were used to estimate the resultant motion of the spicular structure. These perpendicular displacement components were measured using the time-distance (TD) diagrams, generated by applying cross-cuts on intensity (for POS) and Doppler (for LOS) images. The intensity images at Hα line-profile wing positions were used to estimate the POS displacement. The spicular feature was further enhanced with respect to the background emission by the use of a USM procedure. Perpendicular to the spicule axis, cross-cuts were placed at every fourth pixel (∼172 km apart) in height, well above the angular-resolution limit. Similar cross-cuts were also made for the original (non-unsharp-masked) images for error estimation. The oscillatory axis of the spicule structure in the resulting TD diagrams was located as the nearest pixel with the maximum intensity magnitude at each time step.
Spectroscopic measurements in the Hα spectral line at 35 line-scan positions were used to estimate the LOS velocity components. Each image frame with the region-of-interest had pixels with both absorption and emission profiles. This time/ space-dependent behavior of Hα line profiles is best fitted with a single-Gaussian function for pixels near the limb but required a different fitting function for the pixels with asymmetric profiles at heights. This asymmetric behavior in spicular line profiles was explained as a consequence of LOS superposition effects from multiple thread-like structures in spicular waveguides by Skogsrud et al. (2014), generated by K-H vortices (Antolin et al. 2014) formed due to large-amplitude transverse MHD wave propagation (Scullion et al. 2011). All of these factors restrict the use of simple fitting procedures. To address this, a double-Gaussian function was employed, given by is the discrete intensity line profile, "λ" is the wavelength, "b" is the background signal level, "I" is the intensity, "μ" is the mean of the distribution, and "σ" is the standard deviation. The index in subscript refers to the two Gaussians. The estimated LOS Doppler velocity ( ¢ V y ) is then the normalized shift in wavelength, multiplied by the velocity of light.
As the observed structure is inclined to an angle (q =  23 . 6) with the normal, there will be certain contribution from the vertical flow component to the estimated Doppler velocity. This is based on the assumption that the observed structure has similar inclination in the LOS, as measured from the POS intensity images. It is, however, not possible to identify the inclination of the observed structure in the LOS from the current data set. The relationship between vertical flow velocity (V f ), estimated from the longitudinal displacement of the apex position between the two consecutive images (Figure 1), and the apparent LOS velocity was given by Athay & Bessey (1964). According to them, the apparent LOS velocity (V y ) component from the flow velocities can be given as sin .

y f
This contribution from vertical flows to the Doppler velocity measurements is then removed for accurate LOS transverse estimates. The corrected LOS velocity component (V y ) for each time step is given as,

Cross-sectional Width and Intensity Fluctuations
The width of the spicule structure was estimated by fitting a Gaussian function with a linear background to the crosssectional intensity profile (Figures 2(a), (b)). The linear background term in the function incorporates the contribution from the ambient plasma and the geometrical effects of the structure on the radiated optically thick emission. These estimates over the height excluded the regions near the footpoints and at the diffuse apex (Figure 2(d)) of the spicular feature because of flattened/non-Gaussian intensity profiles. Contribution from the background intensity flux near the spicule foot-point and/or from weak contrast in spicule and background intensities at apex can result in a nearly flat intensity profile, which can led to the erroneous estimation of the Gaussian width and hence was not taken into account. The Gaussian function fitted to the cross-sectional intensity profile I(x) at each height of the spicular structure is given by Here, I(x) is the intensity profile across the MFT structure, "a" is the peak flux, "μ" is the central position of the Gaussianfit, "σ" is the Gaussian width while "b" is the linear background contribution. The cross-sectional width ( ¢ W ) of the spicular structure is approximated as the Full-width Halfmaximum (FWHM) of the Gaussian-fit and is given as 2σ ( ) 2 ln 2 1 2 . The mean magnitude of the estimated FWHM for all heights is considered as the unperturbed width (W) of the spicular waveguide. Averaged intensity magnitudes at ±1 pixel to the central position (μ) of the Gaussian-fit are taken for the photometric analysis.

Azimuthal Shear/Torsion
Torsional motion in spicular structures were reported before by De , where they used the relative tilt in the horizontal wavelength directions in the observed spectrograms to measure the magnitude of azimuthal shear/torsion (x t ). Here, the gradient in the Doppler velocities (V y ) at the unperturbed width (W) of the spicular waveguide is used to estimate the angular displacement or azimuthal shear/torsion. In the presence of any shear forces, the magnitude of the Doppler velocities at the two ends would be unequal, with one   edge leading the other. Figure 2(b) shows that the net velocity ( -V V 2 1 ) between the two edges of W provides the displacement (AB) at the leading edge of the spicular cross section. For the given radius (OA) of the unperturbed cross section, the shear-angle (θ) at the tube center is given by The angular displacement ( ¢ AB ) is then taken as q ¢ = · AB O A. The temporal evolution of the estimated azimuthal shear/torsion is given in Figure 3(e).

Multi-taper Spectral Analysis (MSA)
and Phase Relationships To analyze the phase relationships between the oscillatory behavior of the spicular structure with the cross-sectional and photometric variations, we used the MSA method (Thomson 1982;Percival & Walden 1993). This method overcomes the problems (like inefficiency and spectral leakage) posed by discrete Fourier transform methods for small data samples. Other robust techniques available, like, Empirical Mode Decomposition (EMD) and Recurrence quantification analysis (RQA), for nonlinear analysis, were not employed due to the small sample (temporal) size. For data sets with relatively small sample sizes, estimating the amplitude and phase between corresponding component frequencies (spatial or temporal) could suffer from the inefficient statistical approximation of the parameters. Although, this problem could be tackled by the use of a sliding window (single-taper) with the assumption that the input signal is stationary, but that it can further lead to spectral leakage and contribute to the erroneous estimation of the amplitude and phase angle. A possible solution to the above is the use of multi-taper functions. A number of tapers are applied to the data sample, and the power spectrum is estimated for each of those tapers. The resulting averaged spectrum provides the best estimate of any low-amplitude harmonic oscillations in a relatively short sample size with a high degree of statistical significance. This technique has previously been tested (Komm et al. 1998) and used for the solar data analysis in a variety of spatial and temporal scales (Prestes et al. 2006;Kilcik et al. 2010;Mufti & Shah 2011).
The method used in this analysis is adapted from van Hoek et al. (2016). Here, in our case, the spectrum (ĝ) is estimated in the wavenumber (k z ) and frequency ( f ) domains by taking Q tapers to obtain a set of Q eigenspectra, including "q" eigencomponents using the observed parameters (i and j) with N data points and a constant sampling interval, so that = u N 1, 2, 3 ,... . For a given domain (ϝ) in the spatial (k z ) or temporal ( f ) regimes, the set of eigenspectra of i is defined as The cross-spectral density (CSD) provides an estimate of the strength of the coupling between the observed parameters and further the evolution in wavenumber/frequency domain. This is computed using the real valued power-spectral density (PSD) estimate for "i" defined asˆ(ϝ) g ii MT and the complex conjugate of the PSD estimate of j defined as * (ϝ) g jj MT and given as

MT MT
The phase-angle estimate (j (ϝ)) gives an approximation of the relative perturbations between the observed parameters and with that from ambient atmosphere. The estimate is bounded between p  and is the phase difference at each wavenumber/ frequency bin between i and j. It is calculated from the imaginary and the real part of the CSD is the real part of the estimated CSD.

Results and Discussion
Detailed analysis of the different dynamical observables (longitudinal, transverse, cross-sectional width, photometric, and azimuthal shear/torsion variations) are presented here for spicule (SP1), which, along with other similar structures, was studied for wave-mode identification by Sharma et al. (2017). From analysis of the 3D velocity components, it is shown that the concurrent rotational and transverse motions in spicules were most consistent with the kink mode. The observed feature (SP1) had a length of 4.1 Mm with apex-height reaching up to 4.9 Mm from the visible limb at an inclination of  23 . 6 from the normal. The temporal evolution in imaging showed nearly constant inclination of the structure, with apex (marked in Figure 1) undergoing periodic variations along the inclined length. The average longitudinal flow velocity (V f ) is about 41±28 km s −1 , which is consistent with the previous reports for MFT waveguides (De Pontieu et al. 2007a;Pereira et al. 2012;Sekse et al. 2013).
The temporal evolution of the transverse motion of the spicule structure at different heights was studied by TD diagrams, generated by the procedures described in Section 3.2. A visual inspection of a sample TD plot in the imagingspectroscopy (Figure 2(c)) at 3.2 Mm shows the bulk transverse displacement of the spicule with the observed wave-period comparable to the feature's visible lifetime at the line-scan position (±1.204 Å, from the Hα line core). This gives an impression of the presence of a single pulse being responsible for the translation of the structure, with maximum transverse displacement of 180±70 km from its mean axis. This pulselike behavior in other similar features was also evident from the past statistical studies (De Pontieu et al. 2007a;Okamoto & De Pontieu 2011), where the observed lifetimes coincided with the wave-periods estimated from TD cross-cuts. Recently, Martinez-Sykora et al. (2017) modeled the spicule-like events and attributed the role of an impulsive release of the magnetic tension (whiplash-effect), responsible for the observed behavior of spicular structures. TD plots at every four pixels over height in imaging and spectroscopic observations (Figures 3(a), (b)), respectively, were used to study the evolution of MFT dynamics in time and to generate the resultant displacement profiles (Figure 3(c)), which reflects the true motion of the waveguide. Here, it should be noted that previous studies Morton et al. 2012) only considered the MFT motion confined in a particular plane, whereas transverse spicular dynamics is the resultant of the observed motion in the perpendicular observation planes (POS and LOS).
The perturbed resultant displacement (dx r ) along with the cross-sectional width (dx = ¢ -W W W , Figure 3(d)), intensity (dI ), and azimuthal shear/torsion (dx t , Figure 3(e)), with average magnitude of 44±39 km, variations were examined in both spatial and temporal domains in order to understand any possible relationships between these observed parameters. The unperturbed cross-sectional width (W) is found to be 156± 24 km, as compared to 360±120 km for on-disk fibril structures , 260 km for RBE/RREs (Kuridze et al. 2015) and on-disk SCF structures by Gafeira et al. (2017a). The PSD estimate is taken at each time step (t) for = t T 1, 2, 3 ,... , over height to understand the wavenumber (k z ) dependence of these parameters over time. Similar estimates were taken at each height (H) over time, to analyze the evolution of observed parameters in the frequency ( f ) domain. The averaged-PSD is then defined as á ñ á ñ ( )ˆ( ) g k g f , ii z jj and given as  y -axis (on right) are marked here as star. The phase differences are mostly concentrated around 0°and p   with a deviation of   30 . Any further deviation of the mutual phase relations could be indicative of the nonlinear evolution and/or background noise contribution from the small-scale dynamics at the spatial/temporal scales below current observational resolution limit. the peak PSD (l = 2800 z km) of the resultant displacement and the subsequent parameters is highlighted by the shaded region (Figure 4(a)). The coincided peak PSD magnitudes for all four of the observed parameters indicate a strong coherence between these observables at smaller (larger) wavenumbers (wavelengths). The remaining region, with increasing wavenumber (k z ), shows an exponential decay in averaged-PSD estimates for all of the observed quantities. The strong coherent behavior between these observables is also evident in the frequency ( f ) domain. In here, the results are two-fold.
First, the observed parameters (dx dx dI , , r W ) show peaks at 0.013 Hz (77 s) and 0.039 Hz (25.6 s), i.e., a secondary peak at around thrice the frequency of the primary peak, marked by the vertical solid lines in Figure 4 This subsequent period-doubling and -tripling pattern in frequency domain (Linsay 1981;Jiang et al. 1998) for observables could indicate the nonlinear nature of the spicular dynamics. However, the existing nonlinear kink (analytic/ numeric) models have only produced period-doubling (Ziegler & Ulmschneider 1997a, 1997bRuderman et al. 2010;Magyar 2016) and have not yet produced period-tripling. The 3D visualizations of these coupled dynamical components are shown in Figure 5.
Coupled transverse and non-axisymmetric deformation of the MFT cross section could partly be due to the injected photospheric perturbations (as assumed in previous studies) and/or as a consequence of the flux tube to maintain the local pressure equilibrium at the tube surface. The transverse dynamics of the waveguide will be resisted by the inertia of the spicular plasma (ionized, partially ionized and neutrals). This can result in a non-axisymmetric deformation of the MFT cross section with an anti-correlation between plasma and magnetic pressures, in an attempt to maintain the overall pressure equilibrium. While the gradients in plasma pressure due to kink waves were already reported as acceleration at the transverse nodes by Sharma et al. (2017), it is highly reasonable to assume the deviation of spicular cross-sectional geometry to a non-circular case, to compensate for the asymmetry in pressure distribution between forward and backward sides of the moving waveguide (see upper panel, Figure 5).
The physical explanation is well supported by the observed amplitudes for perturbed transverse displacement (dx r ) and the cross-sectional width (dx W ) oscillations, which are comparable (Figure 4(c)). Similar observations were reported by Morton et al. (2012) for on-disk fibrils. As the overall cross-sectional width and the associated photometric variations are directly related to the pressure perturbations, linked with the tube dynamics, it is expected for these to be coupled at similar wavenumbers and frequencies. This aspect is statistically supported by the observations of comparable periods for cross-sectional width and intensity oscillations for on-disk SCF structures, reported recently by Gafeira et al. (2017a).
The coupling in between the transverse (dx r ), cross-sectional width (dx W ), and the azimuthal shear/torsion (dx t ) parameters could be attributed to the magnetic tension forces and/or longitudinal plasma pressure perturbations. When local curvature is formed between the mean axis and the maximum displacement, there will be components, both perpendicular (vertical) and along (horizontal) to the direction of transverse displacement. Over a single wave-period, these vertical components will reverse their sign twice, as compared to the horizontal component. It must, however, be noted that similar vertical components associated with plasma pressure changes would also be present in the spicule due to cross-sectional width variations. The frequency analysis, as shown in Figure 4(b), hints toward the proposed mechanism, with azimuthal shear/torsion component (dx t ) having twice the frequency ( of the primary frequency of the transverse displacement (dx r ) and cross-sectional width (dx W ) parameters. Similar interpretations were proposed by (Ulmschneider et al. 1991;Ziegler & Ulmschneider 1997a) for their numerical analysis of the dynamics of the thin MFT structures. These changes in the vertical component will result in longitudinal compressions and rarefactions, which will further perturb the associated Lorentz forces, visible as the azimuthal shear/torsion (see lower panel, Figure 5). Under the influence of longitudinal compressions and rarefractions, the perturbed azimuthal shear/torsion (dx t ) will further showcase upward and downward motions along the spicular structure (see, supplementary material for Figure 5). This longitudinal motion of the perturbed Lorentz forces was highlighted by Kitiashvili et al. (2013) in their numerical simulations of vortex tube dynamics.
The strong coupling between the observables is also evident from the evolutionary trends estimated at different heights. Figure 4(c) shows the observed parameters at four consecutive heights (1.92, 2.56, 3.20, 3.85 Mm) above the visible limb, with σ error marked as shaded regions. A visible overview of the cross-sectional width (dx W ), photometric (dI ) and azimuthal shear (dx t ) components shows clear signatures of both in-phase and out-of-phase behavior with the resultant displacement (dx r ) of the spicular waveguide. For a better understanding of the coupled behavior, the phase relationships between the paired parameters were estimated using the CSDs in the spatial and temporal domains by the procedure described in Section 3.5. The mutual variations in between the estimated phase angles show the phase difference congregated around 0°and   180 with a deviation of   30 for the spectral bins located in wavenumber (l z in highlighted region, Figures 6(e)-(h)) and frequency ( f=0.013 Hz, Figures 6(a)-(d)). However, the phase angles for cross-sectional width, intensity, and azimuthal shear/torsion pairs in wavenumber (Figure 6(h)) and frequency (Figure 6(d)) domains appears to be strongly concentrated around 0°, indicating collective behavior rather than uncoupled azimuthal shear/torsion and sausage-type motions.

Conclusions
For the first time, the resultant displacement (dx r ) estimated from both the POS and LOS velocity components was analyzed with the cross-sectional width (dx W ), photometric (dI ), and azimuthal shear/torsion (dx t ) variations for off-limb spicules. TD plots show the displacement of the oscillatory axis of the spicular structure in both perpendicular planes of motion, indicating the presence of the kink wave-mode, with periodic longitudinal (field-aligned) flows (V f ). Observations of spicular waveguide and TD analysis further reflected the "pulse-like" behavior of the wave-mode with the dominant period of the order of the flux tube's lifetime.