© ESO 2021

1 Introduction

Chemically peculiar (CP) main-sequence stars of A and B spectral types show abnormal strong or weak line strengths of some elements as compared with normal stars of the same spectral type (Jaschek & Jaschek 1987). These stars present highly chemical structured atmospheres in the vertical and horizontal directions (Ryabchikova et al. 2002; Ryabchikova 2005; Nesvacil et al. 2008).

Non-uniform surface chemical compositions are observed in magnetic and non-magnetic stars, but magnetic stars show the most extreme chemical peculiarities of all (Landstreet 1993).

It is well known that B-type CP stars may show helium abundance anomalies, the so-called helium-weak and helium-strong stars. These stars display photometric and spectral variations that often correlate with the rotation period. In magnetic stars, the amplitude of these variations could be either phased or anti-phased with the longitudinal magnetic field component. Historically, the oblique rotator model (ORM, Stibbs 1950), in which a strong dipole frozen into a rotating stellar plasma is inclined relative to the stellar rotation axis, was proposed to explain these variabilities.

Vallverdú et al. (2014), assuming the ORM, studied the influence of magnetic pressure effects on the atmospheric structure of helium variable stars. These authors found that the Lorentz force can explain variations with the rotation phase up to 3% in the line equivalent widths (EWs) for a dipolar magnetic field of 1000 G. These predicted values are lower than the observed values. These authors stressed that line variations originated by changes in the chemical abundances are more relevant than those expected from an outward-directed (inward-directed) Lorentz force.

Other physical processes that explain the origin of helium inhomogeneities are the presence of stellar winds or chemical diffusion. Vauclair (1975) stated that helium-weak and helium-strong stars are a consequence of the combined effects of atomic diffusion velocities and a radial mass flux (stellar wind). Assuming that the diffusion velocities and wind speed are opposite and depending on whether the mass flow exceeds the diffusion or not, this would cause the removal or sinking of a certain amount of helium from the atmosphere. Thus, the result is a helium enhancing or decreasing with respect to normal stars.

Subsequently, Vauclair et al. (1991) studied the helium inhomogeneities in main-sequence stars considering a magnetic dipole configuration. In these models, the vertical magnetic fields are at the poles, while the horizontal fields are at the equator; thus the wind is expected to flow through the poles originating helium overabundance if the mass loss rate is of the order of that caused by the solar wind. This mechanism could be a good complement to diffusion processes.

Diffusion effects become relevant in slow rotating stars with weak winds, that is late B to early F main-sequence stars (Michaud et al. 1976; Turcotte 2003), white dwarfs (e.g. Fontaine & Michaud 1979; Iben & MacDonald 1985; Althaus & Benvenuto 2000; Panei et al. 2007), and hot horizontal branch stars (Michaud et al. 1983). It is known that in a non-magnetic self-gravity model, gravitational settling and radiative accelerations lead to overabundance of certain atomic species, but not to local inhomogeneities on the stellar surface (Michaud et al. 1976). Later, Michaud et al. (1979) concluded that even considering radiative acceleration, chemical diffusion does not lead to helium-rich stellar atmospheres; these authors showed that the radiative acceleration is always below gravity for a wide range of temperatures. As a consequence, to model stellar atmospheres of CP stars it is necessary to include diffusion effects due to gravitational settling, radiative levitation, and chemical gradients in the presence of magnetic fields (Michaud et al. 2015). It is then expected that the presence of a magnetic field might have a significant impact on the atmospheric structure leading to the formation of non-uniform stellar surface distributions, such as chemical or thermal spots (Kochukhov et al. 2012).

Diffusion in stellar atmospheres of CP stars with magnetic stars was mainly developed in late B- and A-type stars. The interaction between hydrogen ambipolar diffusion and magnetic fields (Lorentz force) could change the structure of the outer layers and produce chemical inhomogeneities showing rings or spots (Babel & Michaud 1991), with a high concentration towards the equatorial regions (LeBlanc et al. 1994). Improved models in 2D and 3D were carried out by Alecian & Stift (2010, 2017). The latter also analysed the effects of a non-axisymmetric magnetic field geometry. Studies on the effects of the radiative acceleration without magnetic fields or weak magnetic fields with mass loss for heavy atomic elements were performed by LeBlanc et al. (2009) and Alecian & Stift (2019).

As diffusion in stellar atmospheres of CP stars is complicated and difficult to model, often some authors proposed to investigate the chemical anomalies by modelling a stratified stellar atmosphere and using empirical abundance profiles derived for a given star. Kochukhov et al. (2009) and Shulyak et al. (2009) suggested performing an iterative self-consistent procedure to determine the stellar parameters, chemical abundances, and stratification. These ad hoc chemical distributions in depth might be obtained by combining high-resolution spectropolarimetry and magnetic Doppler imaging techniques that enable the exploration of the surface chemical patterns of a magnetic star.

In order to go a step further in investigating the formation of non-homogeneous chemical anomalies in helium-peculiar stars, we propose to carry out a self-consistent model in which atomic diffusion is computed in the presence of a magnetic field. To this purpose, we take advantage of the calculation performed in the stellar evolution of white dwarf stars (Iben & MacDonald 1985) as well as theoretical results widely studied in helium and carbon-oxygen cores of white dwarfs (Althaus & Benvenuto 2000; Althaus et al. 2005, 2009a,b; Panei et al. 2007).

Our model assumes non-rotating magnetic stars with plane-parallel photospheric layers consisting of hydrogen, helium, and free electrons. The atomic diffusion is considered in the multi-component treatment of a mix of gases (Burgers 1969), where gravitational settling, thermal and chemical diffusion, and the Lorentz force are included. We study the evolution of the chemical composition in the atmosphere of B-type, helium-peculiar stars, in one evolutionary step, since diffusion timescales are much shorter than evolutionary timescales.

This paper is organised as follows: we present in Sect. 2 the theoretical framework to treat chemical diffusion in stellar atmospheres of early B-type stars; in Sect. 3 the numerical method of solution; in Sect. 4 our results regarding diffusion with magnetic fields and chemical concentration in latitude and depth; and, lastly, in Sect. 5 the discussion and conclusions are outlined. In addition, in Appendices A and B we demonstrate that we deal with N − 1 linear independent equations and describe the algebraic solution, respectively.

2 Input physics

Gravitationally-induced diffusion has shown to be an efficient mechanism in separating elements in stellar interiors. This mechanism leads to sink down all heavy elements and to float up the light elements towards the stellar surface. On the other hand, the inclusion of models with convective dredge-up, mass accretion, and stellar winds have supported the presence of heavy elements in the spectra of white dwarf atmospheres (Fontaine & Michaud 1979; Paquette et al. 1986; MacDonald et al. 1998, amongst others). These studies are based on the equations that describe gravitational settling, and chemical and thermal diffusion for multi-component gas mixtures (Burgers 1969).

To treat diffusion we developed a new code (in Fortran programming language) to model atomic diffusion in the stellar atmospheres of helium-peculiar stars. We assumed a gaseous medium composed by electrons and partially ionised hydrogen and helium, where inertia effects and effects connected with deviatoric stress components have been left out of the account. In our treatment, we considered effective atoms for each species. That is, our mixture consists of electrons and effective atoms of hydrogen and helium (see Sect. 4).

