Magnetohydrodynamic calculation of the temperature and wind velocity profile of the solar transition region. Preliminary results

The sharp almost step like increase the temperature in the transition region (TR) between chromosphere and solar corona is well-known from decades; for first time we are giving a detailed magnetohydrodynamic (MHD) calculation of the height dependence of the temperature. The width of the transition region is evaluated by maximal value of the logarithmic derivative of the temperature. At fixed heating, only MHD can give such a narrow width and in such sense, even the qualitative agreement with the observational data, gives the final verdict what the heating mechanism of the solar corona is. Static profiles of the temperature and wind velocity are calculated for static frequency dependent spectral density of the incoming MHD waves; no time dependent computer simulations. At fixed spectral density of MHD waves, the MHD calculation predicts height dependence of the non-thermal broadening of spectral lines and its angular dependence. For illustration is used one dimensional approximation of completely ionized hydrogen plasma in weak magnetic field, but it is considered that the width of the TR is weakly dependent with respect of further elaboration. The analyzed MHD calculation is a numerical confirmation of the qualitative concept of self-induced opacity of the plasma with respect to MHD waves. The plasma viscosity strongly increases with the temperature. Heated by MHD waves, plasma increases the wave absorption and this positive feedback leads to further heating. The static temperature profile is a result of a self-consistent calculation of propagation of MHD wave through the static background of wind and temperature profile. The numerical method allows consideration of incoming MHD waves with an arbitrary spectral density. Further elaboration of the method are briefly discussed: influence of second viscosity in the chromospheric part of the TR, influence of the magnetic field on the coronal side of the TR and investigation of such type effects on the width of the TR. 1 Alfvén model for corona heating The discovery of the 5303 Å green1 of Fe13+ and 6374 Å red2 of Fe9+ magneto-dipole transition lines in the solar corona spectrum [1] posed an important problem for the fundamental physics – what is the mechanism of the heating of the solar corona and why the temperature of the corona is 100 times e-mail: mishonov@gmail.com 1 Fe XIV, P3/2 → P1/2, 3s23p1 2 Fe X, P1/2 → P3/2, 3s23p5 © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). MATEC Web of Conferences 145, 03009 (2018) https://doi.org/10.1051/matecconf/201814503009 NCTAM 2017

