© N. Gottschling et al. 2021

1. Introduction

Active regions (hereafter ARs) are patches of magnetic field at the surface of the Sun. They are host to sunspots that are visible in white-light images (van Driel-Gesztelyi & Green 2015) and the site of a large variety of phenomena, such as flares and coronal mass ejections. Theories of active region formation include a buoyant rise of magnetic flux tubes from the bottom of the convection zone (see the review by Fan 2009), as well as a formation in the bulk of the convection zone (Nelson et al. 2013) or in the near-surface layers (Brandenburg 2005).

At the surface, active regions are advected on large spatial scales by differential rotation (Snodgrass 1983) and the meridional flow (Duvall 1979). On smaller scales, the magnetic field is diffused by the turbulent convective motions of granulation and supergranulation (Leighton 1964). Methods for inferring the large-scale flows include tracking the motions of small-scale convective or magnetic features, which are displaced by the larger-scale motions, as well as a variety of tools for the analysis of solar oscillations, which are known as local helioseismology (see the review by Gizon & Birch 2005).

In addition, Gizon et al. (2001) found that evolved active regions are surrounded by horizontal converging flows. The authors used time-distance helioseismology (Duvall et al. 1993) to study flows near the surface of the Sun. These surface inflows typically extend to about 10° from the active region and have horizontal velocities of up to 50 m s⁻¹ (Gizon et al. 2001; Haber et al. 2004; Zhao & Kosovichev 2004; Hindman et al. 2009; Löptien et al. 2018; Braun 2019). It is an interesting question to determine what drives these inflows. Enhanced emission of radiation in the magnetic field structures might lead to a horizontal temperature gradient that could drive the flows (Spruit 2003).

The inflows toward active regions are important as they affect the advection and flux cancellation of ARs, and may counteract diffusion (De Rosa & Schrijver 2006). De Rosa & Schrijver (2006) and Martin-Belda & Cameron (2017a,b) included models of the inflows into surface flux transport simulations, showing that they can provide a feedback mechanism that saturates the amplitude of the global magnetic field. This could explain observed variations in solar cycle strengths.

Several studies have applied local helioseismic techniques to investigate different aspects of the properties of these flows. Haber et al. (2004) studied their depth dependence, identifying a transition to outflows in deeper layers for some ARs. Komm et al. (2011, 2012) studied vertical flows and vorticity, finding that emerging flux is associated with weak upflows of less than 1 m s⁻¹ and increasing strength of vorticity, and that decaying flux is associated with downflows and decreasing strength of vorticity.

Recently, Braun (2019), using helioseismic holography (Lindsey & Braun 2000), investigated a sample of 336 ARs and derived flow patterns for subgroups of increasing magnetic flux. This work showed that the amplitude of the flows increases with magnetic flux, and reported a retrograde flow component mainly at the poleward side and, less distinct, the equatorward side of the AR.

Löptien et al. (2017) used local correlation tracking (LCT, November & Simon 1988) of the convective granulation pattern on the solar surface to measure inflows around an average high-flux AR. The regions in their sample have fluxes higher than 5.9 × 10²¹ Mx. Their work confirms inflows of about 20–30 m s⁻¹ that extend about 10° from the AR.

Birch et al. (2019) found a compact (less than 2°) converging flow of about 40 m s⁻¹ in the day before active region emergence, using both LCT and helioseismic holography as independent methods for flow inferences on the surface. Comparisons of flows derived with LCT and local helioseismology showed that the two methods agree well (Švanda et al. 2013; Birch et al. 2016).

In this study, we aim to measure the evolution and structure of the inflows around active regions at the solar surface during the early stages of the active region evolution, that is, before, throughout, and in the days after emergence. We also investigate potential differences in the inflows related to the total unsigned flux and the latitude of the active regions.

We performed a statistical study by averaging over a large sample of ARs. For this, we made use of the Solar Dynamics Observatory (SDO) helioseismic emerging active regions (SDO/HEAR) survey (Schunker et al. 2016). By using this sample rather than performing a case study, we decrease background noise as well as peculiarities of individual emerging active regions (EARs), which can have a dominating effect on the flow structure.

The paper is structured as follows: In Sect. 2 we describe the sample of emerging active regions as well as the flow measurements. Section 3 illustrates the data processing that we applied to derive the evolution of flows around ensemble averages of ARs. Section 4 presents our results, followed by a discussion (Sect. 5).

2. Data

2.1. Sample of emerging active regions

Schunker et al. (2016) described the HEAR survey in detail. We summarize the main aspects below.

The HEAR survey comprises 182 emerging active regions during the time from 2010 to 2014. They were selected from the National Oceanic and Atmospheric Administration (NOAA) solar region reports on the basis of several criteria, requiring that the AR (a) reaches a total sunspot area of at least 10 micro hemispheres (μH; 1μH ≈3 Mm²), (b) is first recorded by NOAA within 50° of the central meridian, and (c) emerges into the quiet Sun, without preexisting large-scale flux within a radius of 18°. These criteria were chosen to minimize contamination by preexisting magnetic field and to reduce projection effects. The emergence time of each active region is defined as the time at which the region has reached 10% of the maximum unsigned flux within the first 36 h after the active region is first recorded by NOAA. The survey also provides the Carrington longitude and latitude, as well as the start and end times of the observations for each region.

In addition to the emerging active regions, the survey includes one quiet-Sun control region for each EAR. These control regions match their corresponding EARs in latitude and distance from the central meridian at a mock-emergence time, at which there are no numbered SHARP regions within a radius of 18° of the control region. These control regions are important as a comparison against which any measurement on the EARs can be tested.

2.2. Flows inferred from local correlation tracking

We used an updated version of the flow maps generated by Löptien et al. (2017), who described their processing procedure in detail. We summarize the important aspects below.

Löptien et al. (2017) applied LCT on full-disk continuum intensity images from SDO/Helioseismic and Magnetic Imager (HMI) (Schou et al. 2012). Specifically, they ran the FLCT code (Welsch et al. 2004; Fisher & Welsch 2008) on hmi.ic_45s series data, starting on 24 April 2010 and ending on 27 April 2016. In order to reduce computation time, only one flow measurement every 30 min was calculated by cross-correlating pairs of images that are 45 s apart. Each measurement consists of a map of velocities v_x in x-direction and a map of velocities v_y in y-direction (with x and y in CCD coordinates).

Local correlation tracking is known to underestimate flow velocities (see, e.g., Löptien et al. 2016a; Verma et al. 2013; Švanda et al. 2007 and references therein). This was addressed by generating calibration data with which the flow velocities were then corrected (see Sect. 2.1 of Löptien et al. 2017).

The noise in LCT calculations increases to the limb due to projection effects. Therefore the flow maps were cropped at 60° from disk center. The maps show systematic effects, such as the orbital motion of HMI and the shrinking-Sun effect (Lisle et al. 2004; Löptien et al. 2016b), a systematic converging flow toward disk center. This was addressed by representing the background signal of the flow maps as a sum of low-order Zernike polynomials (up to radial degree n = 7). The time series of the Zernike polynomial coefficients were then transformed into Fourier space, where frequencies corresponding to known periods of systematic effects were isolated. These are the 24-hour period of the satellite orbit and its 23 observed harmonics, and the 365-day period of the Earth orbit and its first 16 harmonics (down to a period of 22.8 days), which show significant power. The filtered, backtransformed background estimate was then subtracted from the observations.

An exception was made for the time series of Zernike polynomials for the v_x maps and for the v_y maps, which are correlated to each other and were therefore not treated in the way outlined above. For these, the first element of a principle component analysis of the two components was subtracted (as described in Appendix B of Löptien et al. 2017). However, this retained a systematic variation for these terms on the order of several m s⁻¹ with a period of 24 h and its harmonics, which corresponds to large-scale gradients in the flow maps. To correct for this, we repeated the data processing of Löptien et al. (2017), subtracting for the Zernike polynomials for v_x and for v_y the same harmonics as for the other components. In this rerun, the data coverage was increased to 6 March 2020.

