Model of a fluxtube with a twisted magnetic field in the stratified solar atmosphere

TL;DR

This paper models a non-expanding, twisted magnetic fluxtube in the stratified solar atmosphere, providing solutions consistent with observed magnetic field strengths and radii of magnetic bright points.

Contribution

It introduces a new magnetohydrostatic model of a twisted, non-expanding fluxtube with realistic boundary conditions and external pressure, solving the Grad-Shafranov equation analytically.

Findings

Model estimates match observed magnetic field strengths.

Derived radii are consistent with magnetic bright points.

Thermodynamic quantities inside the fluxtube are characterized.

Abstract

We build a single vertical straight magnetic fluxtube spanning the solar photosphere and the transition region which does not expand with height. We assume that the fluxtube containing twisted magnetic fields is in magnetohydrostatic equilibrium within a realistic stratified atmosphere subject to solar gravity. Incorporating specific forms of current density and gas pressure in the Grad--Shafranov equation, we solve the magnetic flux function, and find it to be separable with a Coulomb wave function in radial direction while the vertical part of the solution decreases exponentially. We employ improved fluxtube boundary conditions and take a realistic ambient external pressure for the photosphere to transition region, to derive a family of solutions for reasonable values of the fluxtube radius and magnetic field strength at the base of the axis that are the free parameters in our model. We find that our model estimates are consistent with the magnetic field strength and the radii of Magnetic bright points (MBPs) as estimated from observations. We also derive thermodynamic quantities inside the fluxtube.

Tables

Table 1. Table 1 : A formulary of the derived functions obtained from the solution of GSE. Here, ς ​ ( r ) = C ​ [ F 0 ​ ( − η , a ​ r 2 2 ) + F 0 ∗ ​ ( − η , a ​ r 2 2 ) ] 𝜍 𝑟 𝐶 delimited-[] subscript 𝐹 0 𝜂 𝑎 superscript 𝑟 2 2 subscript superscript 𝐹 0 𝜂 𝑎 superscript 𝑟 2 2 \varsigma(r)=C\big{[}F_{0}\big{(}-\eta,\frac{\sqrt{a}r^{2}}{2}\big{)}+F^{*}_{0}\big{(}-\eta,\frac{\sqrt{a}r^{2}}{2}\big{)}\big{]} and Z ​ ( z ) = e − k ​ z 𝑍 𝑧 superscript 𝑒 𝑘 𝑧 Z(z)=e^{-kz} . The value of the constants are μ ¯ = 1.116 ¯ 𝜇 1.116 \bar{\mu}=1.116 , g = 2.74 × 10 4 𝑔 2.74 superscript 10 4 g=2.74\times 10^{4} cm s -2 , k = 3.248 × 10 − 8 𝑘 3.248 superscript 10 8 k=3.248\times 10^{-8} cm -1 , p 20 = 1.36 × 10 4 subscript 𝑝 20 1.36 superscript 10 4 p_{20}=1.36\times 10^{4} dyne cm − 2 . 2 {}^{-2}. All the quantities in the table are in cgs units.

Functions	$r$–part	$z$–part
$\Psi (r,z)$	$\varsigma \left(r\right)$	$Z\left(z\right)$
$B_r(r,z)$	$\frac{3.248\times {10}^{-8}}{r}\varsigma \left(r\right)$	$Z\left(z\right)$
$B_{\varphi }(r,z)$	$\sqrt{\alpha }\frac{\varsigma \left(r\right)}{r}$	$Z\left(z\right)$
$B_z(r,z)$	$\frac{{\varsigma }^{\prime }\left(r\right)}{r}$	$Z\left(z\right)$
$p(r,z)$	$\frac{a}{8\pi }{\varsigma }^2\left(r\right)+p_{20}$	$Z^2\left(z\right)$
$\rho (r,z)$	$1$	$\left(3.22\times {10}^{-8}\right)Z^2\left(z\right)$
$T(r,z)$	$0.0416\left(\frac{a}{8\pi }{\varsigma }^2\left(r\right)+p_{20}\right)$	1

Table 2. Table 2 : Values of the quantities a , α 𝑎 𝛼 a,\alpha and C 𝐶 C for R = 100 𝑅 100 R=100 km and B 0 = 1 subscript 𝐵 0 1 B_{0}=1 kG for three different modes.