Abstract. The sharp almost step like increase the temperature in the transition region (TR) between chromosphere and solar corona is well-known from decades; for first time we are giving a detailed magnetohydrodynamic (MHD) calculation of the height dependence of the temperature. The width of the transition region is evaluated by maximal value of the logarithmic derivative of the temperature. At fixed heating, only MHD can give such a narrow width and in such sense, even the qualitative agreement with the observational data, gives the final verdict what the heating mechanism of the solar corona is. Static profiles of the temperature and wind velocity are calculated for static frequency dependent spectral density of the incoming MHD waves; no time dependent computer simulations. At fixed spectral density of MHD waves, the MHD calculation predicts height dependence of the non-thermal broadening of spectral lines and its angular dependence. For illustration is used one dimensional approximation of completely ionized hydrogen plasma in weak magnetic field, but it is considered that the width of the TR is weakly dependent with respect of further elaboration. The analyzed MHD calculation is a numerical confirmation of the qualitative concept of self-induced opacity of the plasma with respect to MHD waves. The plasma viscosity strongly increases with the temperature. Heated by MHD waves, plasma increases the wave absorption and this positive feedback leads to further heating. The static temperature profile is a result of a self-consistent calculation of propagation of MHD wave through the static background of wind and temperature profile. The numerical method allows consideration of incoming MHD waves with an arbitrary spectral density. Further elaboration of the method are briefly discussed: influence of second viscosity in the chromospheric part of the TR, influence of the magnetic field on the coronal side of the TR and investigation of such type effects on the width of the TR. larger than the temperature of the photosphere. Those spectroscopic experiments put perhaps the most living unresolved problem in physics. None of the heating mechanisms has been confirmed and none has received an indisputable refutation. The 75 years old problem is still fresh and green.
The first idea by Alfvén [2] was that now called Alfvén waves (AW) [3] are the mechanism for heating the corona. Perhaps AW are generated by the turbulence in the convection zone and propagate along the magnetic field lines as through wave guides. Absorption is proportional to ω 2 and the heating comes from high-frequency AW which are absorbed and not observed in the corona. Alfvén's idea for the viscous heating of plasma by absorption of AW was analyzed in the theoretical work by Heyvaerts [4]. In support of this idea is the work by Chitta [5] (Figures 8, 9 therein). The authors came to the conclusion that the spectral density of AW satisfies a power law with an index of 1.59. This gives a strong hint that this scaling can be extrapolated in the nearest spectral range for times less than 1 s and frequencies in the Hz range. Furthermore in the work by Tomczyk [6] it is stated that there exist very few direct measurements of the strength and orientation of coronal magnetic fields, meaning that the mechanisms responsible for heating the corona, driving the solar wind, and initiating coronal mass ejections remain poorly understood. After the launch of Hinode, Alfvènic type MHD modes were observed [7] and the well-forgotten spatially and temporally ubiquitous waves in the solar corona [8] came again into the limelight and gave strong support for the idea of Alfvén and the frequency of the notion AW as a key worth significantly increased.
A clear presence of outward and inward propagating waves in the corona was noted. k − ω diagnostics revealed coronal wave power spectrum with an exponent of ≈ − 3 2 (cf. Fig. 2 of [6]). The low frequency AW, on the other hand, reach the Earth orbit and thanks to the magnetometers on the various satellites we "hear" the basses of the great symphony of solar turbulence.
The observational data for the temperature profile of the solar corona show that the TR is extremely thin compared to the radius of the Sun [9][10][11][12][13]. All those observational data for the quiet-Sun chromosphere and transition region reveal one dimensional and time independent profiles which theoretically have to be derived by time averaged dynamic consideration. The width of the TR λ might be evaluated by the maximum of logarithmic derivative of the temperature, λ = max( dT T dx ). In order to qualitatively explain this small width, in [14] the idea of self-induced opacity of the plasma for Alfvén waves was introduced. A similar idea was analysed also by Suzuki [15] (see also [16][17][18][19][20][21]), however none of these articles is helpful for comparison with our manuscript because the width of the TR is not calculated in details in any of them. In the current paper we give a numerical realization of this idea by calculating for the first time the width of the TR using the framework of MHD.
First detailed observations by Skylab half a century ago revealed surprisingly that TR is very thin λ-equal in width to metropolitan Los Angeles [10]. As λ is much smaller that the solar radius we have actually a plane one dimensional (1D) problem. In 1D the constant fluxes of mass, energy and momentum give 3 integrals of motion which significantly alleviate the calculations. We can eliminate the density ρ(x) as dynamic variable and for the static profiles of temperature T (x) and solar wind U(x) we have to solve numerically first order differential equations.
Without a doubt the kinetic approach is indispensable for the treatment of low density solar corona but this problem is beyond the purpose of the present work. In our work we take into account the influence of the magnetic field on the viscosity within the hydrodynamic approach.
The purpose of the present work is to examine whether the initial Alfvén idea is correct and to solve the MHD equations which give the dependence of the temperature on the height T (x) and the related velocity of the solar wind U(x) supposing static spectral density of the incoming AW. We illustrate Alfvén's idea by an MHD calculation for completely ionized hydrogen plasma in the beginning of TR where the density of neutral Hydrogen atoms is already negligible. Due to the high density of the TR, MHD is an adequate tool to analyze the beginning of the process.
Our starting point are the MHD equations [22] for the velocity field v and magnetic field B where is the energy density flux, ρ is the mass density, ε is the internal energy density, κ is the thermal conductivity, h is the enthalpy per unit mass; is the Poynting vector and ν m ≡ c 2 ε 0 is the magnetic diffusivity determined by Ohmic resistance and vacuum susceptibility ε 0 ; vacuum permeability is µ 0 . For hot enough plasma is negligible. The total momentum flux is a sum of the inviscid part ρvv + P of the fluid, with pressure P and the viscous part of the stress tensor, with viscosity η and second viscosity ζ, and lastly, the Maxwell stress tensor with δ ik the Kronecker delta. For temperatures above 10 kK we treat coronal plasma as completely ionized hydrogen plasma with the following temperature dependence of kinetic coefficients [23]: where q e is the electron charge, m is the mass of electron, M is the proton mass, T is the temperature and n tot = n e + n p is the total density of electrons and protons; ρ = Mn p . We suppose that µ 0 = 4π and ε 0 = 1/4π, but in the practical system all formulae are the same; as well as in Heaviside-Lorentz units where µ 0 = 1 and ε 0 = 1. As we mentioned above i.e. the hot hydrogen plasma is sticky, dilute, and "superconducting". Here ν k is the kinematic viscosity. We can introduce the magnetic Prandtl number as the ratio between the kinematic and the magnetic viscosity, The coronal heating mechanism can be revealed without additional accessories. Let us mention also the relations κ = 1.5 T/q 2 e leading to and η/κ ≈ 4 9 √ mM. The applicability of our equations is governed by an additional condition: If in the dense plasma the frequency of the proton-proton collisions ν pp is much smaller than the proton cyclotron frequency ω c , the magnetic field can be considered a perturbation and for the kinetic coefficients we can use the zero magnetic field formulas. An elementary gas kinetic consideration of main notions of plasma as proton-proton collision frequency ν pp mean free time τ pp , thermal velocity v T p = √ T/M mean free path l p = v T p τ pp and other is given in the the appendices. Lastly, we need to mention the existence of radiative losses. However, the thin TR has width much smaller not only than Earth's orbit radius but than Earth's radius and thus these losses are neglectable when considering the mechanism for coronal heating; continuum radiation may be neglected for all temperatures fewer than a few million degrees [24]. The hot corona exists in broad ranges where the heating mechanism cannot be effective. For this reason, in the narrow range of the transition region the radiation losses are negligible compared to the intensive heating, no matter what the concrete mechanism is. Quantitatively, this means that radiation power per unit volume P rad nT U/λ, where n and T are the number density and temperature of the corona, λ is the width of the transition region, and U is the velocity of the solar wind. The influence of gravitational field with acceleration g is also negligible for the temperature gradient in the TR Mg λ T ; cf. Ref. [20]