As a last step, we remapped the velocity maps from the CCD coordinate system to Plate Carree projection (Carrington longitude, latitude), with a total size of 3600 × 1200 grid points and a grid spacing of 0.1°. In this projection, averaging over regions at different latitudes and longitudes can be done in a straightforward manner. In the same step, the velocities were transformed from v_x, v_y into the longitudinal and latitudinal velocities v_lon, v_lat. We then applied a Gaussian filter with a width of 0.4° and subsampled to the grid spacing of 0.4° used by Löptien et al. (2017). This processing step differs from Löptien et al. (2017) so that aliasing is removed from the data.

3. Data reduction

3.1. Processing and validation of the flow data

In the flow maps, individual pixels occasionally show unrealistically high values. This affects about 10 pixels per map, mostly close to the limb. We speculate that these arise from the decreased contrast at the limb in combination with the evolution of the granules. To mitigate their effect, we replace outlier pixels that are outside five standard deviations of the mean of the velocity distribution by the mean value of pixels surrounding them in a box of 5 pixel × 5 pixel.

In addition to the systematic effects described above (the orbital motion of SDO and the shrinking-Sun effect), the LCT data exhibit another long-term, large-scale modulation: At the beginning of the time series in 2010, the longitudinal flow component v_lon has a bias toward prograde (retrograde) velocities in the eastern (western) half of the disk, and the latitudinal flow component v_lat has a bias toward southward (northward) velocities in the northern (southern) half of the disk, mostly at latitudes above 30°. Over the years, this modulation changes, such that the hemispheric biases of the two flow components reverse sign in about 2014. The amplitude of this modulation is about 10 m s⁻¹ in a half-year average. The reason for this modulation is not fully understood. Its timescale does not correspond to any known periodic effects, such as the B angle or the orbital motion of HMI.

To correct for this large-scale modulation, we translated all velocity maps from Carrington longitude into Stonyhurst longitude, created running time averages over half a year of data each, centered ten days apart, and subtracted these from every individual map from five days before until five days after the central time of the average. The corrected maps were then mapped back to Carrington longitude. We chose an averaging time of half a year because a shorter timescale might subtract real flows. The spacing of ten days was chosen to avoid drastic changes between two blocks of subtraction. The average maps were smoothed with a Gaussian of σ = 0.8° to reduce small-scale effects.

The above processing was necessary to address the various systematic errors, namely the orbital motion of HMI, the shrinking-Sun effect, and the long-term modulation. The velocities that we measure throughout our analysis are therefore relative to velocities v_lon = 0, v_lat = 0, which are the half-year time averages at each location in Stonyhurst longitude and latitude.

As a test of this processing, we wished to be sure that the Zernike subtraction outlined by Löptien et al. (2017) did not subract the types of flows we wished to study together with the systematic effects. To evaluate by how much this procedure alters the inflows, we created one Carrington rotation of synthetic flow data with inflow features and applied the same filtering on them as on the actual data. Appendix A describes this test. We find that the extent and amplitude of the synthetic flows is well preserved, with changes of less than 5 m s⁻¹ for model inflows with an amplitude of 50 m s⁻¹. Small deviations are introduced to the data, with amplitudes of less than 5 m s⁻¹.

As another test of the LCT flows used in this study, we carried out a control run against an independent data set. In this test, we projected the LCT flows onto the line-of-sight component and compared them to flow measurements obtained from direct Doppler images from SDO/HMI. Appendix B describes this test. The Dopplergrams only provide one velocity component and are therefore not suitable as the main flow measurement of this study. We find that the two methods agree well in general, with correlations exceeding 0.8, except for in close proximity (2°) to sunspots. This might be due to the different depths at which the two measurements are taken, and also to different sensitivities in the presence of strong magnetic field.

To use the LCT flow data on the sample of emerging active regions and control regions, for each region we extracted data cubes of the two flow components v_lon, v_lat, which are centered at the position of the region in Carrington longitude and latitude given in the HEAR survey, and span ±30° in both longitude and latitude. The cubes start (end) at the full or half hour before (after) the first (last) observation time given in the HEAR survey. This is because the HEAR survey defined the emergence time from the HARP vector magnetogram observations that have a cadence of 12 min, whereas the LCT data have a cadence of 30 min, with data at full and half hours. For the same reason, we define the time of emergence as the first full or half hour after the time of emergence t₀ in the HEAR survey. We chose this to ensure that the emergence had consistently occurred for all regions.

The typical observation period of an individual region is about nine days (450 frames). The maximum difference between the start (end) of the observations and the time of emergence is seven days (343 frames) before (after) emergence (leaving about two days of observations after (before) emergence). The averaging over ensembles of active regions was done relative to the time of emergence (see Sect. 4). The total duration of observations is therefore 14 days, with a varying number of active regions at each time relative to the time of emergence. The parts of each frame that are off the visible disk were not included in any averaging.

To further analyze the characteristics of the EARs, we generated magnetogram and continuum intensity data cubes for each region, matching the setup of the flow data described above. We used the same observation start and end times, with the same 30-minute cadence as for the flow maps (described above) on the continuum intensity series hmi.ic_45s to create the intensity data cubes. The data were remapped to Plate Carree projection, with a grid spacing of 0.1°. This is four times higher than that of the flow maps and was chosen for both continuum and magnetogram data to preserve small-scale information. We corrected for limb darkening and normalized the continuum data cubes by convolving each frame with a Gaussian of σ = 5° to retrieve the large-scale variation in the frame (the limb-darkening profile) and then divided the frame by this background. For the magnetograms, we extracted the full-vector magnetic field from hmi.b_720s. From the vector field components, we calculated the radial field B_z using the transformation given by Gary & Hagyard (1990). For pixels in which the field azimuth was not calculated with the minimum-energy algorithm, we used the random field disambiguation (Hoeksema et al. 2014). We then averaged the five frames of the 720 s radial field B_z data together within one hour.

The data of each EAR thus consist of three-dimensional data cubes (time, longitude, and latitude) for the flow components in longitudinal and latitudinal direction, the radial magnetic field, and the continuum intensity. Each cube has frame sizes of 60° ×60° in longitude and latitude, centered at the Carrington longitude and latitude of the AR given in the HEAR survey, and covers about 14 days, centered on the time of emergence, with data coverage of the visible disk for about nine days. The cubes of the velocity maps and the continuum intensity have a cadence of 30 min. The cubes of the radial magnetic field have a cadence of one hour.

3.2. Measuring the location of the magnetic polarities

Over time, the regions move relative to the position given in the HEAR survey. To perform a sensible ensemble average of the active regions relative to their center or their leading or trailing polarity, we measured the positions of the leading and trailing polarities in longitude and latitude in each one-hour radial magnetic field average. To reduce sensitivity to small-scale structure within the AR, we smoothed the magnetograms with a Gaussian of σ = 0.4° for this measurement.

For each of the 182 AR we estimate the position of the polarities at one day after emergence by eye as an initial guess. At this time, all regions show a clear bipolar structure. The procedure for each time frame in each of the 182 data cubes is described below.

We estimated the fluctuation in the radial magnetic field as the standard deviation in a box of 12° ×12° around the position in Carrington longitude and latitude from the HEAR survey. The magnetic field of the AR is part of the estimate. The size of the box was chosen to exclude other AR from the field of view; they would bias the estimate.

We then identified regions as neighboring positive (negative) valued grid points (for brevity, called pixels in this context) above (below) a threshold value, taken as 1.5 times the standard deviation. We used this adaptive threshold because the levels of magnetic activity in the vicinity of the individual ARs differ. The factor of 1.5 was applied to further reduce detection noise and is a result of trial and error. The negative and positive polarities were considered separately. Depending on the AR, the positive (negative) field corresponds to the leading (trailing) polarity, or vice versa.

Regions that are only one pixel in size were excluded. These were considered to be too small for a reliable polarity position estimate. We also excluded regions whose flux-weighted center is farther away from the initial (or first-guess, see below) position than ±2° to avoid selecting unrelated surrounding magnetic field. This is important especially during emergence, when the field of the AR is still weak.

Of the remaining regions, we identified the region that has the highest absolute sum in B_z. We defined the position of the active region polarity as the longitude and latitude of its flux-weighted center.

This method was iterated forward and backward in time, starting from the frame one day after emergence, when we identified a first-guess position by eye. For each frame, the position from the previous frame was used as first guess.

