© ESO, 2017

1. Introduction

Active regions are surrounded by spatially extended near-surface inflows (Gizon et al. 2001; Haber et al. 2002; Gizon 2004; Haber et al. 2004; Hindman et al. 2004; Komm et al. 2007; Hindman et al. 2009; Braun & Wan 2011) with typical flow velocities of 20 − 30 m/s (velocities up to 50 m/s have been reported) that extend up to 10° from the active region. Models suggest that these flows arise due to an enhanced cooling rate in the active regions (Spruit 2003; Gizon & Rempel 2008; Hurlburt & DeRosa 2008). There are indications of outflows at greater depths (see, e.g., Gizon et al. 2010, and references therein).

In flux-transport simulations, the inflows to active regions play an important role in the transport of magnetic flux over the course of the solar cycle. The inflows can affect the dispersal of active regions since they counterbalance the outward diffusion of magnetic flux by convection (De Rosa & Schrijver 2006; Martin-Belda & Cameron 2016). They are also a potential mechanism for modulating the strength of the solar cycle (Cameron & Schüssler 2012; Shetye et al. 2015; Martin-Belda & Cameron 2016).

So far, this inflow has only been studied using local helioseismology. Here we want to measure the inflows to active regions with local correlation tracking (LCT) of solar granulation. Starting from continuum images provided by the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012), we generate a time series of full-disk maps of the horizontal velocity at the solar surface by using the LCT code Fourier local correlation tracking (FLCT; Welsch et al. 2004; Fisher & Welsch 2008).

In order to demonstrate the viability of our technique, we measure a known flow of similar magnitude, the torsional oscillations (Howard & Labonte 1980), which have contemporaneous measurements by a different technique. In addition, we present results for individual active regions and for the average of 243 regions.

2. Generating LCT full-disk flow maps

2.1. Running the FLCT code

We computed maps of the horizontal velocity on the solar surface using the FLCT code (Welsch et al. 2004; Fisher & Welsch 2008). The FLCT code provides the velocity by cross-correlating the granulation pattern between two consecutive intensity images. This code has been used for various applications, such as studying flows connected to solar activity (e.g., Welsch et al. 2009; McKenzie 2013; Su et al. 2013) or systematic errors inherent to local correlation tracking (Löptien et al. 2016b).

We ran the FLCT code directly on full-disk images of the continuum intensity provided by HMI. The time series stretches over six years of HMI data, from 24 April 2010 to 27 April 2016. In order to save computation time, we did not make use of the full cadence of HMI. We applied the FLCT code to pairs of continuum images; the images in the individual pairs were separated by 45 s, but the pairs were separated by 30 min. This results in one flow map every 30 min. This setup slightly increases the noise in the flow maps that originates from the proper motion and evolution of the granules. However, the dominant part of the granulation noise is located at significantly smaller spatial scales than are relevant for this study (see also Fig. 3 in Löptien et al. 2016a). The FLCT code divides the images into several subimages, whose size is given by the parameter σ. Here we used σ = 6 pixel.

Since we did not remap the intensity images before running the FLCT code, the velocities provided by the code do not correspond to physical velocities on the Sun. Instead, they are defined on the grid of the CCD and are in units of pixels per second with v_m denoting the velocity in the horizontal direction (pointing from left to right) and v_n denoting the velocity in the vertical direction (pointing from bottom to top). We converted the velocities from CCD to heliographic coordinates later in the analysis.

Local correlation tracking generally underestimates the true flow velocities on the Sun (see, e.g., Švanda et al. 2007; Verma et al. 2013; Löptien et al. 2016a). This effect depends on the spatial resolution of the data, hence it varies across the disk. We corrected for this effect by generating calibration data (Löptien et al. 2016a) both for v_m and v_n. These calibration data consist of a time series of continuum intensity images where we shifted the individual images using Fourier interpolation in order to add a constant flow. The velocities estimated by the FLCT code depend linearly on the input velocities, with the slope ranging between 0.80 and 0.89.