MHD equations and energy fluxes
First of all, we have to emphasize than MHD is not a model but adequate theory applicable for the dense plasma of TR. Let us now recall the main MHD equations.
The time derivative ∂ t B which implicitly participates in the energy conservation Eq. (2) obeys the equation cf. Ref. [22] Analogously the momentum equation Eq. (3) can be rewritten by the substantial derivative We analyze AW propagating along magnetic field lines B 0 . We focus our attention on the narrow TR, where the static magnetic field is almost homogeneous and the waves are within acceptable accuracy one dimensional. For the velocity and magnetic fields we assume with vertical homogeneous magnetic field B 0 e x perpendicular to the surface of the Sun. The transverse wave amplitudes of the velocity u(t, x) and magnetic field b(t, x) we represent with the Fourier integrals For illustrative purposes it is convenient to consider monochromatic AW with u(t, x) =û(x)e −iωt and b(t, x) = B 0b (x)e −iωt . And later on to restore the summation on different frequencies which finally gives Fourier integration.

Wave equations
For transversal waves the general MHD equations Eq. (15) and Eq. (14) give the following exact system of linear forû(x) andb(x) where is the Alfvén velocity. For illustrative purpose and complete lack of observational data in the present work we analyze only the case of zero transversal to the magnetic field wave vector k y = 0. This system can be Hamiltonized. In our numerical analysis we solve the first order linear system of equations whereγ ≡ d xb ,ŵ ≡ d xû and For homogeneous medium with constant η, ν m , ρ, V A , and U, in short for constant wave-vector matrix K, the exponential substitution Ψ ∝ exp(ikx) in Eq. (22) or equivalently Eq. (19) and Eq. (20) gives the secular equation det (K − k1) = 0, which after some algebra gives the dispersion equation of AW where ω D ≡ ω − kU is the Doppler shifted frequency. Usually either viscosity coefficient is significantly bigger than the other. For example, for hot plasma ν k ν m and with an acceptable approximation we can neglect the smaller one and to work with matrix 3×3. Approximation of ideal fluid is reduced to both negligible viscosity coefficients and corresponding 2×2 matrix gives the wellknown AW dispersion in moving fluid ω 2 D = V 2 A k 2 . A perturbative treatment of dispersion equation ω(k) = ω + iω Eq. (23) gives [22] ω = − 1 2 (ν m + ν k )k 2 . If WKB approximation is applicable for wave amplitude damping we have to use the eigenvalue of K matrix giving smallest damping for right running wave.
In order to make model evaluation of the influence of magnetic field on viscosity for the waves we can use the the Braginskii result [25] ν k,2 = 6 5x wherex ≡ ω Be τ pp .

Wind variables
We solve the wave equation Eq. (22) from some height at which hydrogen plasma is approximately completely ionized x = 0 to some distance large enough x = l, where the short wavelength AW are almost absorbed. As we pointed out, this distance is much bigger than the width of the transition region λ, but much smaller than solar radius. That is why the problem is effectively one dimensional. The horizontal wavelength of the waves responsible for heating is much smaller than the range of area for which we are calculating the temperature profile. The considered one-dimensional 0 < x < l time independent problem has three integrals corresponding to the three conservation laws related to mass, energy and momentum. The mass conservation law Eq. (1) gives the constant flow where ρ 0 = ρ(0), ρ l = ρ(l), U 0 = U(0), and U l = U(l). The energy conservation law reduces to a constant flux along the x-axis Here the first term describes the energy of the ideal wind, i.e. a wind from an ideal (inviscid) fluid. The second termq(x) includes all other energy fluxes; in our notations tilde will denote sum of the non-ideal (dissipative) terms of the wind and wave terms. In detail the non-ideal part of the energy flux q(x) consists of: the wave kinetic energy ∝ |û| 2 , viscosity (wind ∝ 4 3 η + ζ and wave ∝ η components), heat conductivity ∝ κ, and Poynting vector ∝b * , where ξ ≡ 4 3 η + ζ. In order to alleviate the detail analysis depicted in Fig. 3, we have introduced convenient notations for all terms with different nature. Here time averaged energy flux is represented by the amplitudes of the monochromatic oscillations, this is a standard procedure for alternating current processes. In our case we have, for example, The other terms from Eq. (4) are averaged in a similar way like in the equation above. In short, we substituted the formulas for AW Eq. (16) propagating along the magnetic field in x-direction in the general formula for the energy flux Eq. (4) and made the time averaging and summation over waves with different frequencies. We will repeat the same approach to the momentum flux Eq. (6) For constant spectrum of incoming waves the momentum conservation law Eq. (6) gives the sum of the ideal wind fluid and the other terms which take into account the wave part of the Maxwell stress tensor ∝ |b| 2 and viscosity of the wind ∝ ξ. The sum here and in Eq. (27) is actually integral over the spectral density W(ω) of incoming AW Let us mention that all components of the wave Ψ(ω, x) have frequency dependence. One can suppose that the spectral density W of the waves coming from the chromosphere has the same power law dependence, i.e.
The parameter α is between 1.5 and 2: 3 2 [6], 1.59 [5], 2 [26]. C is the unknown parameter of the theory, which we vary for fixed ξ in order to reproduce the temperature increase in the TR. In our earlier works [27] simulations of plasma heating with power law spectral density is given with up to 8 AW with different frequencies.
We have to solve the hydrodynamic problem for calculation of wind velocity and temperature at known energy and momentum fluxes. The problem is formally reduced to analogous one for a jet engine, cf. Ref. [28]. We approximate the corona as completely ionized hydrogen plasma, i.e. electrically neutral mixture of electrons and protons. The experimental data tells us that proton temperature T p is higher than electron one T e . This is an important hint that heating goes through the viscosity determined mainly by protons. However for illustration purpose and simplicity we assume proton and electron temperatures to be equal T e = T p = T.
where, as we mentioned earlier, h is the enthalpy per unit mass and ε is the density of internal energy. Although there are some hints for different values of the adiabatic index γ [29], we will use the traditional value of 5/3 for our calculations since this choice will not change the essence of our presentation. Hydrogen plasma is well theoretically investigated and the 5/3 value of gamma is used in hundreds of works. There are no theoretical hints for significant deviations from this value obtained by ideal gas approximation. Actually a small variation of the value of gamma does not change qualitatively the temperature profile of the TR. In order to alleviate the final formulae we introduce two dimensionless variables χ and τ which represent the non-ideal part of the energy and momentum flux respectively Here we could include gravitational field, Bremsstrahlung or other accessories, which are negligible for the narrow TR. Analogously, for the wind velocity and temperature we have where U 0 = U(0). The energy and momentum constant fluxes Eq. (26) and Eq. (28) in the new notation take the form From the second equation we express the dimensionless temperature Θ and substitute in the first one. Solving the quadratic equation for the wind velocity U we derive where for the discriminant we have If χ = 0 and τ = 0, we get the initial condition U| χ=0,τ=0 = 1. This condition determines the sign in front of √ D in Eq. (36). Here γ is the constant ratio of the heat capacities, and s ≡ c s (0)/U 0 is the ratio of the sound and wind velocity at x = 0. We suppose that initial wind velocity is very small U(0) c s (0). The velocity distribution Eq. (36) can be substituted in Eq. (35) and we derive the dimensionless equation for the temperature distribution expressed by the wind velocity profile The solutions for velocity U(x) Eq. (36) and temperature T (x) Eq. (38) distributions are important ingredients in our analysis and derivation of the self-consistent picture of the solar wind. We use a one dimensional approximation and in addition the constant flux of mass, energy and momentum give 3 integrals of motion. This enables us to solve exactly the nonlinear part the MHD equations for the solar corona heating problem analytically. As the energy and momentum flux density equations, Eq. (26) and Eq. (28) respectively, contain space derivatives, these nonlinear equations are differential equations, in which the variables are explicitly expressed but not the derivatives, as it is usual. The intuitive interpretation is unnecessary, but nevertheless rewriting Eq. (34) and Eq. (35) as one can say that χ(x) is the energy given from all non-ideal terms to the initial energy flux of the ideal wind 1 2 + c p Θ 0 in order ideal wind flux to reach his value 1 2 U 2 + c p Θ at coordinate x. In the same dimensionless units ρ 0 = 1 and U 0 = 1 one can consider the solar wind as a wind of ideal fluid with initial momentum density 1 + Θ 0 in x = 0. In the space interval (0, x) this ideal wind obtains from non-ideal processes momentum flux τ(x) and finally the momentum flux at coordinate x is U + Θ/U. Such a chain of notions could be useful to trace influence of every process to the wind velocity U(x) and temperature T (x) profiles.