All species of particles (electrons and effective atoms) are treated as a binary interaction (named as type s or t) in the same element of volume at the same temperature T, where the interactions are functions of time and coordinates. In addition, we supposed there is no resulting electric charge density and no net mass flow through out the surface enclosing the element volume. Under these considerations, we simultaneously solved the diffusion equation, the heat flow equation, and the mass and charge conservation equations.

2.1 Diffusion equation

According to Burgers (1969), the general expression for the diffusion equation for a given species s states $$\begin{array}{l}\nabla p_s-\frac{{\rho }_s}{\rho }\nabla p-n_s{Z}_{s}eE-(j_s-\frac{{\rho }_s}{\rho }j)\times B/c\hfill \\ ={\sum }_{t\ne s,t=1}^N{K}_{st}(w_t-w_s)+{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}},\hfill \end{array}$$(1)

where p and ρ are the total gas pressure and mass density, respectively. Variables p_s, ρ_s, n_s, e Z_s, and m_s are partial pressure, mass density, number density, mean charge, and mass for species s. The microscopic collisional mechanism is described by binary interactions between species s and t. The sums are over t with t≠ s, and N is equal to the total number of the different types of particles (i.e. ions and neutral atoms) plus electrons, that is N = N_ion,0 + 1. We define m_st = m_s + m_t as the sum of the masses of species s and t. The vectors w_s, r_s, and j_s are the diffusion velocity with respect to the mean mass flow of the gas, the residual heat flow, and the conduction current measured with respect to the mean mass flow velocity for the species s. The vector j_s is defined by $$j_s=n_s{Z}_{s}e{w}_s.$$

The quantities E and B are the electric and magnetic fields, respectively; j is the total induced current (j = ${\sum }_{s=1}^N{j}_s$); e is the electron charge; c is the light speed; and K_st and z_st are the resistance coefficients from Paquette et al. (1986).

2.2 Heat flow equation

Following the Burgers descriptions, the heat flow equation is given by $$\frac{5}{2}{n}_s\text{k}\nabla T-\frac{n_s{Z}_{s}e}{c}(r_s\times B)=-\frac{5}{2}{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t}{m_{st}}(w_t-w_s)-\frac{2}{5}{K}_{ss}z\text{\_}ss\text{^{''}}{r}_s-{\sum }_{t\ne s,t=1}^N\frac{K_{st}}{m_{st}^2}[3m_s{}^2+m_t{}^{2}z\text{\_}st\text{^{'}}+\frac{4}{5}{m}_s{m}_{t}z\text{\_}st\text{^{''}}]r_s+{\sum }_{t\ne s,t=1}^N{K}_{st}\frac{m_s{m}_t}{m_{st}^2}[3+z\text{\_}st\text{^{'}}-\frac{4}{5}z\text{\_}st\text{^{''}}]r_t.$$(2)

In this equation k is the Boltzmann’sconstant and $z\text{\_}st\text{^{'}}$ and $z\text{\_}st\text{^{''}}$ are also resistance coefficients given by Paquette et al. (1986).

2.3 Mass and charge conservation

Assuming there are no net mass flow and no net electrical current, the mass and charge equations satisfy $${\sum }_{s=1}^N{\rho }_s{w}_s=0,$$(3)

and $${\sum }_{s=1}^N{Z}_{s}e{n}_s{w}_s=0.$$(4)

Then, there are 2N + 1 equations: Eq. (1) gives a set of N − 1 independentlinear equations (see demonstration in Appendix A), Eq. (2) has N linear equations, and two additional equations come from Eqs. (3) and (4).

As Eq. (1) is independent of gravity, we must take into account the hydrostatic equilibrium equation $$\nabla p=\rho g.$$(5)

In radial symmetry, the acceleration of gravity g is written as $$g=-\frac{GM(r)}{r^2}\check{r}=-g\check{r},$$

where G is the universal gravitational constant and M(r) is the mass contained within a radius r. In the presence of a magnetic field, Eq. (5) includes the term of Lorentz force per unit volume, $$\nabla p=\rho g+\frac{1}{c}(j\times B).$$(6)

The right member in Eq. (6) can be expressed as $$\nabla p=\rho g_{eff},$$(7)

where g_eff is an effective gravity that takes into account the effects of gravity and Lorentz force. According to Valyavin et al. (2004), $$g_{eff}=g\pm {\sum }_j\frac{{\lambda }_j\sin \theta }{c\rho (1+{\omega }_j{\tau }_j)}{E}_{\text{eq}}{B}_{\theta }.$$(8)

The index j corresponds to all species of charged particles, θ is the magnetic co-latitude angle, E_eq is the equatorial surface electric field induced by the magnetic field evolution, B_θ is the polar component of the magnetic field, λ_j is the non-magnetic conductivity, ω_j is the cyclotron frequency, and τ_j is the mean free path time of the conducting particles (see Valyavin et al. 2004; Vallverdú et al. 2014).

We define g_eff in the same direction of g and opposite to the radial direction, that is $$g_{eff}=-g_{eff}\check{r}.$$(9)

From Eqs. (7) and (8), the radial component of ∇p states that $$\frac{\text{d}p}{\text{d}r}=-\rho g\pm {\sum }_j\frac{{\lambda }_j\sin \theta }{c(1+{\omega }_j{\tau }_j)}{E}_{\text{eq}}{B}_{\theta },$$(10)

where the positive sign refers to the case of an outward directed Lorentz force and the negative sign is for an inward directed magnetic force.

On the other hand, the interactions between photons and particles make particles acquire a radiative acceleration, A_R,s, which can be written in terms of the radiation pressure, ∇p_R,s as $$A_{R,s}=\frac{\nabla p_{R,s}}{{\rho }_s}.$$(11)

From Eqs. (7) and (11), the Eq. (1) becomes $$\begin{array}{lll}\hfill & \hfill & \nabla p_s-{\rho }_s{g}_{eff}+{\rho }_s{A}_{R,s}-n_s{Z}_{s}eE-j_s\times B/c\hfill \\ \hfill & \hfill & ={\sum }_{t\ne s,t=1}^N{K}_{st}(w_t-w_s)+{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}}.\hfill \end{array}$$(12)

Since we want to evaluate the effect of the diffusion in the radial direction, the r-component of cross product j_s ×B is null and the Eq. (12) results, $$\begin{array}{lll}\hfill & \hfill & \frac{\text{d}{p}_s}{\text{d}r}-{\rho }_s(g_{eff}-A_{R,s})-n_s{Z}_{s}eE\hfill \\ \hfill & \hfill & ={\sum }_{t\ne s,t=1}^N{K}_{st}(w_t-w_s)+{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}},\hfill \end{array}$$(13)

where we considered that A_R,s is opposite to gravity, $w_s=w_s\check{r}$ and $r_s=r_s\check{r}$.

In a similar way, the radial component of the heat flow equation (Eq. (2)) can be written as $$\frac{5}{2}{n}_s\text{k}\frac{\text{d}T}{\text{d}r}=\text{\&}-\frac{5}{2}{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t}{m_{st}}(w_t-w_s)-\frac{2}{5}{K}_{ss}z\text{\_}ss\text{^{''}}{r}_s-{\sum }_{t\ne s,t=1}^N\frac{K_{st}}{m_{st}^2}[3m_s{}^2+m_t{}^{2}z\text{\_}st\text{^{'}}+\frac{4}{5}{m}_s{m}_{t}z\text{\_}st\text{^{''}}]r_s+{\sum }_{t\ne s,t=1}^N{K}_{st}\frac{m_s{m}_t}{m_{st}^2}[3+z\text{\_}st\text{^{'}}-\frac{4}{5}z\text{\_}st\text{^{''}}]r_t,$$(14)

