Measurements of Photospheric and Chromospheric Magnetic Field Structures Associated with Chromospheric Heating over a Solar Plage Region

In order to investigate the relation between magnetic structures and the signatures of heating in plage regions, we observed a plage region with the He I 1083.0 nm and Si I 1082.7 nm lines on 2018 October 3 using the integral field unit mode of the GREGOR Infrared Spectrograph (GRIS) installed at the GREGOR telescope. During the GRIS observation, the Interface Region Imaging Spectrograph (IRIS) obtained spectra of the ultraviolet Mg II doublet emitted from the same region. In the periphery of the plage region, within the limited field of view seen by GRIS, we find that the Mg II radiative flux increases with the magnetic field in the chromosphere with a factor of proportionality of 2.38 \times 10^4 erg cm^{-2} s^{-1} G^{-1}. The positive correlation implies that magnetic flux tubes can be heated by Alfv'en wave turbulence or by collisions between ions and neutral atoms relating to Alfv'en waves. Within the plage region itself, the radiative flux was large between patches of strong magnetic field strength in the photosphere, or at the edges of magnetic patches. On the other hand, we do not find any significant spatial correlation between the enhanced radiative flux and the chromospheric magnetic field strength or the electric current. In addition to the Alfv'en wave turbulence or collisions between ions and neutral atoms relating to Alfv'en waves, other heating mechanisms related to magnetic field perturbations produced by interactions of magnetic flux tubes could be at work in the plage chromosphere.


