Changes in the near-surface shear layer of the Sun

We use helioseismic data obtained over two solar cycles to determine whether there are changes in the near-surface shear layer (NSSL). We examine this by determining the radial gradient of the solar rotation rate. The radial gradient itself shows a solar-cycle dependence, and the changes are more pronounced in the active latitudes than at adjoining higher latitudes; results at the highest latitudes (greater than about70 degrees) are unreliable. The pattern changes with depth, even within the NSSL. We find that the near-surface shear layer is deeper at lower latitudes than at high latitudes and that the extent of the layer also shows a small solar-cycle related change.


INTRODUCTION
Helioseismic inversions for solar rotation reveal two distinct shear layers -a very prominent one that is near the base of the convection zone, and one close to the surface (Schou et al. 1998). The former, named the "tachocline" (Spiegel & Zahn 1992) is the main seat of the solar dynamo (see Gilman 2005;Charbonneau 2010, and references therein) in many solar dynamo models. The latter, the near-surface shear layer (NSSL; Thompson et al. 1996), is still somewhat of a mystery. The first evidence for the NSSL came from the observation that emerging active regions rotate faster than the surrounding photosphere (Foukal 1972), helioseismic analyses have since confirmed the presence of the NSSL, and it has been a subject of many helioseismic analyses, whether implicitly in studies of solar zonal flows (e.g., Antia & Basu 2001;Vorontsov et al. 2002;Basu & Antia 2003;Howe et al. 2006;Antia et al. 2008;Howe et al. 2018;Basu & Antia 2019, etc.) or in focused studies of this layer (e.g., Antia & Basu 2010;Barekat et al. 2014Barekat et al. , 2016. It has been argued (Brandenburg 2005;Pipin & Kosovichev 2011) that the NSSL plays a role in a distributed dynamo. Brandenburg (2005) also argued that the radial-shear in the NSSL may explain the migration of the solar-activity belt towards the equator as the solar cycle progresses. Increasingly, the role that this layer can play in solar dynamos is being studied (e.g., Dikpati et al. 2002;Mason et al. 2002;Käpylä et al. 2006;Karak & Cameron 2016;Paradkar et al. 2019;Jha & Choudhuri 2021, etc.).  studied the gradient of the rotation rate in the NSSL and found it to be nearly constant with latitude until about 30 • ; this was subsequently confirmed by Barekat et al. (2014) who found that the gradient is constant to much higher latitudes, up to about 60 • . Both these studies only used fmode data, limiting their results to very shallow depths. Antia et al. (2008) studied changes in the gradients of the rotation rate throughout the convection zone; while they did not explicitly comment on the NSSL, they showed that at 0.95R ⊙ the radial gradient of the rotation rate shows a pattern of migrating band as a function of latitude and time that is similar to that of the zonal flows. Antia & Basu (2010) examined the radial extent of the NSSL and found some solar cycledependent changes, though this study used data covering only one solar cycle. They did not examine any changes to the rotation-rate gradient. Barekat et al. (2016) examined the rotation rate gradient as a function of time and found a well-defined time difference, again with a banded pattern migrating towards the equator with time. They found that sunspots are concentrated at latitudes where the gradient is lower than the average gradient. They claimed that this result contradicts what Antia et al. (2008) found, i.e., sunspots are concentrated at latitudes where the gradient is higher than average; the Barekat et al. (2016) and Antia et al. (2008) results, however, were at different depths inside the Sun, indicating that the change of the rotation-rate gradient is not constant over the NSSL. The NSSL extends up to about 0.95R ⊙ , and thus a time-variation that depends on radius is quite possible. Even the zonal flow pattern at low latitudes show a strong variation in phase around 0.95R ⊙ (Vorontsov et al. 2002;Antia et al. 2008). Antia et al. (2008) compared the variations of the gradient with sunspot locations at 0.95R ⊙ ; the results of Barekat et al. (2016), based only on f-mode data, are for a radius of around 0.99R ⊙ .
In this paper, we analyze helioseismic data obtained over two solar cycles to investigate changes in the NSSL as a function of time at different latitudes. Unlike , Barekat et al. (2014Barekat et al. ( , 2016, we do not limit ourselves to using only the fmodes; this allows us to determine changes throughout the NSSL and also determine whether the extent of the shear-layer changes with the solar cycle.
The rest of the paper is organized as follows. We describe the data used in Section 2, the analysis technique is described in Section 3, we present our results in Section 4, and we discuss the implications of our results in Section 5.

DATA USED
For this work, we use solar oscillation frequencies obtained by the ground-based Global Oscillation Network Group (GONG) (Hill et al. 1996) and the space-based Michelson Doppler Imager (MDI) on board the Solar and Heliospheric Observatory spacecraft (Scherrer et al. 1995) and the Helioseismic and Magnetic Imager (HMI) (Scherrer et al. 2012) on board the Solar Dynamics Observatory.
The data we use from the GONG project cover a period from 1995.05.05 to 2021.01.25. The data are designated by GONG "months", each "month" being 36 days long. Solar oscillation frequencies and splittings of sets starting GONG Month 2 are obtained using 108-day (i.e., 3 GONG months) time series. There is an overlap of 72 days between different data sets, i.e., GONG Month 2 frequencies were obtained from data of GONG Months 1, 2 and 3, those for GONG month 3 from GONG Months 2, 3 and 4, etc. We use data for GONG Months 2 to 260. These data are available from the GONG web-site 1 . GONG datasets are restricted to modes with l ≤ 150. It should be noted that very few modes covered by the GONG data have lower turning point above 0.98R ⊙ and hence it is not possible to infer the rotation rate above this layer using these data sets.
Data from MDI cover the period from 1996.05.01 to 2011.04.24. Solar oscillation frequencies and splittings for these data were obtained from 72-day time series, and the sets have no overlap in time. HMI started obtaining data on 2010.04.30 and these are also obtained from 72-day time series. MDI and HMI sets have p-mode data for l ≤ 200 and f-mode data from l ≈ 120 to l = 300. These data are available to all researchers through the Joint Science Operations Center of the Solar Dynamic Observatory 2 . We have used MDI data from series mdi.vw V sht modes, and HMI series hmi.V sht modes. We combine the MDI and HMI datasets to cover both solar cycles 23 and 24. There is an overlap of 1 year between the MDI and HMI data; we use HMI data for this period. The last HMI dataset that we have used is set 10288, which has an end date of 2021.05.13.
We compare the helioseismic results with the Sunspot butterfly diagram. For this purpose, we use data on the date of emergence and latitude of the emergence of active regions from the National Oceanic and Atmospheric Administration's Solar Region Summary.

ANALYSIS TECHNIQUE
We invert each of the data sets to determine the rotation rate as a function of depth and latitude. We adapt the 2D Regularized Least Squares (RLS) technique described by Antia et al. (1998) to perform the inversions. The rotation rate is represented as a product of cubic B-splines in radius and cos(ϑ), where ϑ is the colatitude. Symmetry between the two hemispheres is ensured through a boundary condition at the equator. We use 20 knots uniformly spaced in cos(ϑ) and 50 knots uniformly spaced in the acoustic radius. The smoothing parameters were chosen to be λ r = 0.0025 and λ θ = 0.04. Iterative refinement of the solution was done as described by Antia et al. (1996), and a fixed number of 9 iterations were performed.
For the GONG data we use the first 8 odd-order splitting coefficients, c 1 , c 3 , . . . , c 15 , while for MDI and HMI data we use 18 splitting coefficients, c 1 , c 3 , . . . , c 35 . We determine the dimensionless radial gradient of the rotation rate ∂ log Ω/∂ log r directly from the B-splines; the use of cubic B-splines obviates the need to calculate nu-merical derivatives. This is different from the approach of Antia et al. (2008) and Antia & Basu (2010), where the derivative was calculated by numerical differentiation of the rotation rate. Analytic derivatives reduce the noise in radial-gradient estimate and allows us to study the changes in the gradient in the NSSL.
We determine temporal changes in the rotation-rate gradient by subtracting out the average gradient for each solar cycle, i.e., where, r, θ and t are respectively radius, latitude and time, Ω is the rotation rate, and the term in the angular brackets represents the time average over a solar cycle. This residual shows large values at high latitudes and to balance out the variation, we multiply the residual by cos θ in all figures. Note that unlike Antia et al. (2008), Antia & Basu (2010), and Basu & Antia (2019), we use the dimensionless gradient for this work. This allows us to compare our results easily with those of  and Barekat et al. (2014Barekat et al. ( , 2016. Guided by our earlier work on zonal flows (Antia & Basu 2013), we average each solar cycle separately and subtract that from data sets of that cycle.
Since we have very little data for solar Cycle 25, we subtract the Cycle 24 average from those sets. To smooth out fluctuations in the time series at each r, θ, we apply running mean over ±6 months from the central point in time.
We define the extent of the near-surface shear layer as the region where the magnitude of the dimensionless gradient is greater than 20% of its maximum value at each epoch. The maximum is usually at the solar surface except at high latitudes; at high latitudes, the maximum gradient sometimes occurs a little below the surface. This allows us to define the position of the base of the shear layer, r s (θ). We have tried different definitions for the extent of NSSL. A smaller value, e.g., 10% causes problems at some epochs at the higher latitudes -the value of r s for those epochs and latitudes suddenly becomes small (i.e., the shear-layer suddenly deepens), perhaps because of noise and also because sometimes the sign of gradient at high latitude does not change until we reach the tachocline region. Since such a sudden change is not expected, the limit of 10% is clearly unsuitable. The 20% value also allows us to be consistent with our previous work (Antia & Basu 2010), except that here we use the dimensionless gradient. The use of a higher limit, e.g., 30% doesn't lead to any sig-nificant difference, except that the average depth of the layer is reduced. It could be argued that the temporal variation in the extent of NSSL is caused by variation of the maximum value. To discount this possibility, we also considered one case where a fixed maximum magnitude of 1 was used to define the extent of the NSSL. Even in this case, the observed temporal variation is similar. Consequently, we retain the definition of Antia & Basu (2010). We also study the temporal variations of r s (θ) by subtracting an average over each solar cycle, just like that for the gradient.

The radial rotation-rate gradient in the NSSL
The time-averaged gradient of the rotation rate in the NSSL is shown in Fig. 1, since GONG data sets have only a few modes with lower turning points above 0.98R ⊙ , we do not show GONG results at 0.99R ⊙ . Note that the magnitude of the gradient decreases as one goes deeper inside the Sun. Like Barekat et al. (2014), we find that the gradient is almost constant up to nearly 60 • latitude, and deviate from the constant at higher latitudes. Our results show a steepening of the gradient with an increase in latitude; it is possible that this is because we do not restrict ourselves to using f-mode data and use p-mode data too. We also find that the MDI and HMI results do not agree at high latitudes implying that there may not be much significance of the deviation of the gradient from −1 at these latitudes. The GONG results show some oscillations in the gradient as a function of the latitude. These are due to the fact that GONG data include only 8 odd-order splitting coefficients and are restricted to modes with l ≤ 150. To confirm this, we repeated the analysis for HMI data using the same restriction on modes and splitting coefficients and find that they also show similar oscillations. These results are also included in the figure. These oscillations are most likely due to the fact that the polynomials in cos ϑ corresponding to the first 8 splitting coefficients are not sufficient to represent the latitudinal variation in the gradient.
Just because the average gradient is almost constant with latitude does not mean that the time-variation of the gradient is the same at all latitudes. The change in the gradient (multiplied by cos θ) at 0.99R ⊙ is shown in Fig. 2; this is the radius that corresponds closest to the depth of the Barekat et al. (2016) results making comparisons easier. As can be seen in the figure, the gradient in the NSSL shows a similar banded pattern as the zonal flows (Howe et al. 2006;Antia et al. 2008;Basu & Antia 2019;Howe 2020, etc.). We also see that the sunspots are confined to regions where the gradient Figure 1. The average gradient in the near-surface shear layer for MDI, HMI and GONG sets at different radii plotted as a function of latitude. The GONG data are averaged over Cycle 23; the results are not significantly different for Cycle 24. We also show HMI results obtained with datasets restricted to l ≤ 150 modes and 8 splitting coefficients. Note that the absolute value of the gradient decreases deeper inside the Sun, and that there are significant differences between the GONG, MDI and HMI results at high latitudes.
is lower than average. This agrees with the result of Barekat et al. (2016). Barekat et al. (2016) had remarked on the disagreement between their results and the results of Antia et al. (2008) with regard to the pattern of concentration of sunspots with respect to that of the radial gradient. The apparent disagreement is most likely because the results were obtained at different depths, and the phase of the pattern has a steep variation with depth in this range. We can see from Fig. 3, where we show our current results at r = 0.95R ⊙ , that the pattern is reversed compared with that at 0.99R ⊙ . At this radius, sunspots coincide with the positive band, as found by Antia et al. (2008) and Basu & Antia (2019). This implies that the phase of temporal variation of the rotation gradient has a radial dependence even within the thin NSSL.
In Fig. 4 we show the change of gradient (again multiplied by cos θ for consistency with the other figures) at 0.97R ⊙ at a few selected latitudes. The figure clearly shows quasi-periodic changes at low latitudes. Although the high-latitude results are dominated by systematic errors in the data sets, there appears to be a hint of periodic behavior even at these latitudes. What we do not completely understand is the disagreement between GONG and MDI+HMI results at 15 • ; it is possible that this is because GONG results have poor resolution at this depth -a similar figure plotted at 0.95R ⊙ does not show such a disagreement at low-and mid-latitudes. However, another reason for the disagreement is probably the slight oscillations in the GONG average gradients (Fig: 1) that feed into the results of the change in the gradient. These oscillations at 0.97R ⊙ are a result of using only 8 splitting coefficients to obtain results with GONG data, compared with 18 coefficients for MDI/HMI. Thus, the MDI and HMI results may be better in terms of examining the time variation. Figures 2 and 3 imply that there is a depth dependence in the change of the gradient, this is shown in Figs. 5 and 6 for latitudes of 20 • and 40 • respectively. At 20 • we see a change in the sign of the residual as one goes deeper in radius.

The extent of the NSSL
The depth of the shear layer will depend on the physical conditions inside the convection zone, and may show variations when the conditions change, which happens with increasing magnetic activity. Hence, a study of the extent of NSSL may give information about the conditions in the near-surface layers of the Sun and provide constraints on theoretical models of NSSL. The extent of the shear layer has a latitude dependence, and this can be seen in Fig. 7. The shear layer is deepest at the equator and becomes shallower with increasing latitude. We also find systematic differences between results obtained with GONG data and those obtained with MDI and HMI data. The differences between GONG and MDI have been found earlier (see e.g., Antia & Basu 2010); there are also differences between GONG and HMI data, as well as MDI and HMI data (Basu & Antia 2019). One possible reason why GONG data give a deeper estimate of the NSSL is GONG datasets have very few modes with lower-turning points shallower than 0.98R ⊙ and consequently, the surface value is essentially determined by extrapolation due to the regularization (smoothing) used in the RLS technique. As a result, the surface value of the gradient, which is used in estimating the depth of the shear layer, may not be correctly determined. The GONG results also show oscillations in latitude. These are caused because GONG sets are limited to 8 splitting coefficients; restricting the HMI sets to 8 splitting coefficients causes the HMI result to be oscillatory, and this can also be seen in Fig. 7. In addition to introducing oscillations as a function of latitude, restricting the mode set and the number of coefficients makes the inferred depth of the shear layer larger. There is also a discrepancy between the MDI and HMI results at high latitude, which produces a shift in r s above 60 • latitude; this is most likely to be a result of the known systematic  differences between MDI and HMI splitting coefficients (see Fig. 1

of Basu & Antia 2019).
The most well-defined solar-cycle related changes in the depth of NSSL are confined to low latitudes. The position of the bottom of the NSSL, r s , shows a change in the active latitudes, as can be seen in Fig. 8; the changes are small just above the active latitudes, but appear to increase again at the highest latitudes. The nature of the shear layer appears to be complicated at high latitudes, though the reliability of the results here is questionable. Antia & Basu (2010) also found a timevariation in r s , but with data available for only one solar cycle and with the numerical differentiation used resulting in larger uncertainties, the changes were not clear. The time variation also shows discrepancies at high latitudes.      The change in r s at different latitudes is clearer if we subtract the time-average of r s and plot the residuals, i.e., r s − < r s >. The results are shown in Fig. 9. As can be seen, the changes in r s also show a migrating pattern, that is similar to, but not as well-defined, as the pattern of the changes in the gradient. The appearance of sunspots coincides with latitudes where the NSSL become shallower than average, i.e., where r s increases.

DISCUSSION
Our investigation of the NSSL layer using both f-and p-modes of solar oscillations show that the gradient of solar rotation is mostly independent of latitude, but depends on radius -the absolute value of the gradient decreases with radius. This is not surprising since the gradient is small below the NSSL, i.e., below about 0.95R ⊙ .
At low latitudes, our results agree with earlier helioseismic analyses, i.e., we find a logarithmic gradient of −1. The behavior of the NSSL at high latitudes is uncertain -different data sets give different results (see Fig. 1). Results obtained with GONG, MDI and HMI differ for latitudes higher than about 60 • . Barekat et al. (2016) had examined systematic differences between MDI and HMI results and had concluded that the 72day high-degree HMI data suffer from systematic errors. However, we find that the HMI and GONG results are more similar -this is what Basu & Antia (2019) also found as seen in Fig. 1 of their paper.
It should be noted that the behavior of the , Barekat et al. (2014) and (Barekat et al. 2016) results are very different from ours at high latitude; their results, obtained with only f modes, show a decrease in the absolute value of the gradient at high latitudes, our results show the gradient becoming steeper around 60 • . It is likely that the inclusion of p-mode data changes the inferred high-latitude behavior of the gradient in the NSSL, and that there may be systematic differences between f-and p-mode data. There could also be systematic differences between the results obtained using p-modes and f-modes because f-modes of high degree are confined very close to the solar surface, where the reference model used to calculate the kernels may not be good enough. The difference between our high-latitude results and those of  and Barekat et al. (2014Barekat et al. ( , 2016 implies that the results cannot be trusted. It may be possible to determine the gradient at the shallow depths at high latitudes with local helioseismic data from off-ecliptic missions like Solar Orbiter or the proposed polar mission, Solaris 3 (Hassler et al. 2020(Hassler et al. , 2021 which aims to observe the solar polar regions for 3 months. The gradient in the NSSL changes with the solar cycle. The pattern of change is reminiscent of the zonal flow pattern, but with one major difference -in the near-surface zonal flow pattern, the low-latitude bands migrate towards the equator, while the high-latitude one move towards the poles, however, the change of gradient in the NSSL shows both low-latitude and high-latitude bands moving towards the equator. This was also seen by Barekat et al. (2016). This is clearer in the deeper layers of the NSSL; in the zonal flows, this feature is seen much deeper, at about 0.8R ⊙ (Basu & Antia 2019). We also find that not just the gradient of the NSSL, but the extent of the NSSL also changes with the solar cycle, and here too, the change shows a migrating pattern with time, with bands moving towards the equator.
Recreating the NSSL correctly in simulations of solar rotation has been challenging. Global MHD simulations are usually done under the anelastic approximation which does not allow simulating rotation and convection correctly in the outer few per cent of the Sun (see, e.g., Miesch 2000), missing out a large part of the NSSL. Such simulations show an NSSL if the density contrast is increased, but even then a shear layer is obtained only at low latitudes (Matilsky et al. 2019). Fully compressible radiative-hydrodynamic simulations in sections of spherical shells are able to form an NSSL close to the surface, but again, only at low latitudes (Robinson & Chan 2001). Kitchatinov & Rüdiger (2005) find that they can form an NSSL using a quasi-linear theory. However, getting the correct magnitude of the gradient at different latitudes has not yet been possible, even with massively parallel, fully compressible codes. There has been some work done to examine the effects of magnetic fields on the NSSL, notably Kitchatinov (2016) found that the shear in this layer is sensitive to the toroidal magnetic fields in these layers, and can serve as a probe for fields of the order of one Kilo Gauss. Our results should, therefore, give constraints on the sub-surface magnetic fields. The change in the extent of the NSSL with activity is also a constraint. The presence of the NSSL has been explained as being due to Reynold stresses (Robinson & Chan 2001) and the change in the gradient because of the suppression of turbulent viscosity by the magnetic field (Kitchatinov 2016), thus the suppression of turbulence by magnetic fields could potentially explain changes in the NSSL, but models of the NSSL need to explain both the change in the gradient and the change in the extent of the NSSL.