## Abstract

Convection in the Sun’s outer envelope generates turbulence and drives differential rotation, meridional circulation, and the global magnetic cycle. We develop a greater understanding of these processes by contrasting observations with simulations of global convection. These comparisons also enhance our comprehension of the physics of distant Sun-like stars. Here, we infer toroidal flow power as a function of wave number, frequency, and depth in the solar interior through helioseismic analyses of space-based observations. The inferred flows grow with spatial wave number and temporal frequency and are confined to low latitudes, supporting the argument that rotation induces systematic differences between the poles and equator. In contrast, the simulations used here show the opposite trends—power diminishing with increasing wave number and frequency while flow amplitudes become weakest at low latitudes. These differences highlight gaps in our understanding of solar convection and point to challenges ahead.

## INTRODUCTION

Despite several decades of space- and ground-based high-resolution observations and substantial theoretical and computational effort, fundamental questions on solar dynamics, such as the emergence of differential rotation, meridional circulation, and solar magnetism, remain unanswered. Convection is thought to be centrally important to these phenomena but is among the least understood aspects of solar dynamics due to the extreme parameter ranges in which it manifests. Gaps in our theoretical appreciation of convection persist. The mixing-length theory (*1*), a widely adopted prescription for modeling stellar convection and evolution (*2*), is a coarse description and does not predict flow details. Numerical simulations and laboratory experiments are also unable to access the relevant parameter regime (*3*, *4*). Under the circumstances, a fruitful avenue to pursue consists of contrasting specific details between simulations and observations, highlighting disagreements and resolving them where possible. Although comparisons are stymied because it is often the case that simulations and observations do not pertain to same quantities, they represent important steps toward an improved understanding. The first such attempt using helioseismology concluded that convective velocity amplitudes were small in comparison to numerical simulations (*5*). This result was subsequently challenged (*6*) by an analysis using a different helioseismic technique, which inferred amplitudes to be consistent with numerical simulations. The disagreement remains unresolved, as a comparison between a multitude of constraints, helioseismic or otherwise, and different numerical simulations (*7*) reveals. A major source of this conflict is the lack of consistency in various parameters in the comparisons: differences in the extents of temporal and spatial averaging, the type of flow (poloidal or toroidal), etc. These efforts, however, were productive, focusing attention on this important question and have moved the field forward, leading to the development of concrete models and new results on the physics of convection (*8*).

Solar convection is a multiscale phenomenon, occurring over broad ranges in spatial and temporal frequencies. The large scales are very interesting since they likely traverse layers deep in the convection zone (*3*, *7*) and are therefore very challenging to model, given that they encounter vast gradients in density and adiabatic index. A major puzzle concerns flows on the largest scales, which are observed at the surface to contain the least power. In contrast, classical turbulence theory suggests that most of the turbulent power lies in spatial scales comparable to the system size. By this definition, the Sun is not a classical turbulent system, since velocity power smoothly decays with decreasing wave number. Supergranulation, another puzzling feature of convection, spans the frequency range from 3 to 11 μHz and spherical harmonic degrees *s* dominantly ranging from 80 to 140, with a power that goes down to very low wave numbers as well. Therefore, the very largest scales of convection are likely to evolve at even lower frequencies, i.e., 3 μHz, and wave numbers of 100. On the basis of these estimates, in the present analysis, we infer toroidal flows within the temporal frequency range of [0,3.05] μHz, for wave numbers of [1,50].

## MATERIALS AND METHODS

Helioseismology provides a window into the internal dynamics of the Sun. Vigorous convective turbulence at the surface generates small-amplitude acoustic waves, which propagate in the solar interior and re-emerge at the surface, forming resonances in cavities (*9*). The modulation of the solar surface by acoustic waves, measured as Doppler shifts in absorption lines, is subjected to seismic analysis. In contrast to other helioseismic techniques, measuring normal-mode coupling from observations is straightforward (*10*–*13*). The additional benefit is that interpretation is tightly connected to well-posed semianalytical models; the price is that these models are nontrivial to evaluate and rely on tedious algebra. Normal modes also provide information about the poloidal-toroidal structure of fluid motions.