The sensitivity of LCT also depends strongly on the magnetic field strength. For strong fields (B ≳ 500 G), the sensitivity of LCT drops significantly (to ~ 0.5) since granulation is suppressed by the magnetic field. This reduction in the sensitivity cannot be easily corrected for because it depends strongly on disk position. This causes the background subtraction (see next section) to not work in regions with strong magnetic field. Hence, our LCT velocities are not reliable within and near sunspots.

2.2. Subtracting the background signal

The LCT full-disk flow maps are dominated by a large-scale background signal (see the mean images in the top panel of Fig. 1), which needs to be removed in order to study the flows around active regions. This background consists of real solar flows (differential rotation and meridional flow) and a systematic error that looks like a flow converging towards disk center, also known as the shrinking-Sun effect (Lisle & Toomre 2004; Löptien et al. 2016b).

The background flow is time-dependent, due to the spacecraft motion of the Solar Dynamics Observatory (SDO). In particular, the way that the flow maps are affected by solar rotation depends strongly on the heliographic latitude of disk center (B₀-angle). The position of a given solar latitude (with its specific rotation velocity) on the CCD changes with the B₀-angle. In addition, in the case of large values of the B₀-angle, differential rotation also affects the flow maps for v_n. The shrinking-Sun effect is also affected by the orbit of SDO. The B₀-angle predominantly exhibits an annual variation due to the Earth’s orbit around the Sun. The orbit of SDO around the Earth causes additional temporal variations of the background flow on a timescale of 24 h. Because of these temporal variations, we were not able to remove the background signal by subtracting a mean image from the flow maps.

The orbit of SDO is known, so in principle, it should be possible to model its influence on the LCT velocities. However, this approach is not able to fully reproduce the 24 h variations (compare with Schuck et al. 2016; Couvidat et al. 2016, who tried to model the influence of the orbit of SDO on the Doppler velocity). In particular, the response of the shrinking-Sun effect to the orbit is only poorly understood as this systematic effect is extremely sensitive to changes in the spatial resolution. Instead, we fit an empirical model to the background signal that varies smoothly in time and space. We modeled the background flow with Zernike polynomials (Zernike 1934), a set of orthogonal basis functions on the area of a circle (see Appendix A). In our model, the coefficients of the individual Zernike polynomials vary smoothly in time in order to account for variations due to the orbit.

Removing the background signal using Zernike polynomials consists of several steps. First, we smoothed the flow maps by convolving them with a Gaussian (σ = 8 pixel) followed by subsampling them 4 × 4 in order to reduce the computation time. Then, we decomposed each individual flow map both for v_m and v_n in Zernike polynomials. Since LCT becomes less reliable when used close to the limb, we fit the Zernike polynomials to the inner part of the disk (up to a heliocentric angle of 60°). We did not use data at heliocentric angles greater than 60° in the further analysis. The center of the disk used for the fitting is the center of the solar disk determined from the limb, not the center of the CCD. The Zernike polynomials are not only sensitive to the background signal, but also to the solar flows we are interested in. Hence, we did not subtract the Zernike polynomials from each flow map individually. Instead, our estimate of the background signal consists of a time series for each Zernike coefficient, stretching over all the data. We filtered these data in Fourier space, selecting only the frequencies that we considered to represent the temporal evolution of the background signal (i.e., signals with a period of 24 h or 1 yr and the corresponding harmonics, see Appendix B for more details). The signal at frequencies that are not affected by the orbit was not included in the time series. After filtering the time series of the Zernike polynomials, the mean and the standard deviation of the individual time series decreased significantly with increasing radial degree n of the Zernike polynomials. We used Zernike polynomials up to n = 7, where for all m both the mean and the standard deviation are ≤ 10^-5 pixel/s (this corresponds to ~ 3.5 m/s at disk center). We then computed the background signal for each flow map from the filtered time series of the Zernike coefficients and subtracted it. Figure 1 shows the mean flow maps with the background subtracted for v_m and for v_n.