since the cross product r_s ×B is contained in the normal planes.

Similarly, Eqs. (3) and (4) can be expressed as $${\sum }_{s=1}^N{\mu }_s{n}_s{w}_s=0,$$(15)

and $${\sum }_{s=1}^N{Z}_s{n}_s{w}_s=0.$$(16)

In Eq. (15) we replaced ρ_s = n_sμ_sM_u. Here μ_s is the molecular weight for s species and M_u is the mass unit. Then, the new system to solve comprises Eqs. (13)–(16) and constitutes a set of 2N + 1 equations with 2N + 1 unknown variables, i.e. N values of w_s, N of r_s and E.

3 Method of solution

3.1 The system of equations

To solve thesystem of Eqs. (13)–(16), we employed the method proposed by Iben & MacDonald (1985) adapted to our problem.

Firstly, we rewrite Eq. (13), in terms of n_s and the ideal gas law, p_s = n_skT, as $$\begin{array}{lll}\hfill & \hfill & \frac{1}{n_s}[{\sum }_{t\ne s,t=1}^N{K}_{st}(w_s-w_t)+{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_s{r}_t-m_t{r}_s}{m_{st}}]-Z_{s}eE\hfill \\ \hfill & \hfill & =-\text{k}T\frac{\text{d}}{\text{d}r}(\ln n_s)+{\theta }_s(g_s,T),\hfill \end{array}$$(17)

where $$g_s=g_{eff}-A_{R,s},$$

and $${\theta }_s(g_s,T)=-{\mu }_s{\text{M}}_{\text{u}}{g}_s-\text{k}T\frac{\text{d}}{\text{d}r}(\ln T).$$(18)

From Eqs. (17) and (18), we see that diffusion velocities, residual heat flows, and electric field have components due to gravity, temperature gradient, and number density gradients; for this reason, we call the diffusion problem gravitational settling and thermal and chemical diffusion. Hereafter, we rename the components due to gravity and temperature gradient as the gravitational settling component, and those due to the gradients in number densities as the chemical diffusion components. According to this, we separate out w_s by setting $$w_s=w_s^g-{\sum }_{t=1}^{N_{\text{ion},\text{0}}}{\xi }_{st}\frac{\text{d}}{\text{d}r}(\ln n_t),$$(19)

where $w_s^g$ is the velocity related to the gravitational settling component. The sum extends over all the particles with the exception of electrons (N_ion,0 = N − 1).

We write a similar expression for the residual heat flows r_s as $$r_s=r_s^g-{\sum }_{t=1}^{N_{\text{ion},\text{0}}}{\eta }_{st}\frac{\text{d}}{\text{d}r}(\ln n_t),$$(20)

and for E, the electric field, $$E=E^g-{\sum }_{t=1}^{N_{\text{ion},\text{0}}}{\zeta }_t\frac{\text{d}}{\text{d}r}(\ln n_t).$$(21)

Similarly, for the heat flow equation, the Eq. (14) can be written as $$\frac{1}{n_s}\{\frac{5}{2}{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t}{m_{st}}(w_t-w_s)+\frac{2}{5}{K}_{ss}z\text{\_}ss\text{^{''}}{r}_s+[{\sum }_{t\ne s,t=1}^N\frac{K_{st}}{m_{st}^2}[3m_s{}^2+m_t{}^{2}z\text{\_}st\text{^{'}}+\frac{4}{5}{m}_s{m}_{t}z\text{\_}st\text{^{''}}]]r_s-{\sum }_{t\ne s,t=1}^N{K}_{st}\frac{m_s{m}_t}{m_{st}^2}[3+z\text{\_}st\text{^{'}}-\frac{4}{5}z\text{\_}st\text{^{''}}]r_t\}=\tilde{\theta }(T),$$(22)

where $$\tilde{\theta }(T)=-\frac{5}{2}\text{k}T\frac{\text{d}}{\text{d}r}(\ln T).$$(23)

We eliminate the electron velocity, w_e, using Eq. (16) $$w_e=\frac{1}{n_e}{\sum }_{s=1}^{N_{\text{ion},\text{0}}}{Z}_s{n}_s{w}_s$$(24)

in Eqs. (15), (17), and (22).

3.2 The equation of continuity

Once we obtained the quantities w_s, r_s and E from Eqs. (19)–(21), we need to resolve the equation of continuity for the s-species, that is $$\frac{\partial {\rho }_s}{\partial t}+\nabla \cdot {\Phi }_s=0,$$(25)

where Φ_s is the mass flow vector defined by $${\Phi }_s={\rho }_s{w}_s=n_s{\mu }_s{M}_u{w}_s.$$(26)

In spherical symmetry this equation becomes $$\frac{\partial n_s}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}(r^2{n}_s{w}_s)=0,$$(27)

where only ions and neutral atoms are taken into account. From Eqs. (19) and (27), we have $$\frac{\partial n_s}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}[r^2{n}_s(w_s^g-{\sum }_{t=1}^{N_{\text{ion},\text{0}}}{\xi }_{st}\frac{\partial }{\partial r}(\ln n_t))]=0,$$(28)

which is solved for any arbitrary number of atoms and ions together with two boundary conditions.

3.3 Boundary conditions

In the outermost part of the atmosphere we assume that the numerical density is close to zero, that is $$n_s{|}_{r_{\text{e}}}=0.$$

Although this value is precisely assigned to the border, a more realistic condition is that $$\frac{\text{d}{x}_s}{\text{d}r}{|}_{r_{\text{e}}}=0,$$(29)

where x_s is the mass fractional abundance, which has null derivatives at the edge to avoid singularities. In this case, r_e is the external radius.

At the base of the atmosphere the boundary condition for the number densities is given by $$\frac{\text{d}{n}_s}{\text{d}r}{|}_{r_{\text{i}}}=0,$$(30)

where r_i is the internal radius. This condition is valid because we expect that the diffusion process that takes place in a stellar atmosphere has a timescale that is much shorter than the timescale for the stellar evolution. In Eqs. (29) and (30) the subscript s takes values from 1 to N_ion,0.

3.4 The solution

To solve the system given by Eqs. (15), (17), and (22), we must first compute the unknown quantities w_s, r_s, and E, that are obtained from Eqs. (19)–(21). The solution matrix scheme for this linear equation system is developed in Appendix B. Once the w_s values are obtained (with s = 1, ..., N − 1), these quantities are replaced in Eq. (27), which are discretised following Iben & MacDonald (1985).

4 Results

We solved the time-dependent flow equations considering gravitational settling as well as thermal and chemical diffusion for atmospheresof B-type stars with and without magnetic fields. Our goal is to investigate the existence of local under- or over-abundances of helium on the stellar surface. To this end, we adopted initial number density distributions that result from atmosphericmodels calculated with a homogeneous chemical profile and then we computed the outcoming density profile when the diffusion process stops. To illustrate the results, we modelled, in particular, stars of spectral types B2 V (T_eff = 22 000 K and log g = 4) and B5 V (T_eff = 15 200 K and log g = 4).