Boundary conditions for the waves
At known background wind variables U(x) and T (x) we can solve the wave equation Eq. (22) for run-away AW at x = l. As we will see later the run-away boundary condition Eq. (55) corresponds to right propagating AW at the right boundary of the interval. The wave equation Eq. (22) is extremely stiff at small viscosity, and numerical solution is possible to be obtained only downstream from x = 0 to x = l. We have to find the linear combination of left and right propagating waves at x = 0, which gives the run-away condition at x = l. The solution of wave equation according to Eq. (26) determines the energy flux related to the propagation of AW where Here j-term represents kinetic energy of the wave, η-term comes from the viscous part of the wave energy flux, and B 0 -terms describe the Poynting vector of the wave. For hot plasmas the Ohmic dissipation is negligible and we can use the approximation ν m = 0. In this case it is not necessary to use dynamic variableγ and we can work with 3×3 matrix. For zero viscosity and Ohmic ressisitivity, equations for AW amplitudes can be reduced to 1D Schrödinger equation [30]. For the illustration of the wave boundary values we will use this approach.
In order to take into account the boundary condition at x = l, we calculate the eigenvectors of the matrix K, which according to Eq. (22) determine the wave propagation in a homogeneous fluid with amplitude ∝ exp(ikx). Then the eigenvalues of K give the complex wave-vectors The three eigenvectors L, D and R are ordered by spatial decrements of their eigenvalues and are normalized by the conditions where the sign corresponds to the direction of wave propagation. Notation L corresponds to left propagating wave, R to right propagating wave, and D for an overdamped at small viscosity mode. For technical purposes we introduce the matrix notations For low enough frequencies ω → 0 and wind velocities the modes describe: 1) right-propagating AW with k R ≈ ω/V A and small k R ≈ ν k ω 2 /2V 3 A k R , 2) left propagating wave k L = −k R , and a diffusion overdamped mode k D ≈ V 2 A /ν k U k D which describes the drag of a static perturbation by the slow wind U V A in a fluid with small viscosity. In this low frequency and long wavelength limit the stiffness ratio of the eigenvalues is very large As we emphasized the wave equations Eq. (22) form a very stiff system and indispensably has to be solved downstream from the chromosphere x = 0 to the corona x = l using algorithms for stiff systems. Let We look for a solution as a linear combination in other words we suppose that from the low viscosity chromosphere plasma do not come overdamped diffusion modes. The strong decay rate make them negligible at x = 0. Physically this means that AW (R-modes) are coming from the Sun and some of them are reflected from the TR (L-modes) Analogously for the configuration of open corona we have to take into account the run-away boundary condition for which we suppose zero amplitude for the wave coming from infinity Written by components This boundary condition gives or by components These systems give the amplitudes of the reflected wave r and transmitted wavet in the solution Eq. (51). For this solution we have the energy flux of transmitted T and reflected R waves T ≡ ψ † (l) g(l) ψ(l) = |t| 2 + |c| 2 + tc * D † (l) g(l) R(l) + c.c. , Then we introduce the absorption coefficient The described solution is normalized by unit energy flux of the R-wave. If we wish to fix energy flux of the right propagating wave to be q wave (0) we have to make the renormalization using a parameter A wave . In this section we have described Absorbing Boundary Conditions (ABC) [31] well known from radar calculations, but realization for AW is more complicated and require eigenvector analysis. Now using Ψ(x) we can calculate the wave part of the energy flux Eq. (42) and the wave part of the momentum flux If we work with 4×4 matrix the generalization is obvious, in this case we have to operate with one additional diffusion mode.