$R$ (km)	$B_0$ (kG)	mode no.	C (${10}^{17}$ Mx)	$\alpha$ (${10}^{-14}$ cm-2)	a (${10}^{-28}$ cm-4)
100	1	1	0.335061	7.50448	20.0417
100	1	2	0.144828	15.6546	107.27
100	1	3	0.0915463	24.0161	268.473

Table 3. Table 3 : Table of the results obtained from our model for r 0 = 80 subscript 𝑟 0 80 r_{0}=80 and 100 100 100 km

$r$	$z$ (Mm)	$B_z$ (G)	$p$ (dyne cm-2)	$\rho$ (g cm-3)	$T$ (K)
0	0	2370	$1.358\times {10}^4$	$3.22\times {10}^{-8}$	5656
0	2	2.19	$3.09\times {10}^{-2}$	$7.33\times {10}^{-14}$	5656
$r_0$	0	0	$1.373\times {10}^4$	$3.22\times {10}^{-8}$	5718
$r_0$	2	0	$3.12\times {10}^{-2}$	$7.33\times {10}^{-14}$	5718

Equations

$$-\nabla p+\frac{1}{4\pi }(\nabla \times \mathbf{B})\times \mathbf{B}+\rho \mathbf{g}=0,$$

$$B_r=-\frac{1}{r}\frac{\partial \Psi }{\partial z};B_z=\frac{1}{r}\frac{\partial \Psi }{\partial r};B_{\varphi }=\frac{I_p}{r}.$$

$$\displaystyle-\frac{\partial p}{\partial r}+\frac{1}{4\pi}\bigg{(}B_{r}\frac{\partial B_{r}}{\partial r}+B_{z}\frac{\partial B_{r}}{\partial z}-\frac{1}{2}\frac{\partial B^{2}}{\partial r}\bigg{)}=0$$

$$\displaystyle-\frac{\partial p}{\partial z}+\frac{1}{4\pi}\bigg{[}\frac{1}{r}\frac{\partial\Psi}{\partial z}\bigg{(}\frac{1}{r^{2}}\frac{\partial\Psi}{\partial r}-\frac{1}{r}\frac{\partial^{2}\Psi}{\partial r^{2}}\bigg{)}-\frac{1}{r^{2}}\frac{\partial\Psi}{\partial z}\frac{\partial^{2}\Psi}{\partial z^{2}}-\frac{1}{2r^{2}}\frac{\partial}{\partial z}(I_{p}^{2})\bigg{]}-\rho g=0.$$

$$\frac{{\partial }^2\Psi }{\partial r^2}-\frac{1}{r}\frac{\partial \Psi }{\partial r}+\frac{{\partial }^2\Psi }{\partial z^2}=-\frac{1}{2}\frac{\partial I_p^2(\Psi )}{\partial \Psi }-4\pi r^2\frac{\partial p_1(\Psi )}{\partial \Psi }.$$

$$\displaystyle-\frac{\partial p_{1}(\Psi)}{\partial z}+\frac{1}{4\pi}\bigg{[}\frac{1}{r^{2}}\frac{\partial\Psi}{\partial z}\bigg{(}\frac{1}{r}\frac{\partial\Psi}{\partial r}-\frac{\partial^{2}\Psi}{\partial r^{2}}\bigg{)}-\frac{1}{r^{2}}\frac{\partial\Psi}{\partial z}\frac{\partial^{2}\Psi}{\partial z^{2}}-\frac{1}{2r^{2}}\frac{\partial I_{p}^{2}(\Psi)}{\partial z}\bigg{]}-\rho g=0.$$

$$p(r,z)=p_1(\Psi )+p_2(z).$$

$$p(r,z)=p_1(\Psi )+p_2(z)$$

$$I_p^2(r,z)=\tilde{\alpha }{\Psi }^2+\tilde{\beta }\Psi +I_0^2$$

$$p_1(\Psi )=\tilde{a}{\Psi }^2+\tilde{b}\Psi ,$$

$$\frac{{\partial }^2\Psi }{\partial r^2}-\frac{1}{r}\frac{\partial \Psi }{\partial r}+\frac{{\partial }^2\Psi }{\partial z^2}=-(a{r}^2+\alpha )\Psi -(b{r}^2+\beta ),$$

$$\frac{{\partial }^2{\Psi }_h}{\partial r^2}-\frac{1}{r}\frac{\partial {\Psi }_h}{\partial r}+\frac{{\partial }^2{\Psi }_h}{\partial z^2}=-(a{r}^2+\alpha ){\Psi }_h$$