With this procedure, we determined the positions in longitude and latitude of each polarity for each EAR at each time step. Before emergence, there is only very little magnetic flux of the active region present by definition, which makes a position measurement unreliable, with the risk of identifying surrounding small-scale field. For consistency, we therefore set the positions in all frames before emergence to the values at emergence. We point out that other choices could be made, for example, extrapolating the post-emergence motion to the times before emergence. This would require additional assumptions on the proportionality of the motion over a certain (flux dependent) time interval, however, and was therefore not attempted here.

In some cases, the inferred positions in a single frame differ largely from those before and after, for example, when a connection between two regions is established by pixels that have magnetic field strengths just above the threshold in the particular frame, and below in the previous and the following frame. These artifacts will be present for any choice of threshold. To mitigate this, we smoothed the time series of positions with a Gaussian of width σ = 3 h. We defined the position of the center as half the distance between the leading and trailing polarity positions. The total motion of an AR over the course of the observations is typically a few degrees. The largest displacements are about 6° for ARs at high latitudes, at which the difference between Carrington rotation rate and solar rotation rate at the surface is largest.

Other methods for inferring the position of AR polarities are possible. To estimate the difference that this might make for our analysis, we applied the method of Schunker et al. (2019). Their algorithm iteratively calculates the field strengths in the surroundings of all grid points and selects the patch as the polarity that lies closest to the patch in the previous frame. Because our method instead selects the strongest patch, this can lead to differences between the two methods because the dominant patch of magnetic field of one polarity in some ARs is later replaced by another patch. Appendix C compares the two methods. We applied the positions of both methods to ensemble averages of the flow data (see the following sections) and found a standard deviation of the difference of the methods of below 2 m s⁻¹ close to emergence, and no systematic difference with respect to the inflows we wished to investigate.

3.3. Weight maps for ensemble averages

When ensemble averaging is performed, we have to take into account that at a given time relative to the time of emergence, the active regions have different disk positions and therefore different background noise levels because the noise increases toward the limb. We accounted for this with weighted ensemble averages. For this, we created data cubes of weight maps for all EARs and control regions as described below.

First, we mapped the full-disk velocity field v_lon, v_lat from the Plate Carree grid with 900 × 300 grid points and a grid spacing of 0.4° (see Sect. 2.2), from Carrington coordinates into heliocentric Cartesian coordinates (Thompson 2006), centered on the Carrington longitude and latitude at disk center. Transforming into this system accounts for the B angle and the P angle. We mapped with a linear interpolation onto a fixed grid with 300 × 300 grid points and a grid spacing of 0.006 solar radii (R_⊙) (corresponding to about 4.2 Mm at disk center). The frame size was 1.8 R_⊙, which is slightly larger than the LCT data coverage (the 60° from disk center of the LCT data corresponds to ).

We did this for the whole ten years of LCT data. We calculated the average flows ⟨v_lon⟩ and ⟨v_lat⟩, and a map of the standard deviation of each flow component (see Fig. 1) on the fixed grid in heliocentric Cartesian coordinates. The velocities in the average maps are about 10 m s⁻¹, which is far lower than the velocities in an individual map (about 1000 m s⁻¹). The reason is that the velocities are equally distributed around zero, with the dominant noise from granulation and supergranulation.

Fig. 1. Temporal averages (left panels) and standard deviations (right panels) over the full-disk flow maps from 24 April 2010 to 6 March 2020. The top and bottom rows show the longitudinal and latitudinal flow component, respectively.

For the longitudinal component, the standard deviation increases mainly toward the east and west limbs. For the latitudinal component, the northern and southern polar regions show the largest standard deviation.

We generated full-disk weight maps by mapping the standard deviation σ of v_lon and v_lat back onto the Plate Carree projection in Carrington coordinates. From these, we created data cubes of the weight maps as for each of the 182 EARs and control regions in the same way as for the data cubes of the velocity maps (see Sect. 3).

3.4. Exclusion of magnetic pixels

Local correlation tracking of solar granules is known to be unreliable within regions of strong magnetic field, where the granulation is attenuated. This is the case for sunspots or plage regions, which appear darker or brighter than the quiet Sun in continuum images. The motion that LCT tracks in these regions is therefore a combination of the motion of the granules and the motion of these features. In Appendix B, we find that in regions within 2° of sunspots, the LCT flows are weaker than the direct SDO/HMI Doppler velocity measurements. This difference can be attributed to physical causes as well, however, that is, the difference in the height of the two measurements.

To ensure that our measurements are not susceptible to this unreliability, we created additional data cubes of normalized continuum intensity for each AR by binning to the same grid spacing as the flow maps (0.4°). We identified pixels with values lower than 0.95 (pores, sunspots) and higher than 1.05 (plage), where the mean value of the quiet Sun is 1. In the flow maps, we excluded these pixels if they were connected to at least one other to avoid excluding random high-valued pixels. The numbers are the result of trial and error and represent a compromise between excluding pixels in the quiet Sun that are above or below the thresholds and not excluding pixels that are systematically lower or higher, for example, near the edge of a sunspot.

3.5. Ensemble averages

To carry out ensemble averages, we first investigated the dependence of the noise in the flow maps on spatial smoothing, temporal averaging, and the number of ARs in the ensemble average to identify parameters that yield appropriate noise levels for our analysis. Previous studies have reported active region inflows with velocities of about several tens of m s⁻¹. We therefore required the noise to be well below 10 m s⁻¹ while still being able to investigate flows on small spatial scales and without loosing much temporal resolution.

For this, we used the sample of control regions of the HEAR survey, where no systematic flows are present. The data cubes of these were prepared in the same way as those of the EARs, with the center of the 60° ×60° maps on the coordinate given in the HEAR survey throughout the observations. For this test, we cropped the frames to the inner 40° ×40°, which yielded a more representative value because the errors increase toward the limb and fewer control regions contribute to the average. Like the active regions, the quiet-Sun regions have different disk positions at a given time relative to the assigned ‘time of emergence’. For a representative estimate, we started at one day before the time of emergence, where the frames from the control regions are close to disk center on average, and computed the average of N subsequent frames in the range of N = 1 (no time average, 30 min resolution) to N = 96 (2 days average). They were additionally smoothed by Gaussians of σ = 0.4, 0.8, 1.2, 1.6, 2.0, and 4.0°. We calculated the median of the standard error over all individual control regions. Second, we carried out the same calculation at a fixed Gaussian smoothing of 0.8° while varying the number of control regions in the average. Figure 2 shows the results for the latitudinal flow component. The results for the longitudinal component are very similar.

Fig. 2. Median of the standard error against averaging time. The lines in the top panel show the cases of no smoothing (yellow) and of Gaussian smoothing of 0.4, 0.8, 1.2, 1.6, 2.0, and 4.0° (green to purple). All lines are from an ensemble average over all 182 ARs. The dotted line shows the least-squares fit to the case of Gaussian smoothing of σ = 1.6°. Bottom: Same as in the top panel for averages over 5, 10, 20, 45, 60, 91, and 182 control regions (yellow to purple). All lines are with a Gaussian smoothing of σ = 0.8°. Both panels show the case of the latitudinal flow component. The scale of the y-axis is different in the two panels.

The top panel of Fig. 2 shows that an ensemble average over all 182 ARs requires a spatial smoothing of 0.8° to obtain a noise level of about 5 m s⁻¹. A time average of 6 or 12 h is sensible because of the timescale on which the flows change. The top panel also shows a least-squares fit of to the case of Gaussian smoothing of σ = 1.6°. The noise decreases more slowly than the fit because subsequent maps are correlated due to the lifetime of supergranules. This holds for all cases. The bottom panel of Fig. 2 shows that for an ensemble average of fewer than 45 regions, temporal averaging of more than 2 days is necessary at this level of spatial smoothing. For an ensemble average of 45 regions, a time average of 12 h is acceptable.

We used these values for spatial smoothing and temporal averaging for the ensemble averages in our analysis. The temporal averaging was aligned relative to the time of emergence, such that the frame of emergence is the first frame in the first time step after emergence. For the ensemble averaging, Hale’s law and Joy’s law are taken care of by reversing the latitudinal direction of all data cubes (magnetic field, continuum intensity, and the v_lon and v_lat flow components) as well as the sign of the magnetogram and the v_lat flow component for regions in the southern hemisphere.