Self-consistent procedure and results
We have to say few words about the numerical methods. For high frequency waves we can use WKB approximation, which for our problem reads As the equations for wind variables are stiff, we can use the implicit Euler method for T (x) and U(x) For the realization of this method we use the explicit formulas Eq. (36) and Eq. (38) for U and T respectively. As dissipative fluxes contains the derivatives −ξUd x U −κ d x T for the energy flux Eq. (27) andΠ ξ (x) ≡ −ξd x U in (29), the explicit expression for the wind variables Eq. (36) and Eq. (38) is actually a system of first order ordinary differential equations. In the WKB approximation complete system for wind and wave variables can be solved simultaneously.
One alternative is to use explicit expressions for derivatives d x U and d x T which age given from the constant energy Eq. (26) and momentum Eq. (28) fluxes. The limit case of ideal fluid with small ξ and κ in the denominators gives the singular perturbation theory. Care is therefore needed for the choice of the numerical method.
For exact wave variables the implicit Euler method applied for Eq. (22) gives For this case the WKB approximation is only a starting point. Then for fixed wind profiles T and U we solve the wave equation with the corresponding boundary conditions. And using wave fluxes of energy and momentum we calculate the next approximation for the wind variables T and U. The procedure has to be repeated to reach self-consistency of wind and wave variables. The simple Euler method can be elaborated by incorporation of predictor-corrector method using for small interpolation the Eitken formula for polynomial interpolation. Sequential polynomial interpolations can be rearranged as Padé interpolation using well-known epsilon algorithm by Wynn and Backer; some details are given in the Appendix II.
Following the described numerical method we have calculated the temperature and wind velocity profiles and the results are depicted in Fig. 1 and Fig. 2.
The energy and momentum flux density of the calculated temperature and wind profiles are respectively given in Fig. 3 and Fig. 4.
Finally, we present a comparison between theory and observation in Fig. 5: the calculated temperature profile from Fig. 1 is placed on the presented observational data from Ref. [13].

Discussions
In spite that iron was suspicious from the very beginning the problem of Coronium was a 70 year standing mystery until unambiguous identification as Fe 13+ by Grotrian and Edlen [33] in 1937. The same 70 year time quantum was repeated. In 1947 Alfvén [2] advocated the idea that absorption of AW is the mechanism of heating of solar corona. Unfortunately the idea by Swedish iconoclast [34] was never realized in its original form: what can be calculated, what is measured, what is explained MHD calculations qualitatively reproduce the well-known observational data [12,13,32]. This step-like behaviour confirms the scenario of self-induced opacity with respect to Alfvén waves. and what is predicted. That is why there is a calamity of ideas still on the arena, for a contemporary review see the SOHO proceedings ( [35]). From qualitative point of view, the narrow width of the TR λ = min |dx/d ln T (x)| is the main property which should be compared against the predictions of other scenarios. For example, in order for the nanoflare hypothesis to be vindicated [36], such reconnections are needed to explain the narrow width of the TR at the same boundary conditions of wind velocity and temperature. Moreover, electric fields of the reconnections heats mainly the electron component of the plasma. How then proton temperature in the corona is higher? Launching of Hinode gave a lot of hints for the existence of AW in the corona [37], see also Ref. [38]. However, most of the those research was in UV region when high frequency AW which heat are already absorbed. All observations are for low frequency (mHz range) AW for which hot corona is transparent. The best can be done is to extract low frequency behavior of the spectral density of AW and to extrapolate to higher frequencies responsible for heating. So observed AW are irrelevant for the heating. In order to identify AW responsible for the heating, it is necessary to investigate high frequency (1 Hz range) AW in the cold chromosphere using optical, not UV spectral lines. We are unaware whether such type of experiments are planned. One of the purposes of the present work is to focus the attention of experimentalists on the 1 Hz range AW in the chromosphere, which we predict on the basis of our MHD analysis. For such purposes we suggest a special attention to be paid to the behavior of oscillations [39] and sunspot waves [40] above active regions. Another possibility is provided by Doppler tomography [41] of Hα or Ca lines. Doppler tomography was successfully used for investigation of rotating objects, such as accretion disks [42] and solar tornados [43]. Here we wish to mark also the Doppler tomography by Coronal Multi-channel Polarimeter build by Tomczyk [6]. For investigation of AW by Doppler tomography we suggest development of frequency dependent Doppler tomography operating as a lock-in voltmeter.
The date from every space pixel should be multiplied by sin ωt and integrated for many wave periods. Finally one can observe time averaged distribution of the AW amplitude. Systematic investigation of such frequency dependent Doppler tomograms will reveal that Swedish iconoclast [34] is again right that AW heat the solar corona, after another 70 years of dramatic launching of vast variety of ideas. In this article we wish to emphasize that the concept of self-induced opacity with respect to MHD waves  (27). Height x dependence of the logarithms of the modulus of different components of the energy flux. Labels of curves correspond to: 0) the total constant energy flux density q, which is integral of motion of MHD equations, 1) the ideal wave energy flux term q ub containing both wave velocityû and magnetic field b, 2) the wind pressure energy flux term q P , 3) the ideal wave energy flux term q b ; for this "aether" term a magnetic energy ∝ |b| 2 is moving with the velocity of the wind U, 4) the ideal wave velocity energy flux term q u , 5) the ideal wind energy flux acceleration term q U , 6) the wind thermal conductivity dissipative energy flux term q κ , 7) the wave dissipative viscous energy flux term q η , 8) the wave dissipative Ohmic resistivity flux term q and 9) the wind viscous dissipative energy flux term q ξ . Both wind dissipative terms q ξ and q κ are negative and describe leeward energy flux. All other terms describe windward energy fluxes. Despite being small, q κ and q ξ contain space derivatives and determine the shape of temperature T (x) and wind velocity U(x) profiles. Eq. (29). Labels of curves correspond to: 0) the total constant momentum flux density Π, which is integral of motion, 1) the wave momentum flux term Π wave , 2) the wind pressure momentum flux term Π P , 5) the ideal wind momentum flux acceleration term Π U and 9) the dissipative viscous momentum flux term Π ξ . As for the energy flux, this small term contains space derivative d/dx and we have singular perturbation theory for a system of differential equations for the wind variables T (x) Eq. (38) and U(x) Eq. (36).
was qualitatively proposed by us several months before launching the Hinode mission; now we are presenting only the detailed MHD analysis confirming our qualitative consideration ten years ago.