$$\frac{{\partial }^2{\Psi }_p}{\partial r^2}-\frac{1}{r}\frac{\partial {\Psi }_p}{\partial r}=-(a{r}^2+\alpha ){\Psi }_p-(b{r}^2+\beta ).$$

$$\frac{S^{\prime \prime }}{S}-\frac{1}{r}\frac{S^{\prime }}{S}+a{r}^2+\alpha =-\frac{Z^{\prime \prime }}{Z}=-k^2,$$

$$Z(z)=C{e}^{-kz},$$

$$\displaystyle\frac{{\rm d}^{2}S}{{\rm d}x^{2}}+\bigg{(}1+\frac{2\eta}{x}\bigg{)}S=0,$$

$$\displaystyle S(r)=C_{1}F_{0}\big{(}-\eta,\frac{\sqrt{a}r^{2}}{2}\big{)}+C_{2}G_{0}\big{(}-\eta,\frac{\sqrt{a}r^{2}}{2}\big{)}.$$

$$F_0(-\eta ,\frac{\sqrt{a}{r}^2}{2})=C_0(\eta )M_{-i\eta ,1/2}(i\sqrt{a}{r}^2),$$

$$G_0(-\eta ,\frac{\sqrt{a}{r}^2}{2})=i{C}_0(\eta )M_{-i\eta ,1/2}(i\sqrt{a}{r}^2)+D_0(\eta )W_{-i\eta ,1/2}(i\sqrt{a}{r}^2),$$

$$C_0(\eta )=\frac{1}{2}|\Gamma (1-i\eta )|e^{-\frac{\pi }{2}(i-\eta )}$$

$$D_0(\eta )=\frac{\Gamma (1-i\eta )}{|\Gamma (1-i\eta )|}{e}^{-\frac{\pi \eta }{2}}.$$

$$\displaystyle S(r)=C_{1}F_{0}\big{(}-\eta,\frac{\sqrt{a}r^{2}}{2}\big{)},$$

$$\displaystyle\Psi_{h}(r,z)=Ce^{-kz}F_{0}\big{(}-\eta,\frac{\sqrt{a}r^{2}}{2}\big{)}.$$

$$B_r(R,z)=0$$

$$B_{\varphi }(R,z)=0$$

$$p(R,0)=p_0$$

$$p(R,z_{tr})=p_{tr}$$

$$p_{ex}(z)=p_0\exp (-2kz),$$

$$p_2(z)=p_{20}\exp (-2kz),$$

$$p(\Psi ,z)$$

$$I_p^2(\Psi )$$

$$\displaystyle-\frac{\partial p}{\partial r}\bigg{|}_{r=R}+\frac{1}{4\pi}\bigg{(}B_{r}\frac{\partial B_{r}}{\partial r}+B_{z}\frac{\partial B_{r}}{\partial z}\bigg{)}\bigg{|}_{r=R}-\frac{\partial}{\partial r}\bigg{(}\frac{B^{2}}{8\pi}\bigg{)}\bigg{|}_{r=R}+j_{\phi}(r)B_{z}(r)\bigg{|}_{r=R}=0.$$

$$\displaystyle-\int_{R-\epsilon}^{R+\epsilon}\frac{\partial p}{\partial r}dr+\frac{1}{4\pi}\bigg{(}\int_{R-\epsilon}^{R+\epsilon}B_{r}\frac{\partial B_{r}}{\partial r}dr+\int_{R-\epsilon}^{R+\epsilon}B_{z}\frac{\partial B_{r}}{\partial z}dr\bigg{)}$$

$$\displaystyle-\int_{R-\epsilon}^{R+\epsilon}\frac{\partial}{\partial r}\bigg{(}\frac{B^{2}}{8\pi}\bigg{)}dr+\int_{R-\epsilon}^{R+\epsilon}j_{\phi}(r)B_{z}(r){\rm d}r=0,$$

Full text

††journal: Journal of LaTeX Templates††journal: Advances in Space Research

1 Introduction

