Reconciliating the Vertical and Horizontal Gradients of the Sunspot Magnetic Field

In the literature, we found 15 references showing that the sunspot photospheric magnetic field vertical gradient is on the order of 3-4G/km, with field strength decreasing with height, whereas the horizontal gradient is nine times weaker on the order of 0.40.5G/km.This is confirmed by our recent THEMIS observations. As a consequence, the vanishing of div ⃗ B is not realized. In other words, a loss of magnetic flux is observed with increasing height, which is not compensated for by an increase of the horizontal flux. We show that the lack of spatial resolution, vertical as well as horizontal, cannot be held responsible for the nonvanishing observed div ⃗ B. The present paper is devoted to the investigation of this problem. We investigate how the magnetic field is influenced by the plasma anisotropy due to the stratification, which is responsible for an “aspect ratio” between horizontal and vertical typical lengths. On the example of our THEMIS observations, made of two spectral lines formed at two different depths, which enables the retrieval of the three components entering div ⃗ B, it is shown that once this aspect ratio is applied, the rescaled div ⃗ B vanishes, which suggests a new methodology for MHDmodeling in the photosphere.


Introduction
In the literature, we found 15 references reporting magnetic field gradient measurements in sunspots, with various lines, instruments, and interpretation methods. We restricted to references later than 1983. A detailed study of the references shows that they may be classified into two types: those where the vertical gradient is measured, which is found on the order of 3-4 G/km with the field strength decreasing with height, and those where alternatively the horizontal gradient is obtained, which is found on the order of 0.4-0.5 G/km. These orders of magnitude were completely retrieved in our recent THEMIS observations. This classification does not suffer any exception. The two gradients differ by a factor of nine. This problem was pointed out by other authors. Solanki [1] investigated the problem and concluded that "no satisfactory solution has been found as yet for the unexpectedly small (horizontal) gradients. " As a consequence, the vanishing of div ⃗ is not realized, and the present paper is devoted to the investigation of this problem.
In Section 2, we report in details about these observations, and we discuss different possible causes on the observational side of the question. We consider the noise effect, the unresolved structure effect, and the line-of-sight and pixel integration effects. As none of these effects has brought the explanation of the phenomenon, we turn to the theoretical side of the question in Section 3. We evaluate in details the different quantities entering the MHD modeling of the solar photosphere. We compute the permeability and the electric conductivity, which we found depth dependent. But we show that the variable conductivity is not yet the answer to our question. We finally find a way out by investigating a plasma effect, the Debye shielding. But this is not yet a sufficient condition. In addition, it has to be considered that the Debye sphere is flattened under the effect of the medium stratification. In this last condition, we show that div ⃗ may mathematically depart from zero. The stratification of the medium is a necessary condition for the effect being possible. We show that the Debye sphere may be flattened under the effect of the anisotropy of the particle velocities. The photosphere being the surface layer of a star is a strongly stratified plasma where such an anisotropy exists. The effect is discussed in Section 4, and finally a new methodology for the MHD modeling of the photosphere is suggested. We conclude in Section 5.
We wrote this paper as we were along the investigation of the different possibilities. During this study, we elaborated a detailed characterization of the photosphere, in particular when we investigated the effect of the variable electric conductivity. The variable conductivity did not finally reveal to be efficient to explain the observed phenomenon, but we found nevertheless interesting to report here the results of our photosphere investigation, summarized in Tables 1 and  2, and to report all the investigations including those that revealed unsuccessful, in order to draw a complete picture of the problem. Reading this paper, it has to be kept in mind that finally we propose the Debye sphere flattening by the stratification as the origin of the observed phenomenon, where div ⃗ may mathematically depart from zero.

The Ensemble of Results.
With THEMIS, we found that the vertical gradient of the magnetic field in the NOAA 10808 sunspot umbra |Δ /Δ | is 3-4 G/km (Figure 1(a)), about nine times larger than the horizontal gradient |Δ /Δ + Δ /Δ |, which we find on the order of 0.4-0.5 G/km (Figure 1(b)). Because the value of | div ⃗ | is obtained from the difference between these two numbers, there is a discrepancy between its theoretical value, zero, and its "observed" value because the difference between these two numbers cannot be zero. In order to investigate this discrepancy, we turned to the literature and we detail below published methods and results. We restricted our investigation to results obtained after the 80 s in order to consider only numerical data acquisition systems. In general, only one of the two gradients was derived in a given paper, and there are only a few papers where the difference was detected and discussed. We report on those papers at the end of the section. The case of the horizontal gradient is particular, because the authors generally present their result as a vertical gradient indeed derived by applying div ⃗ = 0 to the horizontal gradient, which is the true quantity that they measured. It is then necessary to deeply investigate each paper in order to determine the exact nature of the observed gradient. We then divide the papers into two groups, one about the vertical gradient measurements and the other one about the horizontal gradient measurements.