Plasma heating by AW -a historical perspective
What have we learned from the one-dimensional static MHD problem? We have demonstrated that qualitatively predicted self-induced opacity of plasma is an intrinsic property. Absorption of AW causes viscous heating of ions and that is why the proton temperature is higher than the electron one. In this way we have revealed an effective method for ion heating which can be applied to many Figure 5: The thin increasing line is the temperature profile T (x) from MHD calculation presented in Fig. 1. This illustrative MHD example is imposed together with observational data of Ref. [13]. The thick line with mark 'T ' is T (x) from Ref. [13], spectral lines of elements reveal how this observational curve is derived; the thick line with mark 'N' is the particle density n(x) also in logarithmic scale in right ordinate. In the lower part of the TR between 10 4 K and 10 5 K MHD calculation almost perfectly matches the observational data. For higher temperatures the MHD calculation which contains only one monochromatic AW gives lower temperatures than the experimental profile. Optimistically we believe that MHD calculations with several AW could fit to whole observational T (x).
plasma problems. Actually plasma heating by MHD waves is used in the MIT alcator [44]. We suggest however that the toroidal geometry can be replaced by Budker probkotron geometry [45], in which the energy of the AW will be focused in a narrow jet with a hundred times increased temperature. A de Laval nozzle will be realized by strong magnetic fields. We do believe that this will be an effective method for navigation in the Solar System (cf. [46]). Electric power from a nuclear reactor will create a fast electron-proton jet and this will dramatically decrease the initial mass of the rocket. For large-scale Earth-based installations such a jet of high-temperature deuterium will inject a fresh idea in nuclear fusion physics. For this purposes, however, we need of detailed function of temperature and magnetic field dependence of plasma kinetic coefficients. Such investigations started seventy years ago by R. Landshoff [47], E. S. Fradkin in 1951 [48], I. E. Tamm [49] and S. I. Braginskii in 1952 [50]. The main final results are listed in the Kinetics by Lifshitz and Pitaevskii [23]. Unfortunately we are unaware of recent reviews elaborating the old one by Braginskii [25]. Other unresolved problems are the second viscosity, kinetic coefficients and enthalpy of partially ionized hydrogen plasma. Having approximated experimental data with our calculation of temperature T (x) and wind velocity U(x) profiles, our work can be extended to the chromosphere. In such a way, solving the problem of kinetics and magnetohydrodynamics of solar plasma we can concentrate our attention to solve the inverse problem of revealing the spectral density of MHD waves coming from the invisible layers of the sun.
Our MHD analysis confirmed that the narrow transition region is almost a permanent explosion of chromospheric plasma by self-induced opacity with respect of MHD waves. There is no other scenario able to explain the width of the transition region. Except the thermal expansion of completely ionized plasma, a significant contribution for launching of the solar wind gives intensive low frequency MHD waves reflected from the TR as Fresnel reflection by a jump of the mass density ρ(x) = Mn(x) and of the corresponding refraction coefficient . The local Alfvén speed V A (x) = B 0 / µ 0 ρ(x) also has step-like behavior and the jump of the effective refraction coefficient creates strong reflection too, giving momentum to the TR. Unfortunately, up to now there is no reliable observational data for the height profile the solar wind U(x). Using our routines we will continue MHD analysis in this direction but it will be nice to have alternative considerations. The contribution of low-frequency AW to the launching of solar wind is the next problem which has to be set in the agenda of the astrophysics. Existence of intensive Alfénic modes in TR can be confirmed by systematic investigation of nonthermal broadening of spectral lines of Mg, and other elements both in the chromosphere and transition region. Moreover, taking into account that magnetic field is more or less vertical with respect to the Sun surface, the transversal AW create horizontal oscillations of the velocity and the transversal Doppler effect will give maximal broadening close to periphery of the solar disk and minimal at the center.
In the present work we concentrated our attention to the transversal to the magnetic field AW. In the chromosphere however, ionization-recombination processes will create a significant second viscosity ζ and high-frequency slow magnetosonic waves should be absorbed in the chromosphere. For heating of the chromosphere we need to know temperature dependence ζ(T ) and the enthalpy of hydrogen plasma. We are unaware of such papers and will put this task in the agenda for further research.
For our approach, incoming AW are simply boundary condition. AW propagate along the magnetic force lines as through wave-guide. It is now subject of speculations where AW are created by perturbation of the magnetic pipes by turbulent convection or MHD waves are created in the solar tachocline [51] by self-sustained-wave turbulence in shear flow. The last scenario is perhaps working for the jets accompanying accretion [52]. In accretion disks the the effective viscosity and heating is created by self-sustained wave turbulence. Small part of accreting plasma following magnetic force lines is heated and accelerated in jets analogously to the solar wind; this will be a task of space MHD in the next decades.
One dimensional wave-guide approximation can be extended far from the TR to the distances of order of solar radius and even the size of solar system. In this case it is necessary to take into account slow adiabatic increasing of the area of the plasma stream and corresponding decrease of the constant longitudinal magnetic field at approximately constant magnetic flux.
Such type of contribution for the propulsion can be realized in a probkotron used simultaneously as a resonator. In such a way, hydrodynamic approach puts in the agenda many problems which can be easily solved knowing the detailed magnetic and temperature dependence of kinetic coefficients of solar plasma.
The main defect of the present research is that in frequency ω representation we solve static equations which require sophisticated numerical methods. All those results have to be repeated in time t dependent representation using the brute force of Monte Carlo averaging of standard MHD calculations.
A review of various mechanisms for solar corona heating could be easily made -most surely there is no plasma physical process that has not been used as the long-sought mechanism for solar corona heating. Nevertheless, we will point some alternative scenarios, which according to us even qualitatively contradict to the observations.
The longitudinal magnetosonic waves are intensely absorbed in the chromosphere through the large second viscosity and in practice do not reach the corona. And in reality, the intensity of the longitudinal waves in the solar wind is significantly lower than the Alfvén waves intensity.
If nanoflares and magnetic reconnections are the main mechanism for heating, then plasma energy absorption will be by Ohmic heating of the electrons. However, the proton temperature is higher than the electron one in the solar wind [26]. For more thorough review, please see one of the most recent works [32].
According the best we know there are no other attempts for calculation of the temperature profile of TR. In such a way the acceptable agreement between the observational data and our MHD analysis allows to conclude that after 70 years we have reliably returned to Alfvén [2] idea that heating of solar corona is predominantly by AW and the problem of set by Grotrian and Edlen: what heats the solar corona is already solved -80 years perhaps the most long living problem of the science.
Authors are thankful to Yana Maneva [53][54][55][56] and Martin Stoev for the collaboration in the early stages of the present research [14] when the idea of self-induced opacity was advocated and consideration of many problems related to physics of solar corona. Fruitful comments and discussions with Tsvetan Sariisky are highly appreciated.
Authors are also thankful to Ramesh Chandra and Reetika Joshi for discussions of the non-thermal broadening of Mg II lines, Vassil Vassilev for his interest in our work, Ivan Zhelyazkov, Yana Maneva, Steven R. Cranmer, Eckart Marsch, Leon Ofman, Stefaan Poedts, Jaime Araneda, Plamen Angelov, Ramesh Chandra and Yurii Dumin, for short discussions and comments. distance r T between two protons with thermal energy T. The coefficients 0.6 and 0.4 are added to the final exact formulas for respectively the electric conductivity and viscosity coefficients of completely ionized hydrogen plasma.
In a gaseous approximation the reciprocal mean free path for one electron is additive 1 l e = σ ee n e + σ ea n a , where σ ea is the electron-atom cross section and n a is the atom density. The mean free path l e is connected with the mean free time and the collision rate ν e = 1/τ e . The elementary kinetic estimate for the electric conductivity ς and the ohmic resistivity [23] connected with it gives for the magnetic field diffusion coefficient √ mT e e 2 n e n e Λ 0.6 At high temperatures n a → 0 this formula gives the well known result for the Ohmic resistance and electric conductivity of fully ionized hydrogen plasma In such a way, the approximate estimation for the transport section Eq. (64) exactly reproduces the well-known result [23]. Analogously, the heat conductivity coefficient is determined by the heat capacity per one electron c e = For completely ionized hydrogen plasma the heat conductivity coefficient becomes Similarly, for the reciprocal proton mean free path we have 1 l p = n p σ pp + n a σ pa , In the last proportion, it is taken into account that for elastic electron proton collision the exchange of energy is of order of m/M; it can be derived considering the scattering in the system of center of the mass. Our elementary consideration τ ε ep = τ e (M/m) gives where the numerical coefficient 8 √ 2π ≈ 20 in Eq. (42.5) of Ref. [23] requires state-of-the-art consideration. The kinetic of heat exchange between protons and electrons starting from α I = 1/2. We also note that protons determine the mass density n ρ = n p + n a = n p /α I , and then for the pressure we have P = n tot T, n tot = n e + n p + n a = (1 + α I )n ρ .
Ionization degree is important ingredient for the physics of the chromosphere. However, for the temperature of 10 kK and the low density of the upper chromosphere below the TR, the hydrogen is almost completely ionized. Since our calculation starting point is just before the TR, it is justified to assume full ionization α I = 1.
Knowing total cross-section of ionization by electron impact upon atom σ ion we can evaluate the ionization rate and make elementary evaluation of second viscous coefficient of weakly ionized hydrogen plasma in the chromosphere. Large second viscosity in the chromosphere leads to significant damping of the MHD waves having pressure oscillations. That is why in the TR we concentrate our attention on AW.