In addition to the flows around active regions, the flow maps exhibit signatures from the moat flow because sunspots are present in the sample (Sheeley 1972). To better understand their effect on the flow maps, we performed ensemble averages centered on the leading polarity (where we expect the strongest moat flow), in subsamples of the active regions, using a classification scheme to identify whether a sunspot with a clear penumbra is present in an AR at that time. Appendix D describes this test. We find that a moat flow is present in our flow data for regions with a clear sunspot, with typical velocities of about 150 m s⁻¹. The moat flow is not symmetric in this early phase of the active region as the sunspots still grow, in agreement with Vargas Domínguez et al. (2008).

4. Results

4.1. Flows as a function of time and magnetic flux

We investigated the evolution of the flows around active regions of different maximum values of magnetic flux. For each of the 182 ARs, we retrieved the maximum total unsigned flux from the hmi.sharp_720s series keyword USFLUX in the period covered by the HEAR survey. USFLUX records Σ|B_z|dA, with B_z the radial magnetic field component and A the area of the SHARP region (Bobra et al. 2014). From this, we divided the sample into four subsamples, sorted by the maximum unsigned magnetix flux, with two subsamples containing 45 and two subsamples containing 46 EARs.

The left panel of Fig. 3 shows the distribution of maximum USFLUX values of all ARs, together with dotted red lines that indicate the boundaries between the four subsamples. The subsamples have mean total unsigned magnetic flux values of 1.5 × 10²¹ Mx, 3.1 × 10²¹ Mx, 7.0 × 10²¹ Mx, and 1.5 × 10²² Mx, respectively.

Fig. 3. Distributions of the sample of EARs with respect to the unsigned magnetic flux (USFLUX, left panel) and the unsigned latitude |λ| (right panel). The dotted red lines in the two panels separate the distributions into four subsamples, each containing 45 or 46 active regions, sorted by USFLUX and by |λ|, respectively.

Figure 4 shows several 12-hour time steps for all four subsamples. The first row shows the time step of the flows at −18 h (between 24 and 12 h before emergence). At this time, all subsamples exhibit a converging flow toward the center of the AR. The velocities are about 20 to 30 m s⁻¹. The noise level is about 10 m s⁻¹, as discussed in Sect. 3.5. This converging flow is consistent with that reported by Birch et al. (2019). For the two higher-flux subsamples, the converging flow is spatially more extended than for the lower-flux ones. Shortly after emergence, the lowest-flux subsample still shows small-scale inflows, whereas they have ceased for the other subsamples. The three subsamples with higher flux show a prograde flow with velocities of about 40 m s⁻¹ at the position of the leading polarity, similar to the findings by Birch et al. (2019). This feature is not present in the lowest-flux subsample. About two days after emergence (third row in Fig. 4), inflows toward the second-lowest flux subsample set in. For the third subsample, they are clearly visible after about three days, while the highest-flux subsample still shows mostly a strong moat flow at this time. For the highest-flux subsample, the inflows only set in at about four days after emergence. From this, we conclude that the inflows take longer to form for large active regions.

Fig. 4. Evolution (from top to bottom) of the averaged active region magnetic field and flows for the four subsample averages of total unsigned flux (with flux increasing from left to right, cf. the left panel of Fig. 3). Each time step is an average over 12 h, centered on the labeled time. The axis labels are relative to the center of the AR, as defined in Sect. 3.2. For each panel, the total number of active regions that contribute to the average is given at the top. Red (blue) indicates positive (negative) radial magnetic field, saturated at ±150 Gauss. The green circles highlight the time after emergence at which the inflows are first visible in the time steps shown here. The arrows indicate the flows. A reference arrow representing a velocity of 50 m s−1 is given in the upper right corner. The black lines mark the averaging range that was used for the model fit (cf. text and Fig. 6). A moat flow is apparent in the rightmost column at 3 d and 6 h, for instance.

After six days, the EARs of the lowest-flux subsample have already mostly decayed. For these small AR, inflows are seen throughout their lifetime, and they decay along with the magnetic field. Over its lifetime, the lowest-flux subsample never exhibits diverging flows, neither moat flows (due to the lack of sunspots with clear penumbra), nor a prograde flow at the leading polarity of the type observed by Birch et al. (2019).

4.2. Flows as a function of time and latitude

We examined whether the flows associated with the EARs depend on the latitude at which the EARs emerges. To do this, we divided the sample of 182 ARs into four subsamples of 45 or 46 EARs each, sorted by their unsigned latitude |λ| (as provided by the HEAR survey). Each subsample in |λ| therefore contains regions with a variety of maximum total unsigned flux values. The right panel of Fig. 3 shows the distribution of unsigned latitudes of all EARs together with lines indicating the boundaries between the four subsamples.

Figure 5 shows two 12-hour time steps for all four subsamples, one at six hours after emergence, and the other at 4 d and 18 h after emergence. The converging flow toward the center of the active region as well as the prograde flow at the position of the leading polarity around the time of emergence are both present in all four subsamples. Joy’s law is clearly visible after emergence, consistent with Schunker et al. (2020), who showed that active regions on average emerge east-west aligned, and that Joy’s law becomes evident two days after emergence. The tilt angle increases toward the higher-latitude subsamples and in time. There appears to be no striking systematic change in the inflows as a function of latitude.

Fig. 5. Evolution (from top to bottom) of the averaged active region magnetic field and flows for the four subsample averages of unsigned latitude |λ| (with |λ| increasing from left to right, cf. right panel of Fig. 3). The labels, color scales, and arrow scales are the same as in Fig. 4. The diverging flow in the bottom row of the third subsample is due to a strong moat flow in AR 11158, which is the largest AR in this subsample at this time.

4.3. Flows averaged over the full sample as a function of time

While the sections above outline the evolution of the flows in subsamples of active regions with respect to flux and latitude, we wish to set this into the context of the average EAR over the whole sample. Appendix E provides additional plots of the average flows around the time of emergence as well as averages in longitude covering the leading and trailing polarity separately. Overall, the full ensemble average again exhibits converging flows starting one day prior to emergence. The inflows have a maximum extent of about 7° from the center of the AR, in the time between emergence and seven days after emergence, and amplitudes of about 50 m s⁻¹ at a noise level of about 5 m s⁻¹.

4.4. Quantitative model of the inflows as a function of time and magnetic flux

We wish to describe the inflows around active regions quantitatively and to evaluate the differences found in the flux-binned subsamples (Sect. 4.1). To do this, we fit a model to the flows. The main contribution of the inflows is along the latitudinal axis. We therefore averaged the latitudinal flows over the central 5.6° in longitude for each time step and each subsample. This range was chosen to cover a large part of the AR but excludes most of the moat flow signature that is present in the high-flux subsamples. We did not attempt to fit a separate model to the moat flow because the sunspots and with them the moat flow in the individual active regions form at different times relative to the time of emergence, and in some cases break up before the end of the observations (cf. the lower panels for the second-highest flux subsample in Fig. 4).

The model we applied consists of the sum of up to two 1D Gaussians with positive or negative amplitude A, position μ,  and width σ. As an estimate of the background noise, we computed the standard deviation over the same 5.6° in longitude for subsamples of the control regions corresponding to the AR subsamples. If the amplitude of a fitted Gaussian was smaller than 1.3 times the noise, it was discarded. The converging flows before emergence tend to be less spatially extended than the inflows after emergence. We therefore confined the region of the fits before emergence to about 2° in latitude to avoid fitting large supergranules.

To estimate realistic errors on the fit parameters, we added the model fit to the control region frames and reran the model fitting on these frames. The control region flow maps have the same properties (spatial scale, amplitude) as the active region flow maps and provide a background noise that is uncorrelated with the background in the active region maps. Averaging over all realizations of this procedure and taking the standard deviation of each parameter from these realizations yields an error estimate for the model parameters. Of the 30 time steps of the control region subsamples, the first and last five were not used because of their proximity to the limb, which would artificially increase the noise. Thus, for each model fit, there are 4 × 20 = 80 realizations that contribute to the error estimate. Figure 6 shows an example of the model fit. The calculated errors (dashed blue lines) are on the same order as the noise in the data (gray shaded region).