To compute the atomic diffusion throughout the atmosphere, we adopted a mixture of neutral and partially ionised gases composed by hydrogen, helium, and electrons. Hydrogen and helium are assumed as effective atoms. The effective species are calculated using the concentration C_s for s-species. According to Schlattl (2002), this is written as $$C_s=\frac{n_s}{{\tilde{n}}_e},$$(31)

where ñ_e takes into account a possible degeneracy of the electron gas and $${\tilde{n}}_e=\frac{P}{kT}.$$(32)

However, in these stellar atmospheres there is no gas degeneration, then ñ_e = n_e.

Then the concentration C for all species is C = ∑_sC_s, and the weighted average of any variable W is given by $$\overline{W}=\frac{{\sum }_s{C}_s{W}_s}{C}.$$(33)

Thus, using Eqs. (31) and (32) we have $$\overline{W}=\frac{{\sum }_s{n}_s{W}_s}{{\sum }_s{n}_s}.$$(34)

To compute the atomic diffusion throughout the atmosphere we adopted the following components: H I, H II, He I, He II, He III, and electrons. Then, we defined the effective mass and charge of a given element as $\overline{M}={\sum }_s{n}_s{M}_s/{\sum }_s{n}_s$ and $\overline{Z}={\sum }_s{n}_s{Z}_s/{\sum }_s{n}_s$, respectively. Therefore, the charges are written as $$\overline{Z_{\text{ion}H}}=\frac{n_{\text{ion}Hii}}{n_{\text{ion}Hi}+n_{\text{ion}Hii}}\text{and}$$(35) $$\overline{Z_{\text{ion}He}}=\frac{n_{\text{ion}Heii}+2 n_{\text{ion}Heiii}}{n_{\text{ion}Hei}+n_{\text{ion}Heii}+n_{\text{ion}Heiii}}.$$(36)

4.1 Diffusion in absence of magnetic fields

Diffusion in non-magnetic stars was studied by Michaud et al. (1976) and Richer et al. (2000). In this context, we verify the correct operation of the method we implemented.

In all our models, the initial atmosphere structure was calculated using the stellar atmospheric code developed by Rohrmann (2001) and Rohrmann et al. (2002). This code was adapted to include the Lorentz force (Vallverdú et al. 2014).

Firstly, we considered a stellar model in absence of a magnetic field (B = 0). The calculation begins using an atmospheric model with constant effective atoms as initial condition in t₀ = 0. Then, the calculation of diffusion was performed using the Burgers equations and the method of Iben & MacDonald (1985) until convergence of chemical distribution is achieved at t_f (final time). As a criteria, t_f is achieved when the relative variation in the mass fraction is less than 1% for all considered elements. The final time is approximately 1700 yr and 5000 yr for B2 V and B5 V star models, respectively.

In Fig. 1, we show the numerical density and mass fraction abundance changes of effective hydrogen and helium during the diffusion process as a function of the Rosseland optical depth for B2 V and B5 V star models, from time equal zero to t_f at intervals of 400 yr and 1000 yr, respectively. The X_i profiles of each model exhibit opposite behaviours because of different effective charges. In the B2 V star model hydrogen sinks down and the helium rises to the surface, contrary to the case observed in the B5 V star model. In the first star model, the mass fraction does not change with the optical depth (in a logarithmic scale) from − 2.5 to the innermost regions. In the B5 V model the mass fraction does not change between − 6.0 and − 1.5 and between 0.5 and 1.5.

4.2 Diffusion velocities in the presence of magnetic fields

When the magnetic field is considered the atmospheric structure is different. The Lorentz force in combination with a gravitational field can be treated as an effective gravity using the expression given by Eq. (8) with the two possible (inward or outward) directions of the Lorentz force.

To solve the problem of diffusion it is necessary to calculate the speeds of the particles in the mixture (see Eq. (27)). We tested the model by evaluating the temporal evolution of the diffusion velocity with depth. In order to show the speed behaviour, we calculated the velocities for a B2 V star model at the equator (θ = 90°) for a magnetic field of 1 kG and an inward directed Lorentz force.

In Fig. 2, we plot the diffusion velocities for H and He for a B2 V star model. A negative velocity means that a given particle sinks down towards the inner atmosphere, in this case helium atoms. Otherwise, ascending particles (hydrogen atoms) havea positive velocity. In the outer layers of the atmosphere, the velocities for both particles are a few centimetres per second.

Fig. 1 Top left panel: profiles of log Ni at different diffusion times as a function of the Rosseland optical depth for a B2 V star model in absence of magnetic fields (B = 0). Bottom left panel: mass fraction abundance vs. the Rosseland optical depth for the same elements of top left panel. Right panels: same as left panels but for a B5 V star model.

Fig. 2 Left panel: equatorial drift velocity for hydrogen at different times computed for a B2 V star model with an inward Lorentz force and B = 1 kG. Right panel:same as left panel but for helium.

Fig. 3 Initial(solid line) and final (dotted lines except for 10 G in short dashed lines) mass fraction abundances of H and He vs. the optical depth for a B2 V star model with the diffusion processes and different magnetic field strengths (10, 200, 500, 1000, and 2000 G). Left panels: the Lorentz force is inward directed (FL-in). Right panels: the Lorentz force is outward directed (FL-out).

Fig. 4 Same as in Fig. 3 but for a B5 V star model.

4.3 Diffusion with magnetic fields at the equator

For the B2 V star, we compute the diffusion at the magnetic equator with different effective gravities for the following magnetic fields 0, 10, 200, 500, 1000, and 2000 G. Figure 3 shows final abundance distributions for H and He along the atmosphere as a function of the Rosseland optical depth. The abundance profile shows no appreciable change in optical depth between − 2.5 and 1.5, for both Lorentz force directions. When the Lorentz force is inwardly directed (F_L-in), the greatest abundance variation is obtained with the smallest magnetic field. While when the Lorentz force is outwardly directed (F_L-out), the lowest field strength produces the smallest change in abundance.

Figure 4 shows the mass fraction abundance distribution for a B5 V star. The final distributions of H and He are almost similar regardless of the magnetic field strength and the Lorentz force direction.

4.4 Dependence with the co-latitude angle

In Fig. 5 we plot the dependence of the effective gravity with the optical depth and latitude. For the B2 V star, the maximum variation of log g between the equator and the poleis 0.2 dex and 0.4 dex, respectively,when the Lorentz force is inward and outward directed. The slope of the effective gravity in the outermost regions has a relevant effect on the variation of the mass fraction abundance as described below. In the B5 V star model, the maximum variation is 0.1 (for F_L-in) and 0.15 (for F_L-out), and the slope of the effective gravity in the innermost regions rules out the variation of each element.

Figures 6 and 7 illustrate, respectively, the mass fraction distribution X_i and N_He∕N_H ratio for B2 V and B5 V star models with a magnetic field of 1 kG, after diffusion. The figures show the results for three latitude angles: the magnetic pole (θ = 0°), where the effect of the diffusion is the same as in absence of a magnetic field (B = 0); an intermediate latitude (θ = 45°); and the equator (θ = 90°). The diffusion is slightly affected by the strength and the direction of the Lorentz force.

In all cases, mass fraction abundance profiles present an inversion at a given depth. For the B2 V star model with the Lorentz force inward directed, we observe an increase in He in the outer atmospheric layers towards the equator (for optical depths between − 6.0 and − 5.0) but this value decreases between − 4.0 and − 3.0. The opposite occurs when the Lorentz force is outward directed and the difference between the pole and the equator is more marked. The hydrogen have a kind of mirror behaviour. The B5 V star model also presents an inversion in the profile with depth (log (τ) ~−0.25), which is more pronounced when the Lorentz force is inwards directed.