C Bremsstrahlung
In this section we will make a model evaluation of the order of the bremsstrahlung in the TR. We follow the style of model evaluations by Migdal [59] Our staring point is the total radiation energy. The total radiation for the collision of electron and proton is [60], cf. also [61].
In the chromospheric boundary of the TR the temperature is low and electron thermal velocity v Te = √ T/m is much lower than the Bohr velocity v B = e 2 / . For small relative velocities of the charged particles v 0 v B we have to use classical electrodynamics for calculation of the energy ∆E of radiation emission. For the hydrogen plasma this means that temperature should be lower than atomic energy 1 Ry=13.6 eV=158 kK. Let us recall the well-known result for the back scattering of two charged particles [60] where µ = mM/(m + M) is the reduced mass. We are applying this formula for electrons and protons with equal charges and m M. With electron energy and the constant of Coulomb interaction we can construct a parameter with dimension of length r 0 = e 2 /µv 2 0 and to make a model evaluation for the cross-section for energy loss by radiation (effective radiation) where C BS is a dimensionless constant of the order of one. This evaluation agrees with the exact result for the repulsing particles For the hot Corona where v Te v B we have to use the results of the quantum electrodynamics [61]. For two electrons, for example, the dipole emission is zero and the quadrupole emission gives Analogously for the momentum the integrating in one dimensional case the density of gravitational force f g = −ρ(x) d x Φ(x) = −ρ(x)g according to momentum conservation law Eq. (3) gives Similarly, the momentum flux in dimensionless units becomes The negative sign of τ grav means that gravitational forces decrease the momentum density of solar wind. The derived dimensionless energy and momentum gravitational fluxes could be easily added to the calculation, or can be evaluated as a perturbation to the calculated solution. Because both of these fluxes have negligible contribution in our region of interest [20], we have not included them in the numerical example analyzed in the present work.
Moreover, taking into account the TR width λ ∼ 100 km and the temperature in its chromospheric end at T = 10 kK for the worst case scenario, we get Mg λ/k B T ≈ 3 × 10 −2 1, i.e. gravity is a negligible perturbation in the TR.