Fig. 6. Data and model fit of a time step of the second-highest flux subsample (⟨Φ⟩=7.0 × 1021 Mx). The black line shows the latitudinal flow, averaged over longitudes in the range indicated in Fig. 4. The solid blue line shows the model fit of two Gaussians with opposite polarity (see text for details). The dashed blue lines show the 1σ error estimates of the fit. The red markings indicate the fit parameters of the amplitude A (in m s−1), the peak position μ (relative to the center of the active region along latitude, Δ Lat), and the standard deviation σ (in degree) of the two Gaussians. The gray shaded region indicates the rms of the flow data. The fit parameters of the two Gaussians slightly deviate from the fitted curve due to superposition; this is small compared to the error, however.

Figure 7 shows the fit parameters A, μ and σ of v_lat averaged over longitudes for the four flux-binned subsamples in the time range between one day before and six days after emergence (the last time step is left out because the signal-to-noise ratio is low). For the three higher-flux subsamples, the converging flows of the pre-emergence phase vanish shortly after emergence, such that there is no successful fit. After some days, the inflows set in. We defined the onset time of the inflows as the time when the model fits start to be continuously successful. The amplitudes increase from this time on, starting from less than 20 m s⁻¹ and reaching maximum velocities of 50–60 m s⁻¹ in all cases. In the case of the lowest-flux subsample, where there is no gap in the model fits, we defined the onset time as the time when the amplitude of the fits starts to increase, similar to the other subsamples.

Fig. 7. Parameters of the best-fit Gaussians to the longitudinally averaged latitudinal flow (see Fig. 6). The top row shows the amplitudes A+ and A−, the middle row shows the peak positions μ+ and μ−, and the bottom row shows the standard deviations σ+ and σ− of the model fits (cf. Fig. 6). The four columns show the averages over different ranges of flux, with flux increasing from left to right. The error bars are Monte Carlo estimates using quiet-Sun control regions as background noise. No model was fit when there were no inflow signatures above the noise level (see text). The filled dots (circles) and solid (dashed) lines indicate poleward (equatorward) velocity. The dotted black lines indicate the onset time of the latitudinal flow for each subsample. The gray shaded region indicates the rms in the quiet-Sun control regions.

To verify that the moat flow does not affect these measurements, we performed the same analysis again, excluding those pixels from the fit that lie within the area of the moat flow for more than half of the ARs in the respective subsample. We approximated the moat flow area as regions of 2° surrounding positions with reduced continuum intensity (see Sect. 3.4 and Appendix B). This procedure leaves the model fits and onset times unaffected, except for the positive-value Gaussian in the highest-flux subsample, where only that part of the inflow that lies outside of the central region was fit, resulting in reduced amplitude and width of the fit.

The onset times illustrate what we discussed qualitatively in Sect. 4.1: The inflows towards the center of the active region set in at increasingly later times for an increasing amount of flux of the AR, in the range between one and four days after emergence. Figure 8 shows the onset time against the mean flux of each subsample. We point out that the subsamples differ substantially in flux range (cf. left panel of Fig. 4).

Fig. 8. Flow onset times in relation to mean flux. The gray boxes indicate the flux ranges of the subsamples. The plot is truncated at 2 × 1022 Mx; the range for the highest-flux subsample continues to 4.8 × 1022 Mx.

We measured the extent of the inflows from the models as the peak position μ added to the half-width at half-maximum (HWHM) of the fit. With this, we find inflow extents of about 5–7° with an error of about 1°. After five days, the inflows around the lowest-flux subsample decrease in width and amplitude. For the other subsamples, the fits illustrate that the flows have not reached steady state at the end of the observed time period, especially for the highest-flux subsample, where the inflows after emergence only formed two days before the end of the observational period.

5. Discussion

We investigated the evolution of surface flows associated with active regions before, during, and up to six days after their emergence to the surface, using LCT of granulation on the HEAR survey of 182 emerging active regions. We tested the processing of the flow measurements with synthetic input data and by cross-correlation with independent measurements from direct Doppler images.

About half of the active regions in the sample develop a sunspot with a clear penumbra in the observed period. These regions show moat flow structures with velocities on the order of 150 m s⁻¹ in the prograde and equator- and poleward direction.

To study the evolution of the inflows around active regions, we carried out ensemble averages in subsamples ordered by total unsigned magnetic flux and unsigned latitude of the active regions as well as an average over all regions. We then fit a model to the latitudinal flow component in the subsamples ordered by flux.

We find that AR emergence is preceded by converging flows of around 20–30 m s⁻¹, beginning one day before emergence. This agrees with results reported by Birch et al. (2019), and holds for all ARs in our sample, that is, it is independent of the maximum flux or the latitude of the AR. However, the extent of these early inflows is larger for stronger ARs.

We find that these pre-emergence converging flows cease shortly after emergence. Larger-scale latitudinal inflows form in the days following emergence. The time at which they form depends on the flux of the AR, between about one and four days, and increases with flux. There are several possibilities for the cause of this effect: It could be related to the separation speed of the polarities. Schunker et al. (2019) identified different phases of emergence with respect to the separation speed, finding that on average over all ARs in the HEAR sample, the separation stops increasing at 2.5–3 days after emergence. However, the distinction they made between regions with flux higher or lower than the median value indicates that this time also depends on the flux. It is therefore possible that the onset of the inflow coincides with the end of the separation increase. Another possibility is that it could depend on the relative amount of emerged flux, such that the inflows set in at a certain time relative to the peak flux. It could also depend on the ratio of flux within sunspots versus flux within plage because the proposed mechanism for driving the inflows depends on the plage (Spruit 2003). Further investigation is necessary to clarify this.

Our ensemble averages show a maximum extent of the inflows of around 7°, which is somewhat smaller than the 10° reported recently by Braun (2019). Typical velocities during the later stages are about 50 m s⁻¹, which is at the high end of previously reported velocities (Gizon et al. 2001; Haber et al. 2004; Löptien et al. 2017; Braun 2019). This is in part attributable to differences in spatial resolution; with lower resolution (e.g., by spatial smoothing with a broader Gaussian), the inflows would have lower amplitudes and larger extents. Moreover, our sample consists of emerging active regions, whereas previous studies analyzed well-established, long-lived active regions.

Our measurements of the flows show that they are still evolving at the end of the observed time frame. Observations covering the transition to steady state are needed to shed more light on the physical context of the flows. Because the observation of the emergence is necessary for an accurate age determination of an active region, the limitation to seven days before and after emergence is a constraint that at present is difficult to circumvent. One possibility is to add observations from different heliographic longitudes, from the recently launched Solar Orbiter, for example, which increases the time that an active region can be tracked continuously. Another option is to track the ARs from the HEAR sample to the next rotation. This is limited to the very long lived regions, however, and leaves an observational gap of more than two weeks.

Acknowledgments

N. G. is a member of the International Max Planck Research School (IMPRS) for Solar System Science at the University of Göttingen. N. G. conducted the data analysis, contributed to the interpretation of the results, and wrote the manuscript. We thank Jesper Schou, Paul-Louis Poulier and Bastian Proxauf for useful discussions. The HMI data used here are courtesy of NASA/SDO and the HMI Science Team. We acknowledge partial support from the European Research Council Synergy Grant WHOLE SUN #810218. The data were processed at the German Data Center for SDO, funded by the German Aerospace Center under grant DLR50OL1701. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). This work used the NumPy (Oliphant 2006), SciPy (Virtanen et al. 2020), pandas (McKinney 2010) and Matplotlib (Hunter 2007) Python packages.

References

Appendix A: Test of the effect of Zernike subtraction on the flow features

We wished to test the effect that the processing applied to the data (Sect. 3.1), specifically the subtraction of the fitted Zernike polynomials , has on flows resembling the inflows toward active regions that we investigate. To do this, we constructed synthetic flow maps, fit the Zernike polynomials to these maps, and subtracted the Fourier-filtered components of the fitted Zernike time series. This was done in the same way as for the actual LCT data. We created one Carrington rotation of synthetic data, with frames set 30 min apart, by rotating the grid 0.272 degrees from one time step to the next. This corresponds to a rotation period of 27.5735 days.