The study of small scale magnetic structures in the solar photosphere is important because they play a crucial role in the evolution of active regions and sunspots (Muller and Mena, 1987; Centeno et al., 2007). Magnetic bright points (MBPs) are likely to be the fluxtubes observed in the photosphere (Berger et al., 1995; Centeno et al., 2007; Lagg et al., 2010). The topological rearrangement of these magnetic fluxtubes due to the motion of the photospheric foot points or magnetic reconnections, contribute to the coronal heating (Muller et al., 1994; van Ballegooijen, 1986). A three dimensional (3D) single fluxtube model with untwisted magnetic field has been studied by solving linear elliptic partial differential equation by numerical iterative process (Steiner et al., 1986). Schlüter and Temesváry (1958) and Osherovich (1984) studied a 3D fluxtube for sunspots using a self–similar model. The magnetic and thermodynamic structure inside fluxtube with untwisted magnetic field which spans from photosphere to the lower part of the solar corona is studied by Gent et al. (2013). Both 2D and 3D numerical models of fluxtubes with the energy propagation through the torsional Alfven waves have been studied by Murawski et al. (2015a, b), where an empirical form of magnetic flux function was assumed. Vigeesh et al. (2009) assumed an empirical form of gas pressure to investigate the wave propagation and energy transport through a fluxtube. Several interesting results of wave behavior in the solar photosphere and chromosphere have been presented by several authors (Bogdan et al., 2003; Fedun et al., 2009; Shelyag et al., 2010).

In this work, we construct a 3D single cylindrical vertical straight magnetic fluxtube semi-analytically with a twisted magnetic field by obtaining a new solution of poloidal flux function by solving Grad–Shafranov equation (GSE; Grad and Rubin (1958); Shafranov (1958)). We assume a specific form of gas pressure and poloidal current, which has been used to study the equilibrium solution of terrestrial plasma (Atanasiu et al., 2004). An equilibrium solution near the magnetic axis of the plasma torus has been reported previously, using a plasma pressure and poloidal current profile that varies linearly with the poloidal flux function (Solov’ev, 1968). We obtain an analytic solution by assuming a form that is quadratic in the poloidal flux function, and derive the magnetic field structure and thermodynamic quantities inside the fluxtube using the solution that represents an ideal MHS equilibrium. In the future, we will look to explore fully the profile functions that will improve the solution set.

The overview of the paper is as follows. In section 2, the GSE has been derived assuming a specific form of the profile function of gas pressure and poloidal current and the solution of the equation is presented. In section 3, we discuss the boundary condition that is physically acceptable, and can be used for realistic modelling of a fluxtube. In section 4, the mode wise variation of the profile functions are presented and in section 5, we compare the model with the observations. Finally, we conclude with a comparison with other existing models.

2 Solution of Grad-Shafranov equation

We assume an axisymmetric cylindrical geometry, with gas pressure $p$ and take the poloidal current $I_{p}$ constant along a magnetic field line. We express $p(\Psi,z)$ and $I_{p}(\Psi)$ in terms of the poloidal flux function $\Psi(r,z)$ and $z$ and consider a straight vertical axisymmetric fluxtube that spans the altitude from photosphere ($z=0$) to the transition region ($z=2.15$ Mm) that is in equilibrium with the atmosphere outside with the uniform gravity ${\bf g}(=-g\hat{z})$ acting vertically downward. The force balance equation in MHS equilibrium takes the form

[TABLE]

where $\rho$ denotes the mass density and ${\bf B}$ is the magnetic field associated with the poloidal flux function $\displaystyle{\Psi(r,z)=\int_{0}^{r}B_{z}(r^{\prime},z)r^{\prime}{\rm d}r^{\prime}}$ (scaled by the factor $\frac{1}{2\pi}$) in the following form

[TABLE]

This form of $B_{r},B_{\phi}$ and $B_{z}$ ensures the solenoidal condition of magnetic field. Now splitting the MHS force balance equation (1) into $r$ and $z$ directions, we find two different scalar partial differential equations

[TABLE]

If the gas pressure and poloidal current are functions of $\Psi$ alone i.e., $p_{1}(\Psi)$ and $I_{p}(\Psi)$ respectively, then from the (3a, 2) it follows that

[TABLE]

Plugging in $p_{1}$ and $I_{p}^{2}$ in (3b), we find

[TABLE]

By multiplying both sides of (5) by $\displaystyle{4\pi r^{2}\frac{\partial z}{\partial\Psi}}$ and using (4), we obtain $\displaystyle{g\rho\frac{\partial z}{\partial\Psi}=0}$ which implies that $\rho$ is zero, which means that the vertical hydrostatic pressure balance will not be maintained. Therefore, to balance the vertical hydrostatic pressure inside the fluxtube, we introduce a new function, $p_{2}(z)$ such that,

[TABLE]

We assume $p_{1}(\Psi)$ and $I_{p}^{2}(\Psi)$ to be second order polynomials of $\Psi$