The method involves first constructing eigenfunctions for a reference model of the Sun, e.g., model S (*9*), which, for the linear solar wave operator, form an orthonormal and complete basis (*14*). Since the Sun contains convective flows and the reference model does not, the “true” eigenfunctions depart from the reference set. Specifically, convection is treated perturbatively since velocity amplitudes on large scales are substantially smaller than the local sound speed (*3*). Completeness of the reference basis allows us to express the true eigenfunctions as a weighted linear superposition of reference eigenfunctions. Modes of the Sun are thus in a “coupled” state in relation to the reference. The weights, termed as coupling coefficients, may be measured from observations and related to model perturbations through a straightforward inverse problem. The improvement in inferential accuracy obtained by including differential rotation in the background model, which is known to cause changes to the reference basis (*10*), is at most on the order of a few percent for the modes studied here. However, since the computational and technical overhead required to expand around the rotationally distorted eigenfunctions is significantly greater, we expand around the degenerate eigenfunctions of model S.

A useful aspect of this methodology, compared with other techniques of helioseismology, is that it allows us access to time variations of nonaxisymmetric perturbations on a global scale, an appropriate formalism for the application to dynamically evolving phenomena such as convection and Rossby waves. In addition, it provides a rigorous formalism with which to analyze vector flows, enabling a proper comparison with simulations, a notable improvement over prior work. In particular, we apply the Chandrasekhar-Kendall relation, with which vector flows may be decomposed into poloidal (first two terms of Eq. 1 below) and toroidal quantities (last term) (*15*)*Y _{st}* is the spherical harmonic of degree

*s*and azimuthal order

*t*,

*e**is the radially directed unit vector, σ is the temporal frequency associated with the flow, and*

_{r}**∇**

*=*

_{h}**e**

_{θ}∂

_{θ}+

**e**

_{ϕ}(sin θ)

^{−1}∂

_{ϕ}, where θ and ϕ are colatitude and longitude, respectively. Temporal frequency dependence is denoted using a superscript, i.e.,

*f*

^{ω}=

*f*(ω), where

*f*is a function of frequency ω. By construction, toroidal flows are mass conserving, do not have a radial component, and are therefore not involved in radial thermal transport, arising rather as a by-product of convective turbulence. The anelastic approximation (

*16*) may be invoked to express

*13*) in Eq. 1 to

*4*) poloidal and toroidal flow power as a function of solar radius. At every radius, the total flow power is split nearly equally between these flow components. Simulations of the classical Rayleigh-Bénard convection (

*17*) reveal that the ratio of toroidal to poloidal flow energies may be used to classify the regime of turbulence—e.g., buoyancy or rotationally dominated, geostrophic, etc. The simulations suggest that toroidal flows are a critical component of convective turbulence, and imaging them may help shed important insights into solar convection. Regardless of whether these simulations are complete, the measurement of these flows and categorizing their ratios is an important means of testing theories of convection. The similarity in velocity scales suggests that constraints on one might be used to infer the other. This indicates that, toroidal flows, despite not carrying thermal flux, are critical to understanding convection. Although our analysis here is limited to toroidal flows, it is still a significant first step toward the larger goal of measuring the full range of material motions arising from convective transport.

### Data analysis and inversion

The Helioseismic and Magnetic Imager (HMI) (*18*), onboard the Solar Dynamics Observatory, has been observing the Sun since 2010. We analyze 8 years of HMI global-mode data in the ℓ range of 50 to 192 (ℓ is harmonic degree) and select modes with resonant frequencies in the range of 1800 μHz ≤ ω_{nℓ}/(2π) ≤ 3600 μHz, where *n* denotes the radial order (*19*). The upper frequency limit of 3600 μHz is placed to exclude modes with known systematics (*19*). From the lower end, we exclude modes with very small linewidths (that low-frequency modes have), since characterizing the frequency dependence of convection requires sampling modes off resonance. Small linewidths imply that measuring the mode even slightly off resonance can lead to background power seeping into inferences. The global time series data are further divided into eight segments of 360-day-long sets.