To create the data, we imposed synthetic flows on a Plate Carree coordinate grid (longitude, latitude) of size 450 × 450 with a grid spacing of 0.4°. B angle B₀ and P angle Φ₀ were assumed to be zero for the sake of simplicity. As a model for the flows, we used 2D Gaussian derivatives in the longitudinal direction (for the longitudinal flow component) and in the latitudinal direction (for the latitudinal flow component). We created two different setups, a periodic case and a random case. In the periodic case, the flows were equally distributed on a grid along the equator and the central meridian. In the random case, the flows were distributed randomly within ±45° latitude. The peak velocity of the flow features was 50 m s⁻¹ in all cases. We ran two different cases for the periodic and the random setup, with widths of the flow features of approximately 8 and 12° (which is where the velocities drop below 20 m s⁻¹). In addition, we ran the same tests with additional moat flows imposed on the synthetic inflows, which are realized as Gaussians with widths of 1.2° and peak velocities of 500 m s⁻¹. This was done to test whether the small-scale moat flow introduces additional deviations or leaking into the inflows.

First, we transformed the velocities of the flow maps from m s⁻¹ to pixel s⁻¹. This was done by inverting the equations given in the Appendix C of Löptien et al. (2017),

(A.1)

(A.2)

where we followed the definitions and nomenclature of Löptien et al. (2017). We then mapped the Plate Carree coordinate grid (longitude, latitude) to a CCD coordinate grid (x, y) of size 1024 × 1024 grid points, which is the effective size that the processing pipeline uses for the HMI data. The projection was made using the transformation specified by sphere2img() from the Stanford ring-diagram pipeline (Bogart et al. 2011).

We fit the Zernike polynomials, Fourier-filtered them, and subtracted them from the data. Then, we remapped from the CCD back to the Plate Carree projection. This was done using the transformation specified by img2sphere() from the Stanford ring-diagram pipeline (Bogart et al. 2011). The last step was to retransform the velocities from pixel s⁻¹ to m s⁻¹ by applying the equations from Löptien et al. (2017).

Figure A.1 shows the first time step of the random test case with inflow extents of about 12° before and after subtraction. The lower panels show cuts through the maps for v_lon and v_lat. Toward the observation limb at 60°, deviations from the synthetic data are noticeable, with velocities below 5 m s⁻¹. The shape and extent of the flows are unaltered by the processing. The peak velocities are changed by less than 10%. The tests including moat flows showed no artifacts in addition to this.

Fig. A.1. Example time step of the test with synthetic data. Top panel, left column: synthetic flow map in Plate Carree coordinates for vlat (top) and vlon (bottom). Middle: after subtracting the Zernike fits and limiting the field of view to 60° from disk center. Right: sum over the Zernike fits. The red lines indicate the positions of the line plots in the bottom panel. Bottom panel: line plots showing the synthetic data before (black) and after (blue) subtraction of the Zernike fits for vlat (left) and vlon (right).

Appendix B: Comparison to flows from direct Doppler images

We compared the velocities obtained by LCT with those of direct Doppler images on an example active region. Comparisons of flows from helioseismic measurements to Doppler data were carried out by Gizon et al. (2000), Jackiewicz et al. (2008), and Švanda et al. (2013), for example. For the comparison, we transformed the LCT flow maps from (v_lon, v_lat) to the line-of-sight (LOS) component. Roudier et al. (2013) applied similar transformations to calculate the spherical components from x-, y-, and Doppler components. We point out that the conventions used in our derivation below differ from those used by Roudier et al. (2013).

B.1. Projection of LCT to the line of sight

We used a right-handed system in which x points toward the observer and z toward solar north. In this system, latitude λ is defined between [ − 90 : 90] and longitude φ is defined between [0 : 360). Disk center is thus at λ = 0°, φ = 0° (when the B angle B₀ = 0°). The basis vectors of this system are

(B.1)

The transformation between Cartesian and spherical coordinates follows as

(B.2)

In addition, the rotation due to the changing B angle B₀ has to be corrected for. This is a rotation around the y-axis. Thus

(B.3)

With LCT, we only have information on e_λ and e_φ, and the radial component e_r is assumed to be zero. Thus, projecting the LCT velocity vectors onto the line of sight e_x follows the expression

(B.4)

B.2. Data reduction of Doppler and line-of-sight data

To compare the Doppler images to the line-of-sight projected LCT velocity images, large-scale systematics depending on disk position need to be removed. In addition, the data sets need to share a common frame of reference as well as temporal and spatial resolution.

First, we applied the same mapping procedure as was used for the magnetograms and the intensity maps to the full-disk Dopplergrams from hmi.v_720s. This creates a Plate Carree mapped cube for each EAR, with frames of 60° ×60° and a grid spacing of 0.1°. The cubes were centered on the coordinate specified in Schunker et al. (2016) and cover the time period between the start and end time provided in the HEAR survey. In this step, we also subtracted the background signal that stems from solar rotation, which was taken from averages over one-third of a Carrington rotation (hmi.v_avg120). Each image was subtracted by the first average frame whose central time as specified by the MidTime keyword lies after the observation time as specified by the T_REC keyword.

In addition to solar rotation, the observer motion in radial direction has to be corrected for. This was done by subtracting the velocity specified by the frame OBS_VR keyword for each frame.

Next, we averaged in time over one day relative to the time of emergence and applied Gaussian smoothing of σ = 0.4°. We then resampled the grid spacing of the map from 0.1° to the 0.4° of the LCT velocity maps.

The convention for the sign of the HMI Dopplergrams is positive for movement away from the observer (redshift) and negative for movement toward the observer (blueshift). This is opposite to the convention used in the transformation in Sect. B.1, where the x-axis points toward the observer. The sign of the Dopplergrams is therefore flipped for the further analysis.

B.3. Results

Figure B.1 shows the comparison of a Doppler map and the corresponding line-of-sight projected LCT map for AR 11158. The maps are averages over 24 h, starting four days after the time of emergence.

Fig. B.1. Comparison of Doppler velocity (left panel) and line-of-sight projected LCT velocity (right panel). The green contours outline the absolute radial magnetic field at 10 Gauss. The dark and bright red contours indicate the umbral and penumbral boundary, respectively. The orange shaded region in the right panel indicates the region that is excluded in the following (see text for details). The red line in the left panel indicates the position of the line plot in Fig. B.2. The axis labels are in Carrington coordinates. The center of the maps is approximately 23° from disk center, which is in the direction toward the upper left corner.

Figure B.2 shows a cut along the red line in Fig. B.1. Figure B.3 shows a density plot of the two velocities. We fit a line through the data for the relation between the Doppler and the LOS data. Because both are subject to error, we performed the fit with orthogonal distance regression (ODR), with the standard error of the 24-hour time averages as errors. The errors close to disk center are about 15 m s⁻¹ for the Doppler frames and 5 m s⁻¹ for the LOS frames, and they increase outward to 40 m s⁻¹ and 60 m s⁻¹, respectively. These are rough estimates because the individual frames are correlated. We also computed the correlation between the two flow measurements.

Fig. B.2. Plot along the red line indicated in Fig. B.1. The black and blue lines indicate the velocity measurements from Doppler and line of sight, respectively. The error bars indicate the standard error of the average. The dashed green line indicates the radial magnetic field strength. The red and orange shaded regions correspond to the red and orange regions in Fig. B.1.

The comparison (Figs. B.1, B.3, and B.2) illustrates that the Doppler velocity and the LOS velocity inferred from LCT agree well, except in the presence of a strong magnetic field. Here, velocities from LCT are lower than the Doppler velocities. This has several reasons: LCT is known to underestimate velocities at a strong magnetic field (Fisher & Welsch 2008; Löptien et al. 2017). A strong magnetic field suppresses convective blueshift, which results in a redshift signature in the Doppler maps (Stix 2002). Furthermore, the two measurements are taken at different depths (Jackiewicz et al. 2008).

Fig. B.3. Density plot of the velocity data in Fig. B.1. The red shaded dots indicate all pixels, the blue shaded dots only those outside of the orange shaded region in Fig. B.1. The continuous (dashed) line indicates the fit to the full (magnetic-field-filtered) sample, respectively. The legend provides the fit parameters for both fits as well as the correlations between the measurements.

When regions of 2° are excluded around positions where the LCT flows are excluded due to reduced or increased continuum intensity (see Sect. 3.4), the systematic outliers are ruled out and the slope of the fit between the Doppler and the inferred line-of-sight velocities is closer to 1 (cf. Fig. B.3). The correlations between the Doppler and the LOS velocities are between 0.7 and 0.9 for all observed one-day averages.