Vertical Gradient.
Let us consider first those authors who directly measure the vertical gradient of the vertical component of the magnetic field Δ /Δ . The determination of the -variation results from a line formation modeling. It may be remarked that a modern inversion provides the total magnetic field strength and not the line-of-sight (l.o.s.) component one. But the spot under study is generally located close to disk center, and the field is mainly vertical in the spot umbra, so that the vertical gradient and the l.o.s. one may be confused in the discussion. We point out the cases where the observed spot is far from disk center.
Westendorp Plaza et al. [2] used the SIR inversion code to invert spectropolarimetric data of the Fe I 6301.5 and 6302.5Å photospheric lines observed with the Advanced Spectropolarimeter (Sacramento Peak). The SIR inversion code [3] computes the line response functions so that the line formation depth was also derived. The temperature gradient was also retrieved and LTE was assumed, which is convenient for photospheric lines. Westendorp Plaza et al. [2] reported Table 1: Some parameters in the solar quiet photosphere from the VAL C atmosphere model: height ℎ above the 5000 = 1 level, temperature , neutral hydrogen H and electron densities, Debye length , Coulomb electron-ion ] ei , H − formation by electron attachment ] H − and elastic but deviating e − + H collision ] el frequencies, electric conductivity , magnetic diffusivity = 1/ 0 , magnetic Reynolds number , and magnetic susceptibility .  1.5-2 G/km (decreasing field with increasing altitude), and 4 G/km in a previous analysis. Balthasar and Schmidt [4] interpreted spectropolarimetric observations made at the VTT (the German Vacuum Tower Telescope at the Observatorio del Teide, Tenerife) in three lines Fe I 6302.5 and 6842.7Å and Fe II 6149.2Å. The magnetic field was retrieved by comparing the observed profiles with theoretical ones computed following the method developed by Grossmann-Doerth et al. [5] and Grossmann-Doerth [6] where the Unno-Rachkovsky equations were integrated in a model atmosphere (three different empirical models were used), and the difference in formation depth was derived from the line contribution functions. They obtained 2.5-3 G/km for the vertical gradient |Δ /Δ |. Schmidt and Balthasar [7] analyzed the umbral dots of the same observation and found the same gradient value of 2.5-3 G/km, although the field strengths were weaker in the umbral dots.
With the Gregory-Coudé Telescope (GCT) at Tenerife, Pahlke and Wiehr [8] observed the circular polarization of six lines in a spot umbra, namely, Si I 6142.9, Zr I 6143.2, Fe II 6149.2, Ti I 6149.7, Fe I 6151.6, and Na I 6154.2Å. In the spot umbra the Zeeman splitting was visible in the circular polarization so that the direct measurement of the field strength was possible. A refined model of the circular polarization was applied with a magnetic field vertical gradient and the best agreement between the six lines resulted in a gradient of 2 G/km. Also with the GCT, Collados et al. [9] observed two spots in three Fe I lines at 6297.8, 6301.5, and 6302.5Å. The SIR code was run for the inversion. The authors found 4.7 and 3 G/km as vertical gradients in the two spots, respectively.
Bruls et al. [10], reanalyzing FTS (Kitt Peak) observations by Hewagama et al. [11] in two infrared Mg I lines at 12.22 and 12.32 m, obtain a vertical gradient of 2-3 G/km in the sunspot umbra. These lines are formed in the upper photosphere. The theoretical profiles are issued from the non-LTE multilevel MULTI code [12], based on several model atmospheres, and followed by the DELO Stokes profile synthesis code [13][14][15].
With the Tenerife Infrared Polarimeter (TIP) mounted on the VTT, Mathew et al. [16] observed two infrared Fe I lines at 15,648.5 and 15,652.8Å that are formed deeply in the photosphere. The spectropolarimetric observations were inverted by using the SPINOR code [17], based on the Unno-Rachkovsky solution of the radiative transfer, LTE hypothesis, and response functions for an efficient computation of derivatives of the observation/modeling merit function with respect to the free parameters, similarly to the SIR code. These response functions enable the derivation of the line formation depth. Mathew et al. [16] obtained a vertical gradient of 4 G/km in the umbra.
Finally, Berlicki et al. [18] observed a sunspot far from disk center, at cos = 0.72, in different points along the Na I D 1 line profile with THEMIS/MSDP. These wavelengths along the profile are associated to different heights in the atmosphere via response functions computed with the MULTI code [12]. The derived value of the gradient was given in relative value, that is, divided by the magnetic field strength itself. Assuming the typical spot field strength of 2500 G, and taking into account the heliocentric angle given by cos = 0.72, the corresponding vertical gradient is 2 G/km. The analysis would not be complete if we do not cite one measurement by Rüedi et al. [19] that is made between a photospheric and a chromospheric line. Their observations were taken at the Mc Math Pierce telescope at Kitt Peak in umbra and penumbra of a spot not too far away from disk center. The photospheric line was Si I 10,827.14Å, unfortunately not an ideal tool for measuring photospheric fields, since it is insufficiently Zeeman sensitive. The chromospheric line was He I 10,830Å. The authors obtained a vertical gradient of 0.4-0.6 G/km in the umbra, weaker than the purely photospheric results, but the difference in depth formation between a photospheric line and a chromospheric line is probably too large to derive a precise value of the gradient, or the gradient is not the same in the photosphere and the chromosphere.
If we discard this result that involves a chromospheric line, it can be concluded that without any exception, in the photosphere the observed vertical gradient is on the order of 3-4 G/km.

Horizontal Gradient.
The second type of methods is as follows: discarding modeled determination of thedependence, the transverse field is measured in a sunspot and the ambiguity is intuitively solved from the spot polarity, which enables the determination of −(Δ /Δ + Δ /Δ ) that is claimed to be Δ /Δ via div ⃗ = 0.
With the Tenerife Infrared Polarimeter (TIP) mounted on the VTT, Balthasar [20] observed three Fe I lines at 15,648, 15,452, and 10,896Å in 8 spots. The spots were not all close to disk center, but the results were corrected from 3D effects. He applied the SIR inversion to determine the horizontal magnetic field (but not the formation heights) and derived 0.5 G/km in the umbra.
With the MSFC magnetograph (Hunstville, Alabama), Hagyard et al. [21] observed the Fe I very magnetically sensitive 5250.2Å line in one sunspot located near disk center. From the transverse field components measured by this magnetograph, they report a horizontal gradient of 0.1-0.3 G/km in the umbra and internal penumbra.
Liu et al. [22] applied the same method to observations of a sunspot in the Fe I 5324.2Å line with the Solar Magnetic Field Telescope (SMFT) at Huairou. The measurement technique was also the magnetograph one. The horizontal gradient was found to range between 0.1 and 0.7 G/km, with average value on the order of 0.2 G/km. The spot was not so close to disk center but the average value was unchanged once the geometrical correction for disk center departure was introduced.
In their observations of the two infrared Fe I lines at 15,648.5 and 15,652.8Å, Mathew et al. [16] derived also the radial (or horizontal) gradient and found it ranging between 0.05 and 0.18 G/km in the umbra and 0.12 G/km in the penumbra.
Balthasar and Collados [23] derived also 0.18 G/km in the umbra and internal penumbra of an isolated sunspot 24 ∘ away from the disk center (the 3D effect was corrected), observed with the TIP mounted on the VTT, and by applying the SIR inversion code to the Fe I 10,896.3Å line.
Beside all these measurements, the typical order of magnitude can be considered: assuming a penumbral field = +1500 G that switches to = −1500 G at the opposite side of a typical umbra of 10,000 km diameter, the resulting average gradient would be 0.3 G/km.
It can be concluded that without any exception in the photosphere the observed horizontal gradient is at most on the order of 0.4-0.5 G/km.