Our estimate of the background signal only includes temporal variations on timescales that are affected by the orbit of SDO. The flows around active regions evolve on different timescales (days to weeks), meaning that they are not affected by the background subtraction. In principle, it would also be possible to mask the active regions before decomposing the flow maps into Zernike polynomials. However, in this case, the Zernike polynomials would no longer form an orthogonal basis. The leakage between the different polynomials would depend on what area on the disk is masked. This would be time-dependent since the active regions evolve and move across the disk due to solar rotation. In addition, masking a significant fraction of the disk would cause the Zernike polynomials to be more sensitive to convective noise, which increases the background noise in the power spectrum.

The mean image for v_n with the background subtracted (bottom right part of Fig. 1) exhibits a streakiness that is not present in the image for v_m. These features are caused by supergranulation. Since we did not track the data with the solar rotation when computing the mean images, individual supergranules propagate across the visible disk.

In the mean image for v_m in Fig. 1, a part of one of the columns (m = 2727) has a lower velocity than the neighboring pixels. This is a bad column of the CCD. We did not include these data in the analysis.

2.3. Remapping the LCT velocities

In the next step, we tracked and remapped the flow maps from the CCD coordinates to heliographic coordinates in order to generate Carrington maps of the flow velocities. The resulting flow maps have a spatial resolution of 0.4°. We also converted the flow velocities from the reference frame of the CCD (coordinate vectors ê_m and ê_n) to spherical coordinates (coordinate vectors ê_φ and ê_θ, see Appendix C). In addition, we generated Carrington magnetograms using the same tracking rate and the same spatial resolution as for the LCT flow maps.

3. Flows derived using LCT

3.1. Carrington maps

Figures 2 and 3 show examples of Carrington maps of the line-of-sight component of the magnetic field (averages of the magnetic field over the entire disk passage). We present results for three Carrington rotations at solar minimum (CR 2097, CR 2098, and CR 2099 in Fig. 2) with only a few active regions and three examples of Carrington maps at solar maximum (CR 2147, CR 2148, and CR 2149 in Fig. 3) with a high coverage by active regions. Since the velocities provided by LCT are dominated by velocities on small spatial scales (proper motion of granulation, supergranulation) we smoothed the flow maps by convolving them with a Gaussian (σ = 4°). In addition, we subtracted for each latitude the longitudinally and temporally averaged flow velocity (both for v_φ and v_θ).

The Carrington maps exhibit a complex flow pattern, both in the quiet Sun and near active regions. In the quiet Sun, there are several examples of cyclonic inflows and anticyclonic outflows. A prominent feature is visible in the Carrington rotations 2097−2099 at ~ 45° latitude and ~ 150° longitude. In Hathaway et al. (2013), this flow cell was also detected and identified as a giant cell. Most active regions in Figs. 2 and 3 are surrounded by converging inflows with typical flow velocities of 20−30 m/s. The inflows are far from being azimuthally symmetric, extending locally up to 10° from the active region and reaching flow velocities of up to 40 m/s.

3.2. Torsional oscillations

In order to test the quality of the LCT flow maps, we compared the torsional oscillations seen in our data and the flows seen by global mode helioseismology. These are inversions of the f-mode frequency splittings using HMI full-disk data (Larson & Schou 2017) and were processed as described in Schou (1999). Figures 4 and 5 show the flows in the φ-direction after subtracting the mean flow at each latitude from both data sets. The LCT data exhibits signal at smaller spatial scales than the f-mode data, and also contains the component of the torsional oscillations that is antisymmetric with respect to the equator (global helioseismology is only sensitive to the symmetric component of the flow). However, after symmetrizing the LCT data around the equator and smoothing them with the averaging kernels of the global helioseismology data, the results derived from the two methods are in good agreement and the annual averages shown in Fig. 5 have a comparable rms over time (σ = 0.3 m/s).

3.3. Significance of the inflows: case study of AR 12007