[TABLE]

where

[TABLE]

and the parameters $\tilde{a},\tilde{b},\tilde{\alpha},\tilde{\beta}$ and $I_{0}^{2}$ are to be determined by appropriate boundary conditions. The function $p_{2}(z)$ is to be evaluated later. The substitution of $p$ given by (6a) in (3a) gives (4) and we obtain the following second order scalar partial linear inhomogeneous differential equation

[TABLE]

with the rescaled parameters, $\displaystyle{a=8\pi\tilde{a};\alpha=\tilde{\alpha};b=4\pi\tilde{b};\beta=\tilde{\beta}/2}$. To solve (7), we split $\Psi$ in two parts: a homogeneous part, $\Psi_{h}(r,z)$ and an inhomogeneous part $\Psi_{p}(r)$, i.e. $\Psi(r,z)=\Psi_{h}(r,z)+\Psi_{p}(r)$. Using this form in (7), we separate the homogeneous and the inhomogeneous parts to obtain the following expressions:

[TABLE]

To solve the homogeneous part, we seek a solution of the form $\Psi_{h}(r,z)=S(r)Z(z)$. Then we separate out the $r$ and $z$ part in (8a) as follows

[TABLE]

where $k$ is an arbitrary real constant. Motivated by the fact that the poloidal flux function $\Psi(r,z)$ decreases with $z$, we assume that the solution of the $z$-part of (9) takes the form

[TABLE]

where $C$ is an arbitrary constant. To solve the $r$–part of (9), we substitute $x=\frac{\sqrt{a}r^{2}}{2}$ (where $a>0$) and insert it in (9) to find

[TABLE]

whose solutions are given by Coulomb wave functions $F_{L}(-\eta,x)$ and $G_{L}(-\eta,x)$ (Abramowitz and Stegun, 1972) (page $537-544$) with $L=0$ and $\displaystyle{\eta=\frac{\alpha+k^{2}}{4\sqrt{a}}}$. The solution of (11) takes the following form

[TABLE]

Here F_{0}\big{(}-\eta,\frac{\sqrt{a}r^{2}}{2}\big{)} and G_{0}\big{(}-\eta,\frac{\sqrt{a}r^{2}}{2}\big{)} are called the regular and irregular Coulomb wave functions respectively which are complex quantities with real arguments (Boersma, 1968), given by

[TABLE]

where $M_{-i\eta,1/2}(i\sqrt{a}r^{2})$ and $W_{-i\eta,1/2}(i\sqrt{a}r^{2})$ are called the Whittaker-M and Whittaker-W function (see Figure 3 below) and the constants $C_{0}(\eta)$ and $D_{0}(\eta)$ are defined by

[TABLE]

and

[TABLE]

The Whittaker function has been used by several authors in their models of the solar atmosphere albeit in different physical problems (eg. Tsinganos (1979) in the context of inviscid flows, and also in the context of MHD waves by Hindman and Jain (2008) and Erdélyi and Fedun (2010)). Now $B_{z}(r,z)$ has to be a finite quantity that varies linearly with the term $\displaystyle{\frac{1}{r}\frac{\mathrm{d}S}{\mathrm{d}r}}$ but $\displaystyle{\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}[G_{0}(-\eta,\frac{\sqrt{a}r^{2}}{2})]}$ blows up at $r=0$; therefore for $B_{z}$ to be finite on the axis of the fluxtube $C_{2}$ in (12) must vanish. As a result $S(r)$ takes the form

[TABLE]

and the homogeneous part of the solution is given by

[TABLE]

A similar but a different solution, which is oscillatory in $z$ is used for laboratory plasma for both a $D$-shaped plasma and toroidally diverted plasma (Atanasiu et al., 2004). The general solution of (7) is given by the sum of the homogeneous part $\Psi_{h}(r,z)$ given above and an inhomogeneous part $\Psi_{p}(r)$ which will be presented in a paper in preparation. We have found that the presence of $\Psi_{p}(r)$ term in the poloidal flux function $\Psi(r,z)$, implies that $p$ and $I_{p}^{2}$ cannot be simultaneously positive for any combination of $b$ and $\beta$ in the physical parameter domain space for all $r$ and $z$. For avoiding these unphysical effects we present the case of $\Psi=\Psi_{h}$ and an exploration of the general solution $\Psi=\Psi_{h}+\Psi_{p}$ will be studied later. Since $\Psi(r,z)$ and its complex conjugate function, $\Psi^{*}(r,z)$ are the valid solutions of (7), we construct a solution of (7) by redefining $\displaystyle{\frac{\Psi(r,z)+\Psi^{*}(r,z)}{2}\rightarrow\Psi(r,z)}\equiv\varsigma(r)Z(z)$.