THEMIS Observation Results.
With the polarizationfree telescope THEMIS, we observed the double sunspot (spot) of the active region NOAA 10808 on 13 September 2005 between 14 : 25 and 15 : 25 UT. THEMIS was operating in MTR (MulTiRaies) mode where the solar image is scanned over the spectrograph slit. We analyzed the Stokes spectra of the pair of lines Fe I 6301.5 and Fe I 6302.5Å. For deriving the magnetic field, we performed the UNNOFIT inversion [24,25] to Fe I 6302.5Å. This is Milne-Eddington inversion, based on the Unno theory [26], further generalized for magneto-optical effects by Rachkovsky [27,28]. For Fe I 6301.5Å, which is not a normal Zeeman triplet line, we applied the more general (but more time consuming) UNNOFIT2 code also developed by Landi Degl'Innocenti and coworkers, which is based on the more general definition of the elementary absorption coefficients with = −1, 0, +1 for the , , and components of the Unno-Rachkovsky theory, respectively [29] (p. 386, for notations).
( , V) is the Voigt function.
In order to determine div ⃗ , the formation depth difference between the two lines has to be determined. In the quiet sun, this difference was recently derived from HINODE observations by Faurobert et al. [30] who applied a phaseshift analysis. They obtained as observed value 63.2 ± 0.9 km, Physics Research International 5 corroborated by the value of 69 km derived by the same phase-shift technique applied to theoretical profiles computed with the non-LTE Uitenbroek's code [31]. These values are in full agreement with the 3D magnetoconvection simulation of Khomenko and Collados [32]. Besides, Khomenko and Collados [32] show that although the difference is well on the order of 60-70 km in the quiet sun, it is rather of 100 km in the active region. As we observed an active region, we applied the value of 100 km to our computations.
In addition, Khomenko and Collados [32] obtain that the formation depth difference of these two lines remain constant even if each formation depth varies. We characterize this feature by saying that these two lines behave in a parallel manner as for their formation depth. Such a parallelism probably originates in the fact that these two lines Fe I 6302.5 and 6301.5Å belong to the same multiplet (n. 816) and have different values ( = 0.180 for 6301.5 and = 0.0627 for 6302.5-data taken from the Kurucz databases). Since differential non-LTE effects within multiplets are thought to be very small (as proven by detailed, multilevel, non-LTE computations) this implies that the absorption coefficient of 6301.5 is 3 times larger than the absorption coefficient of 6302.5. There is then no doubt that 6301.5 forms higher than 6302.5.
Given the and pixel sizes, which are Δ = 1160 km and Δ = 581 km after a 2 × 2 pixel binning for noise reduction, we were able to determine div ⃗ . We obtained 3-4 G/km for the vertical gradient |Δ /Δ | in the umbra (Figure 1(a)), with field strength decreasing with height, and 0.4-0.5 G/km for the horizontal gradient |Δ /Δ + Δ /Δ | (Figure 1(b)), after having rotated the results into the solar frame. These two values are, respectively, in full agreement with the other published results as described above.

Some Modeling Results.
As for the sunspot magnetic modeling, the modeled vertical gradient is also found about one order of magnitude weaker than the observed vertical gradient. Such a behavior is reported by Eibe et al. [33] where the observed vertical gradient is derived from longitudinal magnetic field measurement at different wavelengths inside the Na I D 1 line profile observed with THEMIS/MSDP. These wavelengths along the profile are associated to different heights in the atmosphere via response functions computed with the MULTI code [12]. The modeled vertical gradient is obtained by force-free extrapolation [34]. The authors report a difference of one order of magnitude between the observed and modeled vertical gradients.
In a magnetostatic equilibrium modeling (not forcefree), Pizzo [35] reports a modeled vertical gradient of 0.2-0.4 G/km (absolute value, for large tube radii that model sunspots), which is also of the same order of magnitude as the observed horizontal gradient, but one order of magnitude weaker than the observed vertical one.

The Particular Case of the Recent HINODE/SOT/SP
Observations. Puschmann et al. [36] inverted an HIN-ODE/SOT/SP observation in a portion of penumbra where the ambiguity solution was simple, so that the full magnetic field vector was derived. The inversion code was SIR. They plotted the histogram of the measured div ⃗ values and it was obtained that the histogram is clearly not centered on zero, but on 0.2 G/km instead. As this 0.2 value corresponds to the shift of the histogram but not to its width, Puschmann et al. [36] pointed out that the 0.2 does not result from noise. In order to investigate the problem, Puschmann et al. [36] plotted individual histograms of the three contributions to div ⃗ , namely, / , / , and / . It is obtained that the shift of 0.2 G/km remains in the / histogram, whereas the histograms of / and / are both centered on zero. This result shows us that in these observations also, the contribution of the vertical gradient is prevailing in div ⃗ , as in the ensemble of results cited above. The fact that the 0.2 G/km value of Puschmann et al. [36] is smaller than the difference between the generally observed vertical and horizontal gradients comes from the fact that they studied a portion only of penumbra, whereas we reported results for the whole umbra and inner penumbra. But the same inexplicable feature is present in all the results. Puschmann et al. [36] finally ascribed this unexpected value of nonzero div ⃗ to unresolved magnetic structures, as suggested by Sánchez Almeida [37], who also detected the problem of inconsistency between different observation results when div ⃗ = 0 is applied in between.

Research of Cause of the Observed Discrepancy.
We can then conclude that without any exception, in the sunspot photosphere the observed vertical gradient is on the order of 3-4 G/km and the horizontal one is on the order of 0.4-0.5 G/km. As a consequence, there is a problem to ensure the vanishing of |div ⃗ |, whose value is obtained by subtracting these two values. One value is about ten times larger than the other one, so that their difference may not be zero. We analyze below different suggested causes on the observational side, without getting any positive conclusion.