These time series are transformed to the temporal Fourier domain and denoted by *m* is azimuthal order. We model the wavefield as a sum over weighted mode eigenfunctions, with each mode excited stochastically (the amplitude and phase are random variables). We therefore calculate the ensemble-averaged correlation *13*). Correlations are measured from 360-day-long data, which are therefore averages over a large number of source ensembles (∼10^{5} wave periods). Because convective cells are advected by the rotation of the Sun (Ω_{⊙}), we track the data in a corotating frame at the rate Ω/(2π) = 453 nHz, associated with the equatorial band, i.e., at latitudes roughly less than 30^{∘}. Super-rotating (*t* > 0) and subrotating (*t* < 0) features are identified with respect to this rate. Note that the Sun exhibits differential rotation in both radius and latitude and, thus, a variety of different tracking rates may equally be applied (see fig. S3 of section S5).

Correlations thus computed contain several independent parameters ω, σ, *t*, *m*, and ℓ and are also noisy, rendering interpretation difficult. Consequently, we introduce *B* coefficients, which condense and average the information present in raw correlations. The coefficient *W* is a theoretically derived weight function [see section S2 and (*20*) for details], is the coupling coefficient between modes (ℓ, *m*, and *n*) and (ℓ, *m* + *t*, and *n*). It is sensitive to the properties of the medium (*13*) through

The function *B* coefficient is taken over the discrete temporal frequency grid (see section S7 and the Fourier transform convention). We may thus infer toroidal flow *21*) through the measurement of Rossby waves in the Sun (*22*). Directly interpreting *B* coefficients is challenging since it is a more abstracted quantity than conventional helioseismic measurements such as travel times or mode frequencies. Essentially, it captures the extent of eigenfunction mixing between (ℓ, *m*, and *n*) and (ℓ, *m* + *t*, and *n*): Translating the *B* coefficient, a complex number, into a flow requires tedious algebra, much of which is summarized in the Supplementary Materials [see section S2 and (*13*)].

Modes in the Sun are stochastically driven, and the excitation process is assumed to be statistically independent for each wave number and temporal frequency channel because, at any given time, a very large number of isotropically distributed (roughly) granules across the surface of the Sun are each exciting waves with random phases and amplitudes. This implies that, for reference modes, *m*′) to (ℓ, *m*) results in nonzero values of *B* coefficient arises from the combination of leakage and realization noise and must be subtracted to obtain estimates of the flow power (Fig. 2). Toroidal power is thus given by *17*). Figure 2 shows *B*-coefficient power and the noise background (computed using eqs. S15 and S16 in section S3) as a function of *s* − ∣*t*∣ for each of the 8 years of HMI data that we have analyzed; the fact that the theoretically derived noise background and *B*-coefficient power are very close to each other, especially at higher values of *s* − ∣*t*∣, suggests that our noise background model is likely reasonable. The coefficients α_{nℓ}(*r*_{0}) are weights determined using the technique of subtractive optimally localized averages (*20*, *23*). The spatial mask arising from our partial view of the Sun induces systematical bias, resulting in broadening in the spectral domain and making it impossible to isolate individual spherical harmonics [called leakage (*24*)]. The consequent amplitude aberrations in inferred velocities are modeled using leakage matrices (*24*) and removed from α_{nℓ}(*r*_{0}) (*20*) (see eq. S22). These leakage matrices take into account spatial windowing and line-of-sight projection. Synthetic tests have also been performed to verify this method in the context of Rossby modes (*25*).

Last, we comment on the trustworthiness of the result. We have experimented with different sets of modes; specifically, we studied the results using the entire set of modes, i.e., those that have been successfully identified by the HMI pipeline from analyses of 360-day datasets (*26*). We found that the inferred amplitudes of prograde and retrograde flows are somewhat sensitive to the high-frequency modes (>3600 μHz), which are known to be problematic (*19*). However, one feature that is robust to all choices of mode sets is the variation of toroidal flow power with *s* − ∣*t*∣, i.e., sectoral modes are dominant and there is very little power in tesseral and zonal modes (for the spatiotemporal band under consideration here).

### Simulation setup and importance of toroidal flows

At the outset, it is important to note that this simulation is merely one of innumerable possibilities, with many possible choices for turbulence modeling, boundary conditions, resolutions, numerical schemes, etc. The comparisons with simulations here are therefore solely a commentary on the setups chosen for this work.