Fig. 5 Behaviour of effective gravity with the co-latitude angle.

Fig. 6 Variations of H and He abundances with the optical depth and the co-latitude angles (0° at the pole, 45°, and 90° at the equator in green, blue, and red, respectively) for the two directions of the Lorentz force. The models correspond to B2 V (top panels) and B5 V (bottom panels) stars with B = 1 kG. The left (right) plots are for an inward (outward) directed Lorentz force.

Fig. 7 Variations of He numerical density relative to H with the optical depth and the co-latitude angles (0° at the pole, 45°, and 90° at the equator in green, blue, and red, respectively) for B2 V and B5 V star models in both directions of Lorentz force with B = 1 kG.

4.5 Chemical concentration with latitude at different depth

We now analyse the profiles of number density of atoms or ions of i species, obtained with Saha ionisation equation, with latitude at different optical depths. In Fig. 8, we plot log N_i versus log τ for a B2 V star model, where N_i is the numerical chemical abundance of the species i (e.g. H II, He I, and He II). The magnetic field strength is 1 kG and the two directions of the Lorentz force are considered. As we can see, proton abundance increases monotonically inwards the atmosphere at all latitudes. Something similar occurs for He II, except at the base of atmosphere (about log τ = 0.8), where the abundance falls inside. For He I, there is a maximum about log τ = −1.2. In all cases, we plot curves for different co-latitudinal angles with strong increasing or decreasing gradient in latitude at a given depth. This latitudinal dependence could show different behaviours at different optical depths for a given chemical species.

In Figs. 9–11, we illustrate in false colours log N_i vs. cos θ for a B2 V star model for H II, He I, and He II at two different optical depths: log τ = 0.0 and − 3.0, showing the effects of the abundance gradients. In these cases, the Lorentz force is inward directed. The abundances decrease towards the magnetic poles (cos θ = 1) for the three chemical elements at the mentioned optical depths.

The gradients in false colour depend on the slope of the curve in the diagram log N_i versus cos θ. In the successive stellar disc graphs (Figs. 12, 14 and 15), we employ the same false colours to represent approximately the same variation interval of number density (Δ log N_i).

When the Lorentz force is outward-directed the behaviour of each element is opposite to the previous case (inward-directed), as expected. As an example we illustrate the gradient of He I with latitude (for log τ = 0.0 and − 3.0), which leads to a concentration of this element at the poles (see Fig. 12).

Hereunder, we analyse the variations of log N_i with log τ and cos θ for H II and He I, corresponding to a B5 V star model with a 1 kG magnetic field and both directions of the Lorentz force f (see Figs. 13). We observe that the H II number density increases inward in the atmosphere for both directions of the Lorentz force, except for about log τ = −0.4. However, the abundance of He I grows inward until log τ ≈−0.4 and then falls down (see Fig. 13, bottom panels), independently of the direction of the Lorentz force.

With an inward Lorentz force, the number densities of H II and He I increase from the magnetic pole to the equator for all the optical depths. At log τ = −1.0 the number density contrast between the equator and pole N_i, e∕N_i, p ~ 2.5. This value is greater than that of the B2 V model, which is ~1.5. Figure 14 illustrates these behaviours for H II and He I at log τ = 0.5 and − 1.0, respectively. Instead, the number densities decrease from the magnetic pole to the equator for the same elements when the Lorentz force is in the opposite direction.

In both star models when the Lorentz force is directed inward, the hydrogen and helium distributions resemble a ring structure. Opposite behaviours occur when the Lorentz force is directed outward. In this case we would observe polar spots.

As a consequence, depending on the orientation of the magnetic axis with respect to the rotation axis, helium lines that form at different regions would show different behaviours during a rotation period. In Fig. 16, we depict six different phases for a ring-shaped density distribution. These 3D plots show density changes that depend on the orientation of the magnetic axis with respect to the rotation axis.

Fig. 8 Profiles of log Ni vs. log τ for following species H II, He I, and He II in two atmospheric models of a B2 V star with a magnetic field of B = 1 kG. Left panels: inward directed Lorentz force; right panels: outward force. Each panel shows different latitudes.

Fig. 9 For a B2 V star model with a magnetic field of B = 1 kG for inward directed Lorentz force. Right panels: log Ni vs. cos (θ) (θ, co-latitudinal angle) for H II are plotted at two different profundities, log τ = 0.0 and log τ = -3.0. This ion shows a concentration towards the equatorial plane. Left panels: the stellar discs are represented in the corresponding depths in false colours.

Fig. 10 Same as in Fig. 9 for He I at the same depths. There is an increment of the He I number density towards the equatorial plane.

Fig. 11 Same as in Fig. 9 for He II at the same depths.

Fig. 12 Same as in Fig. 10 but for an outward Lorentz force.

Fig. 13 Profiles of log Ni vs. log τ for the following species: H II and He I, for a B5 V star model with a magnetic field of B = 1 kG for an inward and outward directed Lorentz force (left and right panels, respectively).

5 Discussion and conclusions

It is known that helium deficiency in absence of magnetic fields is naturally explained by diffusion processes. In addition, Michaud et al. (1979) also confirmed that radiation force is not able to support helium atoms that would sink into the atmosphere in all the types of stars. As a consequence, the helium depletion modifies the helium ionised regions, thereby also affecting the pulsation activity of the star. On the other hand, one would expect that magnetic fields could reduce the sinking, and the mass loss could also prevent helium from sinking.

In this work we explore in more detail the diffusion process in helium-peculiar stars with magnetic fields. To this end we considered the ORM with a dipolar magnetic field and transport by atomic diffusion throughout the stellar atmosphere. We analysed the influence of the magnetic field and diffusion in the atmospheric structure of B-type stars and its impact in the formation regions of helium lines.

We assumed that the centre of the star and the bulk composition of the entire star have normal chemical abundance mixtures that reflect the composition of the gas clouds from which they are formed. On the other hand, we neglected both convective mixing and meridional mass flux, hence it is assumed that the stability is provided by the magnetic field. As a consequence, our model assumes a magnetospheres in hydrostatic equilibrium, in which the combined effect of a gravitational field and the Lorentz force generates an effective gravity. This effective gravity could significantly increase the efficiency of gravitational settling at different magnetic latitudes, leading to vertical and horizontal stratifications.

To model the atomic diffusion we wrote a code based on the scheme developed by Iben & MacDonald (1985) that solves the multi-component flow equations, which take into account gravitational settling as well as chemical and thermal diffusion, following the theory established by Burgers (1969). In this work, we compute atomic diffusion for effective hydrogen, effective helium, and electrons. After diffusion process the different ionisation stages are obtained from both effective elements.

These models were calculated for dwarf stars with temperatures of 15 000 and 22 000 K and magnetic fields between 0 and 2000 G. We find that the combined effect of gravitational settling, thermal and chemical diffusion with a magnetic field allows for local effective helium over-abundance at surface for a B2 V and a helium depletion in a B5 V. This result is in agreement with what it is expected in He-strong and He-weak stars. Although the models explain variations of chemical abundances with depth and latitude, the contribution of the diffusion to helium line profiles is small. In a B2 V star model the resulting line EW increases around 1%, as shown in Fig. 17. This value is lower than the effect produced by the orientation of the Lorentz force (Vallverdú et al. 2014). Therefore, we conclude that the present model approximation does not seem to give a satisfactory explanation to the abnormal strong or weak line strengths observed in these stars. For instance, typical variations in the helium line are larger than 20% (Rivinius et al. 2011; Shultz et al. 2015).