Noise Effect.
What would be the observation noise effect? If the difference was observation noise, as the difference is about 90% of the largest value |Δ /Δ |, the noise would be as large as it so that all the values would be within noise. In that case, one would not obtain that the difference is always of the same sign (decreasing field strength with height) and order of magnitude. All the results would be dispersed, which is not the case. One always observes a flux absolute decrease with increasing height in the atmosphere.
It is however necessary to investigate in a deeper manner the question of the respective accuracy of the determination of the different magnetic field components, because it is generally believed that with the Zeeman effect the longitudinal field is better determined than the transverse one. This belief is based on the weak field laws of the Zeeman effect, where the Stokes parameters of the emitted line are expressed as functions of the derivatives of the intensity profile. In these laws the circular polarization Stokes parameter is found to depend linearly on the longitudinal component of the magnetic field cos , where is the angle between the magnetic field vector and the line of sight, and on the first derivative of the intensity / . Alternatively, the linear polarization and Stokes parameters are found to depend quadratically on the transverse component, in 2 sin 2 associated to the second derivative of the intensity profile 2 / 2 [29] (p. 397 for a detailed derivation). Two ideas could follow: the first one that given a polarimetric noise, the magnetic noise will be larger for the transverse component than for the longitudinal one, and the second one that the Zeeman effect is not so linear due to the presence of the quadratic term. The first idea is verified when the weak field laws are used for interpreting the polarization measurements in terms of magnetic field vectors. This was the case of the observation studied by Bommier et al. [38], where the noise is found to be 5-10 G for the longitudinal field and 50-100 G for the transverse field in the weakest field regions. But Stenflo [39] derived from observations that the magnetic field is highly inhomogeneous, formed by strong field flux tubes filling a few percent of the space, which is the so-called magnetic filling factor. If the weak field laws were obeyed by the solar magnetic field outside active regions, the quadratic versus linear dependence of the Stokes parameters would lead to / and / of the order of 10 −2 when / is of the order of 10 −1 , for instance, so that the linear polarization Stokes parameters would be notably weaker than the circular polarization ones, in the weak field regions. This is not what observed, confirming thus the result by Stenflo [39] of a strong local field so that , , and are all of the same order of magnitude, and the weak field laws have not to be applied. One then turns towards inversion methods and codes where the four Stokes parameters are entered and treated in a simultaneous manner and not in a separate circular versus linear manner. Moreover, the magnetic filling factor has to be introduced in the inversion. Using the test data of Bommier et al. [24], we have verified that in such inversions the magnetic noise is the same for the three components of the magnetic field. In these test data the introduced polarimetric noise is 1.5 × 10 −3 and the resulting magnetic noise is found to be 5 G for the longitudinal component and 10 G for the transverse component, each averaged by the filling factor, which is quite the same.
As for the linearity of the Zeeman effect, it is ensured because either the weak field laws are valid and then / and / , which are quadratic, are much smaller than / , which is linear so that the largest term has a linear behavior, or the weak field laws do not apply, which means that the magnetic field is strong even if it does not fill the whole space, and when the magnetic field is strong it manifests in the same linear manner in the four Stokes parameters, depending on the atomic sublevels magnetic splitting, which is linear. The linear character of the Zeeman effect may thus be claimed in all the cases.
Anyway, the sunspot magnetic field studied in this paper is not a weak field, and its filling factor is found to be very close to unity, so that the result holds even for inversions without filling factors, for the sunspot magnetic field.
Pure noise effect or differential noise effect cannot then be retained as the cause of the different behavior of |Δ /Δ | and |Δ /Δ + Δ /Δ |.

Unresolved Structure Effect.
Unresolved structure effect is frequently suggested as the solution. It is obviously assumed that each unresolved structure obeys div ⃗ = 0. But div ⃗ is a linear function so that the div ⃗ of the average structure is the average of the div ⃗ of the individual structures and remains zero. In other words, if the magnetic flux is conserved in elementary structures, one cannot explain why it would not be conserved at a wider scale.
In the case of the sunspot magnetic field, an unresolved structure could be made of small magnetic loops inserted between outgoing vertical field lines, in the umbra. The observed phenomenon is that the magnetic flux at higher level is weaker than the flux at lower level and that this difference is not offset by the flux through the side surfaces. If small loops are considered, with summit lower than the higher altitude considered, they return flux at the lower level that is not evacuated at the higher level. But if the two feet of each loop are contained in the lower considered surface, the brought flux is also evacuated at the lower level so that the balance remains zero. Conversely, if the loop does not evacuate the flux through the lower surface, it does that through the side surface having its second foot outside the box, so that this would be detected also in the measurement of the side flux and the balance would also remain zero. Thus, it has not been possible to imagine an unresolved structure capable of generating nonzero div ⃗ .
A particular kind of unresolved structure is offered by turbulent microstructures. In hydrodynamics where it is div ⃗ V = 0 that is expected in incompressible media, Dubrulle and Frisch [40] studied the superposition of a small-scale turbulent perturbation to a large-scale basic structure. They obtain that "in general the [large-scale velocity field] does not [locally] have a vanishing divergence" (due to a contribution of the small-scale perturbation), but they add "in spite of the vanishing of the divergence of [the large-scale velocity field averaged on the small-scales]. " In other words the divergence of the large-scale velocity may locally depart from zero due to a contribution of the small-scale turbulent perturbation, but when the large-scale velocity is averaged over the small scales, div ⃗ V = 0 as usual. Returning to the solar observations, the average is taken over the unresolved structures so that departure from a vanishing divergence cannot be expected also due to the unresolved structures.
To resume, due to the linearity of div ⃗ and of the Zeeman effect and of the measurement technique, we estimate that unresolved structure effects cannot be retained as the searched for cause.
Physics Research International 7

Line-of-Sight and Pixel
Integrations. Due to the linearity of the Zeeman effect as previously discussed, the line-of-sight integration can be modeled as where is the contribution function, and where we have assigned a "depth of formation" to the final result. This depth of formation is generally close to the maximum of the contribution function, as visible in the examples computed by Bruls et al. [41]. The contribution function acts as a filter in the above equation.
Analogously, the pixel integration in or can be modeled with a convolution by a crenel function.
If now we compute the divergence of the observed magnetic field, we can apply the derivation of a convolution product as recalled below in Section 3.4. In the case of the line-of-sight integration, this is so that the divergence of the observed magnetic field is zero because div ⃗ is locally zero. Analogous derivation can be made for the filtering with a crenel function for the pixel integration. It results from the derivation property of the convolution product that any filter would lead to the same result, which is that the divergence of the filtered quantity is the filtering applied to the divergence of the local quantity, so that if the last one is zero, the first one is also zero.
One cannot then ascribe the departure from zero of the divergence of the observed field to the line-of-sight and/or pixel integration effects, as well as to any filtering effect.
As we were unable to solve the problem by investigating the observational side of the question, we turn now to the theoretical side.