There is simply no means by which simulations will be able to reach solar convection parameter regimes in the foreseeable future. However, it is within the realm of possibility that there exists an effective theory of convection, i.e., where all the ingredients required to properly recover large-scale properties of these flows are captured in the calculations. An instance of this is the simulation of surface granulation, which is able to reproduce line formation, granule sizes, shapes, correlation functions, etc. [e.g., (*27*)] despite the very large range of scales prevalent in near-surface turbulence. Thus, the question is as follows: Are simulations of global convection at a stage where we can trust them to accurately reproduce large-scale flows? And if not, then can we place observational constraints to eliminate or differentiate one calculation from another?

In this effort, we have tried two other numerical setups (not reported here)—the first being a nonrotating high-resolution simulation (1536 × 3072 × 512), with the upper boundary at *r*/*R*_{⊙} = 0.99 and another with the same boundary and resolution, but with rotation included. Neither of these appeared satisfactory to us because the first had no rotation and was therefore isotropic across the 2-sphere; the second showed antisolar rotation (i.e., with poles spun up relative to the equator), and we therefore discarded both. The reason is that high-resolution calculations tend to have high Rossby numbers, lower viscosities and thermal diffusivities; the latter two effects introduce large (relatively speaking) convection velocities. Further, to obtain solar-like (fast equator) differential rotation, the Rossby number must be kept as low as is possible (*28*). Consequently, we settled on a low-resolution case with the upper boundary at *r*/*R*_{⊙} = 0.94 and solar-like rotation, i.e., equatorial spin up. To reiterate, this is just a small subset of the range of simulations that are possible with a variety of different codes.

The numerical simulation we use here (*29*) resolves rotating (Fig. 1, A and B) convection using 384 × 192 × 64 points, along the longitude, latitude, and radius, respectively, with *r*/*R*_{⊙} ≈ 0.96 serving as the upper boundary of the computational domain. The hydrodynamic equations in three-dimensional spherical geometry (with the entropy equation instead of the total energy equation) are solved using a fourth-order central derivative scheme and stepped forward in time with the Runge-Kutta method. We also apply a slope-limited artificial diffusion to all variables to stabilize the calculation. Explicit diffusion is applied for velocity and entropy, whose values at the top boundary are 1 × 10^{13} cm^{2} s^{−1} and 2 × 10^{13} cm^{2} s^{−1}, respectively, and proportional to ρ^{−1/2} (where ρ is the density), to reproduce a rotation profile with a relatively fast rotating equator in comparison to the pole. In comparisons presented here, the static *t* = 0, σ = 0 component of the numerical data is removed.

## RESULTS

Figure 3 (A and B) shows that, in contrast to simulations where power is largest at low spatial wave numbers and temporal frequencies, observationally inferred velocity monotonically decreases with decreasing wave number and frequency. This puzzling trend, in contrast to classical turbulence, where small wave numbers are energy containing, has been the focus, albeit inconclusive, of a variety of theoretical and computational efforts (*4*, *30*–*32*). While the disparity in amplitude is not as significant as suggested by prior time-distance methods (*5*), the two behaviors are markedly different. To benchmark the seismic results against independent measurements of the surface toroidal-flow spectrum obtained from Doppler data, we compare the inferred power spectrum at *r*/*R*_{⊙} = 0.995 with surface-Doppler estimates (*33*) in Fig. 3C. Specifically, we use figure 12 of (*33*), multiplying it by *33*)]. One important difference is that the temporal frequency ranges over which the sums of flow power are taken are different, ranging from 0.3 to 11 μHz for the surface Doppler curve and 0 to 3.05 μHz in the present seismic analysis. Nevertheless, the curve provides an independent confirmation that our seismic inferences are possibly in the correct ballpark. Figure 3D compares the theoretical and inferred root mean square (RMS) velocities over this frequency range, finding the discrepancy to be much smaller than previously reported, for this frequency and wave number range and for these depths. Figure 3E graphs this ratio as a function of harmonic degree, indicating that simulations significantly overestimate the low wave numbers.