3 Boundary conditions and the reduced form of $p$ and $I_{p}$

The ideal magnetic fluxtube is embedded in a magnetic field free region with no current outside the fluxtube boundary. We make the following standard assumptions $B_{r}(r=R,z)=0$ and $B_{\phi}(r=R,z)=0$ to ensure that there is no net current $I_{p}$ at the fluxtube boundary. The pressure at the photosphere ($z=0$) outside the fluxtube is $p_{0}=1.228\times 10^{5}$ dyne cm*-2* and at the transition region ($z_{tr}=2.15$ Mm) is $p_{tr}=0.1058$ dyne cm*-2* and is taken from Avrett-Loeser model (Avrett and Loeser, 2008). We summarize the boundary conditions below

[TABLE]

Assuming that pressure decreases exponentially from photosphere to transition region, we use the following expression for the external pressure

[TABLE]

where \displaystyle{k=\frac{1}{2\times 2.15}\ln\big{(}\frac{p_{0}}{p_{tr}}\big{)}} Mm${}^{-1}=3.248$ Mm${}^{-1}.$ Matching the pressure scale heights, we see that $p_{2}(z)$ (6a) also decreases exponentially with $z$ as

[TABLE]

where $p_{20}$ will need to be calculated. Taking $\Psi(r,z)=\Psi_{h}(r,z),$ the reduced forms of $p$ and $I_{p}^{2}$ are given by

[TABLE]

and

[TABLE]

Taking the radial component of the MHS force balance equation (1) and adding the contribution of the radial force due to the presence of sheet current $j_{\phi}$ at the boundary we write the force balance equation

[TABLE]

Now the sheet current $j_{\phi}$ can be expressed as a delta function $j_{\phi}(r)=j_{\phi s}\delta(r-R)$ which is non zero only at the boundary. Integrating (24) w.r.t. $r$ from $r=R-\epsilon$ to $r=R+\epsilon$ where $\epsilon$ is an infinitesimal positive quantity we obtain

[TABLE]

which leads to the MHS force balance at the boundary to be given by

[TABLE]

where $[...]_{R}$ denotes the jump condition at the boundary and $B_{i}$ and $B_{e}$ are the magnetic fields inside and outside the fluxtube boundary. This is an improved boundary pressure condition for a magnetic fluxtube, as previous studies have ignored the sheet current. Now to calculate $j_{\phi s}$, we assume an infinitesimal current loop at the boundary which has a vertical height of length $L$ and radial extent from $R-\epsilon$ to $R+\epsilon$. Using Stokes line integral theorem along the closed loop (see Fig.1), we see that

[TABLE]

which implies $\displaystyle{j_{\phi s}=\frac{B_{z}}{4\pi}}$. Since $B_{r}(R)=B_{\phi}(R)=B_{e}(R)=0$ for any height $z$ and $p_{in}(R,z)=p_{2}(z)$, from (26), the total pressure at the boundary inside the fluxtube is $\displaystyle{p_{2}(z)+\frac{3B_{z}^{2}(R)}{8\pi}}$ and matching the pressures gives

[TABLE]

The mass density inside the fluxtube obtained from (3b) is

[TABLE]

and the density inside the fluxtube varies only with $z$ and at the transition region ($z_{tr}=2.15$ Mm), which should match with the external density which is typically $\rho_{tr}=2.77\times 10^{-14}$ g cm*-3* (Avrett and Loeser, 2008). From (28) and (20) we see that

[TABLE]

In our model $g$ is assumed not to vary much from photosphere to the transition region and its value is taken to be $g=274$ m s*-2*, the value at the solar surface and that determines $p_{20}=1.36\times 10^{4}$ dyne cm*-2*. Using the forms of $B_{z}(r,z)$ from (2, 18) we obtain $B_{z}(0,0)=4\sqrt{a}C\equiv B_{z0}$, and the relation between $C,B_{z0}$ and $a$ is derived in the following.

The real component of $\Psi$ is given by,

[TABLE]

and by using (13) and (18) we can express the flux function in the form of Whittaker–M functions as

[TABLE]