The flows around active regions are a superposition of the inflows and convective noise (mostly supergranulation). Even after smoothing the flow maps with the Gaussian (σ = 4°), the rms level (~ 8 m/s) is comparable to the magnitude of the inflows. This requires a test to determine whether the inflow signal is real or just caused by noise. We tested this for the active region AR 12007 (highlighted by the blue boxes in Fig. 3). This active region emerges on the far side of the Sun between CR 2147 and CR 2148 and persists over the entire disk passage in CR 2148 without significant evolution. Figure 6 shows the flows around this active region one Carrington rotation after emergence (CR 2148) and the flows in the same region on the quiet Sun one Carrington rotation before the emergence (CR 2147). After emergence, there is an inflow to the active region, while there is no sign of a converging flow in the Carrington map before emergence. However, the magnitude of the flow velocities is similar in both Carrington maps. In order to test whether the inflow signal is real and not just noise, we show in Fig. 7 the daily variation of v_φ and v_θ over the course of one Carrington rotation after averaging them along strips in latitude and longitude with a width of 15° that are centered around the center of gravity of the unsigned magnetic field of the active region (between the dashed lines in Fig. 6). In the quiet Sun (CR 2147), the flow velocities are dominated by short-lived noise (predominantly supergranulation). After the emergence of the active region, the noise level has not changed, but an inflow towards the active region has developed. The temporal variation of the flows is consistent with the noise level in the quiet Sun, indicating that the inflow does not undergo substantial temporal evolution and is indeed a real signal.

3.4. Flows around the averaged active region

As discussed in the previous section, the inflows to individual active regions are affected by supergranulation noise. This means that studying individual active regions is not sufficient for making general statements about the inflows. We addressed this issue by computing the flow field around an averaged active region. We automatically identified active regions with high magnetic flux and then added up the flow field around the selected active regions. This process consists of several steps. We started by convolving Carrington maps of the unsigned magnetic field with a Gaussian (σ = 2°). In these maps, we then identified all regions where the unsigned magnetic field exceeds 30 Gauss as active regions. In total, we detected 971 active regions. Since we are only interested in active regions with a high total magnetic flux, we then integrated the unsigned magnetic field over the area of the individual active regions. In the following, we consider only those active regions, where the total flux is within the top 25% (larger than 5.9 × 10²¹ Mx). This results in 243 active regions. We then average the magnetic field and the flow field for these selected active regions by aligning the center of gravity of the unsigned magnetic field (active regions on the southern hemisphere are flipped in latitude and the sign of the magnetic field is switched). We do not consider the flow velocities within the individual active regions in order to avoid contamination of the flow maps by the magnetic field. Also, when averaging, we do not take into account that the individual active regions have different sizes, tilt angles, or shapes. Figure 8 shows the magnetic field and the flow field for the resulting averaged active region. The averaged active region clearly exhibits an inflow that extends up to 200 Mm from the center of the active region and reaches magnitudes of ~ 15 m/s. This is lower than for individual active regions and is caused by averaging over active regions with different sizes and shapes. The inflow around the averaged active region is far from being azimuthally symmetric. The inflow converges predominantly to the trailing polarity. Also, most of the inflow occurs in the θ-direction, with higher latitudes contributing more than lower latitudes. The shape of this mean inflow depends on the tracking rate and on the subtraction of the background flows. As stated in Sect. 3.1, we subtract the longitudinally and temporally averaged flow velocities from each latitude before averaging the flow velocities.

The magnetic field also exhibits an asymmetry. In the leading polarity, the magnetic field reaches up to 146 Gauss, while it is only 97 Gauss in the trailing polarity. Morphological asymmetries of active regions have been reported before (see, e.g., review by Fan 2014).

4. Discussion