Theoretical Investigation
In the MHD modeling, the magnetic field results from the local density current via the Maxwell-Ampère's law where the displacement current is neglected as usual at low frequencies. And the density current results from the fields via the Ohm's law where is the electric conductivity, ⃗ and ⃗ are the electric field and magnetic induction, respectively, and ⃗ V is the local fluid velocity. The reason for the existence of the ⃗ V × ⃗ term in the above equation is the frame change from the local frame moving at ⃗ V velocity, where Ohm's law is ⃗ = ⃗ , to the local rest frame where ⃗ = ⃗ + ⃗ V × ⃗ . Here we use the MKSA unit system. Introducing the Maxwell-Faraday law of induction and combining these equations, one derives the induction equation At stationary state ⃗ / = 0, this equation describes the equilibrium between the field created by the current and the current created by the field and the fluid motion. This is the situation that we are going to study, but before developing this equation we have first to investigate two questions: (i) the relation between the magnetic field ⃗ and the magnetic induction ⃗ ; (ii) the variation of the electric conductivity with height in the solar photosphere.
For calculating these relations or parameters, we consider the average quiet sun VAL C atmosphere model [43] and we report the results in Tables 1 and 2. Indeed, the series of models of Vernazza et al. [43] from the darkest region A to the brightest region F do not give very different results.

About ⃗ and ⃗ :
The Permeability of the Photosphere. The difference between the magnetic field ⃗ and the magnetic induction ⃗ lies in the fact that charged particle motions may exist at small scale, or spins, that are not free to participate to the current that creates ⃗ , but as charged particle motion they contribute to the ⃗ creation via the magnetization ⃗ where 0 and are the vacuum and medium permeabilities, respectively. As visible in Table 1, the solar photosphere is mainly made of a gas of neutral hydrogen atoms, which are paramagnetic because the total quantum kinetic momentum in the fundamental level is nonzero being = 1/2. Thus, ⃗ is parallel to ⃗ and the susceptibility defined by is given by the Curie's law where is the material-specific Curie constant and is the temperature. The Curie constant for the gas of neutral hydrogen atoms is where is the hydrogen atom density, is the Boltzmann constant, and eff is given by where is the Landé factor, which is 2 for the hydrogen fundamental level, is the Bohr magneton ℏ/2 , where and are the electron charge and mass, respectively.

Physics Research International
We have computed the susceptibility of the neutral hydrogen atom gas in the quiet photosphere conditions. The results are given in Table 1. It can be seen that in the density and temperature conditions of the photosphere, is very small, on the order of 10 −10 , so that in the photosphere ⃗ and ⃗ may be considered as identical with = 0 .

The Electric
Conductivity of the Photosphere. Following for instance Meyer-Vernet [44], the electric conductivity in a plasma is given by where is the electron density and ] is the frequency of the collision that deviates electrons. The Coulomb interaction with the plasma ions is the first candidate for these collisions, but in the photosphere where the neutral hydrogen atom density is 10 4 times larger than the electron density, we have also to consider the collisions between the electrons and the neutral hydrogen atoms. Let us recall that the photosphere is a plasma formed from a gas of neutral hydrogen atoms essentially. The plasma is globally neutral so that negative free charges circulate, which are electrons, and positive free charge circulate in equal density, which are essentially protons. Because their mass is much larger, the proton thermal velocity is much smaller than the electron thermal velocity. In the photosphere the charge density remains much smaller than the neutral particle density (see Table 1). We consider below (i) the formation of H − by electron attachment to H; (ii) the elastic collisions − + H but with electron deviation; (iii) the inelastic collisions − + H.

Electron-Ion Collisions. Following Meyer-Vernet [44], the frequency ] ei of the electron-ion collisions is given by
where V is the electron thermal velocity √3 / , Le is the Landau radius and is the Debye length We have computed the frequency ] ei of the electron-ion collisions in the quiet photospheric conditions and reported the results in Table 1.

Formation of − by Electron Attachment.
The electron may form H − by attaching to an hydrogen atom, a photon being emitted. For the collision rate per atom and electron, we have considered the average value between Callaway [45] and Janev and Van Regemorter [46], which is about H − = 4× 10 −15 cm 3 s −1 at the photospheric temperature. Considering this first approximation value, we have computed the collision frequency ] H − = H − H . We have reported the results in Table 1, where it can be seen that this frequency is about 10 7 times smaller than ] ei , so that its contribution to the conductivity can be considered as completely negligible and it is not necessary to use more refined rate values.

Elastic Collisions
Most of the e − + H elastic collisions do not deviate the electron, because the hydrogen atom is neutral. But deviating collisions may however occur. Differential cross-sections were observed by Gilbody et al. [47] and computed by Burke and Schey [48], excluding the forward scattering which is very strong. Considering the average value of these differential cross-sections of 1× 2 0 per unit of solid angle, and considering 4 as the total solid angle, and the electron thermal velocity, we obtained the values given in Table 1, which are about 10 6 times smaller than ] ei and also completely negligible so that it is not necessary to use more refined rate values.

Inelastic Collisions
The impacting electron may also excite the = 2 level of neutral hydrogen. The corresponding cross-section was computed by Kingston et al. [49]. It can be seen in their Tables 7 and 8 that, at the electron thermal energy in the photosphere, which is on the order of 1 eV maximum, the cross-section of the 1 -2 and 1 -2 excitation is smaller than the elastic collision cross-sections considered above, so that their effect is also completely negligible.

Conclusion about the Electric Conductivity.
As a result, the contribution to the electric conductivity is essentially the one of the electron-ion Coulomb interaction. The conductivity in the photosphere is reported in Table 1. It can be seen in Figure 2 that the conductivity behavior follows exactly the one of the temperature as suggested by Meyer-Vernet [44]. The conductivity is found decreasing with height in the photosphere.