The present model assumes the diffusion of effective particles, as stated from Eqs. (34) to (36), and this hypothesis might not be completely accurate. A more real model would be to treat of diffusion of real charges and consider that the term on the right-hand side of Eq. (27) is non-null. This implies the description of collisional and radiative interactions or charge exchange mechanisms that it would be important when studying diffusion processes that involve silicon and helium particles. On the other hand, it is important to study this process in combination with the contribution of stellar winds. However, this study is out of scope of the present work.

Our main results can be summarised as follows:

• The combined effect of gravitational settling as well as thermal and chemical diffusion with a magnetic field allows for a local effective helium over-abundance at surface for a B2 V (prototype of a He-strong star), while in a B5 V (prototype of a He-weak) the effective helium sinks down.

• Furthermore, the variations in the number density of atoms (or ions) between the magnetic pole and equator depend mainly on the direction of the Lorentz force. This effect leads to under- or over-abundances giving the appearance of rings (equator) or spots (pole).

• We find that the latitudinal abundance gradient (Δlog N_i) depends on the optical depth. However, the change in the helium line EW between the pole and the equator that results from diffusion processes is around of 1%.

• The abundance distribution for a dipolar magnetic configuration inclined a certain angle with the rotation axis might produce changes with the rotation phase due to local enhancement or depletion of helium with latitude.

Our next step is to explore the diffusion process in other atomic species. To this goal, we must improve the model by including radiative accelerations, as they cannot be neglected in anyway.

Fig. 14 For aB5 V star model with a magnetic field of B = 1 kG for inward directed Lorentz force. Right panels: log Ni vs. cos (θ) (θ, co-latitudinal angle) are plotted for H II at log τ = 0.5 and for He I at log τ = − 1.0. Left panels: the stellar discs are represented by the corresponding elements at the same depths in false colors.

Fig. 15 Same as in Fig. 14 but for an outward directed Lorentz force. This configuration produces polar spots.

Fig. 16 Different phases for a ring-shaped density distribution in 3D. In all cases, the Z-axis emerges, the direction of X-axis goes from left to right (equal to direction of ω-vector), and Y -axis goes from down to up. The rotational axis ω is indicated with red arrows and the magnetic axis is indicated with blue arrows. With (α, β, γ) the rotation in degrees around axis X, Y, and Z, respectively, is shown. Top panel: (270, 0, 0), (300, 0, 0), and (330, 0, 0). Bottom panel: (0, 0, 0), (30, 0, 0), and (60, 0, 0).

Fig. 17 B2 V star model: EWs of He I λ5875 vs. latitude (0 is the magnetic pole and π∕2 is the equator) calculated with (blue circles) and without diffusion (red circles). The filled (open) circles denote the case when the Lorentz force is inward (outward) directed.

Acknowledgements

We would like to thank our anonymous referee for his or her valuable comments and suggestions. We would also like to thank René Rohrmann for his suggestions. This research was funded by grants from CONICET (PIP N° 0177), Agencia Nacional de Promoción Científica y Técnica (préstamo BID, PICT 2016-1971) and the Universidad Nacional de La Plata (Proyecto I+D 11/G162).

Appendix A: Diffusion equation

Here we demonstrate that Eq. (1) represents a set of N − 1 linearly independent equations and not N. The N^th equation is a linear combination of the N − 1 remaining equations. We assumed that we have N-equations and if we add all together, we reach an identity of form zero equal zero.

Let us sum from s = 1 to N in Eq. (1) as follows: $$\begin{array}{lll}\hfill & \hfill & {\sum }_{s=1}^N\nabla p_s-{\sum }_{s=1}^N\frac{{\rho }_s}{\rho }\nabla p-{\sum }_{s=1}^N{n}_s{Z}_{s}eE-{\sum }_{s=1}^N(j_s-\frac{{\rho }_s}{\rho }j)\times B/c\hfill \\ \hfill & \hfill & ={\sum }_{s=1}^N{\sum }_{t\ne s,t=1}^N{K}_{st}(w_t-w_s)+{\sum }_{s=1}^N{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}}, \hfill \end{array}$$(A.1)

then $$\begin{array}{lll}\hfill & \hfill & {\sum }_{s=1}^N\nabla p_s-[{\sum }_{s=1}^N{\rho }_s]\frac{1}{\rho }\nabla p-[{\sum }_{s=1}^N{n}_s{Z}_{s}e]E\hfill \\ \hfill & \hfill & -[{\sum }_{s=1}^N{j}_s]\times B/c+[{\sum }_{s=1}^N{\rho }_s]\frac{1}{\rho }j\times B/c\hfill \\ \hfill & \hfill & ={\sum }_{s=1}^N{\sum }_{t\ne s,t=1}^N{K}_{st}(w_t-w_s)+{\sum }_{s=1}^N{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}}, \hfill \end{array}$$(A.2)

where $$\nabla p={\sum }_{s=1}^N\nabla p_s,\rho ={\sum }_{s=1}^N{\rho }_s,j={\sum }_{s=1}^N{j}_s,$$

and the net charge is null $${\sum }_{s=1}^N{n}_s{Z}_{s}e=0.$$

Hence, $${\sum }_{s=1}^N{\sum }_{t\ne s,t=1}^N{K}_{st}(w_t-w_s)+{\sum }_{s=1}^N{\sum }_{t\ne s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}}=0.$$(A.3)

In this equation, each double sum gives a number of terms equal to V (N, 2), that is variations of N-elements taken 2; since the order in which thesubscripts s and t are written matters (that is demonstrated by inductionmethod). Then, V (N, 2) = N(N − 1), where N(N − 1) is an even number. Thus, we can write each double sum in two parts of N(N − 1)∕2 terms. Then, $$\begin{array}{lll}\hfill & \hfill & {\sum }_{s=1}^N{\sum }_{t>s,t=1}^N{K}_{st}(w_t-w_s)+{\sum }_{s=1}^N{\sum }_{t<s,t=1}^N{K}_{st}(w_t-w_s)\hfill \\ \hfill & \hfill & +{\sum }_{s=1}^N{\sum }_{t>s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}}\hfill \\ \hfill & \hfill & +{\sum }_{s=1}^N{\sum }_{t<s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}}=0.\hfill \end{array}$$(A.4)

Because, s y t are dummy indices, then we can swap them, and K_st = K_ts, z_st = z_ts and m_st = m_ts are symmetrical. Then, $$\begin{array}{lll}\hfill & \hfill & {\sum }_{s=1}^N{\sum }_{t>s,t=1}^N{K}_{st}(w_t-w_s)+{\sum }_{t=1}^N{\sum }_{s<t,s=1}^N{K}_{ts}(w_s-w_t)\hfill \\ \hfill & \hfill & +{\sum }_{s=1}^N{\sum }_{t>s,t=1}^N{K}_{st}{z}_{st}\frac{m_t{r}_s-m_s{r}_t}{m_{st}}\hfill \\ \hfill & \hfill & +{\sum }_{t=1}^N{\sum }_{s<t,s=1}^N{K}_{ts}{z}_{ts}\frac{m_s{r}_t-m_t{r}_s}{m_{ts}}=0,\hfill \end{array}$$(A.5)