We generated a six-year time series of full-disk maps of horizontal flows at the solar surface generated using LCT. This includes a new approach for removing systematic errors. When systematic errors exhibit a strong time dependence, for example due to variations of the B₀-angle, the commonly used approach of subtracting a mean image does not work very well. This method requires averaging the data over a certain amount of time in order to reduce the noise level to a satisfactory level. During this time, however, the systematic error might change. Our approach of describing the spatial dependency with Zernike polynomials (or any other suitable set of functions) and then filtering the resulting time series of the model parameters in Fourier space reduces this problem. The approach takes into account the time evolution of the systematic errors and is not affected by noise.

Our data set covers a long time series and a large fraction of the solar surface. In addition, LCT allows us to resolve a broad range of spatial scales, from supergranulation up to large-scale flow patterns. This means that these data can be used to study a broad range of science topics. As an example, we use the data to study the inflows to active regions. We can, for the first time, confirm the existence of the inflows with a method that is independent from local helioseismology. In addition, by averaging over a large sample of active regions, we reveal an asymmetry of the inflows. The flows converge predominantly to the trailing polarity of the active region and are predominantly in the θ-direction (with respect to the longitudinally and temporally averaged flow velocities). This effect has not been taken into account by the simple parametrized models used for flux-transport dynamo simulations (De Rosa & Schrijver 2006). Studying the statistics of the inflows over a large number of active regions with different parameters (magnetic flux, latitude, time since emergence, etc.) could also help us test different physical models for the origin of the inflows.

Acknowledgments

We acknowledge support from Deutsche Forschungsgemeinschaft (DFG) through SFB 963/1 “Astrophysical Flow Instabilities and Turbulence” (Project A1). Support was also provided by European Union FP7 projects SPACEINN and SOLARNET. The German Data Center for SDO, funded by the German Aerospace Center (DLR), provided the IT infrastructure for this project. L.G. acknowledges support from the Center for Space Science, NYU Abu Dhabi Institute, UAE.

References

Appendix A: Definition of the Zernike polynomials

Zernike polynomials (Zernike 1934) form a set of orthogonal functions on a disk of radius one. Here, n (radial degree, ≥0) and m (azimuthal degree) are integers with n ≥ | m |. The variables ρ and ϕ are the radial distance and the azimuthal angle, respectively. The Zernike polynomials are defined in the following way:

The radial component is defined as

The Zernike polynomials are only defined if n − m is even.

We choose to normalize the Zernike polynomials as follows:

Appendix B: Time evolution of the Zernike coefficients

The top right part of Fig. B.1 shows the time evolution of the Zernike coefficients for the Zernike polynomial with n = 3 and m = − 1 for v_m. The time series is dominated by a yearly oscillation due to variations of the B₀-angle. In addition, the orbit of SDO causes a series of peaks in the power spectrum at a period of 24 h and its harmonics. Since the orbit of SDO around the Earth is modulated by the orbit of the Earth around the Sun, each of these harmonics is a multiplet consisting of several individual peaks. These peaks are separated by a frequency of ~ 2.8 × 10^-3 1/day, which corresponds to a period of one year. Apart from that, the power is almost constant, suggesting that it originates from real solar flows or white noise. The time series for the individual Zernike coefficients do not all look the same. Not all Zernike coefficients exhibit a yearly period and which harmonics of the 24 h period can be seen is also not always the same.

Here, we account for the influence of the orbit using a purely empirical approach. For each Zernike coefficient we generate a time series consisting of the signal measured at the frequencies that we consider to be affected by the orbit. First, we remove the extreme outliers from the time series. In some cases, the HMI data are corrupted (mostly due to eclipses), which affects the flow maps. If the coefficient for n = 0 and m = 0 (the mean of the flow map) for v_m or for v_n deviates by more than 5σ from the mean, we assume this time step to be corrupted and replace all Zernike coefficients for that time step by their respective mean values. We do not include these data in our analysis. We perform the outlier correction two times consecutively. In the next step, we compute the Fourier transform of the time series and filter for only those frequencies that are affected by the orbit of SDO plus the mean of the time series. These are the daily and yearly variations plus the corresponding harmonics. In the case of the yearly variations, the main frequency (2.8 × 10^-3 1/day) plus the first 16 harmonics (up to a frequency of 4.4 × 10^-3 1/day) exhibit significant power and are taken into account. In the case of the daily variations, we use the main frequency (1/day) and all 23 harmonics that can be resolved with the cadence of our data (30 min, frequency of the last harmonic: 24/day). Since the daily variations and their harmonics are not isolated peaks but multiplets with the individual peaks being separated by a frequency of 1/year, we extract the signal in a range of frequencies around the central frequency (± 7.5 × 10^-3 1/day) in order to cover the entire multiplet. Taking account only of the signal at the frequencies selected in this manner, we finally generate a time series for each Zernike polynomial (see Fig. B.1 for an example).