The Induction Equation under the Differential Form.
Because the photosphere paramagnetism is negligible, the induction equation may be rewritten as where is the magnetic diffusivity and introducing explicitly thedependence of these quantities. Returning to the Ohm's law equation (6) and assuming for the moment that ⃗ = 0, the density current is Besides, under the assumption ⃗ = 0, and at the stationary state ⃗ / = 0, (18) may be integrated as Alternatively, this equation is obtained by equating the current that creates the magnetic field following the Maxwell-Ampère's law given in (5), with the current created by the fields following the Ohm's law given in (6). However, for investigating div ⃗ , we found it preferable to use the integral form of this equation that we will introduce below.

The Induction Equation under the Integral Form.
The magnetic field ⃗ created at point ⃗ by the electric current described by the density current ⃗ at the various points ⃗ is given by the Biot and Savart's law The differential form of the Biot and Savart's law is the Ampère's law given in (5) which we will rewrite now as because we have obtained that in the photosphere = 0 .
This describes the magnetic field created by the electric current. But alternatively the electric current is simultaneously created by the fields combined with the fluid motion following the Ohm's law given in (6). Inserting the Ohm's law in the Biot and Savart's law, one obtains the integral form of the induction equation at the stationary state where we let be spatially variable, and we neglect the electric field ⃗ ( ⃗ ) for the moment. In the following we present a detailed calculation of div ⃗ from this equation. A simpler and more direct calculation exists anyway, but some details of the present calculation will be later useful.
The objective is to apply the ⃗ ∇ operator to (24). In order to determine how this operator has to be applied under the integral, we first remark that (22) is the sum of elementary convolution products of the form In order to evaluate the derivative, let us evaluate ( + d ).
One has by changing the variable = − d . One has then for the derivative ( ) which is that the derivative of the convolution product is the convolution (by the same function) of the derivative. Applying this result to (22), where it must be noted that the ⃗ ∇ operator has then to be applied to ⃗ ( ⃗ ) but not to Replacing now ⃗ ( ⃗ ) by its value from the Ohm's law 10 Using integral by parts of the form it can be shown that the contribution of (30) multiplied by ( ⃗ ) to (28) reduces to which is exactly balanced in (28) by the contribution of the first term of the right-hand side of (29). The first term of (31) is 0 at large distances, and the last term can be summed over the components into which is also 0. Finally (28) reduces to The result would be the same if ⃗ ( ⃗ ) was replaced by any quantity of the type ( ⃗ ) ⃗ ( ⃗ ), whatever ⃗ ( ⃗ ) be. This result shows that even with a variable conductivity, the Biot and Savart's law leads to div ⃗ = 0. This result can be obtained more simply and directly from but in the following we investigate the limitation of the Biot and Savart's integral by the Debye shielding and we will use some intermediate results of the calculation presented above.

Effect of the Debye Shielding.
In a plasma, it is well known that a given charge is shielded by the surrounding ones, so that at a certain distance called the Debye length, the charge has no effect. Following, for instance, Meyer-Vernet [50], "[in a plasma] any charge has a dress of size the Debye length which makes it 'invisible' from larger distances". It is also well known that the Debye shielding applies to the calculation of the electric field created by the charge, and the Debye shielding in itself is an effect of the electrostatic potential created by each charge. But what about the magnetic field?
Consider the electrostatic potential created by an elementary charge located at origin This charge repels the charges of same sign and attracts the charges of opposite sign, so that at distances > , the total potential vanishes. Consider now that this charge has a velocity ⃗ V. It then creates the magnetic vector potential provided that V ≪ , the speed of light. The form of the vector potential is quite the same as the one of the electrostatic potential, in particular as for the -dependence. The velocity ⃗ V to be considered here is not the thermal velocity, but the fluid or current velocity. It can be seen in Table 1 that the order of magnitude of the Debye length is the micron in the photosphere. Thus it can be assumed that the fluid or current velocity ⃗ V is constant over the Debye volume, so that the velocity can be factorized in the screening formalism. The form of the screening becomes thus the same as for the electrostatic potential.
More precisely, the screened electrostatic potential of the charge is where ( ) and ( ) are, respectively, the electron and ion density excess or depletion in due to the presence of the charge at origin. Here is the absolute value of the electron charge. Following Meyer-Vernet [50], this reduces to As for the vector potential, assuming that all the particles surrounding the charge undergo the same velocity ⃗ V, the total vector potential in is which analogously results in This shows that the magnetic vector potential is also screened, provided that the fluid or current velocity is the same in all the volume under interest. As a consequence, the Debye shielding applies also to the magnetic field in the photosphere, and the Biot and Savart's law of (28) has indeed to be limited to the Debye volume. Indeed, the particle velocity ⃗ V is the sum of the thermal velocity ⃗ V and of the fluid velocity ⃗ V The magnetic vector potential is then the sum of both thermal and fluid contributions, for each particle. As for the potential created by the particle under interested and recalled in (37), the thermal contribution is compensated for by the similar contributions of other particles having different and random thermal velocities. The thermal velocity is not responsible for a magnetic field. As for the screening particles whose contribution to the vector potential is given by the triple integral in (40), the thermal contribution has also to be examined. The thermal velocity being different for each elementary particle of the screen, the thermal velocity does not factorize out of the integral, but on the contrary, the random character of the thermal velocity makes the corresponding triple integral vanishing. Thus, there is no shielding of the vector potential under the effect of thermal velocities. But as for the fluid velocity, it can be considered as constant over the volume integral in its main contributing region, which is on the order of the Debye sphere whose typical radius in the photosphere is 1-2 microns (see Table 1), so that it can be factorized as done in (40). The fluid velocities are then responsible for shielding the magnetic vector potential. The typical convolution product as in (25) entering the Biot and Savart's law becomes With the same variable change as before, it can be shown that the derivative of ( ) with respect do is given by In the detailed calculation of div ⃗ , the only difference is that the first term of (31) has to be reconsidered because its limits are now − and + . By summing over all the terms in and taking the triple integral, it can be shown that the contribution of the terms of this kind is of the type so that by summing over the three components, the contribution of these terms is finally which represents the flux through the Debye volume surface of the elementary field created at the point of interest ⃗ by the current located on this surface. If the Debye volume is a sphere as usual, ⃗ − ⃗ and d ⃗ are parallel so that ⃗ ∇ ⋅ ⃗ ( ⃗ ) = 0 whatever the current is. But if the Debye volume is not a sphere, ⃗ ∇ ⋅ ⃗ ( ⃗ ) may not be zero.