E Momentum equation and distribution of the temperature
After some algebra following Eq.
For an ideal fluid with negligible kinetic coefficients η, ζ and all terms of this equation are zero.
It is possible to check that upper equation for entropy transport and production is mathematically equivalent to the equation for the conservation of energy Eq. (2). In this sense Eq. (2) can be considered as a proof the evolution equation of the magnetic field Eq. (14). That is why one can say the the complete set of MHD equations Eq. (14) and Eq. (15) can be considered as consequence of mass Eq. (1), energy Eq. (2), and momentum Eq. (3) conservation equations, and the equation of the temperature field T (r, t) is is also consequence of the conservation equations.

F Numerical methods
The differential equations we related to MHD analysis of solar corona are stiff and for them we have to apply the implicit Euler equation in the simplest case The wind variables y = (T, U), and the conservation laws Eq. (26) and Eq. (28) give the possibility to express their derivatives explicitly. We however used the explicit formulas for wind variables (T (x), U(x)): Eq. (36) and Eq.
A few words we have to say for the numerical determination of limes. Let s 0 , s 1 , s 2 , . . . are the values of y(x+h) calculated by steps h, h/2 1 , h/2 2 , h/2 3 , . . . . We can apply ε-algorithm for the calculation of y(x + h) as the limes of the sequence {s 0 , s 1 , s 2 , s 3 , . . . }. In such a way, from the Euler method we can build a high-order method using Padé approximants. Analogously, let the sequence y 0 , y 1 , y 2 , . . . is the Aitken polynomial extrapolation of y(x + h) calculated using 0,1,2,3, . . . interpolation points. Applying again the epsilon algorithm we arrive at N-point Padé approximant predictor method for solving system of differential equation; according to the best we know such an opportunity is not used up to now. The ε-algorithm is fast convergent and the optimal order is determined by the error of truncation, i.e. our method has an adaptive order. In such a way, in spite of time consuming, we have a precision method for solving system of ordinary differential equations, which utilizes all numerical resources at a fixed machine precision. Details of the numerical realization of ε-algorithm and will be given in a separate publication.

G Influence of the magnetic field on the kinetic coefficients
The influence of magnetic field on the kinetic coefficients is significant for the solar corona, but for the TR it is also important to be taken into account. Unfortunately it is not standard for the MHD considerations of the solar corona and that is why the full set of equations is given only in this appendix. Let us start our consideration with the tensor of magnetic viscosity in the Drude τapproximation for the electron drift velocity V dr and corresponding tensors for the conductivityσ, resistivityˆ Ω and magnetic diffusivityν m j = n e q e V dr =σE, E =ˆ Ω j,ν m ≡ ε 0 c 2ˆ Ω .