Observations and data reduction
In this work we present observations acquired on 14 April 2025 between 21:38 ut and 21:41 ut using the world’s largest telescope for optical and infrared observations of the Sun, DKIST21. During the data acquisition, the atmospheric conditions were variable, with Fried parameter values during the best periods measured around \({r}_{0}\approx 12\,\mathrm{cm}\). The DKIST wavefront correction system was operating in diffraction-limited mode at the time of the observations, as it was measuring and correcting all optical modes that it is capable of detecting. In this work, we focus on the first target of the day, a region containing pores near an active region (NOAA 14060) close to the disk centre with helioprojective Cartesian coordinates (X = −162 arcsec, Y = 168 arcsec), corresponding to \(\mu \approx 0.97\) (Fig. 1).
A dataset was acquired with a diagnostic (FastCam) set-up at DKIST, which was built as a collaborative effort between the National Solar Observatory and the Max Planck Institute for Solar System Research. The set-up was installed in front of the spectrograph entry slit of the visible spectro-polarimeter and removed at the end of the data acquisition for the diagnostic test. The optical set-up enhanced the operation of the feed telescope of the visible spectro-polarimeter, which did not have any other powered optics, resulting—with the pixel cell size of 5.5 μm for the Max Planck Institute camera—in a pixel size of ∆s ≈ 0.00825 arcsec px−1, corresponding to approximately 6 km px−1 on the solar surface. The wavelength chosen for the observation was 416 nm, and we used a narrow band-pass with a full-width at half-maximum of 0.5 nm. The region of interest read out by the camera was set to 2K × 1K px2, which meant that the sensor could be read out at a rate of 740 frames per second while maintaining an exposure time of 100 μs. The camera was equipped with a phase-diversity beam-splitter assembly, which imaged two fields, one in-focus and the other with a de-focus of approximately 0.4 waves root mean square. The two fields were imaged side-by-side onto the sensor in a phase-diversity set-up38,39, which allowed us to determine the residual wavefront aberrations that remained after the correction by the adaptive optics system40. The solar FoV visible in both the focused and the de-focused images corresponds to approximately 5,800 × 4,350 km2 on the Sun (8 × 6 arcsec2). The time span covered by the dataset analysed in this work is about 3 min. The frames observed with this camera were image-reconstructed using multi-frame blind deconvolution41, a technique that can be used to numerically model and remove residual atmospheric aberrations from the data. The algorithm used 2,000 dark- and gain-calibrated camera frames to compute one reconstructed science frame. After image reconstruction, the effective cadence of the imaging data was 2.7 s. We estimated that most frames achieved a spatial resolution of 19 km, the Rayleigh diffraction limit at 416 nm, as inferred from the full-width at half-maximum of the finest structures in the data.
Context imaging was provided by the visible broadband imager (VBI)42, an imaging instrument capable of acquiring large-FoV short-exposure data sequences at a single wavelength within 3 s. In this experiment, the VBI was configured to acquire sequences of 80 frames alternating in wavelength between the red continuum (668.4 nm) and the Hα line (656.3 nm; not shown), with sampling in the images corresponding to 12.3 km px−1 (0.017 arcsec px−1). The frames of each consecutive sequence were image-reconstructed using speckle reconstruction algorithms43 to remove residual seeing effects over the large FoV of approximately 50 × 50 Mm2 (69 × 69 arcsec2). Even though the VBI red channel was started with a 3-min delay with respect to the FastCam data that are presented here, the FoV clearly still contains the (evolved) structure of the higher spatial resolution dataset (Fig. 1).
Full disk observations were provided by the helioseismic and magnetic imager onboard the Solar Dynamics Observatory (SDO/HMI)44. Figure 1a shows a photospheric continuum image from the HMI at 21:40 ut. The HMI image was coaligned with DKIST/VBI/FastCam data using the SolarSoft auto_align_images function through a cross-correlation (Fig. 1).
Numerical simulations and spectral synthesis
Here we compare the observations with MHD simulations performed with the MURaM code22. So that we could perform the simulations at a very high spatial resolution, we focused on a small domain of extent 6.144 × 6.144 × 2.048 Mm3 (with the third dimension, z, being vertical), which is comparable with the high-resolution FastCam FoV. In the vertical direction, the domain extended approximately 1.367 Mm below z = 0 Mm and 0.681 Mm into the overlying atmosphere. A solar-like configuration with granulation driven by radiative losses could not be achieved in the numerical simulation without the overlying stable atmosphere. A realistic photosphere was essential for making the connection between numerical models and observations. We started from a relaxed hydrodynamic simulation to which we added the vertical magnetic field component extracted from the corresponding HMI magnetogram for the observed region (Supplementary Fig. 1). We initialized the domain with a strictly vertical magnetic field of the form \(n{B}_{z}-(n-1)\langle {B}_{z}\rangle \), where \({B}_{z}\) was taken from an HMI magnetogram and was enhanced to compensate for the field dispersal that occurred due to granular motions when the simulation was further relaxed for about 1.8 h (n is the enhancement factor). After the relaxation, the individual flux concentrations were no longer impacted by the initial state. To achieve a realistic simulation of magneto-convection, the critical parameter is the net flux imbalance, which was unchanged by the enhancement. Therefore, the enhanced HMI magnetogram was used solely to create a set-up that looks more closely like the observed distribution, and it was not needed for the development and study of KHIs. We found that a choice of \(n=2.5\) led to a large-scale magnetic field distribution comparable with the HMI magnetogram and the formation of pores that resemble those observed (Supplementary Fig. 1). After initialization with the enhanced magnetogram, the simulation ran for 4,980 s with 12.8-km grid spacing, followed by 1,200 s with 6.4 km and then another 1,200 s with 3.2-km grid spacing. The sequence we analysed also covers 480–720 s in the 3.2-km sequence or 5,460–5,700 s since the initialization with the HMI magnetogram. The simulation was computed with 12 opacity bins and used the Asplund 2009 opacities45.
The grid scale of 3.2 km was greater than the diffusive length scale expected in the photosphere of approximately 1.4 km based on a Spitzer diffusivity \(\eta \) = 2 × 108 cm2 s−1 (ref. 46) and a timescale \(\tau \) = 100 s using \(l\approx \sqrt{\eta \tau }\). The use of the numerical diffusivity22 is, thus, justified for the simulations presented here. However, a spatially dependent Spitzer diffusivity would need to be used for higher resolution simulations.
Developed snapshots of the simulation were used to generate the synthetic spectrum for a wavelength interval between 415.32 nm and 416.7 nm, using 500 spectral points with the one-dimensional version of the Rybicki–Hummer radiative transfer code47 in local thermodynamic equilibrium (Extended Data Fig. 4c). The syntheses included molecular (CH and CN) and atomic lines and the continuum within the band-pass. Intensities were also synthesized at a single continuum wavelength of 500 nm. To save computational time, the vertical grid space was reduced to 9.6 km by using only every third grid point along the vertical axes. We synthesized spectra at disk centre and towards the solar limb at \(\mu =1\) and \(\mu =0.97\), respectively. To synthesize 416 nm intensities at \(\mu =0.97\), the MURaM cube was modified by shifting each horizontal layer relative to the layer below by \(\Delta z\,\tan (\theta )\), where \(\Delta z\) represents the vertical grid spacing and θ is the heliocentric angle, defined as the angle between the line of sight and the solar surface normal. This adjustment aligned the slanted line-of-sight direction vertically. Pixel sampling was foreshortened along the y direction by a factor of \(\mu \) and increased along the z direction by a factor of \(1/\mu \) for the \(\mu =0.97\) data. To closely replicate the DKIST 416-nm filtergrams, the synthetic data were multiplied by the transmission profile of the interference filter employed in the diagnostic set-up (Extended Data Fig. 4c) and integrated over the wavelength range covered by the filter.
Extended Data Fig. 4a,b demonstrates that all common features observed in the 416-nm filtergrams were reproduced in the synthetic images, including KHI vortices and striations. Extended Data Fig. 4c shows the synthetic spectrum averaged over the full FoV of the simulation, overplotted with a solar atlas48 with the same wavelength range as observed by a Fourier transform spectrometer. The two spectra agree to a high level of detail, demonstrating that the simulation accurately reproduced the observed spectral features. The close agreement between the synthetic and observed spectra and images confirms that the simulated data serve as a powerful diagnostic tool for interpreting 416-nm intensity observations.
From the MURaM snapshot, we calculated the height of the surface where the optical depth reached unity at different wavelength positions along the synthesized spectral interval (Extended Data Fig. 9k,l). Extended Data Fig. 9l shows that the heights of \(\tau =1\) and \({\tau }_{415.91}=1\) (as well as other 416-nm continuum wavelengths) are very similar with only marginal differences, indicating that these two continuum wavebands formed at nearly the same geometrical heights. We used the optical depth at 500 nm to select the heights at which KHI parameters were directly studied in the MURaM cube and compared with the synthetic 416-nm filtergrams.
Historically, striations have been observed and studied primarily at inclined viewing angles27,28,29,30. Our simulations show that they are also present in the disk centre. Extended Data Fig. 9k confirms that the structures seen in the disk-centre filtergrams arose from opacity variations. Their appearance in the emergent intensity is closely correlated with spatial variations in \({B}_{z}\) and anticorrelated with density in the layers above the KHI formation heights (Supplementary Fig. 9). These variations shift the geometrical height at which the emergent intensity forms, giving rise to the appearance of striations in the disk-centre synthetic filtergrams. This is the same physical mechanism that produces striations at inclined viewing angles28,29,30,31. An inclined viewing geometry makes these features appear longer and more extended towards the limb and, therefore, easier to detect. Our results, however, indicate that they should also be observable at disk centre, given a sufficiently high spatial resolution.
KHI growth rate in linear MHD theory
The linear theory of the KHI in MHD systems shows that the instability grows fastest when the shear flow is perpendicular to the magnetic field and the wavevector is aligned with the shear velocity23. For finite and continuous velocity profiles, the fastest-growing wavelengths are determined by the thickness of the velocity shear layer.
Chandrasekhar23 investigated the KHI in a slab geometry consisting of two adjacent incompressible fluids with different densities separated by a finite-width shear layer of intermediate density in which the velocity varies linearly from −\({U}_{0}\) to +\({U}_{0}\) across the layer. Specifically, the velocity and density profiles are given by equation (1) (Supplementary Fig. 10):
$$\begin{array}{lll}x < -d/2, & {\rho }_{1}={\rho }_{0}(1+{\epsilon }), & U=-{U}_{0},\\ -d/2 < x < d/2, & \rho ={\rho }_{0}, & U=2{U}_{0}x/d,\\ x > d/2, & {\rho }_{2}={\rho }_{0}(1-{\epsilon }), & U={U}_{0}.\end{array}$$
(1)
Here \({\rho }_{0}\) is the density in the shear layer, \({\epsilon }\) is a factor that defines the density contrast between two layers, \(x\) is the distance across the KHI surface and \(d\) is the thickness of the velocity shear layer. The continuity of the perturbation equations results in the dispersion relation23:
$${{\rm{e}}}^{-2\kappa }=\left[1-\frac{\kappa {(\nu +1)}^{2}}{(\nu +1)+1/2{\epsilon }\kappa {(\nu +1)}^{2}}\right]\cdot \left[1+\frac{\kappa {(\nu -1)}^{2}}{(\nu -1)+1/2{\epsilon }\kappa {(\nu -1)}^{2}}\right],$$
(2)
where \(\nu =\gamma /k{U}_{0}\), \(\kappa ={kd}\), and \(\gamma \) is the frequency and \(k\) is the wavenumber of the KHI. In this equation, the magnetic field contribution and Richardson number are neglected because the shear flow is perpendicular to the magnetic field and the wavevector is perpendicular to the gravity vector.
The dispersion relation (equation (2)) is a transcendental equation for the complex frequency \(\gamma ={\gamma }_{{\rm{r}}}-{\rm{i}}{\gamma }_{{\rm{i}}{\rm{m}}}\), whose imaginary part \({\gamma }_{\mathrm{im}}\) indicates an instability process in the system, in particular, the exponential growth rate of unstable modes. The real part \({\gamma }_{{\rm{r}}}\) is the frequency related to the propagation speed νph (phase speed) of KHI modes through \({\gamma }_{{\rm{r}}}=k{v}_{\mathrm{ph}}\). Extended Data Fig. 6 shows the sensitivity of KHI wavelength, growth rate and phase speed on the parameters \({U}_{0}\), \(d\) and \({\epsilon }\) for values within the ranges found in our data (Extended Data Table 1).
Estimating the KHI parameters
To estimate the KHI growth rate and phase speed, we measured the relevant parameters in the simulation and observation for several prominent cases, which are summarized in Extended Data Table 1.
The simulation shows that the KHI occurred coherently across several horizontal layers. Figure 4, Extended Data Fig. 9 and Supplementary Video 5 demonstrate that KHI vortices are three-dimensional structures with different amplitudes at different heights. The KHI patterns observed in the 416-nm filtergrams are primarily defined by the layer where the continuum optical depth at 416 nm (or 500 nm) approached unity. Owing to the Wilson depression, the \(\tau =1\) layer formed at different geometrical heights inside and outside MFCs (Extended Data Fig. 9). We measured the KHI parameters at the height where the average continuum optical depth at 500 nm reached unity.
First, we identified the KHI interface at this height, calculated the component of the two-dimensional velocity vector parallel to it \({v}_{\parallel }\), and estimated its average one-dimensional profile perpendicular to the interface (Extended Data Fig. 7f and Supplementary Figs. 3f–6f). The average thickness of the velocity shear layer was defined as the spatial extent over which the parallel velocity component \({v}_{\parallel }\) transitions between the two layers. The shear velocity \(\Delta U\) was estimated as the difference in \({v}_{\parallel }\) across the shear layer, determined from a linear fit to the velocity \({v}_{\parallel }\) profile (Extended Data Fig. 7f and Supplementary Figs. 3f–6f). Average densities on either side of the shear layer, along cuts perpendicular to the KHI surface, were also determined. The wavelength \(\lambda \) was measured as the spatial separation between adjacent developing vortices during the initial (linear) growth phase of the instability (Extended Data Figs. 2c and 5c).
The growth rate was measured from the temporal change in KHI amplitude. In the linear phase, the perturbation amplitude of KHI increased exponentially, \(A(t)=A({t}_{0})\,\exp ({\gamma }_{{\rm{i}}{\rm{m}}}t)\), where \(A({t}_{0})\) and \(A(t)\) are the amplitudes at \({t}_{0}\) and \(t\), respectively, and \({\gamma }_{\mathrm{im}}\) is the growth rate of the instability, which corresponds to the imaginary part of the frequency in the theoretical dispersion relation presented and analysed above (equation (2)). We measured the instability amplitude directly from the corrugated boundary of the unstable interface. To do this, we drew at different time steps two lines parallel to the mean orientation of the shear layer to mark the full transverse extent of the deformed boundary (Extended Data Figs. 3b and 7d). The curve that follows the separation (red curves in Extended Data Figs. 3b and 7d) between these two lines was taken as a proxy for the KHI amplitude as it captures the degree of boundary corrugation and the transverse displacement of the boundary. The temporal evolution of \(A(t)\) was then fitted with an exponential function (Supplementary Fig. 11). The amplitudes and growth rates were derived from the early linear phase of the instability, before the onset of nonlinear saturation.
To determine the phase speed, we used time–distance diagrams and calculated the average speed at which different vortices were displaced along the KHI interface in the image plane. In equation (1), the reference frame is defined by the velocity at the centre of the shear layer (\(x=0\)). For symmetric velocity profiles, this corresponds to \(v=0\) km s−1. The time–distance diagrams measure the apparent phase speed in the image frame along the shear layer. To maintain consistency with equations (1) and (2), the absolute phase speed was computed as \({v}_{\mathrm{ph}}={v}_{\mathrm{app}}-{v}_{{\rm{c}}}\), where \({v}_{\mathrm{app}}\) and \({v}_{{\rm{c}}}\) denote the apparent velocity and the velocity at the centre of the shear layer, respectively. Extended Data Table 1 summarizes the derived parameters.
Analogous procedures were applied to the observational data to determine the growth rate, average apparent phase speed and wavelength of the KHI (Extended Data Figs. 1–3 and Supplementary Fig. 11).
The theoretical growth rate was calculated by solving equation (2) using the shear layer thickness, shear velocity and density contrast derived from the simulations for the five cases listed in Extended Data Table 1. The results, shown in Extended Data Fig. 6, demonstrate good agreement between Chandrasekhar’s theoretical predictions and our measurements.
Statistics of KHI wavelengths
We manually identified KHI interfaces in the observed and synthetic images (by visual inspection) and studied the distribution of Kelvin–Helmholtz wavelengths. The selection criterion was the presence of at least two KHI-like vortices at the interface between the MFC and the granules. To quantify the Kelvin–Helmholtz wavelength, we measured the median distance between the dips appearing in the intensity profiles computed from cuts through the vortices along the KHI interfaces. Extended Data Fig. 2 displays the locations of KHI interfaces and an example of the intensity profile along one selected case.
Similar analyses were carried out for synthetic DKIST filtergrams generated from MURaM simulations and a spectral synthesis. The results are presented in Extended Data Fig. 5. KHI wavelengths were also measured directly from the MURaM \({B}_{z}\) map at the height corresponding to the average \(\tau =1\) (Supplementary Fig. 2).
Although the KHI wavelength histograms peak at values above the spatial resolution limit of DKIST for both observations and numerical simulations, we refrain at this point from concluding that smaller characteristic scales do not exist on the Sun, given the remaining differences in the two distributions.