Discussion
We have shown above that div ⃗ may not be zero when the Debye volume is not a sphere. We think that this is realized in the photosphere, because, as the surface of a star, the photosphere is a stratified medium. Such a medium may be characterized by an aspect ratio (see for instance [42]), which is the ratio between the vertical characteristic length and the horizontal one. What are these lengths in the photosphere? The horizontal characteristic length could be the granule typical diameter, 1000 km. The vertical characteristic length could be the pressure scale height, 110 km [51]. The aspect ratio in the photosphere results in being = 1/9 = 0.11. Let us cite other aspect ratios: 0.1-0.01 in the deep ocean, 0.01-0.001 in the terrestrial atmosphere close to the ground and in the thermocline [52], which is the zone of rapid thermic transition between surface and deep water in the oceans. This is only a trivial remark. Returning to the photosphere case, we interestingly remark that this ratio is the same as the one observed between the magnetic field gradients.
The degree of stratification of the medium can be evaluated by computing the horizontal Froude number. Referring to Brethouwer et al. [42] and adopting their notations of for the typical horizontal velocity, ℎ and V for the characteristic horizontal and vertical lengths, respectively, so that we introduce the Brunt-Väisälä frequency , which is the oscillation frequency of a particle undergoing a vertical motion in the fluid where the gravity acceleration is and the density is decreasing with height, so that The horizontal Froude number is then given by and analogously for the vertical Froude number. In Table 2, we give the Brunt-Väisälä frequency evaluated in the photosphere following the VAL C atmosphere model and considering the gas of neutral hydrogen which is the most dense with respect to the free electrons and ions. Considering for the small scale velocity of the particles as given by the VAL C model and listed in Table 2, which is in good agreement with the values observed by Bommier et al. [53], and taking for ℎ the typical granule diameter stated above, we derive the horizontal Froude number given in Table 2. It can be seen that we have ℎ ≪ 1 in the photosphere, so that the medium may be considered as strongly stratified.
On the basis of direct numerical simulations of stratified flows in which the vertical characteristic length V is let as a free parameter, Brethouwer et al. [42] show that two regimes may be distinguished following the value of the buoyancy Reynolds number given by where is the Reynolds number of the fluid. We have evaluated this number for the gas of neutral hydrogen of the photosphere where ] HH is the collision frequency of the gas particles, which are neutral hydrogen atoms. Following Meyer-Vernet [44] we simply adopt the billiard ball model for evaluating ] HH for the gas of neutral hydrogen, so that we consider the cross-section HH = 10 −19 m 2 . For the typical length entering , we have taken the horizontal length scale ℎ to be coherent with the horizontal Froude number. We find ≫ 1 in the photosphere, but ≪ 1, a regime that Brethouwer et al. [42] depict as the "viscosity-affected stratified flow regime, " more precisely they say that "if ≪ 1 the large-scale dynamics is determined by a balance between inertial and viscous forces due to vertical shearing. Therefore no inertial cascade can develop and the dissipation occurs predominantly at the largest scales. . . Kelvin-Helmholtz type of disturbances cannot be present when ≪ 1". This result was already derived and experimentally confirmed by Godoy-Diana et al. [54]. Such a conclusion can interestingly be brought together with the result obtained by Bommier [55] who has shown that the structure of the quiet sun magnetic field as derived from recent inversions is not turbulent as generally admitted, but organized in connected vertical opening flux tubes.
Moreover, Brethouwer et al. [42] show that in the case of the viscosity-affected stratified flow regime, the vertical characteristic length V that was let free in the simulation results is so that the theoretical value of the aspect ratio th is We have computed th for the photosphere and it can be seen in Table 2 that we obtain a theoretical value in very good agreement with the empirical one of 1/9 that we determined from the granule typical diameter and the pressure scale height, and that we observed in div ⃗ . Concerning the anisotropy of the velocities, Brethouwer et al. [42] show that when the horizontal velocity scales as , the vertical one scales as 2 ℎ / , so that the theoretical aspect ratio for the velocities is We give the inverse of this number in Table 2, where it can be seen that the ratio horizontal/vertical is even higher for the velocities than for the typical lengths, in the photosphere. It can then be concluded that the velocities are highly anisotropic, in this ratio. Meyer-Vernet [50] has recomputed the Debye length discarding the usual hypothesis of the thermodynamical equilibrium, and basing the calculation on the charge velocity. The model assumes isotropic velocities, but it can be generalized to anisotropic ones, by performing an affinity along one of the reference axes, which leads to different Debye lengths ℎ and V for the different vertical and horizontal typical velocities V ℎ and V V : so long as the hydrogen atom velocity ratio is transferred to the charges by their collisions between themselves. Indeed, the particle velocity is the sum of the thermal velocity and of the fluid velocity, as in (42). The anisotropy described above is the one of the fluid velocity. The thermal velocity is isotropically distributed, but the total velocity remains anisotropic, leading to the different horizontal and vertical Debye lengths derived above. As a result the Debye volume is not a sphere and div ⃗ may mathematically depart from zero. Thus div ⃗ is not zero as a consequence of the stratification, as observed. The aspect ratio is transferred to div ⃗ .
The value of the isotropical and equilibrium Debye length for the VAL C model is given in Table 1 as an order of magnitude of Debye lengths.
In Table 1 we give also the magnetic Reynolds number, computed by taking as typical length the horizontal length ℎ in coherence with the Froude and Reynolds number calculation. The magnetic Reynolds number makes intervene the Coulombian electron-ion collisions via the electric conductivity (see (14)), whereas the Reynolds number is concerned by the collisions between the neutral hydrogen atoms (see (51)). The photosphere plasma is a complex medium, having the neutral particle density four orders of magnitude higher than the charged particle density. All these particle collide, the neutral particles with the neutral particles, the charged particles with the charged particles,  but also the charged particles with the neutral ones. All the collision frequencies are given in Tables 1 and 2. Due to the ensemble of collisions, the fluid properties should be common for the neutral and charged particles. The relative effect of the magnetic forces is scaled by the magnetic Prandtl number, which is simply the ratio between the magnetic Reynolds number and the Reynolds number = / , but is also the ratio between the kinematic viscosity and the magnetic diffusivity. We give the value of the magnetic Prandtl number in the photosphere in Table 2. It can be seen that this number is large with respect to unity. This means that the dynamical viscosity dominates the magnetic diffusivity, so that the photosphere can effectively be considered as a highly stratified medium, even in the presence of magnetic field. In Table 2 we give also the plasma coefficient, which is the ratio between the gas pressure and the magnetic pressure: We have computed this number assuming a magnetic field = 1000 G, as a typical value for sunspot umbra and penumbra. We obtain that is of the order of unity, which also shows that the medium structuration is not dominated by the magnetic field and that the stratification acts besides.