The intercept of the fit is in the range of −11 to −16 m s⁻¹ for all one-day averages, meaning that the line-of-sight velocities are systematically lower than the Doppler velocities (i.e., the Doppler velocities have a preference toward upflows at disk center). Because the projection onto the line of sight only takes into account v_lon and v_lat, the line-of-sight velocity is necessarily zero at disk center, whereas the Doppler velocity is sensitive to any motion along the line of sight, including radial motions. Several studies have investigated the rms velocity of Doppler images at disk center and the radial velocity of supergranulation, and reported velocities between 5 m s⁻¹ and 15 m s⁻¹ (Giovanelli 1980; Hathaway et al. 2002; Duvall & Birch 2010). Our findings are in accordance with this.

Appendix C: Comparison of methods for measuring the positions of the AR

We compared the two different methods that were used to measure the positions of the polarities of the active regions (cf. Sect. 3.2). Figure C.1 shows two example active regions.

Fig. C.1. Comparisons of the two position-finding methods on AR 11310 (left) and AR 11640 (right). The blue, orange, and green symbols mark the position of the leading and trailing polarity and the center of the AR, respectively. The crosses and diamonds refer to coordinates derived with the method described in Sect. 3.2 and the method described by Schunker et al. (2019), respectively. The axis labels are relative coordinates provided by the HEAR survey. The background images show the radial field Bz. The gray scale saturates at ±500 Gauss.

Appendix D: Effect of moat flows

We wished to identify the presence of moat flows in the flow data on our sample of active regions. To do this, we classified time series of six-hour averages of each active region in the sample according to the presence of a sunspot with a clear penumbra, which would potentially host a moat flow. We then performed an ensemble average on the resulting subsamples to identify the moat flow signature. The classification was made with the procedure described below.

We cropped the normalized continuum intensity frames and the magnetograms to the inner 12° ×12° around the position of the leading polarity. We then created a mask of all pixels where the absolute radial field |B_z| is below 10 Gauss and applied that mask to the continuum frame. This was used as the total number of pixels corresponding to the active region. We then counted the number of pixels in the masked normalized continuum intensity frame that fall in bins of below 0.6 (corresponding to umbra), between 0.6 and 0.95 (penumbra), and above 0.95 (no intensity darkening). The threshold at 0.95 is the same as was used to exclude pixels in Sect. 3.4, the threshold at 0.6 is a result of testing.

Connected regions of pixels with normalized intensity below 0.6 are identified as umbrae. For the largest region, the total pixel count of umbra and the percentage of penumbra in a box of 2.2° around the center of the umbra were calculated. Each time step was then classified into one of five categories with the following scheme:

A classification of 4 (no intensity darkening) was assigned if the total number of pixels in the umbra is zero and the total percentage of pixels in the penumbra, relative to the total number of active region pixels, is lower than 1. A classification of 3 (small pore) was assigned if the total number of pixels in the umbra is zero and the total percentage of pixels in the penumbra is higher than 1. A classification of 2 (large pore, no penumbra) was assigned if the number of pixels in the umbra of the largest spot is lower than 12 and the percentage of penumbra around this umbra is lower than 20. A classification of 1 (large pore, with some penumbra) was assigned if either the number of pixels in the umbra of the largest spot is lower than 12 and the percentage of penumbra around this umbra is higher than 20, or if the number of pixels in the umbra of the largest spot is larger than 12 and the percentage of penumbra around this umbra is lower than 20. A classification of 0 (clear penumbra) was assigned if the number of pixels in the umbra of the largest spot is larger than 12 and the percentage of penumbra around this umbra is higher than 20.

This classification scheme was derived and checked by inspecting several ARs. A table of all time steps between −2 and +7 days for all 182 active regions is available online¹.

Figure D.1 shows an example of the classification. The colored regions in the right panel indicate pixels with an absolute radial field above 10 Gauss, which we associate in this context with the AR. The small features at large distances from the sunspots do not contribute significantly.

Fig. D.1. Example of the sunspot classification. Left: normalized six-hour average continuum map of AR 11158. The axis labels are relative to the center of the leading polarity, as defined in Sect. 3.2. Right: same as left, with a colored overlay of pixels for which |Bz| is higher than 10 Gauss. Colors indicate no intensity darkening (green), penumbra (yellow), and umbra (red; see text for definition). The blue cross indicates the center of the largest umbral region. The time step shown here has the sunspot classification SSQ = 0.

Figure D.2 shows an example time step for averages over three different sets of data: All 161 active regions for which data are available at this time relative to emergence, only those regions without a sunspot with clear penumbra (classifications of 1, 2, 3, and 4), and only those regions with a sunspot with clear penumbra (classification of 0). It is apparent that the diverging flow at the leading polarity is connected to the presence of developed sunspots with a clear penumbra, which is consistent with the moat flow. The velocities are about 150 m s⁻¹, which is lower than typical moat flow velocities observed in evolved sunspots. A reason for this could be that the sunspots in our sample are just forming, together with the (in this context) comparatively large grid spacing of 0.4° and additional Gaussian smoothing of σ = 0.8° of our flow data, which reduces localized peak velocities.

Fig. D.2. Six hour average flow maps, centered at 2 d and 21 h after emergence. The left panel shows the average of 161 active regions, the middle panel shows the average over 119 active regions with no or only partial sunspot presence, and the right panel shows the average over 42 active regions with a sunspot with clear penumbra (corresponding to classification 0, see text). The active regions are averaged relative to the position of their leading polarity. Red (blue) indicates positive (negative) radial magnetic field saturated at ±150 Gauss. The arrows indicate the flows.

Furthermore, the moat flow seen here is not symmetric: The westward (i.e., prograde) component is stronger than the eastward (retrograde) component. The components in the north- and southward directions have approximately the same strength as the westward component. Švanda et al. (2014) also found an asymmetry between the east- and westward moat flows.

Our findings agree with those of Vargas Domínguez et al. (2008), who showed that moat flows are predominantly present in the directions where there is penumbra. The sunspots in our sample are only just forming, with additional flux joining the sunspot from the emergence site, that is, the retrograde direction. When the accumulation of flux is complete and the penumbra surrounds the umbra in all directions, the moat flows become more symmetric. At later times, however, the number of ARs that host a spot decreases, leading to high background noise.

Appendix E: Average over all EARs

E.1. Flows during the time of emergence

We performed an ensemble average over all 182 EARs in the sample to connect to the work of previous studies. In this case, we averaged in time over six hours and applied Gaussian smoothing of σ = 0.8° to the flow maps. Figure E.1 shows the early evolution of flows in several time steps. Starting one day before emergence, a converging flow toward the center of the average EAR is clearly visible, becoming stronger toward the time of emergence. This converging flow is noticeable in all four flux-sorted subsamples (see Sect. 4.1) and was found by Birch et al. (2019). During emergence, these flows vanish. At the same time, a prograde flow at the position of the leading polarity forms. This is again consistent with Birch et al. (2019).

Fig. E.1. Evolution (from left to right) during the emergence phase of the average active region for the ensemble average over all 182 ARs in the sample. Each time step is an average over six hours, centered on the labeled time. The flows are smoothed with a Gaussian of σ = 0.8°. Red (blue) indicates positive (negative) radial magnetic field, saturated at ±300 Gauss.

E.2. Flow variation in time, per polarity

Following the method of Braun (2019), we calculated longitudinal averages of the flows. For each polarity, we averaged over 6° from the center of the active region toward east (west) for the trailing (leading) polarity.

Figure E.2 shows the results for the flows outside of a strong magnetic field. The longitudinal flow components corresponding to time steps from three days onward show the retrograde flow that Braun (2019) detected. For the leading polarity, they span out to a maximum of approximately 7°, with velocities of about 20 m s⁻¹. The velocities here are lower than for the flux-binned subsamples in Sect. 4.1 because of the different averaging ranges.

Fig. E.2. Evolution of the flow in a longitudinal band of 6° width for each polarity. Left (right) column: leading (trailing) polarity. First (second) row: latitudinal (longitudinal) flow component, and the third row shows the absolute value of the radial magnetic field. The fourth row shows the normalized continuum intensity. All data are time-averaged over six hours. The flows are smoothed with a Gaussian of σ = 0.8°. The curves represent different time steps (see legend).