INTRODUCTION
Since the temperature in the solar chromosphere and corona is higher than that expected for a state of radiative equilibrium, some mechanical heating is required to maintain the atmospheric temperature (Athay 1966;Withbroe & Noyes 1977;Vernazza et al. 1981;Anderson 1989). The strongest quasi-steady heating occurs in the active chromosphere and in particular within plage regions (Lockyer 1869;Shine & Linsky 1974;Fontenla et al. 1991), which emit brightly in chromospheric lines such as Hα (Chintzoglou et al. 2021 further discuss the definition of plage). To balance the large energy loss due to radiative cooling from the spectral lines, a heating energy flux of ∼ 6 × 10 7 erg cm −2 s −1 is required to maintain the tem-perature in these regions (Avrett 1981). Howard (1959) and Leighton (1959) found that spatial patterns of strong magnetic field in the photosphere correspond to regions with enhanced line emission. In addition, the Poynting flux in the photosphere is large enough to heat the chromosphere (Welsch 2015). Hence, chromospheric heating in plage regions should be related to the magnetic field; however, the details of the energy transport and dissipation are less clear (see, e.g., Carlsson et al. 2019).
One goal of this work is to observationally discriminate between the various heating mechanisms at work in plage. This, however, is a difficult task as both the dynamical and dissipation length scales are unresolved. High-resolution Interface Region Imaging Spectrograph (IRIS, De Pontieu et al. 2014) observations of the Mg II h & k spectral profiles within plage are only currently reproduced using models with high micro-turbulence De Pontieu et al. 2015;de la Cruz Rodríguez et al. 2016;da Silva Santos et al. 2020). This implies that shocks, torsional motions, or Alfvén wave turbulence can occur along the line of sight, or on a spatial scale smaller than 250 km, and may play significant roles in the energy dissipation. As a result, many potential interrelated processes have been identified, which we will now review in order to identify possible observational discriminants. Table 1 lists the various heating mechanisms proposed for the chromosphere above regions that consist of magnetic flux tubes with strong magnetic field, i.e., plage or network. In addition, the table includes how each mechanism's heating rate is expected to scale with the magnetic field.
Shock wave heating was first proposed by Schatzman (1949) as a mechanism of chromospheric heating. In association with the characteristic pattern of shock waves found in temporal evolution of spectra, Skogsrud et al. (2016) observed small bright patches occurring in the transition region over a plage, and Chintzoglou et al. (2021) measured electron temperature enhancements in the chromosphere. Having compared acoustic energy flux with radiative cooling energy flux using spectroscopic methods, Sobotka et al. (2016) and Abbasvand et al. (2020a) claimed that the acoustic wave energy expected to be dissipated in the chromosphere by the shock waves is the dominant energy source to heat plage regions. In contrast, Abbasvand et al. (2020bAbbasvand et al. ( , 2021 used the same technique and concluded that the acoustic energy deposited in plage regions is too low to balance the radiative loss, while the acoustic energy deposited in quiet Sun regions is a major contributor in balancing the radiative loss for the mid-chromosphere. A definitive conclusion about shock wave heating, therefore, has not been reached observationally. Shock waves may be generated by several mechanisms (Osterbrock 1961;Hollweg 1971;Suematsu et al. 1982;Kudoh & Shibata 1999;De Pontieu et al. 2004;Shibata et al. 2007;Nakamura et al. 2012;Takasao et al. 2013). The amplitude of upward propagating Alfvén waves excited by granule motions increases with the expansion of the flux tubes and becomes large enough to generate compressive shock waves that heat the chromosphere inside magnetic flux tubes (Matsumoto & Suzuki 2012Wang & Yokoyama 2020). Assuming that the magnetic field strength of the magnetic flux tubes is a constant, i.e., about 1 kG in the photosphere (Stenflo 1973;Parker 1978;Webb & Roberts 1978;Spruit 1979;Nagata et al. 2008;Lagg et al. 2010), the strength of the shock waves is inversely proportional to the magnetic field strength in the chromosphere (Table 1, Cranmer & Woolsey 2015). In a region consisting of multiple magnetic flux tubes, numerical simulations for weak shock waves repeatedly occurring with a short time scale less than 100 seconds yield chromospheric heating near the central part of the magnetic flux tube (Table 1, Hasan & van Ballegooijen 2008).
Torsional motions have been observed in the chromosphere resulting from vortex flows in the photosphere (Brandt et al. 1988; Wang et al. 1995;Bonet et al. 2008;Wedemeyer-Böhm & Rouppe van der Voort 2009;Wedemeyer-Böhm et al. 2012). Some numerical simulations show that torsional vortex motions can generate electric currents and shear flows inside flux tubes, which in turn heat the chromosphere via Ohmic heating and viscous dissipation (Table 1, Moll et al. 2012;Yadav et al. 2021). In one of the simulation results, the horizontal structure of the chromospheric temperature does not correlate with the vertical component of the magnetic field in the chromosphere (see Figure 6 of Yadav et al. 2021). In addition, Snow et al. (2018) simulated two flux tubes independently perturbed by vortex motions in the photosphere. At the location where the flux tubes merge in the chromosphere, many small-scale transient magnetic substructures are generated from the interaction of the twisted or perturbed flux tubes. Moreover, the interaction also amplifies the perturbations in the velocity and the magnetic field and produces upward propagating shocks that heat the chromosphere (Table 1).
Heating due to Alfvén wave turbulence is another candidate mechanism suggested for chromospheric heating (Matsumoto & Suzuki 2012;Liu et al. 2014;Ragot 2019). van Ballegooijen et al. (2011) presented in their numerical test that the heating energy rate in the chromosphere increases with increasing magnetic field strength in the corona. A steeper Alfvén wave velocity gradient enhances the reflection of Alfvén waves, creating stronger counter-propagating waves and correspondingly stronger turbulence in the chromosphere (Table 1).
Magnetic flux concentrations expanding with height are associated with strong gradients in the magnetic field as well as a large current density (Table 1, e.g., de la Cruz Rodríguez et al. 2013;Morosin et al. 2020).
Collisions between ions and neutral atoms enhance the dissipation of magnetic energy through various processes (Osterbrock 1961;Goodman 1997;Liperovsky et al. 2000;Krasnoselskikh et al. 2010;Khomenko & Collados 2012;Martínez-Sykora et al. 2020). In general, the heating efficiency increases with the diffusion coefficient of the magnetic field due to collisions between ions and neutral atoms and the current density perpendicular to the magnetic field (e.g., Goodman 2004;Khomenko et al. 2018). As a first example, torsional Alfvén waves can be damped by collisions between ions and neutral atoms (e.g., De Pontieu et al. 2001;Tu & Song 2013;Soler et al. 2015a). Soler et al. (2017) performed numerical simulations and found that the chromospheric heating due to torsional Alfvén waves damped by collisions decreases with the increase of the expansion rate of the magnetic flux tube. Assuming the magnetic field strength at the roots of the magnetic flux tubes is a constant in the photosphere (e.g., Stenflo 1973), the expansion rate should be inversely proportional to the magnetic field in the chromosphere (Table 1). Second, torsional or transverse Alfvén waves can also be damped by resonant absorption or by the Kelvin-Helmholtz instability Antolin et al. 2015;Zaqarashvili et al.  2015; Giagkiozis et al. 2016;Terradas et al. 2018). In chromospheric jets, collisions between ions and neutral atoms can play a dominant role in the diffusion of vortices produced at the boundary of the flux tubes by the Kelvin-Helmholtz instability and also in the dissipation of the magnetic energy (Table 1, Kuridze et al. 2016). Third, collisions dissipate the magnetic energy, not only in Alfvén waves, but also for fast and slow magnetohydrodynamic waves and their multifluid counterparts (Soler et al. 2013(Soler et al. , 2015bShelyag et al. 2016). In this case, these additions open a multitude of pathways for heating dependent on the magnetic geometry. Shelyag et al. (2016) demonstrated that magnetic diffusion due to collisions between ions and neutral atoms leads to large and continuous dissipation of magnetic waves along the radius of a magnetic flux tube (Table 1). In realistic threedimensional radiative-magneto-hydrodynamic simulations, strong heating frequently occurs along inclined canopies in the photosphere or on fibril-like structures connecting opposite magnetic polarities (Table 1, Khomenko et al. 2018). Finally, magnetic reconnection associated with magnetic flux cancellation or magnetic reconnection associated with horizontal magnetic fields has been suggested as sources of chromospheric heating (e.g., Isobe et al. 2008;Ishikawa & Tsuneta 2009;Priest et al. 2018). In an analytical solution derived by Priest et al. (2018), the heating energy rate is proportional to the magnetic field strength overlying the magnetic reconnection region (Table 1).
As Table 1 summarizes, important discriminators of the possible heating mechanisms include the location of the heating and the correlation between the magnetic field properties in the chromosphere and the local heating rate.
A few observations already address the heating discriminants described in Table 1. de la Cruz Rodríguez et al. (2013) concluded that, in a plage region, the chromospheric temperature at the edges of magnetic flux tubes is similar to that near the flux tube axis from comparing observed intensity profiles of a chromospheric line, Ca II 854 nm, with synthetic line profile computed from a numerical simulation (Gudiksen et al. 2011). However, they did not derive which part is heated more strongly quantitatively. Pietrow et al. (2020) and Morosin et al. (2020) measured vertically oriented fields with a field strength ∼ 400 G in chromospheric plage. Ishikawa et al. (2021) inferred the line-of-sight component of the magnetic field strength at the top of the chromosphere in a plage region from a Mg II spectro-polarimetric observation, and they presented a strong positive correlation between the field strength and the Mg II line core intensity, which is a signature of chromospheric heating. However, due to the low spatial resolution of these observations, it remains to be determined where the heating occurs on spatial scales comparable to individual magnetic flux tubes with a diameter of a few arcseconds in the chromosphere (e.g., see Figure 16 of Morosin et al. (2020)). The positive correlation might also result from the global structure of the plage region, rather than individual small heating events.
Our aim is to investigate the relationship between the signatures of heating in a plage region and the magnetic field in the photosphere and the chromosphere using spectropolarimetric observations to identify what mechanisms are at work. In order to determine the heating mechanisms, we must be able to infer the spatial distribution of the magnetic field strength and electric current density on a sub-arcsecond spatial scale (∼ 500 km), which is fine enough to clarify where heating occurs locally and in relation to individual magnetic flux tubes (Morosin et al. 2020) in the plage chromosphere. We also have strong temporal constraints due to rapid dynamic events that affect the formation of chromospheric lines; e.g., spicules evolve on a temporal sclae of ∼10 seconds (Hansteen et al. 2006;De Pontieu et al. 2007;Anan et al. 2010). Therefore, we need to measure Stokes spectra in two spatial dimensions covering magnetic flux tubes within ∼10 seconds. Integral field spectro-polarimeters are the most suitable instruments to meet these requirements (e.g., Kleint & Gandorfer 2017).  (Berkefeld et al. 2012;Berkefeld et al. 2018) improves image quality, which is blurred and distorted by the Earth's atmosphere. During the observation, we measured a Fried's seeing parameter (Fried 1966) of 6.2 ± 2.1 cm. Figure 1 shows images of the active region taken at 9:38 UT. The region was bipolar with two pores which had opposite polarities (panels a and c). Around the pores, there were plage regions, which are bright in chromospheric lines, e.g., the H I 656 nm (Hα) and the Mg II k lines (panels b and d).
Stokes spectra in a near infrared spectral window containing the He I triplet and Si I 1082.7 nm were obtained by the GREGOR Infrared Spectrograph (GRIS, Collados et al. 2012), which uses an image-slicer type Integral Field Unit (IFU) (Dominguez-Tagle et al. 2018, ;Dominguez-Tagle et al. in prep). The IFU enables measurement of Stokes spectra in two spatial dimensions with a field of view of 6.075 × 3.0 arcsec 2 . The rectangular field of view is sliced by 8 slitlets each having dimensions 6.075 × 0.375 arcsec 2 . In the raw data, the spectral images from each slitlet are clearly separated by dark bands. The first step of the reduction process identifies their position and removes the dark bands from all files (map, dark, flatfield and polarimetric calibration files), creating a continuous spectral image in the spatial direction as if it were coming from a continuous long slit. From there on, the standard procedure for dark current and flat field correction was applied as described in Schlichenmaier & Collados (2002) for the Tenerife Infrared Polarimeter (TIP, Collados 1999;Mártinez Pillet et al. 1999) data. Calibration optics installed at the secondary focus, before any oblique reflection, were used to measure the demodulation matrix at the time of the calibration. A model of the telescope was used to update the modulation matrix to the time of the observations, updating it at every individual scan position. The polarization model of the telescope uses as free parameters the refractive indices of the mirrors, which are obtained after calibration data taken with the telescope pointing at around 10 different orientations on the sky over the span of multiple days. As the last reduction step, the spectral images of the dual-beam setup are combined, and a cross-talk correction from I to Q , U and V is calculated and applied to force the continuum polarization to zero for each individual Stokes spectrum. Finally, the calibrated spectral images corresponding to each slice are separated and put in their corresponding position and time in the measured 2 dimensional field of view, creating four-dimensional data cubes of I , Q , U , V (λ, x, y, t).
During this observation, the solar image was scanned across the IFU using 2 steps in the direction parallel to the short side of the IFU resulting in a total field of view of 6.075 × 6.0 arcsec 2 . The two-dimensional spatial distributions of the full Stokes spectra were sequentially obtained with a temporal cadence of 26 seconds. In total we observed three regions marked by oblique boxes in Figure 1 panels (b) -(d): 11 maps of Plage A (black box), 13 maps of Plage B (red box) and 48 maps of Plage C (blue box). To increase the sensitivity, we binned the data in both the spatial (i.e. along the slit direction) and spectral dimensions to 0.276 arcsec pixel −1 and 5.4 pm pixel −1 , respectively. The final polarization sensitivity, defined here as the standard deviation of the locally normalized continuum intensity, are ∼ 0.6%, ∼ 0.03%, ∼ 0.04%, and ∼ 0.03%, for I, Q, U, and V, respectively. Figure 2 shows enlarged images of the plage region inside the black square in panel (a) of Figure 1, which includes the GRIS field of view for Plage C. GRIS and IRIS maps were co-aligned with HMI observations by comparing features in the continuum images (see the black contour lines in Figure  2). A two-channel slit-jaw imager provides context imaging for the GRIS IFU with filters to cover the 770 nm continuum and Hα line core. A reflective field stop placed at the entrance of the IFU directs light to the slit-jaw system, while the light passing through a rectangular hole in the field stop reaches the image slicer. Figure 1(b) shows an image of the active region taken by the Hα imager.
While GRIS observed the plage region, the IRIS spacecraft scanned the active region including the plage region nine times (OBSID: 3600504455). The scan step was 0.35 arcsec, and the raster scan cadence was 204 seconds. Spectra of the ultraviolet Mg II h & k doublet as well as the ultraviolet continuum near 283 nm were obtained by IRIS at the same time. They were calibrated with dark subtraction, flatfield correction, and geometric and wavelength calibration. The Mg II k intensity map and the ultraviolet continuum intensity map are shown in panel (d) of Figure  of the Mg II h & k emissions in order to investigate the association between the magnetic field structure and the chromospheric heating. In the chromosphere, the heating can be assumed to be balanced by the radiative cooling. We analyzed maps taken when both GRIS and IRIS scanned the same region at the same time. The number of the maps that satisfy this condition for Plage A, B, and C are 1, 2, and 6, respectively. As the Hα and the Mg II k images show (Figure 1), Plage A and B fully overlap the bright Cycle I Q U V I Q U V 1 3 0 0 0 2 0 0 0 2 0 1 1 1 0 1 1 1 3 3 2 2 2 2 2 2 2 plage region, while Plage C was at the edge (or periphery) of the plage region.

Inversion
An example full Stokes spectrum obtained by GRIS is shown in Figure 3. In order to fit all the observed spectra and infer physical quantities, we applied the HAnle and ZEeman Light v2.0 inversion code (HAZEL2, Asensio Ramos et al. 2008) to the Stokes spectra of the Si I 1082.7 nm and the He I 1083.0 nm (the red lines in Figure 3). HAZEL2 carried out three optimization cycles using weights summarized in Table  2. At the same time, a telluric line of H 2 O at 1083.211 nm was fitted using a Voigt function.
Inversions of the Si I Stokes spectra, which are formed in the upper photosphere (Bard & Carlsson 2008; Shchukina et al.  2017), were performed by HAZEL2 using Stokes Inversion based on Response functions (SIR, Ruiz Cobo & del Toro Iniesta 1992) in local thermodynamic equilibrium. The atmospheric model consists of physical variables specified at a number of fixed-height nodes. Interpolation between nodes provides the variation of each quantity with height, and perturbing the value at each node leads to the best-fit model. We used 5 nodes for temperature, 2 nodes for bulk velocity, and assumed that both the magnetic field and micro-turbulent velocity are uniform with height. With this assumption, the inferred magnetic field can be interpreted as an average over the line of sight in the photosphere. We also assume the filling factor is equal to 1.
All magnetic fields inferred from the inversion for the Si I Stokes profiles, whose signals are dominated by the Zeeman effect, suffer from a 180 • azimuthal ambiguity in a line-of-sight reference frame, Φ B , (e.g., Landi Degl'Innocenti & Landolfi 2004). We derived the other candidate solution from the one inferred from the inversion by changing Φ B by 180 • . We then classified the two candidates into the vertical and the horizontal ones in the reference frame of the local solar vertical. In Figure 4, we compare the histograms of the two classes of solutions, across the field-of-view and in the local vertical reference frame, with that found in the associated Space-weather HMI Active Region Patches (SHARPs) dataset (Bobra et al. 2014). The azimuthal ambiguity has been resolved in the SHARP data using a minimum energy method (Metcalf 1994;Leka et al. 2009). We select the vertical class of solutions as the true solution based on their favorable overall agreement with the SHARP data. The finding that most of the magnetic field vector in the photosphere was almost perpendicular to the solar surface over the plage region is consistent with that reported by Martínez Pillet et al. (1997) and Buehler et al. (2015).
The He I 1083.0 nm Stokes profiles were fitted taking into account all relevant physical mechanisms producing or modifying polarization signals in the spectral line: the Zeeman effect, the Paschen-Back effect, scattering polarization, and the Hanle effect. The atmospheric model consists of a constant property slab located above the solar surface at a height of 2 arcsec. The slab is characterized by the magnetic field vector, an optical depth of the line, the Doppler broadening, the damping parameter of the Voigt absorption profile, and the bulk velocity.
Magnetic fields inferred from the He I triplet inversion in the Hanle saturated regime are subject to 90 • and 180 • azimuth ambiguities described by Merenda et al. (2006), Asensio Ramos et al. (2008), and Schad et al. (2013). Which ambiguities influence an inversion result depends on geometrical configurations of the actual magnetic field orientation and the location of the observed region on the solar disk. For example, Figure 3 of Schad et al. (2015) shows regions of inversion ambiguities over an active region. In the south of the sunspot that they observed, the inferred magnetic fields are affected by all +90 • , −90 • and +180 • ambiguities, meanwhile some regions are unaffected by ambiguities. We generated four solution candidates for each inversion and we use the following algorithm to narrow down the solution space: 1. run the first inversion described above, 2. run three more inversions using fixed thermal parameters determined in the inversion at point 1 as well as where Φ B,1 is given by the first inversion, 3. for each of the four inversions, calculate the sum of squared residuals in Stokes Q and U in a wavelength range between the spectral line center ± the inferred line width, 4. identify the best fit solution that has the smallest sum of squared residuals, and 5. select candidates where the difference of the sum of squared residuals from the smallest one are less than the square of the observational noise, in addition to the best fit solution. Figure 5 shows the inferred magnetic field strength, B, its inclination angle with respect to the local solar vertical, θ B , and its azimuth angle with respect to the direction to the solar disk center, φ B . For the chromosphere, most of the Stokes spectral profiles have a unique solution without ambiguities. It might correspond to the case reported by Schad et al. (2015), in which there are some regions where the inferred magnetic field does not have any ambiguity. To confirm the uniqueness of the solution, we performed 2000 inversions for a Stokes profile observed at the location marked by the cross in Figure 5 with chromospheric angles for Θ B and Φ B fixed to randomly distributed values in the ranges shows the mean square of residuals normalized by the noise level. It demonstrates that the observed Stokes V excludes a magnetic field vector that has Θ B > 90 • from the solution. Moreover, for Stokes Q and U , the residuals of the best-fit solution are significantly smaller than those in the other local minima.
The inclination angle and the azimuth angle maps for the chromosphere display a relatively uniformly distributed horizontal magnetic field primarily directed towards the solar disk center (the direction indicated by the black arrow in panel c of Figure 1), while for the photosphere the magnetic field was almost vertical to the solar surface ( Figure 4 and 5). The measured change in direction of the magnetic field between the photosphere and chromosphere is broadly consistent with our expectations for this region, where the field should rapidly turn over with height to connect to the nearby opposite polarity across the large scale diffuse polarity inversion line around the head of the black arrow in Figure 1c. The dark threads connecting across the polarity inversion line observed in both Hα and the Mg II k support this interpretation (panels b and d of Figure 1). Note that for a comparable set of observations both Pietrow et al. (2020) and Morosin et al. (2020) found that the field remained mostly vertical, but their observations were situated very close to pores.
Finally, we calculated the electric current density J = (4π/c) ∇ × B using the inferred magnetic field vectors in the photosphere and chromosphere. For the derivative in the X-Y plane, which is horizontal to the solar surface, a central difference scheme with the truncation error of O(∆x 2 ) was used. For the calculation of ∂B Z /∂X or ∂B Z /∂Y , we took the average of B Z in the photosphere and B Z in the chromosphere taking into account the apparent displacement derived in section 3.3 before evaluating the derivative. For the derivative in the vertical direction, a height difference between the photosphere and the chromosphere of 2 arcsec (= 1450 km) was assumed by considering the apparent displacement.

Inversion errors
We investigated inversion errors by applying HAZEL2 to 10,000 synthetic Stokes profiles, which were calculated with randomly distributed parameters of the magnetic field and the line-of-sight velocity, V LOS , and to which random noise with variances matching the observational noise were added. For the Si I inversions, the considered range of B, Θ B , Φ B , and V LOS are, respectively, 0 G < B < 1000 G, 0 • < Θ B < 180 • , −180 • < Φ B < 180 • , and −10 km s −1 < V LOS < 10 km s −1 , where Θ B is the inclination angle of the magnetic field with respect to the line of sight. Figure 7 shows the root mean square of the difference between the inferred values and the values used for synthesizing Stokes profiles. The error in the magnetic field strength increases with the field strength. For the range of field strengths inferred in the observed plage region (∼ 100 to 800 G in Figure 5), the error in the photospheric inclination and azimuth angles are approximately equal to 5 • . The error in the lineof-sight velocity is less than 0.5 km s −1 .
For the He I inversions, since the sensitivity of the linear polarization to changes of the magnetic field strongly depends on the magnetic field, we narrowed down the range for the parameters of the synthetic Stokes spectra. Figure 8 shows probability density functions for chromospheric B, Θ B and Φ B inferred from the GRIS data. The narrowed ranges were determined to be 0 G < B < 600 G, (c) Figure 9. PDF of (a) the J Z,ph error, (b) the J Z,ch error, and (c) the JH error for plage region (blue) and its periphery (black).
such that 98% of the parameters vary within these ranges for the plage region and 83% for its periphery. The range of V LOS are the same as those for the Si I inversions above. The errors in the magnetic field strength, the inclination, the azimuth angles, and the line-of-sight velocity are almost same as those for the Si I inversions (Figure 7). In order to evaluate errors of the electric current density, we calculated the current density 100 times after first adding random noise to the inferred magnetic field vectors with a standard deviation equal to the errors derived above. Because the formation height of the He line can spatially change with a peak-to-valley of 2 Mm , we also add this variance in the height difference between the photosphere and chromosphere as a source of error. We set the variance to be 0.3 Mm in the calculations. Figure 9 shows histograms of the standard deviation of the 100 calculation results for each pixel for vertical components of the current density in the photosphere, J Z,ph , and in the chromosphere, J Z,ch , and its horizontal component, J H = J 2 X + J 2 Y . We found typical current densities of 5 × 10 3 esu s −1 cm −2 , while typical errors of J Z,ph , J Z,ch and J H for the plage region are 2 × 10 3 , 3 × 10 3 , and 2 × 10 3 esu s −1 cm −2 , respectively, and for the periphery region 1 × 10 3 , 2 × 10 3 , and 1 × 10 3 esu s −1 cm −2 , respectively.

Mg II h & k flux
The ultraviolet Mg II doublet typically appears as two broad absorption lines with self-reversed emission cores in solar spectra emitted from outside of sunspots (e.g. Durand et al. 1949;Lemaire & Skumanich 1973). Figure  10 shows a spatially averaged Mg II h & k flux profile over the regions of interest. The two self-reversed emission cores of the Mg II k and h appear at wavelengths of 279.6 nm and 280.3 nm, respectively. The flux was estimated using the iris calib.pro routine, part of the IRIS SolarSoftware tree (Freeland & Handy 1998). The routine derives the flux using the observed counts (in the data number units) with some information at the time of the observation, e.g., the number of photons per data number, the effective area, the size of the spatial pixels, the size of the spectral pixel, the exposure time, and the slit width.
We obtained the total radiative cooling energy flux of the Mg II h & k emissions by integrating the flux with respect to the wavelength in two ranges between 279.571 nm and 279.698 nm and between 280.296 nm and 279.698 nm (the red vertical lines in Figure 10). The error of the total radiative cooling energy flux was estimated to be about 0.6 % from the photon noise and the camera readout noise (De Pontieu et al. 2014).

Alignment of chromospheric maps with photospheric ones
We investigate the structure of the magnetic field not only in the chromosphere but also in the underlying photosphere. However, because the plage region was observed close to the solar limb (µ ∼ 0.65), the projection effect is significant. In addition, the horizontal magnetic field in the chromosphere implies that the magnetic flux tubes were bent or inclined toward the disk center. Therefore, we have to consider apparent spatial displacements in order to investigate the spatial correspondence between photospheric structures and strong chromospheric heating. For the spatial correspondence between structures and heating in the chromosphere, we assume that the formation height of the He I line is the same as those of the Mg II h & k lines. To derive the apparent displacements, we used spatial distributions of Mg II k intensity formed at different heights in the atmosphere. The formation of the Mg II h & k lines has been studied theoretically (Athay & Skumanich 1968;Milkey & Mihalas 1974;Gouttebroze 1977;Lemaire & Gouttebroze 1983;Uitenbroek 1997;Leenaarts et al. 2013a). The intensity in a central depression of the self-reversed profile, which is often referred to as Mg II h3 or k3, is formed in the upper chromosphere or transition region (e.g., Vernazza et al. 1981). In addition, it corresponds to a geometrical height at which the optical depth is equal to unity at the line core wavelength (Leenaarts et al. 2013b). The intensity of the two emission peaks in the line core, which are often referred to as Mg II h2 or k2, corresponds to the temperature in the middle chromosphere (Leenaarts et al. 2013b). Because departure from LTE is small in the atmosphere below the formation heights of Mg h2 and k2, the intensity profile from the peak to the line wing displays temperature variation along the line of sight from the middle chromosphere to the photosphere. The side dips of the self-reversed profiles, often referred to as Mg II h1 or k1, are formed near the temperature minimum (e.g., Ayres & Linsky 1976).
Because thermal and magnetic structures in the chromosphere are too different from those in the photosphere to compare directly, we calculated apparent spatial displacements between intensity maps in successive wavelength bands, moving from the Mg II k line wing to line core, marked by vertical lines in Figure 11(a), and the continuum image at 283 nm, which is relatively free of spectral lines. Each wavelength approximately corresponds to thermal structures at successively greater heights between the photosphere and chromosphere. At first, a rigid displacement was derived us- ing a cross correlation technique between the intensity maps of the continuum and Mg II k line wing at the location indicated by index 1 in Figure 11. Figure 12 shows that the spatial alignment of the line-wing-intensity map with the continuum intensity map shifted by the derived offset is better than that with the original continuum intensity map in the center portion of the field of view. Next, we iteratively derived rigid displacements between intensity maps with indexes i and i + 1 as i ranged from 1 to 11. Figure 11(b) and (c) show the inferred apparent displacements of an intensity map from the continuum image in each wavelength. The field of view of the analyzed intensity maps is 6 × 6 arcsec 2 and the center of the field is the same as that of the IFU observation. In order to increase the signal-to-noise ratio, intensities were integrated with respect to wavelength with a spectral range of 0.013 nm, which corresponds to the interval of the vertical lines in Figure 11. Variation in the lines of the same color, e.g., red or blue, reflect the sensitivity of this method, because the spatial samplings of IRIS in X and Y are 0.17 and 0.35 arcsec/pixel, respectively.
The total apparent displacements from the continuum images to Mg II k2 tend to be positive in X (0.0 -0.3 arcsec, i.e., 0 -200 km), which is consistent with the direction of the displacement caused by the projection effect. On the other hand, there are no significant displacements in Y except for one indicated by the black line (−0.6 arcsec, i.e., −450 km).
The coalignment strategy described above allows us to directly compare Mg II h & k emissions with the photospheric magnetic field using the cumulative displacement between maps of the continuum intensity and peak Mg II intensity (see Figure 13 rows 3-5 and subsequent analysis, below). This procedure essentially assumes that all misalignment between thermal structures at different heights acts like projection effects, but in reality there is degeneracy between thermal structure, magnetic structure, and projection effects. The observed magnetic field has both varying inclination, and changes in inclination with height, throughout the FOV. In addition, thermal structures do not necessarily represent magnetic structures (de la Cruz Rodríguez & Socas-Navarro 2011;Schad et al. 2013;Leenaarts et al. 2015;Martínez-Sykora et al. 2016;Asensio Ramos et al. 2017). Therefore, interpreted in the most conservative fashion, we can derive general relationships between the magnetic field in the photosphere and chromospheric heating, but not individual relations. At the same time, because the chromospheric field has a relatively uniform orientation, it is likely that our alignment follows magnetic structures reasonably well, so that the spatial distribution of chromospheric heating can be related to the spatial distribution photospheric magnetic flux.

RESULTS
In order to distinguish the heating mechanisms operating in the chromosphere over plage regions, we investigate relations between the heating energy flux and the magnetic field in the photosphere and the chromosphere. Figure 13 shows maps of the total integrated radiative cooling energy flux of the Mg II h & k lines, the chromospheric magnetic field strength, B ch , the photospheric magnetic field strength, B ph , the photospheric line-of-sight velocity, the solar vertical component of the electric current density in the photosphere, J Z,ph , and the horizontal component of the electric current density, J H . For the velocity, blue and red indicate velocities of the Si I moving toward and away from the observer, respectively. Since the thermal conduction is negligible in the chromosphere, the heating energy flux should balance with the radiative cooling energy flux in a steady state (Priest 1982). The Mg II h & k radiative cooling energy flux maps display where chromospheric heating occurs because about 15% of the radiative energy is released as emissions of the Mg II lines (Avrett 1981;Vernazza et al. 1981). The rest of the radiative energy is mainly released as H I and Ca II lines. In the photospheric maps (the three bottom rows), the red contour lines indicate the place where the chromospheric heating energy flux is large when shifted by the displacement derived in section 3.3 to compare the Mg II flux with photospheric structures.
From left to right, the columns of Figure 13 display the mapped regions for Plage A at 9:27 UT, Plage B at 9:30 UT, Plage B at 9:34 UT, and Plage C, i.e., the periphery of the plage region, at 9:38 UT. In Plage B, the magnetic field structure in the photosphere and the chromosphere, and both components of the current density structure were stable, while the area with large heating energy flux changed in 4 minutes.

Chromospheric radiative flux vs. magnetic field strengths
We compare chromospheric heating to the photospheric and chromospheric field strengths for both the plage interior and plage periphery data sets. We find that the plage periphery exhibits a roughly linear trend between heating and magnetic field strength, while the trend breaks down when for the plage-interior regions. As we argue below, this change in trend is not due to noise, and therefore may reflect a true physical relationship. The photospheric magnetic field is strong in down flow regions (red in the Figure 13 maps V ph ) where granule-scale convection gathers the magnetic flux to the inter-granular lanes and forms strong magnetic field concentrations. This effect is more pronounced in the plage interior (Plage A and B) compared to the periphery (Plage C).
Visual comparison between the Mg II heating at the magnetic field gives somewhat mixed results. Clearly the strongest heating is often not co-spatial with the strongest magnetic field strength: compare the red contours of Mg II flux to the field strengths in Figure 13 rows 2 and 3. In the plage interior there appears to be a preference for strong Mg II h & k flux at the edges or between patches of strong photospheric magnetic field (the blue contours in Figure 13 row 3), but this was not always the case. The spatial correspondence between the Mg II h & k flux and B ch is less conclusive than that to B ph , due to the more uniform field strengths found in the chromosphere.
To better tease out real correlations we generated 2D histograms of the Mg II flux versus B ph and B ch in the plage interior and periphery. Figure 14 shows the results, where the color scale indicates the number of pixels in each bin out of 4150 and 8300 total pixels for the plage and periphery regions, respectively. The histograms involving B ph include the displacement derived in section 3.3. As reflected in the figure, we recover the general trend supporting a non-linear increase in chromospheric energy flux versus the photospheric field strength, as reported by Skumanich et al. (1975), Schrijver et al. (1989), and Barczynski et al. (2018).
Here that trend is also reflected in the trend between the radiated flux and B ch . We note that in some cases, the magnetic field in the chromosphere is stronger than in the photosphere, i.e. at a given spatial pixel. This is particularly the case in the periphery, and we expect it can be explained by the horizontal expansion of the magnetic flux tubes in the chromosphere.
In the periphery of the plage region, the Mg II h & k flux had strong positive correlation with B ph and B ch with linear Pearson correlation coefficients, r, of 0.72 and 0.70, respectively. We fitted the relation between the Mg II flux and B ch with a linear function with a proportionality factor equal to (2.38 ± 0.03) × 10 4 erg cm −2 s −1 G −1 (the white solid line in Figure 14c). The dotted white line shows the same relation but shifted in +Y direction from the linear fit by two times of the standard deviation of the residual. The white solid and dotted lines in Figure 14d are the same as those in Figure 14c.
The Mg II h & k flux and the magnetic field strengths are continuously distributed from the periphery to the plage on the density plot (see yellow contour lines in panels b and d), since a part of the field-of-view for Plage B overlapped that for Plage C (Figure 1). In the plage region, the magnetic field strengths and the Mg II h & k flux were stronger than those in the periphery of the plage region, and they did not have significant correlation in contrast to the periphery (Figure 14 b and d). Since the errors of B ph and B ch are 30 -40 G at B = 400 G (Figure 7a) and the error of Mg II h & k flux is 0.6 % (Section 3.2), the lack of correlation cannot be explained by the observational noise on these spatial scales; rather, it may be reflective of the true physical relationship. Moreover, the high density peak above the 2-sigma fit in Figure 14d (dotted line) shows that the Mg II h & k flux in the plage interior is significantly stronger than the linear relation derived from the plage periphery region.
To isolate the regions with significantly enhanced radiative flux we selected regions that are above the dotted line in Figure 14d. In Figure 15 we overplot contours of these highflux regions on the maps of B ph (after accounting for spatial displacements described in section 3.3). This figure demonstrates that the strongest heating within the plage interior regions appears to come from either the edges or in-between the strong photospheric concentrations.

Chromospheric radiative flux vs. electric currents
In regions where the Mg II h & k flux was large, J Z,ph was weak and had positive or negative signs, i.e., was consistent with zero, while the strong photospheric patches correlate relatively well with strong positive J Z,ph (Figure 13). In particular, strong Mg II h & k flux seems to be emitted from areas around strong J Z,ph located at (707", −224")  Figure 14). Otherwise, there was no significant correlation. Since the 1σ error of J Z,ch for the periphery and the plage regions are 2 × 10 3 and 3 × 10 3 esu s −1 cm −2 , respectively (Figure 9), we might not be able to find significant relations in panels (a) and (b) of Figure 16 due to sensitivity of the measurements.

Temporal variations of the magnetic field and the line-of-sight velocity
We investigate not only spatial relations between the Mg II flux and the magnetic field structures but also relations of the Mg II flux to amplitudes of temporal variations in the chromospheric magnetic field and the velocity. If there is a significant correlation among them, waves can be at work. Figure  17 shows two-dimensional histograms of the Mg II h & k flux versus the standard deviation of temporal variation of physical quantities in 104 seconds before IRIS scanned the region. We do not find any significant correlation between the Mg II h & k flux and the standard deviations of the magnetic field or the line-of-sight component of the velocity in the chromosphere.

SUMMARY & DISCUSSION
Basal heating of chromospheric plage regions requires the largest quasi-steady heating rate of any region of the solar atmosphere from the photosphere through the corona. Table 1 lists many proposed plage heating mechanisms and how heating relates to the magnetic field. High spatial resolution observations of chromospheric magnetic fields and heating are crucial for distinguishing between the proposed mechanisms. We inferred the magnetic field vector in the photosphere and the chromosphere over a plage region using spectro-polarimeteric data taken with the integral field unit mode of GRIS installed at the GREGOR telescope. The integrated radiative cooling energy flux of the Mg II h & k emissions, which is used as a proxy for the total radiative losses that balance the deposited energy, was calculated from the spectra obtained by IRIS. Figure 18 provides a schematic summarizing our findings. In the periphery of the plage region within the limited field of view seen by GRIS we found that the radiative cooling energy flux strongly correlates with B ph and B ch with linear Pearson correlation coefficients of 0.72 and 0.70, respectively. The linear relation holds over an extensive range of chromospheric field strengths (50 -400 G) with a proportionality factor of (2.38 ± 0.03) × 10 4 erg cm −2 s −1 G −1 .
We found somewhat different behavior in the plage interior regions (Plage A and B): 1. The radiative cooling energy flux, i.e, heating energy flux, did not have significant correlation with B ch .
2. The heating energy flux was typically larger between magnetic patches in the photosphere or at the edges of the patches, but this was not always the case.
3. Between the magnetic patches or at their edges, the heating energy flux is significantly larger than expected from the linear function of B ch , which is derived in the plage periphery region.
4. Strong heating seems to occur around a magnetic patch that has strong electric current density in the photosphere.
5. Regions where strong heating occurs were observed to change their shapes on timescales of approximately 4 minutes.
6. The heating energy flux and the magnetic field strengths were stronger than those in the periphery of the plage region.
7. The heating energy flux did not have significant correlation with J H or J Z,ch .
8. The heating energy flux also did not correlate with the standard deviation over temporal variations of physical quantities in the chromosphere.
9. We verified that B ph is strong in the down flow regions, which are inter-granular regions. We also verified that B ch was stronger than B ph in some cases.
Based on these findings, we discuss what mechanisms may heat the chromosphere in the plage region and its periphery. For the periphery of the plage region, the heating energy positively correlates with B ch (Figure 14). Positive correlation between B ph and the radiation from Ca II H & K have been also reported at network regions where the magnetic field has an important role for heating in the chromosphere (e.g., Rezaei et al. 2007;Beck et al. 2013). Assuming the magnetic field in the corona is the same as that in the upper chromosphere, magnetic flux tubes in the plage periphery can be heated by Alfvén wave turbulence or by collisions between ions and neutral atoms relating to propagating torsional Alfvén waves (Table 1, van Ballegooijen et al. 2011;Soler et al. 2017). Interestingly, both mechanisms relate to Alfvén waves driven at or below the photosphere. Although magnetic reconnection could also be responsible for the heating, it would have to be in the strong-guide field regime (Katsukawa et al. 2007;Sakai & Smith 2008;Nishizuka et al. 2012): the fraction of horizontal or opposite polarity field is too small in our observations to produce large-scale anti-parallel magnetic fields (Figure 4). It is still possible that small and weak opposite polarity patches were Figure 18. Relation between the radiative cooling in the chromosphere and the magnetic structures in the photosphere and the chromosphere (a) in the plage region and (b) in the plage periphery. The thick blue lines labeled B signify magnetic field lines. The thick black and gray arrows indicate J Z,ph and V LOS,ph , respectively. The colored ovals indicate heated chromosphere. The redder the color in the ovals, the stronger the heating. hidden beneath the expanding flux tubes (Buehler et al. 2015;Wang et al. 2019;Chitta et al. 2019).
The two-dimensional histogram of the Mg II h & k radiative energy flux versus B ch for the plage region lies approximately along the linear trend derived from the radiative energy flux and B ch in the periphery region (panels c and d of Figure 14). These findings suggest that the same mechanisms can heat both the plage interior and periphery regions. However, we find no significant correlation between the radiative energy flux and B ch within the plage region (finding 1). In addition, radiative energy flux that was significantly larger than expected from the periphery region was emitted from between magnetic patches or at edges of the magnetic patches in the photosphere (findings 2 and 3). This suggests that other additional mechanisms can heat the plage chromosphere between the patches or at their edges.
Based on the heating location near or between flux tubes, possible heating mechanisms include shock wave heating due to interactions between magnetic flux tubes (Snow et al. 2018), ion-neutral collisions at vortices produced by resonant absorption or Kelvin-Helmholtz instability (Kuridze et al. 2016), or the ion-neutral collisions related to electric current frequently located along edges of flux tubes (Khomenko et al. 2018) (see Table 1). We propose the following mechanism in this paper putting together these three processes. The observed magnetic flux tubes can be perturbed or twisted by the granule motions. Such perturbations will propagate upward along the tubes (e.g., Shelyag et al. 2016) and produce vortices in the chromosphere (Yadav et al. 2021). Resonant absorption or Kelvin-Helmholtz instability might also contribute to the vortex formation at the boundaries of the tubes (e.g., Kuridze et al. 2016). When multiple flux tubes interact or merge at the chromosphere, velocity and magnetic field structures of the vortices present in each tube become superposed and generate complex substructures between flux tubes (Snow et al. 2018). Such substructures would accompany electric current density perpendicular to the magnetic field, where the magnetic energy strongly dissipates due to the ion and neutral collisions (Khomenko et al. 2018;Yadav et al. 2021).
The above scenario is supported by the lack of a significant correlation between the heating energy flux and B ch in the plage region (finding 1). This is consistent with what Yadav et al. (2021) found in their simulation for chromospheric heating due to electric currents associated with vortex interactions, although they studied vortex interactions inside a magnetic flux tube. Snow et al. (2018) proposed a chromospheric heating mechanism due to upward propagating shock waves with a flow speed of ∼ 50 km s −1 , which is generated by superposition of vortices between magnetic flux tubes. However, no such high speed upward movement in the plage chromosphere was observed by GRIS for 104 seconds before IRIS scanned the region. Because the high-speed plasma can pass the formation layer of the He I line in a short time, our time cadence of 26 seconds might be too long to detect the signature of the shock wave.
Since our observation showed that the heating changed in time, and that the heating energy was not always large in all regions between the magnetic patches and at their edges (findings 2 and 5), determining the location where the additional heating occurs is a necessary but not sufficient condition for distinguishing between different heating mechanism. In the plage, oscillations in intensity of spectral lines formed in the chromosphere occur with a period of ∼ 5 minutes (Bhatnagar & Tanaka 1972). Kostik & Khomenko (2016) reported that the brightness of Ca II H correlates with the power of 5 minute oscillation in the photosphere. As Tarr et al. (2017) and Tarr & Linton (2019) described in detail, magneto acoustic energy can dissipate at topological features of the magnetic field, e.g., magnetic nulls, magnetic minima, or separatrix surfaces. This allows another pathway to spatially and temporally localize heating that would act in combination with the mechanisms above but modify them based on the global magnetic connectivity.
One requirements for the formation of a plage region is a closed field region that connects opposite polarity patches through a hot dense corona (Carlsson et al. 2019). Indeed, it is implied by our observation that the plage region was connected to the opposite polarity region (section 3.1.1). It is possible that some aspects of the overlying corona may play an important role in heating of the plage regions.
We found that all strong photospheric flux regions have the same sign vertical photospheric current density ( Figure  13 and 18). The expansion of same-twist photospheric flux tubes will result in anti-parallel (plus guide-field) magnetic fields at the boundaries of the flux domains in the chromosphere. This can lead to magnetic energy dissipation due to ion-neutral collisions or magnetic reconnection (e.g., Leake et al. 2013;Ni et al. 2020).
We do not find a correlation between the radiative energy flux and the electric current density (finding 7, Figure 13 and 16). Because the width of an electric current sheet is only a few 10s of kms and their polarity is mixed at interfaces of the vortices (see Figure 6 of Yadav et al. (2021)), the spatial resolution of our observations is not enough to resolve the current sheets. We also did not find perturbations in the velocity or the magnetic field on magnetic flux tubes, which is expected in our proposed process ( Figure 17). Shelyag et al. (2016) derived amplitudes for magnetic field fluctuations due to Alfvén and fast-mode magnetoacoustic waves up to 5 G, which is too small to be detected using the present measurements (Figure 7). Upcoming observations with the National Science Foundation's Daniel K. Inouye Solar Telescope will enable us to search for and resolve electric current sheets or the magnetic field fluctuations (DKIST, Rimmele et al. 2020).
The authors thank Ms. K. Gerber for operating the GRE-GOR telescope during the observation. We also thank Dr. X. Sun, Dr. M. Kramar, and Dr. S. Mahajan, who discussed with us every week. The authors are grateful to Dr. A. Asensio Ramos for his guidance in the use of HAZEL v2.0. TA extends thanks to Dr. S. Nagata for discussing various science objectives during tea time at Hida observatory. The concept in this paper developed from those discussions. The National Solar Observatory (NSO) is operated by the Association of Universities for Research in Astronomy, Inc. (AURA), under cooperative agreement with the National Science Foundation. The GREGOR solar telescope used in this study was built by a German consortium under the leadership of the Leibniz-Institute for Solar Physics (KIS) in Freiburg with the Leibniz Institute for Astrophysics Potsdam, the Institute for Astrophysics Göttingen, and the Max Planck Institute for Solar System Research in Göttingen as partners, and with contributions by the Instituto de Astrofísica de Canarias and the Astronomical Institute of the Academy of Sciences of the Czech Republic. The SDO data are provided courtesy of NASA/SDO and the AIA and HMI science teams. This work was supported by JSPS KAKENHI Grant Number 20KK0072.