hence we obtain an identity zero equal zero, QED.

Appendix B: The solution vector

As was stated, the number of different components: electrons plus N_ion,0 (ions and neutral atoms) is N = N_ion,0 + 1.

Then, we have a system of linear equations represented by the following matrix equation: $$M^{2N\times 2N}\cdot Z^{2N\times 1}=T^{2N\times 1},$$

where the matrix M and the vector T have the structures shown in Fig. B.1 and both can be partitioned in blocks. The solution vector Z is given by $$Z={(w_i,\mathrm{...},w_{N-1},r_i,\mathrm{...},r_N,E)}^t.$$

Hereinafter, we denote N_ion,0 = N − 1 and use the subscript e to represent the electrons (the element N).

The different matrix coefficients are written as follows:

Block V: $$V_{ss}=\frac{1}{e{n}_s}{\sum }_{t\ne s,t=1}^{N-1}{K}_{st}+\frac{K_{se}}{e{n}_s}-\frac{K_{se}{Z}_s}{e{n}_e},$$

for s = 1, ..., N − 1,

and $$V_{st}=-\frac{1}{e{n}_s}{K}_{st}-\frac{K_{se}{Z}_t}{e{n}_e}\frac{n_t}{n_s},$$

for t≠s, and s, t = 1, ..., N − 1.

Block Q: $$Q_{s,N-1+s}=-\frac{1}{e{n}_s}{\sum }_{t\ne s,t=1}^{N-1}{K}_{st}{Z}_{st}\frac{m_t}{m_{st}},$$

for s = 1, ..., N − 1,

and $$Q_{s,N-1+t}=\frac{1}{e{n}_s}{K}_{st}{Z}_{st}\frac{m_s}{m_{st}},$$

for t≠s, s = 1, ..., N − 1 and t = 1, ..., N.

Block F: $$F_{s,2N}=-Z_s,$$

for s = 1, ..., N − 1.

Block W: $$W_{N-1+s,s}=-\frac{5}{2}\frac{1}{n_s}{\sum }_{t\ne s,t=1}^{N-1}{K}_{st}{Z}_{st}\frac{m_t}{m_{st}}-\frac{5}{2}\frac{1}{n_s}{K}_{se}{Z}_{se}\frac{m_e}{m_{se}},$$

for s = 1, ..., N − 1, with m_se = m_s + m_e,

and $$W_{N-1+s,t}=\frac{5}{2}\frac{1}{n_s}{K}_{st}{Z}_{st}\frac{m_t}{m_{st}}+\frac{5}{2}\frac{1}{n_s}{K}_{se}{Z}_{se}\frac{m_e}{m_{se}}\frac{1}{n_e}{Z}_t{n}_t,$$

for t≠s, and s, t = 1, ..., N − 1.

Block R: $$R_{N-1+s,N-1+s}=\frac{2}{5}\frac{1}{n_s}{K}_{ss}{Z^{″}}_{ss}$$ $$+\frac{1}{n_s}{\sum }_{t\ne s,t=1}^N{K}_{st}[3\frac{m_s^2}{m_{st}^2}+\frac{m_t^2}{m_{st}^2}{Z^{\prime }}_{st}+\frac{4}{5}\frac{m_s{m}_t}{m_{st}^2}{Z^{″}}_{st}],$$

for s = 1, ..., N − 1,

and $$R_{N-1+s,N-1+t}=-\frac{1}{n_s}{K}_{st}\frac{m_s{m}_t}{m_{st}^2}[3+{Z^{\prime }}_{st}-\frac{4}{5}{Z^{″}}_{st}],$$

for t≠s, and s = 1, ..., N − 1 and t = 1, ..., N.

Block X: $$X_{2N-1,t}=-\frac{5}{2}\frac{1}{n_e}\{Z_t\frac{n_t}{n_e}[{\sum }_{u=1}^{N-1}{K}_{eu}{Z}_{eu}\frac{m_u}{m_{ue}}]-K_{et}{Z}_{et}\frac{m_t}{m_{et}}\},$$

for t = 1, ..., N − 1.

Block S: $$S_{N,t}=-\frac{1}{n_e}{K}_{et}\frac{m_e{m}_t}{m_{et}^2}[3+{Z^{\prime }}_{et}-\frac{4}{5}{Z^{″}}_{et}],$$

for t = 1, ..., N − 1.

and $$S_{N,N}=\frac{2}{5}\frac{1}{n_e}{K}_{ee}{Z^{″}}_{ee}+\frac{1}{n_e}{\sum }_{t=1}^{N-1}{K}_{et}[3\frac{m_e^2}{m_{et}^2}+\frac{m_t^2}{m_{et}^2}{Z^{\prime }}_{et}+\frac{4}{5}\frac{m_e{m}_t}{m_{et}^2}{Z^{″}}_{et}],$$

for t = 1, ..., N − 1.

Block Y: $$Y_{2N,t}=({\mu }_t+{\mu }_e{Z}_t)n_t,$$

for t = 1, ..., N − 1.

Block A: $$A_s=-{\mu }_s\frac{M_{\text{u}}{g}_s}{e}-\frac{\text{k}T}{e}\frac{\text{d}}{\text{d}r}(\ln T),$$

for s = 1, ..., N − 1.

Block B: $$B_{N-1+s}=-\frac{5}{2}\text{k}T\frac{\text{d}}{\text{d}r}(\ln T),$$

for s = 1, ..., N − 1.

Block C: $$C_{2N-1}=-\frac{5}{2}\text{k}T\frac{\text{d}}{\text{d}r}(\ln T).$$

Block D: $$D_{2N}=0.$$

Calculation of the vector Z

To obtain the solution of system, Eqs. (19)–(21) can be written as follows: $$w_s=w_s^g+{\tilde{w}}_s,$$ $$r_s=r_s^g+{\tilde{r}}_s,$$

and $$E=E^g+\tilde{E},$$

where $${\tilde{w}}_s=-{\sum }_{t=1}^{N_{\text{ion},\text{0}}}{\xi }_{st}\frac{\text{d}}{\text{d}r}(\ln n_t),$$(B.1) $${\tilde{r}}_s=-{\sum }_{t=1}^{N_{\text{ion},\text{0}}}{\eta }_{st}\frac{\text{d}}{\text{d}r}(\ln n_t),$$(B.2)

and $$\tilde{E}=-{\sum }_{t=1}^{N_{\text{ion},\text{0}}}{\zeta }_t\frac{\text{d}}{\text{d}r}(\ln n_t).$$(B.3)

Fig. B.1 Linear system of equations, defined by M2N×2N ⋅ Z2N×1 = T2N×1.

Firstly, to calculate $w_s^g$, $r_s^g$, and E^g, we set $\frac{\text{d}}{\text{d}r}(\ln n_t)=0$ and obtain $w_s=w_s^g$, $r_s=r_s^g$ and E = E^g by solving the linear system.

Then, toderive the coefficients ξ_st, η_st, and ζ_st, we must set the gravitational terms equal zero: $w_s^g=0$, $r_s^g=0$, and E^g = 0. Furthermore, we must vanish the thermal terms to obtain particular solutions for the diffusion processes (see Eqs. (18) and (23)) as follows: $${\theta }_s(g_s,T)=-{\mu }_s{\text{M}}_{\text{u}}{g}_s-\text{k}T\frac{\text{d}}{\text{d}r}(\ln T)=0,$$(B.4)