Whittaker–M function can be expressed in terms of hypergeometric function by the standard relation (Dixit and Moll, 2015)

[TABLE]

Here $F_{1}^{1}$ represents the hypergeometric functions with the parameters $t,m$ and argument $z$. Therefore (31) takes the form

[TABLE]

It follows from (2) that $B_{z}$ takes the form

[TABLE]

So, $B_{z}(r,z)$ at $r=z=0$ is by definition $B_{z0}$ takes the form from (34) as,

[TABLE]

But F_{1}^{1}\bigg{(}\frac{4\sqrt{a}+i(k^{2}+\alpha)}{4\sqrt{a}},2,0\bigg{)}=1 is an identity and therefore we finally have

[TABLE]

We define a physical observable $B_{0}$ which is the average magnetic field strength at the base within the fluxtube as

[TABLE]

where $R$ is the radius of the fluxtube. Therefore we have two free parameters $R$ and $B_{0}$ that we can tune to fit our model with the observations. At $z=0$, from (28), we get

[TABLE]

Therefore, from (19a, 19b, 19c) we determine $a$, $\alpha$ and $C$ in terms of the free parameters $R$ and $B_{0}$ and hence the thermodynamic quantities within fluxtube. The temperature within the fluxtube is calculated by the ideal gas law according to the following form

[TABLE]

where $R_{g}=8.314$ J mol*-1* K*-1* is the universal gas constant and

[TABLE]

is the average value of the mean effective molar mass from photosphere to transition region given by an empirical formula \displaystyle{\mu_{eff}(z)=1.288\big{[}1-0.535\big{(}\frac{z}{2.152}\big{)}^{3}\big{]}} (Solov’ev and Kirichek, 2015) in the domain of $0<z<2.152$ Mm. A formulary of the different quantities are listed in Table 1.

4 Mode analysis of different profile functions

The quantities $a,\alpha$ and $C$ are functions of the free parameters $R$, $B_{0}$ and the mode number $n$ whose values are given in Table 2 for a sample set of the free parameters.

The solutions to $\Psi,B,p$ and $T$ are shown for different mode numbers, in Figs. 2–5 respectively.

All the profile functions are normalized to their peak values and the radial distance to the total radius of the fluxtube $R$. As per the boundary conditions, the flux function, $\Psi$ vanishes both at the axis and at the boundary of the fluxtube, where the the total gas pressure is $p_{2}(z)$. The solutions of higher modes have, the profile functions with higher frequency along the radial direction and realistically we may not have such reversible fields, as they are unstable. Therefore, we use the fundamental mode ($n=1$) for further analysis in the paper. The $3$D topology of the magnetic field lines inside the fluxtube for the fundamental mode is shown in Fig. 9.

5 Comparing the model with observations

Now we compare our model with the observations reported from high resolution and high cadence instruments. Small scale magnetic structures, i.e., MBPs, are the best candidates for comparison because such structures can be assumed to consist of fluxtubes. MBPs are seen in G-band filtergrams or are identified by making spectro-polarimetric measurements (Utz et al., 2009, 2013; Yang et al., 2016). The radial variation of the profile functions $p,\rho$ and $T$ and the magnetic components $B_{r},B_{\phi}$ and $B_{z}$ are independent of $z$, but the amplitude decreases exponentially with $z$ except for $T$. In the following, we validate the model by comparing the observed magnetic field strengths and radius of MBPs with those calculated in our model, and estimate the magnetic field strength and thermodynamic quantities at the transition region which may be verified by future observations.

The MBPs number distribution, magnetic field strength and size distribution has been reported by Utz et al. (2009, 2013) at photosphere. The size distribution of MBPs peaks around $200$ km and $160$ km for low and high spatial sampling rates, respectively (Utz et al., 2009). The magnetic field strength distribution is bimodal with two peaks at $\sim 1400$ G and $\sim 200$ G (Utz et al., 2013). Since MBPs are observed as the region of unipolar flux concentrations, we construct a cylindrical boundary inside the simulation domain where the vertical magnetic field $B_{z}$ is positive. We call this cut–off radius as $r_{0}$, where the line of sight magnetic field $B_{z}$ vanishes. The value of $B_{z}$ after this grid line becomes negative. In Fig. 6, the vertical grid line denotes the boundary radius $r_{0}$. We study two different cases for $r_{0}=80$ and $100$ km which corresponds to the peak values for the MBP size distribution, for which $R$ is found to be $127$ and $159$ km respectively.