All Figures

	Fig. 1. Temporal averages (left panels) and standard deviations (right panels) over the full-disk flow maps from 24 April 2010 to 6 March 2020. The top and bottom rows show the longitudinal and latitudinal flow component, respectively.
In the text

Fig. 2. Median of the standard error against averaging time. The lines in the top panel show the cases of no smoothing (yellow) and of Gaussian smoothing of 0.4, 0.8, 1.2, 1.6, 2.0, and 4.0° (green to purple). All lines are from an ensemble average over all 182 ARs. The dotted line shows the least-squares fit to the case of Gaussian smoothing of σ = 1.6°. Bottom: Same as in the top panel for averages over 5, 10, 20, 45, 60, 91, and 182 control regions (yellow to purple). All lines are with a Gaussian smoothing of σ = 0.8°. Both panels show the case of the latitudinal flow component. The scale of the y-axis is different in the two panels.

In the text

	Fig. 3. Distributions of the sample of EARs with respect to the unsigned magnetic flux (USFLUX, left panel) and the unsigned latitude |λ| (right panel). The dotted red lines in the two panels separate the distributions into four subsamples, each containing 45 or 46 active regions, sorted by USFLUX and by |λ|, respectively.
In the text

Fig. 4. Evolution (from top to bottom) of the averaged active region magnetic field and flows for the four subsample averages of total unsigned flux (with flux increasing from left to right, cf. the left panel of Fig. 3). Each time step is an average over 12 h, centered on the labeled time. The axis labels are relative to the center of the AR, as defined in Sect. 3.2. For each panel, the total number of active regions that contribute to the average is given at the top. Red (blue) indicates positive (negative) radial magnetic field, saturated at ±150 Gauss. The green circles highlight the time after emergence at which the inflows are first visible in the time steps shown here. The arrows indicate the flows. A reference arrow representing a velocity of 50 m s−1 is given in the upper right corner. The black lines mark the averaging range that was used for the model fit (cf. text and Fig. 6). A moat flow is apparent in the rightmost column at 3 d and 6 h, for instance.

In the text

Fig. 5. Evolution (from top to bottom) of the averaged active region magnetic field and flows for the four subsample averages of unsigned latitude |λ| (with |λ| increasing from left to right, cf. right panel of Fig. 3). The labels, color scales, and arrow scales are the same as in Fig. 4. The diverging flow in the bottom row of the third subsample is due to a strong moat flow in AR 11158, which is the largest AR in this subsample at this time.

In the text

Fig. 6. Data and model fit of a time step of the second-highest flux subsample (⟨Φ⟩=7.0 × 1021 Mx). The black line shows the latitudinal flow, averaged over longitudes in the range indicated in Fig. 4. The solid blue line shows the model fit of two Gaussians with opposite polarity (see text for details). The dashed blue lines show the 1σ error estimates of the fit. The red markings indicate the fit parameters of the amplitude A (in m s−1), the peak position μ (relative to the center of the active region along latitude, Δ Lat), and the standard deviation σ (in degree) of the two Gaussians. The gray shaded region indicates the rms of the flow data. The fit parameters of the two Gaussians slightly deviate from the fitted curve due to superposition; this is small compared to the error, however.

In the text

Fig. 7. Parameters of the best-fit Gaussians to the longitudinally averaged latitudinal flow (see Fig. 6). The top row shows the amplitudes A+ and A−, the middle row shows the peak positions μ+ and μ−, and the bottom row shows the standard deviations σ+ and σ− of the model fits (cf. Fig. 6). The four columns show the averages over different ranges of flux, with flux increasing from left to right. The error bars are Monte Carlo estimates using quiet-Sun control regions as background noise. No model was fit when there were no inflow signatures above the noise level (see text). The filled dots (circles) and solid (dashed) lines indicate poleward (equatorward) velocity. The dotted black lines indicate the onset time of the latitudinal flow for each subsample. The gray shaded region indicates the rms in the quiet-Sun control regions.

In the text

	Fig. 8. Flow onset times in relation to mean flux. The gray boxes indicate the flux ranges of the subsamples. The plot is truncated at 2 × 1022 Mx; the range for the highest-flux subsample continues to 4.8 × 1022 Mx.
In the text

Fig. A.1. Example time step of the test with synthetic data. Top panel, left column: synthetic flow map in Plate Carree coordinates for vlat (top) and vlon (bottom). Middle: after subtracting the Zernike fits and limiting the field of view to 60° from disk center. Right: sum over the Zernike fits. The red lines indicate the positions of the line plots in the bottom panel. Bottom panel: line plots showing the synthetic data before (black) and after (blue) subtraction of the Zernike fits for vlat (left) and vlon (right).

In the text

Fig. B.1. Comparison of Doppler velocity (left panel) and line-of-sight projected LCT velocity (right panel). The green contours outline the absolute radial magnetic field at 10 Gauss. The dark and bright red contours indicate the umbral and penumbral boundary, respectively. The orange shaded region in the right panel indicates the region that is excluded in the following (see text for details). The red line in the left panel indicates the position of the line plot in Fig. B.2. The axis labels are in Carrington coordinates. The center of the maps is approximately 23° from disk center, which is in the direction toward the upper left corner.

In the text

	Fig. B.2. Plot along the red line indicated in Fig. B.1. The black and blue lines indicate the velocity measurements from Doppler and line of sight, respectively. The error bars indicate the standard error of the average. The dashed green line indicates the radial magnetic field strength. The red and orange shaded regions correspond to the red and orange regions in Fig. B.1.
In the text

	Fig. B.3. Density plot of the velocity data in Fig. B.1. The red shaded dots indicate all pixels, the blue shaded dots only those outside of the orange shaded region in Fig. B.1. The continuous (dashed) line indicates the fit to the full (magnetic-field-filtered) sample, respectively. The legend provides the fit parameters for both fits as well as the correlations between the measurements.
In the text

Fig. C.1. Comparisons of the two position-finding methods on AR 11310 (left) and AR 11640 (right). The blue, orange, and green symbols mark the position of the leading and trailing polarity and the center of the AR, respectively. The crosses and diamonds refer to coordinates derived with the method described in Sect. 3.2 and the method described by Schunker et al. (2019), respectively. The axis labels are relative coordinates provided by the HEAR survey. The background images show the radial field Bz. The gray scale saturates at ±500 Gauss.

In the text

Fig. D.1. Example of the sunspot classification. Left: normalized six-hour average continuum map of AR 11158. The axis labels are relative to the center of the leading polarity, as defined in Sect. 3.2. Right: same as left, with a colored overlay of pixels for which |Bz| is higher than 10 Gauss. Colors indicate no intensity darkening (green), penumbra (yellow), and umbra (red; see text for definition). The blue cross indicates the center of the largest umbral region. The time step shown here has the sunspot classification SSQ = 0.

In the text

Fig. D.2. Six hour average flow maps, centered at 2 d and 21 h after emergence. The left panel shows the average of 161 active regions, the middle panel shows the average over 119 active regions with no or only partial sunspot presence, and the right panel shows the average over 42 active regions with a sunspot with clear penumbra (corresponding to classification 0, see text). The active regions are averaged relative to the position of their leading polarity. Red (blue) indicates positive (negative) radial magnetic field saturated at ±150 Gauss. The arrows indicate the flows.

In the text

	Fig. E.1. Evolution (from left to right) during the emergence phase of the average active region for the ensemble average over all 182 ARs in the sample. Each time step is an average over six hours, centered on the labeled time. The flows are smoothed with a Gaussian of σ = 0.8°. Red (blue) indicates positive (negative) radial magnetic field, saturated at ±300 Gauss.
In the text

Fig. E.2. Evolution of the flow in a longitudinal band of 6° width for each polarity. Left (right) column: leading (trailing) polarity. First (second) row: latitudinal (longitudinal) flow component, and the third row shows the absolute value of the radial magnetic field. The fourth row shows the normalized continuum intensity. All data are time-averaged over six hours. The flows are smoothed with a Gaussian of σ = 0.8°. The curves represent different time steps (see legend).

In the text