Spectra as a function of *s* − ∣*t*∣, which capture the geometry of convection, show a dominance of sectoral modes in the Sun, i.e., small *s* − ∣*t*∣ (Fig. 4, A and B), with prograde (*t* > 0) and retrograde (*t* < 0) components indistinguishable to within error bars (for the tracking frequency of 453 nHz). In contrast, simulations contain a much broader set of convective modes. This is more markedly seen in Fig. 4 (C and D), which highlight the differences between the latitudinal profiles of simulated and observed flow, where we plot power as a function latitude for simulations (Fig. 4C) and seismology (Fig. 4D). The error bars on the observed curves are invisible on this plot, whereas the simulated curve is very noisy: This occurs because the low spatial and temporal wave numbers are preferentially powered in the simulations, thereby requiring a much longer integration time to achieve convergence. It is interesting that the most vigorous toroidal flows are cospatial with the fastest rotating layers, i.e., differential rotation and toroidal power peak at the equator. Whether anisotropic toroidal flow power induces cospatial rotational shear, or the other way around, remains to be investigated. The similarity with recent measurements of Rossby waves in purely sectoral channels [e.g. (*21*, *22*)] tempts one to speculate if there are more general underlying processes at play. The absence of nonsectoral modes in toroidal flow may be due to multiple reasons: low-amplitude excitation, exclusion of tesseral/zonal modes due to specific boundary conditions, or rapid damping (possibly by rotational shear). The most vigorous toroidal flows appear to lie within the spherical wedge delineated by the intersection of the cylinder tangent to the radiative interior with the sphere, especially so in the near-surface layers.

Figure 5, a reconstruction of flow images from observations, shows that equatorially localized north-south–aligned flows are preferentially powered—these are the high-amplitude sectoral modes, i.e., for *s* − ∣*t*∣ ≲ 12 in Fig. 4 (A and B). Simulated velocity images qualitatively differ from the inferences, with large-scale features distributed across all latitudes (Fig. 5, C and D). In contrast, the small scales appear to be preferred in the Sun (Fig. 5, A and B), due to the inferred power spectrum increasing continuously as a function of wave number. This spectral dependence results in the highest spatial frequencies having the most power, whereas the simulated spectra peak at low spatial frequencies and decrease with increasing wave number.

In Fig. 6 (A and B), we show the dependence of toroidal flow power on temporal frequency—the observed velocity scale increases continuously with frequency, broadly in line with the idea that convection at higher frequencies is more vigorous than at lower frequencies, in contrast with computation (Fig. 6C). Figures 3 (A to C) and 6 (A and B) together indicate that solar flow power is continuously increasing as a function of wave number and frequency (over the ranges considered here), eventually overtaking and exceeding the simulations. Prior comparisons [e.g., (*5*)], which relied solely on amplitudes, thus miss the point that velocity-amplitude comparisons are sensitive to the spatial and temporal frequency bands under consideration and that, even if the comparisons were favorable, simulations and inferences could still be capturing intrinsically different physics (e.g., classical versus nonclassical turbulent systems, as discussed early on in this section). In the present work, we highlight the importance of making detailed comparisons with observations to improve faith in theoretical and computational models.