A Proposal for MHD Modeling in the Photosphere.
It could be considered that in reduced spatial scales, where the horizontal scales are multiplied by , the Debye volume becomes spherical so that div ⃗ ≈ 0. This is true for the observations reported at the beginning of this paper, having the vertical gradient as 3-4 G/km and the horizontal one as 0.4-0.5 G/km. If the vertical gradient is multiplied by , it results in 0.33-0.44 G/km, of the same order of magnitude as the horizontal gradient so that div ⃗ ≈ 0. We have applied this aspect ratio scaling to our THEMIS observations of NOAA 10808. The result is reported in Figure 3(b), where it can be seen that this rescaled div ⃗ resembles now to noise, which is in other words div ⃗ ≈ 0.
Thus we are led to suggest that the usual methods of the MHD could nevertheless be applied for calculations in a stratified medium as the photosphere, by considering that (i) the div operators have to be computed in the reduced length scale; (ii) the curl operators have to be computed in the normal length scale. The magnetic field remains in any case given by the Biot and Savart's law (limited to the Debye volume) that reduces locally to Maxwell-Ampère's law of (23) where the curl is in the normal scale.
As an example, one can consider the magnetostatic sunspot model by Pizzo [35]. In this paper, the prototype spot model has a diameter of about 2000 km, whereas the average diameter is 10,000 km in the observations [51]. If the reduced length scale had been used in the computation, the spot model diameter would be 1/ times larger so that it would become of the same order of magnitude as the average observed spot diameter.

Conclusion
We have shown that, without any exception, the observed value of the so-called vertical gradient of the sunspot photospheric magnetic field |Δ /Δ | is 3-4 G/km, with field strength always decreasing with height, whereas the horizontal gradient |Δ /Δ + Δ /Δ | is 0.4-0.5 G/km. These observations involve various lines, visible or IR, various instruments, ground-based or spaceborn, various interpretation methods. In some cases the same inversion code performs the magnetic field and depth difference determination, in other cases they result from two different codes. The fact that no exception was found in spite of the variety of lines, instruments, and methods argues in favor of a true difference and not a result of noise. Also the fact that it is always a flux loss and never a flux gain that is observed with increasing height. Moreover their ratio is on a factor of 9, which is rather high. As, also, the difference is large and nearly as large as the largest gradient, ascribing the difference to noise would imply that all the measurements would be not larger than noise. This is not compatible with the fact that the difference is always of the same sign.
This difference implies that div ⃗ derived from these two gradients would not be zero. This is usually ascribed to unresolved magnetic field structures. But we are not convinced, because if each elementary unresolved structure has div ⃗ = 0, we do not see how the field divergence of the whole region could be nonzero, because the divergence operation is linear. In other words, how the field could be flux non-conservative at large scale together with resulting of the combination of smaller flux conservative structures?
On the theoretical side of the problem, div ⃗ = 0 is a direct consequence of the Biot and Savart's law. In this paper, we propose to solve the discrepancy by considering a plasma effect, the Debye shielding of the charges. Without challenging the Biot and Savart's law, we remark that in the plasma the charge effect is shielded so that the integral of the Biot and Savart's law has to be limited to the Debye volume. But when the volume is a sphere, we obtain that nevertheless div ⃗ = 0. This is only when the volume is not a sphere that div ⃗ may depart from zero. We show that this could be the case of the solar photosphere, because this is the surface layer of a star, and thus it is a stratified medium. We verify this point by computing the horizontal Froude number. We evaluate its aspect ratio as 1/9 by considering the granule typical radius as the horizontal characteristic length, and the pressure scale height as the vertical characteristic length. Their ratio is 1/9, interestingly the same as the one of the observed horizontal and vertical magnetic field gradients. This observed value is interestingly in very good agreement with the theoretical one derived following Brethouwer et al. [42]. Following these authors we obtain also that the photosphere ranges in the "viscosity-affected stratified flow regime, " where the inertial cascade cannot develop. As a result of the stratification, the horizontal velocity is highly different from the vertical one. It is possible to generalize the result by Meyer-Vernet [50] that the Debye length is proportional to the charge velocity, to such anisotropic velocities by performing an affinity of one of the reference axes. Thus, the Debye volume is a flattened sphere and div ⃗ may depart from zero in the photosphere, as observed.
A magnetic monopole can be the source of nonzero div ⃗ . But magnetic monopoles have not been observed. The logical assertion is that the presence of a monopole implies nonzero div ⃗ . The negative assertion is that div ⃗ = 0 implies the absence of any magnetic monopole. But the negative assertion is not that the absence of detection of any monopole would imply div ⃗ = 0, as generally believed. In other words, the presence of a magnetic monopole is a sufficient condition, but not a necessary condition, for div ⃗ ̸ = 0. Other sufficient conditions may exist. Anisotropic shielding constitutes an alternative to the existence of monopoles for being responsible for nonzero div ⃗ in a medium.
It can be concluded that such an effect of div ⃗ departing from zero is very specific of the photosphere, because it is a combined plasma and stratification effect.
Finally, this leads us to introduce a spatial scaling by the aspect ratio, where div ⃗ ≈ 0. We then propose this scaling method to nevertheless perform MHD modeling in the photosphere, where the div operators have to be computed in the reduced length scale, and the curl operators have to be computed in the normal length scale.
We have also evaluated the Froude number at the basis of the solar corona and found 3.6, increasing with height, which shows that the corona is conversely not at all a stratified plasma, so that the conditions are very different there.