and $$\tilde{\theta }(T)=-\frac{5}{2}\text{k}T\frac{\text{d}}{\text{d}r}(\ln T)=0.$$(B.5)

This way, we keep the same coefficient matrix M^2N×2N but with a different single column matrix. Then the coefficients ξ_st, η_st, and ζ_st can be derived one by one by setting for each element t, $\frac{\text{d}}{\text{d}r}(\ln n_1)=1$, and the rest $\frac{\text{d}}{\text{d}r}(\ln n_t)=0$, for t = 2, ..., N − 1. From Eq. (B.1) we now have $${\tilde{w}}_s=-{\xi }_{s1}={\tilde{\xi }}_{s1},$$

for s = 1, ..., N − 1.

Proceeding in the same way, we obtain $${\tilde{r}}_s=-{\eta }_{s1}={\tilde{\eta }}_{s1},$$

for s = 1, ..., N − 1, and $$\tilde{E}=-{\zeta }_1={\tilde{\zeta }}_1.$$

Hence, we solve the system $$M^{2N\times 2N}\cdot {\tilde{Z}}_1^{2N\times 1}={\tilde{T}}_1^{2N\times 1},$$

where $${\tilde{Z}}_1={({\tilde{\xi }}_{11},{\tilde{\xi }}_{21},\mathrm{...},{\tilde{\xi }}_{N-1,1},{\tilde{\eta }}_{11},{\tilde{\eta }}_{21},\mathrm{...},{\tilde{\eta }}_{N1},{\tilde{\zeta }}_1)}^t,$$

and $${\tilde{T}}_1={(-\frac{kT}{e},0,\mathrm{...},0)}^t.$$

The vector ${\tilde{T}}_1$ has 2N − 1 zeros.

Similarly, for the coefficients of the second element (t = 2) ξ_s2, η_s2, and ζ₂ we set $\frac{\text{d}}{\text{d}r}(\ln n_2)=1$ and $\frac{\text{d}}{\text{d}r}(\ln n_t)=0$ for t = 1, ..., N − 1 and t≠2. Thus, $${\tilde{w}}_s=-{\xi }_{s2}={\tilde{\xi }}_{s2},$$

for s = 1, ..., N − 1, $${\tilde{r}}_s=-{\eta }_{s2}={\tilde{\eta }}_{s2},$$

for s = 1, ..., N − 1, and $$\tilde{E}=-{\zeta }_2={\tilde{\zeta }}_2.$$

The system to solve is given by $$M^{2N\times 2N}\cdot {\tilde{Z}}_2^{2N\times 1}={\tilde{T}}_2^{2N\times 1},$$

where $${\tilde{Z}}_2={({\tilde{\xi }}_{12},{\tilde{\xi }}_{22},\mathrm{...},{\tilde{\xi }}_{N-1,2},{\tilde{\eta }}_{12},{\tilde{\eta }}_{22},\mathrm{...},{\tilde{\eta }}_{N2},{\tilde{\zeta }}_2)}^t$$

and $${\tilde{T}}_2={(0,-\frac{kT}{e},0,\mathrm{...},0)}^t.$$

Again the vector ${\tilde{T}}_2$ has 2N − 1 zeros.

This procedure is repeated to calculate the coefficients for the rest of the elements.

References

All Figures

	Fig. 1 Top left panel: profiles of log Ni at different diffusion times as a function of the Rosseland optical depth for a B2 V star model in absence of magnetic fields (B = 0). Bottom left panel: mass fraction abundance vs. the Rosseland optical depth for the same elements of top left panel. Right panels: same as left panels but for a B5 V star model.
In the text

	Fig. 2 Left panel: equatorial drift velocity for hydrogen at different times computed for a B2 V star model with an inward Lorentz force and B = 1 kG. Right panel:same as left panel but for helium.
In the text

	Fig. 3 Initial(solid line) and final (dotted lines except for 10 G in short dashed lines) mass fraction abundances of H and He vs. the optical depth for a B2 V star model with the diffusion processes and different magnetic field strengths (10, 200, 500, 1000, and 2000 G). Left panels: the Lorentz force is inward directed (FL-in). Right panels: the Lorentz force is outward directed (FL-out).
In the text

	Fig. 4 Same as in Fig. 3 but for a B5 V star model.
In the text

	Fig. 5 Behaviour of effective gravity with the co-latitude angle.
In the text

	Fig. 6 Variations of H and He abundances with the optical depth and the co-latitude angles (0° at the pole, 45°, and 90° at the equator in green, blue, and red, respectively) for the two directions of the Lorentz force. The models correspond to B2 V (top panels) and B5 V (bottom panels) stars with B = 1 kG. The left (right) plots are for an inward (outward) directed Lorentz force.
In the text

	Fig. 7 Variations of He numerical density relative to H with the optical depth and the co-latitude angles (0° at the pole, 45°, and 90° at the equator in green, blue, and red, respectively) for B2 V and B5 V star models in both directions of Lorentz force with B = 1 kG.
In the text

	Fig. 8 Profiles of log Ni vs. log τ for following species H II, He I, and He II in two atmospheric models of a B2 V star with a magnetic field of B = 1 kG. Left panels: inward directed Lorentz force; right panels: outward force. Each panel shows different latitudes.
In the text

	Fig. 9 For a B2 V star model with a magnetic field of B = 1 kG for inward directed Lorentz force. Right panels: log Ni vs. cos (θ) (θ, co-latitudinal angle) for H II are plotted at two different profundities, log τ = 0.0 and log τ = -3.0. This ion shows a concentration towards the equatorial plane. Left panels: the stellar discs are represented in the corresponding depths in false colours.
In the text

	Fig. 10 Same as in Fig. 9 for He I at the same depths. There is an increment of the He I number density towards the equatorial plane.
In the text

	Fig. 11 Same as in Fig. 9 for He II at the same depths.
In the text

	Fig. 12 Same as in Fig. 10 but for an outward Lorentz force.
In the text

	Fig. 13 Profiles of log Ni vs. log τ for the following species: H II and He I, for a B5 V star model with a magnetic field of B = 1 kG for an inward and outward directed Lorentz force (left and right panels, respectively).
In the text

	Fig. 14 For aB5 V star model with a magnetic field of B = 1 kG for inward directed Lorentz force. Right panels: log Ni vs. cos (θ) (θ, co-latitudinal angle) are plotted for H II at log τ = 0.5 and for He I at log τ = − 1.0. Left panels: the stellar discs are represented by the corresponding elements at the same depths in false colors.
In the text

	Fig. 15 Same as in Fig. 14 but for an outward directed Lorentz force. This configuration produces polar spots.
In the text

Fig. 16 Different phases for a ring-shaped density distribution in 3D. In all cases, the Z-axis emerges, the direction of X-axis goes from left to right (equal to direction of ω-vector), and Y -axis goes from down to up. The rotational axis ω is indicated with red arrows and the magnetic axis is indicated with blue arrows. With (α, β, γ) the rotation in degrees around axis X, Y, and Z, respectively, is shown. Top panel: (270, 0, 0), (300, 0, 0), and (330, 0, 0). Bottom panel: (0, 0, 0), (30, 0, 0), and (60, 0, 0).

In the text