Figure 6D tracks “deep” and shallow convective velocities, i.e., the inferences at radii *r*/*R*_{⊙} = 0.90,0.99, as a function of the magnetic cycle. It is conceivable that these turbulent flows may be weakened during periods of high solar activity because magnetism suppresses free movement of plasma (*34*). However, Fig. 6D shows that there is no apparent correlation between surface and interior velocities with the sunspot cycle (*35*).

## DISCUSSION

In summary, the state-of-the-art numerical simulation used here (*29*) reproduces certain features such as overall amplitudes—to within a factor of a few—of the toroidal flow but are not good models of solar turbulence, as evidenced by their inherent difference from seismic observations: (i) Large-scale flows are preferentially powered in simulations, whereas they are suppressed in the Sun, and (ii) simulated flows are weak at the equator and vigorous at high latitudes, in contrast to inferences, which show that they are confined to the equatorial regions and weak at high latitudes.

Differential rotation is thought to be sustained by a combination of a thermal wind and Reynolds stresses (*3*). A thermal wind could be induced by a latitudinal entropy gradient, in turn caused by differences in convection between the poles and equator (*36*), resulting in changes to the turbulence. Our inferences here of latitudinally nonuniform toroidal flow appear to match these symptoms, although they do not directly point to the cause. It would be very useful to theoretically study the impact of thermal winds on toroidal flows in future investigations. Inferences of Reynolds stresses would go a substantial distance toward discerning the origin of differential rotation, both in the Sun and stars (*37*), which state-of-the-art simulations are not successful at identifying. However, these are much more challenging to measure since they require the complete set of phases and amplitudes of toroidal and poloidal flows at all spatial and temporal scales (see the Supplementary Materials). Since inferences of toroidal flows here are obtained from studies of same-ℓ couplings, issues arising from leakage are significantly mitigated. Modes for ℓ < 150 or so are typically well separated in frequency (on the order of several to tens or more of microhertz), and attempts to infer perturbations varying on temporal frequencies much smaller than the nominal frequency separation between pairs of modes, as used in the analysis here, are hindered by the influence of spatial leakage and background power on the measurement. As a consequence, cross-ℓ couplings must be interpreted much more carefully and are deferred to a future study. The penalty for considering only same-ℓ couplings is that we are only able to infer odd harmonic degree toroidal flows and have practically no sensitivity to poloidal components.

Convection comprises spatially complex and temporally evolving structures; testing different theories of this process requires building a framework with which to make consistent comparisons. The theory and computation of convection are sensitive to a large variety of necessary inputs, e.g., parameters and boundary conditions, and eventual solutions may vary substantially. Determining whether a simulation is “correct” is for instance a difficult question—would the successful recovery of differential rotation and meridional circulation be sufficient cause to declare success or should we seek a greater degree of agreement, i.e., turbulence structure functions, Reynolds stresses, etc.? We have shown, for example, that convective amplitudes alone, while providing some broad indications of the success of a theory, are insufficient for making accurate and meaningful comparisons. To eliminate or validate a theory may thus demand agreement with both large-scale circulations (e.g., differential rotation) and other constraints, such as provided here. Thus, the present disagreement should be viewed as a significant step toward building a framework with which to make comparisons of theory and computation with observation and adding to the range of available measurements of convection. One lesson from the direct comparison with simulations here (*29*) is that high resolution alone is insufficient to achieve agreement and that a great distance still remains to be traversed between computation and observation.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/30/eaba9639/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**S.M.H. wishes to thank T. Larson for retrieving MDI and HMI observational data from the JSOC pipeline and is grateful to M. Woodard for many conversations on the topic that helped motivate and clarify several basic issues. S.M.H. also thanks A.C. Birch for useful discussions and P. Subramanian for a careful reading of the manuscript, pointing out some errors and providing a useful plot in the Supplementary Materials. All calculations in this analysis were carried out at the SEISMO cluster at TIFR and the Stanford University solar group computational facility. The numerical simulation was carried out on the FX100 supercomputer at Nagoya University. Last, we thank two anonymous referees whose careful and critical reading of the manuscript and constructive comments helped to substantially tighten the arguments presented here.

**Funding:**S.M.H. acknowledges support from Ramanujan Fellowship SB/S2/RJN-73 and the Max-Planck Partner group program. This research was supported, in part, by the International Centre for Theoretical Sciences (ICTS) during a visit for participating in the program—Turbulence from Angstroms to Light Years (ICTS/Prog-taly2018/01). H.H. was supported by MEXT as “Exploratory Challenge on Post-K Computer” (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar System), the Project for Solar-Terrestrial Environment Prediction, toward a unified view of the universe: From large scale structures to planets, and MEXT/JSPS KAKENHI grant number JP18H04436 and JP20K14510.

**Author contributions:**S.M.H. designed the helioseismic analysis and the framework for comparing with simulations, interpreted results, and wrote the manuscript. K.R.S. helped design the research, interpret results, and wrote the paper. H.H. performed large-scale numerical computation to simulate convection, helped interpret results, and wrote the paper.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are available publicly for download from Stanford University’s Joint Science Operations Center (JSOC) pipeline. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).