This approach works for most Zernike polynomials, but not for all of them. The time series for the polynomial with n = 0 and m = 0 (the mean of the velocity) both for v_m and for v_n are affected by the orbit on all temporal scales. In these two cases, we do not select specific frequencies but use the full frequency spectrum for the further analysis. Issues also arise for the polynomials with n = 1 and m = − 1 for v_m, and n = 1 and m = 1 for v_n. The power spectra of the time series for these polynomials exhibit a significantly higher background power than the other Zernike polynomials. The time series for these two Zernike coefficients are highly correlated, suggesting that this background noise originates from rotating the images. It is known that the p-angle of the HMI data exhibits oscillations whose origin is not fully understood (Löptien et al. 2016a). We account for this by performing a principal component analysis of the two time series. We assume that the correlated component originates from the p-angle oscillations, while the uncorrelated part represents real solar flows.

Appendix C: Converting velocities from CCD to spherical coordinates

The velocities derived using LCT are defined on a 2D Cartesian coordinate system that is aligned with the plane of the CCD. The m-axis points to the right, and the n-axis points upwards. Tracking and remapping these velocities to heliographic coordinates results in v_m(λ,φ) and v_n(λ,φ), with λ and φ being the heliographic latitude and longitude, respectively. The velocities are now given as a function of the heliographic coordinates, but their coordinate vectors, ê_m and ê_n, are still defined on the 2D Cartesian coordinate system of the CCD. In addition, the velocities are still given in units of pixel/s and need to be converted to m/s. In order to transform the coordinate vectors from CCD to spherical coordinates, we need to make the assumption that we only measure velocities on the surface of the Sun (i.e., r = R_⊙) and that they are purely horizontal (i.e., v_r = 0).

In the first step, we take into account the finite distance of SDO to the Sun and convert the velocities from pixel/s to m/s. We define a new coordinate system centered on the Sun given by the basis vectors , , and , with the x′-axis pointing towards the observer, the -axis pointing to the right (-direction), and the z′-axis pointing upwards (-direction). The coordinates in the system are related to the CCD coordinates by Here, s is the plate scale of HMI (in rad/pixel) and D is the distance from SDO to the center of the Sun. In order to convert the velocities, we take the partial derivative of these equations with respect to time and solve for and . Since we assume our velocities to be purely horizontal, we can express as a function of and . This leads to the following equations: \newpageHere, the elements of the matrix J originate from expressing as a function of and (in order to enforce v_r = 0) and are given by In the next step, we rotate the coordinate system given by the basis vectors , , and around the y′-axis by an angle − B₀. In the resulting coordinate system with unit vectors , ŷ, and ẑ, the z-axis is aligned with the rotation axis of the Sun. Then, we transform from the Cartesian coordinate system, , ŷ, and ẑ, to spherical coordinates , (direction of increasing heliocentric longitude), and (direction of increasing co-latitude). Again, we use the assumption that that the velocities on the Sun are purely horizontal (i.e., v_r = 0). This results in the following equations for v_φ and v_θ: The elements of the matrix a and the factor A are given by Inserting Eqs. (C.3) and (C.4) in (C.8) and (C.9) gives the final result The coefficients are given by

All Figures