For both cases, we calculate $B_{0}$ and the mean value of $B_{z}$, $\bar{B_{z}}$ in the radial direction up to $r_{0}$ and find that for realistic values of the thermodynamic quantities inside the fluxtube, the upper limit of the vertical magnetic field strength $B_{z0}$ is $2.37$ kG. Beyond this value of $B_{z0}$, the viable solutions will shift to the higher modes. The temperature inside the fluxtube increases as the value of $B_{z0}$ decreases and temperature inside the fluxtube becomes greater than the typical photospheric temperature when $B_{z0}<2.31$ kG. Thus it can be considered as the lower cut off limit of the magnetic field strength. The $\bar{B_{z}}$ value is only sensitive to $B_{z0}$ but not on $r_{0}$, and the thermodynamic quantities inside the fluxtube remain the same for both $r_{0}=80$ and $100$ km. We found $\bar{B_{z}}=1.42$ and $1.4$ kG for $B_{z0}=2.37$ and $2.31$ kG respectively. The radial and vertical variations of the vertical magnetic field strength, gas pressure, density and temperature inside the fluxtube are shown in Figs. 7 and 8. We see that the vertical magnetic field strength decreases from $2.37$ kG (on the axis at $z=0$) to zero at the MBP boundary ($r_{0}$, see Fig. 6). The variation of gas pressure and temperature from axis to the MBP boundary is very small; at the photosphere the gas pressure changes $1.358\times 10^{4}$ (on the axis at $z=0$) to $1.373\times 10^{4}$ dyne cm*-2* (at MBP boundary) and it decreases with $z$ to $3.12\times 10^{-2}$ dyne cm*-2* (at MBP boundary) at the transition region ($z=2$ Mm). The temperature changes from $5656$ K (on the axis) to $5718$ K (at MBP boundary) which is small compared to the outside photosphere temperature ($6583$ K) (Avrett and Loeser, 2008). The average temperature inside fluxtube has been calculated by integrating the temperature from axis to the MBP boundary and is found to be $5679$ K. The density distribution is constant along the radius of the fluxtube at a given height which decreases with height from $3.22\times 10^{-8}$ g cm*-3* at the photosphere to $7.33\times 10^{-14}$ g cm*-3* at the transition region. The values of the quantities estimated from our modelled are summarized in Table 3.

6 Conclusions and discussion

In this paper, we constructed a single fluxtube with twisted magnetic field by solving GSE analytically. We summarize our results below:

We have an improved boundary condition by incorporating the sheet current as compared to the previous studies e.g. Solov’ev and Kirichek (2016). 2. 2.

Our model depends on the form of the external pressure distribution which is assumed as an exponentially decreasing function with $z$. Future observations leading to the more accurate form of pressure distribution from photosphere to transition region can be used to improve our model. The plasma $\beta$ parameter inside the fluxtube remains constant with $z$ but it varies along $r$; $\beta<1$ is obeyed from the chromosphere to the transition region but not in the photosphere and lower atmosphere. Therefore the magnetic effects will dominate the gas dynamics throughout the simulation domain. 3. 3.

In our model, the temperature varies along the radial direction, but it is constant along the vertical direction $z$. In other models e.g. Gent et al. (2013) the temperature rises with height from photosphere to the transition region and in Vigeesh et al. (2011) the temperature decreases from $6300$ K at surface to $4000$ K (at $z=600$ km) and then it remains the same up to $1200$ km. 4. 4.

The effects of shock wave dissipation and magnetic reconnection starts to dominate in the corona which causes the coronal heating. We have not considered these mechanisms in our model and therefore, we have not model the region in the corona or higher and have restricted our simulation domain to end at the transition region. 5. 5.

Recently Hewitt et al. (2014); Uitenbroek and Criscuoli (2013); Riethmüller and Solanki (2016), have simulated bright points using MuRAM and Copenhagen-Stagger code. We find that the magnitude of magnetic field strengths, pressure and densities reported in these studies are in fair agreement with our predictions but the temperature distribution along $z$ is in variance with the results of the numerical simulations.

The fluxtube model gives values of magnetic field and thermodynamic quantities consistent with observations and compares well with simulations. We plan to present a more detailed study of the model in a paper in preparation.

We thank Prof. P. Venkatakrishnan for useful discussions and the anonymous reviewers for insightful comments and helpful suggestions. We also thank the support staff of the IIA HPC facility and Sandra Rajiva for proofreading the manuscript.