1 Introduction

Dark matter is one of the key elements of modern scenaria of the Universe evolution, of the galactic formation and rotation [1,2,3]. Dark matter forms halos and subhalos, filaments and walls; their structures give a significant part of information about the gravitational field of the central bodies, which are surrounded by these galactic units. It is impossible to observe dark matter directly, since it neither emits, nor scatters light, that is why theoretical predictions based on the structure modeling of the dark matter halos, subhalos, filaments and walls play the important role in the investigations of these enigmatic cosmological units.

We follow the line that the axions, massive pseudo-Goldstone bosons, predicted in the works [4,5,6], are the particles, which are responsible for the phenomena, associated with the dark matter. During the last forty years the basic aspects of the axion theory, the cosmological and astrophysical applications of the axion models were well documented in many works (see, e.g., [7,8,9,10,11,12,13,14,15,16]). During the last decade new trends in the theory of the axionic dark matter appeared. First of all, we mean the non-minimal extensions of the axion models, which are based on the F(R) modifications of the theory of gravity (see, e.g., [17,18,19,20,21]), on the non-minimal extensions of the theory of the axion-photon coupling [22], and on the axionically extended models of the dynamic aether [23]. The second trend is associated with the recent surge of interest to the black hole mergers and supermassive black holes: the axionic dark matter can play the role of a marker revealing specific features of the strong gravitational field (see, e.g., [24, 25] as illustrations). The third trend is connected with the modeling of axionic structures of different scales, associated with the galactic halos, miniclusters, halos surrounding the magnetars, dyons, etc. (see, e.g., [26,27,28,29]). We have to emphasize, that one of the most influential idea in these new models is the idea of internal self-interaction in the axionic systems. The axionic dark matter, apparently, is not a simple collisionless gas without pressure, which interacts by the gravitational field only; probably, the axions form correlated systems. The most known model of this type is the model of axionic Bose–Einstein condensate, presented in [30]. The alternative description of correlated axionic systems is based on the models containing the potential of the pseudoscalar field \(V(\phi ^2)\), which can have the polynomial form \(\phi ^n\), the \(\phi ^4\) Higgs-type structure, or the nonlinear periodic form (see, e.g., [31]).

We focus on the description of the axionic systems based on the analysis of the equation for the pseudoscalar field \(\phi \), associated with axions [32, 33]. New aspects of this theory are connected with the structure of the pseudoscalar field potential V. The standard potential is considered as an even function of the argument \(\phi \) only, i.e., \(V=V(\phi ^2)\); we extend the theory by introducing the function of \(s+1\) arguments \(V(\phi ^2, \xi _{(1)},\xi _{(2)}, \ldots \xi _{(s)})\), where \(\left\{ \xi _{(a)} \right\} \) is the set of some scalars (\(a=1,2, \ldots s\)). What can be the origin of these scalars?

(i) Scalars associated with the Einstein–Aether theory.

The extension of the potential was prompted by the analogy with the Einstein–Aether theory [34,35,36], which is based on the introduction of a timelike dynamic vector field \(U^i\) normalized by unity (\(g_{ik}U^iU^k=1\)). In the framework of the Einstein–Aether theory one can introduce into the potential V four auxiliary scalars based on the decomposition of the covariant derivative \(\nabla _mU_n\) (see, e.g., [37,38,39,40] for details and references). The first scalar is the divergence of the vector field (\(\xi _{(1)} \rightarrow \Theta \equiv \nabla _m U^m\)). The second scalar is the square of the acceleration four-vector \(a^i = U^k \nabla _k U^i\) (\(\xi _{(2)} \rightarrow a^2 \equiv g_{ik} a^ia^k\)). The third scalar is the square of the symmetric shear tensor \(\sigma _{pq}\) (\(\xi _{(3)} \rightarrow \sigma ^2 \equiv \sigma _{ik} \sigma ^{ik}\)), and the fourth scalar is constructed using the square of the skew-symmetric vorticity tensor \(\omega _{ik}\) (\(\xi _{(4)} \rightarrow \omega ^2 \equiv \omega _{ik}\omega ^{ik}\)). The scalar coinciding with the modulus of the vector \(U^i\) is equal to one and thus it can not be considered as a guiding function. In other words, when the theory possesses intrinsic unit vector field, we can add four new arguments to the potential of the pseudoscalar field V, and can use this covariant extension for detailing the structure of the dark matter configuration.

There are three illustrations of this idea. The first one relates to the Friedmann type cosmology, for which \(a^i =0\), \(\sigma _{ik}=0\), \(\omega _{ik}=0\), and the only scalar \(\Theta \) proportional to the Hubble function H(t) (\(\Theta = 3H\)) is the nonvanishing guiding function (see, e.g., [41, 42]). The second illustration is connected with the pp-wave symmetric spacetimes with two auxiliary nonvanishing scalars: \(\Theta \) and \(\sigma ^2\) [43, 44]. The third example appeared in the Gödel type model, in the framework of which the nonvanishing parameter \(\omega ^2\) can be considered as the guiding function [45].

(ii) Scalars associated with the Killing vectors.

If the model does not contain vector field \(U^i\) associated with the dynamic aether, but possesses some specific symmetries, described by the set of Killing vectors and/or conformal Killing vectors \(\xi ^i_{(a)}\) (\(a=1,2, \ldots s\)), we suggest to use the new scalar quantities of two types. First, we can consider the moduli of the Killing vectors, \(\xi _{(a)}= \sqrt{|g_{mn} \xi ^m_{(a)}\xi ^n_{(a)}|}\) (or their combinations) as the auxiliary arguments of the modified axionic potential. It was impossible in the Einstein–Aether theory, since the modulus of the velocity field \(U^i\) is equal to one. The scalar quantities of the second type contain nonvanishing convolutions of the derivative \(\nabla _m \xi _{n}\). When we deal with the standard Killing vector, which satisfies the equation \(\nabla _m \xi _{n}{+}\nabla _n \xi _{m}=0\), the divergence \(\nabla _m \xi ^{m}\) is equal to zero, but for the so-called conformal Killing vector, which satisfies the equation \(\nabla _m \xi _{n}{+}\nabla _n \xi _{m}= \frac{1}{2} g_{mn} \nabla _s \xi ^{s}\), the scalar \(\nabla _s \xi ^{s}\) is nonvanishing and can be used for the physical modeling. In its turn, for the standard Killing vector the skew-symmetric quantity \(\frac{1}{2} \left[ \nabla _m \xi _{n}{-}\nabla _n \xi _{m}\right] = \nabla _m \xi _{n}\) is nonvanishing, and its square can be used in analogy with the square of vorticity tensor \(\omega _{mn}\) appeared in the Einstein–Aether theory. Illustration of this idea can be found in [28], where the model with the static spherically symmetric spacetime of the Reissner–Nordström type is analyzed.

(iii) Why the new guiding scalar functions are necessary for the extended analysis?

In order to describe the new units in the dark matter halos (filaments, walls, etc.), we need of covariant scheme to fix some lines and surfaces. For instance, when we deal with a static spherically symmetric gravitational field and try to separate one spatial domain from the other, we could use the well-known Heaviside function h; however, we can not insert into the Lagrangian the term \(h(r{-}r_*)\), since the difference of radial coordinates is non-invariant quantity. Nevertheless, we can use the term \(h(\xi -\xi _*)\), where \(\xi \) is the appropriate additional scalar, and \(\xi _*\) is its critical value; this scheme is the covariant one.

(iv) How the paper is organized?

In order to provide the self-consistency of the approach, which is based on the extended formalism of Killing vectors, we have to introduce the Killing vector field as the dynamic field, i.e., we have to add into the Lagrangian the scalar terms quadratic in the covariant derivative \(\nabla _m \xi _n\), and to develop the variation formalism in analogy with the formalism used in the Einstein–Aether theory. In Sect. 2.1 of this work we establish the strict covariant formalism justifying this idea in general case.

Since the additional scalars, based on the Killing vector fields, are designed to be used for description of the pseudoscalar (axion) field, we have to modify the axionic periodic potential and to derive the correspondingly modified master equations for the axion, electromagnetic and gravitational fields; in Sect. 2.2 we present these modified master equations and discuss the concept of the equilibrium state of the axionic subsystem.

In Sect. 3 we apply the formalism for description of two-level distribution of the axionic dark matter near the static spherically symmetric objects with magnetic and electric fields; the analysis of the obtained key equation is presented in Sect. 4. In Sect. 5 we discuss the features of the profiles of the axionic dark matter obtained in the framework of the suggested model.

2 Extended formalism, which includes Killing vector fields

2.1 Analogy with the Einstein–Aether theory

2.1.1 Action functional

The action functional of the Einstein–Aether theory is known to have the form

$$\begin{aligned} S_{(\mathrm{EA})}= & {} \int d^4 x \sqrt{-g} \ \frac{1}{2\kappa }\left\{ (R+2\Lambda ) + \lambda (U^mU_m-1) \right. \nonumber \\&\left. + K^{abmn}\nabla _a U_{m} \nabla _b U_n \right\} , \end{aligned}$$
(1)

where R is the Ricci scalar; \(\Lambda \) is the cosmological constant; \(U^i\) is the vector field, \(\nabla _k\) is the covariant derivative; \(K^{abmn}\) is the Jacobson’s constitutive tensor:

$$\begin{aligned} K^{abmn}= & {} C_1 g^{ab}g^{mn} {+} C_2 g^{am}g^{bn} \nonumber \\&{+} C_3 g^{an}g^{bm} {+} C_4 U^a U^b g^{mn}, \end{aligned}$$
(2)

which contains four phenomenological constants \(C_1\), \(C_2\), \(C_3\) and \(C_4\). The scalar quantity \(\lambda \) is the Lagrange multiplier; variation with respect to \(\lambda \) gives the constraint \(g_{mn}U^mU^n=1\); this quadratic relationship is, in fact, the normalization condition for the vector field.

When the vector field \(U^i\) is absent, but there exists the Killing vector field \(\xi ^i\), we suggest to work with the action functional similar to (1):

$$\begin{aligned} S_{(\mathrm{EK})}= & {} \int d^4 x \sqrt{-g} \ \frac{1}{2\kappa }\left\{ {(R+2\Lambda ) }\right. \nonumber \\&{+} {\tilde{\lambda }} \left[ \nabla _{(k}\xi _{m)} \nabla ^{(k}\xi ^{m)}{-} \frac{1}{4} (\nabla _n \xi ^n)^2 \right] \nonumber \\&\left. {+} {{\mathcal {K}}}^{abmn}\nabla _a \xi _{m} \nabla _b \xi _n \right\} , \end{aligned}$$
(3)

where the parentheses denote symmetrization: \(\nabla _{(i}\xi _{m)} = \frac{1}{2} \left( \nabla _{i}\xi _{m} {+} \nabla _{m}\xi _{i} \right) \), and the tensor \({{\mathcal {K}}}^{abmn}\)

$$\begin{aligned} {{\mathcal {K}}}^{abmn} = {{\mathcal {K}}}_1 g^{ab}g^{mn} + {{\mathcal {K}}}_2 g^{am}g^{bn} + {{\mathcal {K}}}_3 g^{an}g^{bm} \end{aligned}$$
(4)

contains only three phenomenological constants \({{\mathcal {K}}}_1\), \({{\mathcal {K}}}_2\) and \({{\mathcal {K}}}_3\) (we assume that the constitutive tensor \({{\mathcal {K}}}^{abmn}\) contains the metric tensor, but does not include the Killing vector itself). In other words, we consider the vector field \(\xi ^i\) as a dynamic quantity in analogy with the aether velocity field, however, instead of the algebraic constraint we use the differential one.

2.1.2 Variation with respect to \({\tilde{\lambda }}\) and Killing equations

Variation of the action functional (3) with respect to the Lagrange multiplier \({\tilde{\lambda }}\) gives the constraint

$$\begin{aligned} \nabla _{(k}\xi _{m)} \nabla ^{(k}\xi ^{m)} - \frac{1}{4} (\nabla _n \xi ^n)^2 = 0 , \end{aligned}$$
(5)

which can be rewritten as follows:

$$\begin{aligned} \left[ \nabla _{(k}\xi _{m)} {-} \frac{1}{4} g_{km} \nabla _n \xi ^n \right] \left[ \nabla ^{(k}\xi ^{m)} {-} \frac{1}{4} g^{km} \nabla _s \xi ^s \right] = 0 . \end{aligned}$$
(6)

This equation is satisfied, if

$$\begin{aligned} \nabla _{k}\xi _{m} {+} \nabla _{m}\xi _{k} = \frac{1}{2} g_{km} \nabla _n \xi ^n . \end{aligned}$$
(7)

Generally speaking, (7) can be classified as the sufficient but not the necessary condition for the fulfillment of the quadratic relationship (6), nevertheless, here we do not discuss this fine mathematical detail.

Clearly, (7) are the equations for the so-called conformal Killing vector, if \(\nabla _n \xi ^n \ne 0\), and for the standard Killing vector, if \(\nabla _n \xi ^n =0\) (see, e.g., [46] for details). Calculation of the divergence of the left and right-hand sides of (7) yields

$$\begin{aligned} \nabla _m \nabla ^m \xi _k + R_{kj} \xi ^j = - \frac{1}{2} \nabla _k (\nabla _n \xi ^n) . \end{aligned}$$
(8)

Also, when we deal with the standard Killing vector field (i.e., \(\nabla _k \xi ^k=0\)) we can use the relationships

$$\begin{aligned} \nabla _s \nabla _m \xi _k + R_{jskm} \xi ^j =0 , \end{aligned}$$
(9)

as the integrability condition of the first order of the Killing equation. Here, as usual, \(R_{kj}= R^m_{\ \ kmj}\) is the Ricci tensor, and \(R^i_{\ kmj}\) is the Riemann tensor.

2.1.3 Variation with respect to \(\xi ^i\)

When one deals with the Einstein–Aether theory, the variation with respect to the velocity four-vector \(U^i\) gives the master equation

$$\begin{aligned} \nabla _a \left[ K^{abjn} \nabla _b U_n \right] = \lambda U^j + C_4 U^b \nabla _bU_n \nabla ^j U^n . \end{aligned}$$
(10)

Thus, the normalization condition yields

$$\begin{aligned} \lambda = U_m \nabla _a \left[ K^{abmn} \nabla _b U_n \right] - C_4 U^b \nabla _bU_n \ U^a \nabla _a U^n . \end{aligned}$$
(11)

When we deal with the extended formalism, which includes the Killing vector field, variation of the functional (3) with respect to \(\xi ^j\) gives the equations

$$\begin{aligned} \nabla _a \left\{ {\tilde{\lambda }} \left[ \nabla ^{(a} \xi ^{m)}- \frac{1}{4} \nabla _n \xi ^n \right] + 2{{\mathcal {K}}}^{abjn} \nabla _b \xi _n \right\} =0 . \end{aligned}$$
(12)

Keeping in mind the relationships (4), (7) and (8) we can reduce these equations to the relations

$$\begin{aligned} \left( {{\mathcal {K}}}_3 {-} {{\mathcal {K}}}_1 \right) R_{js} \xi ^s {+} \left( {{\mathcal {K}}}_2 {+} {{\mathcal {K}}}_3 {-} \frac{1}{2}{{\mathcal {K}}}_1 \right) \nabla _j \nabla _m \xi ^m =0 , \end{aligned}$$
(13)

which does not include \({\tilde{\lambda }}\). The Eq. (13) is satisfied identically, when

$$\begin{aligned} {{\mathcal {K}}}_2 = - \frac{1}{2} {{\mathcal {K}}}_1, \quad {{\mathcal {K}}}_3 = {{\mathcal {K}}}_1 \end{aligned}$$
(14)

for arbitrary conformal Killing vector with \(\nabla _n \xi ^n \ne 0\). When we deal with the standard Killing vector, i.e., \(\nabla _n \xi ^n = 0\), the parameter \({{\mathcal {K}}}_2\) remains arbitrary and hidden.

2.1.4 Variation with respect to metric

When one deals with the Einstein–Aether theory the contribution of the vector field into the total stress-energy tensor is known to have the following form

$$\begin{aligned} T^{(\mathrm{U})}_{ik}= & {} \frac{1}{2} g_{ik} {{\mathcal {J}}}^{am} \nabla _a U_m\nonumber \\&{+}\nabla ^m \left[ U_{(i}{{\mathcal {J}}}_{k)m}\right] {-} \nabla ^m \left[ {{\mathcal {J}}}_{m(i} U_{k)} \right] {-} \nabla _m \left[ {{\mathcal {J}}}_{(ik)} U^m\right] \nonumber \\&+C_1\left[ (\nabla _mU_i)(\nabla ^m U_k) {-} (\nabla _i U_m \nabla _k U^m) \right] \nonumber \\&{+}C_4 (U^a \nabla _a U_i)(U^b \nabla _b U_k), \end{aligned}$$
(15)

where the following definition is used:

$$\begin{aligned} {{\mathcal {J}}}^{am} = K^{abmn} \nabla _b U_n. \end{aligned}$$
(16)

In order to obtain the similar formula for our extended model, let us find the corresponding analogies. We assume now that \({{\mathcal {K}}}_1 = {{\mathcal {K}}}_3\), \({{\mathcal {K}}}_2 = - \frac{1}{2} {{\mathcal {K}}}_1\), and obtain that

$$\begin{aligned}&{{\mathcal {K}}}^{abmn} = {{\mathcal {K}}}_1 \left( g^{ab}g^{mn} + g^{an}g^{bm} - \frac{1}{2} g^{am}g^{bn} \right) , \end{aligned}$$
(17)
$$\begin{aligned}&{{\mathcal {K}}}^{abmn} \nabla _b \xi _n = 2 {{\mathcal {K}}}_1 \left[ \nabla ^{(a} \xi ^{m)} - \frac{1}{4}g^{am} \nabla _n \xi ^n \right] , \end{aligned}$$
(18)
$$\begin{aligned}&{{\mathcal {K}}}^{abmn} \nabla _a \xi _m \nabla _b \xi _n = 2 {{\mathcal {K}}}_1 \nonumber \\&\quad \left[ \nabla ^{(a} \xi ^{m)} \nabla _{(a} \xi _{m)}{-} \frac{1}{4} \left( \nabla _n \xi ^n \right) ^2 \right] . \end{aligned}$$
(19)

This means that the corresponding term in the integral (3) can be rewritten as follows

$$\begin{aligned} \lambda ^{*} \left[ \nabla ^{(a} \xi ^{m)} \nabla _{(a} \xi _{m)} - \frac{1}{4} \nabla _n \xi ^n \right] , \quad \lambda ^{*} \equiv ({\tilde{\lambda }}+ 2{{\mathcal {K}}}_1) . \end{aligned}$$
(20)

Thus, the analogies are the following: first, we have to omit \(C_4\); second, the analog of the term \({{\mathcal {J}}}^{am}\) is symmetric and it vanishes on the solutions to the Killing equations. In other words, following this analogy we obtain the vanishing stress-energy tensor associated with the Killing vector field.

If we fulfil the direct variation of the action functional (3) with respect to metric, the corresponding effective stress-energy tensor \(T_{pq}^{(\mathrm{K})}\) associated with the contribution of the Killing vector field can be written as follows:

$$\begin{aligned} T_{pq}^{(\mathrm{K})}= & {} \frac{(-2)}{\sqrt{-g}}\frac{\delta }{\delta g^{pq}} \left\{ \sqrt{{-}g}\lambda ^{*} g^{ab} g^{mn} \right. \nonumber \\&\times \left[ \nabla _{(a} \xi _{m)} {-} \frac{1}{4} g_{am}(\nabla _s \xi ^s) \right] \left. \left[ \nabla _{(b} \xi _{n)} {-} \frac{1}{4} g_{bn}(\nabla _l \xi ^l) \right] \right\} .\nonumber \\ \end{aligned}$$
(21)

Using (7) and (8), we can show that this tensor takes zero value, \(T_{pq}^{(\mathrm{K})} = 0\). This means that in the suggested scheme there are no additional contributions to the equations of the gravitational field associated with conformal and/ot standard Killing vector fields.

2.1.5 Short summary

  1. 1.

    If we consider the proposed scheme of using of the Killing vector field as a dynamic field based on the action functional (3) with coupling parameters (14), constraint (7) and its differential consequence (8), we obtain non-violated equations for the gravitational field.

  2. 2.

    If we consider the s-parameter group and deal with the set of Killing vector fields \(\left\{ \xi _{(a)} \right\} \) (\(a=1,2,..s\)), we can extend the model by inserting the term \(\sum _{a} \lambda ^{*}_{(a)}\left[ \nabla _{(m} \xi _{(a)n)} \nabla ^{(m} \xi _{(a)}^{n)} - \frac{1}{4} (\nabla _n \xi ^n_{(a)})^2\right] \) instead of \(\lambda ^{*}\left[ \nabla _{(m} \xi _{n)} \nabla ^{(m} \xi ^{n)}-\frac{1}{4} (\nabla _n \xi ^n)^2 \right] \) into the Lagrangian.

  3. 3.

    For the spatially homogeneous cosmological models of the Friedmann type with the scale factor a(t) there are one conformal time-like Killing vector \(\xi ^i_{(0)} = a(t) \delta ^i_0\), and three space-like divergence-free Killing vectors \(\xi ^i_{(\alpha )} = \delta ^i_{\alpha }\), where \(\alpha = 1,2,3\). Clearly, using the moduli of these Killing vectors we can construct only one additional scalar \(\xi = a(t)\), or can choose more convenient quantity \(x=\frac{a(t)}{a(t_0)}\) (see, e.g., [47,48,49]).

  4. 4.

    When we consider a static model, we can use the time-like Killing vector \(\xi ^i_{(0)} = \delta ^i_0\), so that its modulus \(\xi _{(0)} = \sqrt{g_{00}}\) can play the role of the guiding function (see [28]).

  5. 5.

    For spherically symmetric models we can take the azimuthal Killing vector \(\xi ^i_{(\varphi )} = \delta ^{i}_{\varphi } \) to obtain the additional scalar \(\xi _{(\varphi )} = r \sin {\theta } \), where, as usual, r is the radial variable, \(\theta \) is the meridional angle, \(\varphi \) is the azimuthal angle.

  6. 6.

    When we deal with the spacetimes with the so-called pp-wave symmetry, we obtain one covariant constant null Killing vector \(\xi ^i_{(v)} = \delta _0^i {-} \delta _1^i\) and two space-like Killing vectors \(\xi ^i_{(2)} = \delta _2^i\) and \(\xi ^i_{(3)} = \delta _3^i\). Thus, we can use two moduli \(\xi _{(2)} = \sqrt{|g_{22}|}\) and \(\xi _{(3)} = \sqrt{|g_{33}|}\), as well as, the scalar product \(g_{ik}\xi ^i_{(2)}\xi ^k_{(3)} = g_{23}\), as additional guiding functions.

2.2 Total action functional and extended master equations for interacting fields

2.2.1 Extended action functional

Let us consider now the total action functional

$$\begin{aligned} S_{(\mathrm{total})}= & {} \int d^4 x \sqrt{-g} \left\{ L^{(\mathrm{matter})} + \frac{1}{2\kappa }\left( R+2\Lambda \right) \right. \nonumber \\&\left. + \frac{1}{2\kappa } \sum _{a} \lambda ^{*}_{(a)}\left[ \nabla _{(m} \xi _{(a)n)} \nabla ^{(m} \xi _{(a)}^{n)} - \frac{1}{4} (\nabla _n \xi ^n_{(a)})^2 \right] \right. \nonumber \\&+ \frac{1}{4} \left( F_{mn} {+} \phi F^{*}_{mn}\right) F^{mn} \nonumber \\&\left. - \frac{1}{2} \Psi _{0}^2\left( \nabla _m \phi \nabla ^m \phi {-} V \right) \right\} . \end{aligned}$$
(22)

It describes the electromagnetic field, represented in terms of the Maxwell tensor \(F_{mn}\) and its dual \(F^{*}_{mn}\), which interacts with the pseudoscalar (axion) field \(\phi \) with the potential V; \(L^{(\mathrm{matter})}\) is the Lagrangian of a baryonic matter; the quantity \(\Psi _0\) is reciprocal to the coupling constant of the axion-photon interaction \(\frac{1}{\Psi _0}= g_{A \gamma \gamma }\) (see, e.g., [50] for its observational constraints).

2.2.2 Ansatz about an equilibrium state of the axion system

We assume now that the potential of the pseudoscalar field has the periodic structure

$$\begin{aligned} V(\phi , \xi _{(1)}, \ldots \xi _{(s)}) = \frac{m^2_A \Phi ^2_{*}}{2\pi ^2} \left[ 1- \cos {\left( \frac{2 \pi \phi }{\Phi _*}\right) } \right] , \end{aligned}$$
(23)

where \(\Phi _*=\Phi _*(\xi _{(1)}, \ldots \xi _{(s)})\) is a function of the moduli of the Killing vectors \(\xi _{(a)}\). When \(\phi = n \Phi _*\) with arbitrary integer n, the potential and its derivative take zero values

$$\begin{aligned} V_{|\phi =n\Phi _*} = 0 , \quad \left( \frac{\partial V}{\partial \phi }\right) _{|\phi =n\Phi _*} =0 , \end{aligned}$$
(24)

the values \(\phi =n\Phi _*\) correspond to minima of the potential. The coefficient in front of the periodic function in (23) is chosen so that, when \(\phi \) has a small deviation from the minimum value \(\phi =n\Phi _*\) (i.e., when \(\phi =n\Phi _* {+} \psi \) and \(|\psi |<<1\)), the potential converts into \(V \rightarrow m^2_A \psi ^2\). Keeping in mind the mechanical analogy that equilibrium states of dynamic systems can be realized just in the minimum of the corresponding potential, we use below the special term equilibrium state of the axion field, \(\phi _{(\mathrm{eq})}\), if the potential of the pseudoscalar field \(\phi \) and its derivative with respect to \(\phi \) takes zero value at \(\phi = \phi _{(\mathrm{eq})}\).

2.2.3 Extended master equations

The variation procedure gives the system of master equations of the model, which contains four sub-sets.

(i):

The first sub-set appears as the result of variation of the total functional (22) with respect to the electromagnetic potential \(A_i\); it includes the extended Maxwell equations

$$\begin{aligned}&\nabla _k \left[ F^{ik} {+} \phi F^{*ik} \right] = 0 , \end{aligned}$$
(25)
$$\begin{aligned}&F_{ik} = \nabla _i A_k - \nabla _k A_i \ \Rightarrow \nabla _k F^{ik*} =0 . \end{aligned}$$
(26)

When the axion field is in the equilibrium state, we have to replace \(\phi \) with \(n \Phi _{*}\) in (25).

(ii):

The second sub-set consists of one equation for the pseudoscalar field:

$$\begin{aligned} \nabla _k \nabla ^k \phi + \frac{1}{2} \frac{\partial V}{\partial \phi } = - \frac{1}{4\Psi _{0}^2} F^{*}_{mn}F^{mn} , \end{aligned}$$
(27)

it is the result of variation of the total action functional with respect to \(\phi \). When the axion field is in the equilibrium state this equation converts into

$$\begin{aligned} \nabla _k \nabla ^k \Phi _{*} = - \frac{1}{4 n \Psi _{0}^2} F^{*}_{mn}F^{mn} . \end{aligned}$$
(28)
(iii):

The third subset appears as the result of variation with respect to the vector fields \(\xi ^i_{(a)}\); it has the form

$$\begin{aligned}&\nabla _k \left\{ \lambda ^{*}_{(a)} \left[ \nabla ^{(k} \xi ^{m)}_{(a)} - \frac{1}{4} g^{km} (\nabla _n \xi ^n_{(a)})^2 \right] \right\} \nonumber \\&\quad = - \kappa \Psi ^2_0 \left( \frac{\partial V}{\partial \Phi _{*}}\right) \left( \frac{\partial \Phi _{*}}{\partial \xi _{(a)}}\right) \left( \frac{\xi ^m_{(a)}}{\xi _{(a)}}\right) . \end{aligned}$$
(29)

When we deal with the standard and/or conformal Killing vectors \(\xi _{(a)}\), the left-hand side of this equation vanishes, as it was shown in Section IIA. When the axion field is in the equilibrium state, the quantity \(\frac{\partial V}{\partial \Phi _{*}}\) [see (23)] also is equal to zero. In this sense, our scheme of extension of the potential of the pseudoscalar (axion) field using the moduli of the Killing vectors is self-consistent, when we assume that the axion field takes one of the equilibrium values \(\pm \Phi _*, \pm 2\Phi _*, \ldots , \pm n \Phi _*, \ldots \) In other words, the equations for the fields \(\xi _{(a)}\) coincide with the Killing equations (7).

(iv):

Variation with respect to metric gives the equations for the gravitational field in the following form:

$$\begin{aligned}&\frac{1}{\kappa }\left[ R_{ik} {-} \frac{1}{2} R g_{ik} {-} \Lambda g_{ik} \right] \nonumber \\&\quad = T^{(\mathrm{m})}_{ik} {+} T^{(\mathrm{axion})}_{ik}{+} T^{(\mathrm{em})}_{ik} {+} T^{(\mathrm{V})}_{ik}. \end{aligned}$$
(30)

Here the stress-energy tensors of the matter, of the pseudoscalar (axion) field and of the electromagnetic field are given, respectively, by

$$\begin{aligned}&T^{(\mathrm{m})}_{ik} = -\frac{2}{\sqrt{-g}} \frac{\delta }{\delta g^{ik}}\left[ \sqrt{-g} L^{{(\mathrm{matter})}} \right] , \end{aligned}$$
(31)
$$\begin{aligned}&T^{(\mathrm{axion})}_{ik}= \Psi ^2_0 \left[ \nabla _i\phi \nabla _k\phi -\frac{1}{2} g_{ik} \left( \nabla _n \phi \nabla ^n \phi -V \right) \right] ,\nonumber \\ \end{aligned}$$
(32)
$$\begin{aligned}&T^{(\mathrm{em})}_{ik}=\frac{1}{4}g_{ik}F_{mn}F^{mn}-F_{in}{F_{k}}^n . \end{aligned}$$
(33)

The last term relates to the contribution of the interaction between the axion field and Killing vector field; it can be written as follows:

$$\begin{aligned} T^{(\mathrm{V})}_{ik} = \frac{1}{2} \Psi ^2_0 \left( \frac{\partial V}{\partial \Phi _*}\right) \sum _{a}\frac{\xi _{i(a)}\xi _{k(a)}}{\xi _{(a)}} \left( \frac{\partial \Phi _*}{\partial \xi _{(a)}} \right) . \end{aligned}$$
(34)

The term \(\frac{\partial V}{\partial \Phi _*}\) can be presented in the form:

$$\begin{aligned} \frac{\partial V}{\partial \Phi _*} = \frac{2}{\Phi _{*}} V - \frac{m^2_A \phi }{\pi } \ \sin {\left( \frac{2 \pi \phi }{\Phi _*}\right) } . \end{aligned}$$
(35)

Clearly, when the axion field is in the equilibrium state, i.e., \(\phi = n \Phi _*\), the term \(T^{(\mathrm{V})}_{ik}\) vanishes, and the gravitational field equations remain non-violated.

3 Master equations for the model of static spherically symmetric gravitational field

3.1 Geometry of the outer zone and representation of the guiding functions

Let us consider now an outer zone of a static spherically symmetric object, which possesses a magnetic charge \(\mu \) and an electric charge q; the term \(|Q|\equiv \sqrt{\mu ^2{+} q^2}\) describes the hybrid charge. The mentioned object can be presented, for instance, by a charged black hole; in this case the outer zone covers the spacetime domain from the outer black hole horizon to the infinity. When we consider a dyon, we deal with the zone, which covers the domain from the boundary of the solid body of this object to the infinity. We assume that in both cases the metric in the outer zone can be presented as follows:

$$\begin{aligned} ds^2 = N(r)dt^2 -\frac{dr^2}{N(r)} - r^2 \left( d\theta ^2 + \sin ^2{\theta } d\varphi ^2 \right) , \end{aligned}$$
(36)

and the metric coefficient N(r) is effectively described by the Reissner–Nordström function

$$\begin{aligned} N=1-\frac{2M}{r} + \frac{{Q^2}}{r^2} . \end{aligned}$$
(37)

3.1.1 First guiding function, \(\xi \)

The static spacetime with the metric (36) is known to admit the existence of the time-like Killing vector field \(\xi ^i=\delta ^i_0\), and the modulus of this four-vector is \(\xi \equiv \sqrt{N(r)}\). The value of the modulus \(\xi \) on the outer Reissner–Nordström horizon \(r_{+}= M {+} \sqrt{M^2{-}Q^2}\) is equal to zero, \(\xi (r_+) =0\), and its value at infinity is equal to one, \(\xi (\infty ) =1\). In other words, the scalar \(\xi \) takes values in the interval (0, 1) and can be used as the invariant delimiter of the first type. This means that if one needs to distinguish the specific sphere \(r=r_{*}\), one can use this first scalar \(\xi \) and the appropriate delimiter \(\xi =\xi _{*}=\sqrt{1-\frac{2M}{r_*}{+}\frac{Q^2}{r^2_{*}}}\).

3.1.2 Second guiding function, \(\eta \)

The spacetime with the metric (36) admits the so-called azimuthal Killing vector \(\xi ^i_{(\varphi )} = \delta ^i_{\varphi }\), whose modulus is \(\xi _{(\varphi )}= r \sin {\theta }\). Based on this fact we introduce the transformed scalar quantity

$$\begin{aligned} \eta \equiv \arcsin {\frac{\xi _{(\varphi )} \left( 1-\xi ^2 \right) }{\left[ M+\sqrt{M^2-Q^2\left( 1-\xi ^2 \right) }\right] }}, \end{aligned}$$
(38)

which coincides with the meridional angle \(\theta = \eta \); it can be used as a delimiter to distinguish some special value of the angle \( \theta = \theta _{*}\).

3.1.3 Third guiding function, \(\zeta \)

Also, the spacetime with the metric (36) admits two Killing vectors

$$\begin{aligned} \xi ^i_{(1)}= & {} \sin {\varphi } \ \delta ^i_{\theta } {+} \mathrm{ctg}{\theta } \cos {\varphi } \ \delta ^i_{\varphi } ,\nonumber \\ \xi ^i_{(2)}= & {} \cos {\varphi } \ \delta ^i_{\theta } {-} \mathrm{ctg}{\theta } \sin {\varphi } \ \delta ^i_{\varphi } . \end{aligned}$$
(39)

The difference of their squares

$$\begin{aligned} \xi ^2_{(1)}-\xi ^2_{(2)} = r^2 \sin ^2{\theta } \cos {2\varphi } \end{aligned}$$
(40)

gives the hint for reconstruction of the third guiding function \(\zeta \):

$$\begin{aligned} \zeta \equiv \frac{1}{2} \mathrm{arccos}\left[ \frac{\xi ^2_{(1)}-\xi ^2_{(2)}}{\xi ^2_{(\varphi )}} \right] , \end{aligned}$$
(41)

which coincides with the azimuthal angle \(\varphi = \zeta \). This scalar can be used as a delimiter to distinguish some special value of the angle \( \varphi = \varphi _{*}\).

To conclude, one can say the following: first, when we intend to introduce on the invariant level the special sphere (e.g., to mark the dark matter wall) we have to use the first invariant \(\xi \) and the delimiter \(\xi =\xi _{*}\); second, when we intend to describe the straight dark matter filament, we can use two requirements \(\eta = \theta _{*}\) and \(\zeta = \varphi _{*}\) to fix the corresponding line. In this paper we will illustrate the idea to use the scalar \(\xi \) only; we plan to discuss the problem of filaments in the next work.

3.2 Master equations of the axionic magneto-electro-statics

3.2.1 The key equation

The Eq. (25) for the electromagnetic field near the static spherically symmetric dyon are known to be reduced to one equation

$$\begin{aligned} \frac{d}{dr}\left( r^2 \frac{d A_0}{dr} + \mu \phi \right) = 0 , \end{aligned}$$
(42)

where \(A_0(r)\) is the electrostatic potential, and the magnetic charge \(\mu \) is associated with the magnetic potential \(A_{\varphi }= \mu (1{-}\cos {\theta })\) (see, e.g., [28] for references). Integration of (42) yields

$$\begin{aligned} \frac{d A_0}{dr} = \frac{{{\mathcal {Q}}}(r)}{r^2} , \quad {{\mathcal {Q}}}(r) \equiv Q_* - \mu \phi (r) , \end{aligned}$$
(43)

where \({{\mathcal {Q}}}(r)\) is the so-called effective charge, which depends on the axion field, and \(Q_*\) is the integration constant. When we search for solutions continuous on the interval \(r_{+}<r<\infty \), we can link the parameter \({{\mathcal {Q}}}(\infty ) =Q_{*}{-}\mu \phi {(\infty )}\) with a total electric charge q of the object, which could be found by a distant observer. This means that the quantity \(-\mu \phi {(\infty )}\) plays the role of an effective axionically induced electric charge, and it is predetermined by the behavior of the pseudoscalar field at infinity. When we search for solutions continuous on the interval \(r_{+}<r< r_{*}\), we have to link the constant \(Q_{*}\) with the value of the pseudoscalar field \(\phi (r_{*})\) on the delimiting surface \(r=r_{*}\). Clearly, the electric potential can be found in quadratures, when the profile \(\phi (r)\) is known:

$$\begin{aligned} A_0(r) = A_0(r_0) + Q_* \left( \frac{1}{r_0} - \frac{1}{r}\right) - \mu \int ^r_{r_{0}} \frac{d\rho }{\rho ^2} \ \phi (\rho ) . \end{aligned}$$
(44)

The equation for the axion field (27) with the potential (23) can be reduced to the form:

$$\begin{aligned} \frac{d}{dr} \left( r^2 N \frac{d \phi }{dr} \right) - \frac{m^2_{A} r^2 \Phi _*}{2\pi } \sin {\left( \frac{2\pi \phi }{\Phi _*}\right) }= - \frac{\mu }{\Psi ^2_0} \frac{d A_0}{dr} . \end{aligned}$$
(45)

Excluding \(A_0\) from this equation using (43), and applying the ansatz concerning the equilibrium state (i.e., \(\phi = n \Phi _*\), \(n \ne 0\)), we obtain the key equation of the model in terms of the variable r:

$$\begin{aligned} r^2 \frac{d}{dr}\left( r^2 N \frac{d \Phi _{*}}{dr}\right) = \frac{\mu ^2}{\Psi ^2_0}\left( \Phi _{*} - \frac{Q_*}{n \mu }\right) . \end{aligned}$$
(46)

Since \(\xi = \sqrt{N(r)}\), we can rewrite the Eq. (46) in terms of \(\xi \) using the relationship

$$\begin{aligned} r = \frac{M + \sqrt{\xi ^2 Q^2 + (M^2-Q^2)}}{1-\xi ^2} . \end{aligned}$$
(47)

The sign in front of square root is chosen so that the outer horizon \(r=r_{+}= M{+}\sqrt{M^2{-}Q^2}\) corresponds to the value \(\xi =0\). We denote the derivative with respect to \(\xi \) as a prime, and obtain the following differential equation:

$$\begin{aligned} \Phi _{*}^{\prime \prime }(\xi ) {+} \left( \frac{1}{\xi } {+} \frac{\xi }{\xi ^2 {+} \nu } \right) \Phi _{*}^{\prime }(\xi ) - \frac{\mu ^2 \left( \Phi _{*}{-} \frac{Q_*}{n\mu } \right) }{\Psi ^2_0 Q^2 \left( \xi ^2 {+} \nu \right) } =0 , \end{aligned}$$
(48)

with the guiding parameter \(\nu \) given by

$$\begin{aligned} \nu = \frac{M^2{-}Q^2}{Q^2} . \end{aligned}$$
(49)

Below we indicate (48) as the key equation.

3.2.2 Stepwise equilibrium functions

The ansatz about the equilibrium function \(\phi = n \Phi _*(\xi )\) can be extended as follows. We assume now, that there are two domains \(0<\xi <\xi _{*}\) and \(\xi _{*}<\xi <1\), divided by the spherical surface indicated by the value \(\xi _{*}\) of the scalar \(\xi \). We assume that the integer n takes the values \(n_1\) and \(n_2\) in the first and second domains, respectively. The equilibrium function describing the pseudoscalar field can be now presented by the stepwise function

$$\begin{aligned} \phi = \phi _{(\mathrm{eq})}= \Phi _*(\xi ) \left[ n_1 h(\xi _{*} -\xi ) + n_2 h(\xi - \xi _{*}) \right] , \end{aligned}$$
(50)

where h(z) is the Heaviside function, equal to one, when \(z \ge 0\) and to zero, when \(z<0\). Why this extension is considered to be interesting? From the mathematical point of view, this extension keeps the fundamental properties of the potential (23), i.e., the potential and its first derivative take zero values, when the axion field is in the equilibrium state. Indeed, we can easily check the following formulas:

$$\begin{aligned}&\frac{2\pi ^2 V(\phi _{(\mathrm{eq})})}{m^2_A \Phi ^2_{*}} = 1{-} \cos {\left\{ 2\pi \left[ n_1 h(\xi _{*} {-}\xi ) {+} n_2 h(\xi {-} \xi _{*}) \right] \right\} }\nonumber \\&\quad = 1{-} \cos {\left[ 2\pi n_1 h(\xi _{*} -\xi )\right] } \cos {\left[ 2\pi n_2 h(\xi -\xi _{*})\right] }\nonumber \\&\qquad {+} \sin {\left[ 2\pi n_1 h(\xi _{*} {-}\xi ) \right] } \sin {\left[ 2\pi n_2 h(\xi {-}\xi _{*})\right] } =0 , \end{aligned}$$
(51)
$$\begin{aligned}&\frac{\pi }{m^2_A \Phi _{*}}\frac{dV}{d\phi }(\phi _{(\mathrm{eq})}) = \sin {\left\{ 2\pi \left[ n_1 h(\xi _{*} {-}\xi ) {+} n_2 h(\xi {-} \xi _{*}) \right] \right\} } \nonumber \\&\quad =\sin {\left[ 2\pi n_1 h(\xi _{*} {-}\xi )\right] } \cos {\left[ 2\pi n_2 h(\xi -\xi _{*})\right] } \nonumber \\&\qquad + \sin {\left[ 2\pi n_2 h(\xi -\xi _{*})\right] } \cos {\left[ 2\pi n_1 h(\xi _{*} {-}\xi ) \right] } = 0 . \end{aligned}$$
(52)

In principle, the mentioned unique properties of the periodic potential allow us to extend the described procedure from the example of the two-level profile to a multi-level profiles. In other words, one can consider three, four, etc. levels in the axionic dark matter profiles, using the extension of the formula (50) for the set of integers \(n_1\), \(n_2\), \(n_3\), \(n_4\), etc. As for the function \(\Phi _{*}(\xi )\), its profile plays the role of a common envelope.

3.2.3 Conditions on the boundary \(\xi =\xi _{*}\)

For the two-level profiles one has to solve the key equation (48) in two domains \(0<\xi <\xi _{*}\) and \(\xi _{*}<\xi <1\) separately, and to assume, that the parameter \(Q_{*}\) takes different values in these domains \(Q_{*}^{(1)} \ne Q_{*}^{(2)}\). The envelope function \(\Phi _{*}(\xi )\) is considered to be continuous near the surface \(\xi =\xi _{*}\),

$$\begin{aligned}{}[\Phi _{*}] \equiv \lim _{\epsilon \rightarrow 0} \left\{ \Phi _{*}(\xi _{*}+ \epsilon )- \Phi _{*}(\xi _{*}- \epsilon ) \right\} = 0 . \end{aligned}$$
(53)

This is possible, when the parameters \(Q_{*}^{(1)}\) and \(Q_{*}^{(2)}\) are linked as follows:

$$\begin{aligned} \frac{Q_{*}^{(1)}}{n_1} = \frac{Q_{*}^{(2)}}{n_2} \equiv \frac{Q_{*}}{n_*} , \end{aligned}$$
(54)

providing the key equation (48) to be the same in both spatial domains. When \(n_1 \ne n_2\), we deal with the jump of the pseudoscalar field (50) on the boundary \(\xi =\xi _{*}\)

$$\begin{aligned}{}[\phi ] \equiv \lim _{\epsilon \rightarrow 0} \left\{ \phi (\xi _{*}{+} \epsilon ){-} \phi (\xi _{*}{-} \epsilon ) \right\} = (n_2{-}n_1) \Phi _{*}(\xi _{*}) . \end{aligned}$$
(55)

The value of this jump \([\phi ]\) is equal to zero, if the delimiting value \(\xi _{*}\) coincides with the null of the envelope function, i.e., \(\Phi _{*}(\xi _{*})=0\). Similarly, the jump of the derivative \([\phi ^{\prime }]\) vanishes, when two conditions are satisfied: \(\Phi _{*}(\xi _{*})=0\) and \(\Phi ^{\prime }_{*}(\xi _{*})=0\). Since the integer \(n_1\) appears in the denominator [see (54)], the case \(n_1=0\) should be considered separately; now the scheme of analysis is consistent, if we put \(Q^{(1)}_{*}=0\).

From the physical point of view, the condition \([\phi ] \ne 0\) means that we deal with a wall with a non-zero surface density of axions. Modeling of the two-level distributions of the axionic dark matter is also interesting for description of profiles near black holes. We mean that a radius can exist, say \(r_{*}\), which marks the specific boundary: \(\xi =\xi _{*}\). When \(\xi <\xi _{*}\) we find that \(n_1=0\) and thus \(\phi _{(\mathrm{eq})} = 0\), i.e., all the dark matter particles are absorbed by the black hole. When \(\xi >\xi _{*}\) the axionic dark matter profile is not empty, i.e. \(n_2 \ne 0\), and particles rotating around the black hole can resist to the gravitational attraction.

4 Analysis of the key equation

4.1 The Heun equation

The key equation (48) is the particular case of the known Heun equation

$$\begin{aligned} Y^{\prime \prime } + Y^{\prime } \left[ \frac{\epsilon }{x{-}a} {+} \frac{\delta }{x{-}1} {+} \frac{\gamma }{x} \right] + Y \frac{\alpha \beta x {-} \rho }{(x{-}a)(x{-}1)x} = 0 , \end{aligned}$$
(56)

(see, e.g., [51, 52]), which is in its turn the particular case of the Fuchs equation [53, 54]. The solution of this equation is regular at infinity, when \(\epsilon {+}\gamma {+}\delta =\alpha {+}\beta {+}1\). The equation (48) can be transformed to (56) using the relationship \(x = \frac{i \xi }{\sqrt{\nu }}\), if we put

$$\begin{aligned} \epsilon= & {} \delta = \frac{1}{2} , \quad \gamma =1 , \quad a=-1 , \quad \rho =0 ,\nonumber \\ \alpha= & {} \frac{1}{2} + \sqrt{\frac{1}{4} + \frac{\mu ^2}{\Psi ^2_0 Q^2}} , \quad \beta = \frac{1}{2} - \sqrt{\frac{1}{4} + \frac{\mu ^2}{\Psi ^2_0 Q^2}} . \end{aligned}$$
(57)

The dimensionless guiding parameter \(\nu = \frac{M^2{-}Q^2}{Q^2}\) plays essential role in the analysis of the key equation.

  1. 1.

    The standard model relates to the positive value \(\nu >0\) (or \(M^2>Q^2\)); in this case we deal with the object (e.g., the charged black hole), which has the inner and outer horizons. The key equation (48) is characterized by one real singular point only, \(\xi =0\), which is situated on the left boundary of the admissible interval for \(\xi \).

  2. 2.

    When \(\nu \) is negative, we deal with the so-called naked central singularity in the spacetime (singularity without horizons). In this case in the key equation there are two real singular points: \(\xi =0\) and \(\xi = \sqrt{|\nu |}\) (we consider \(\xi \) to be positive).

  3. 3.

    The intermediate case \(\nu =0\) (or \(M^2=Q^2\)) describes the so-called extremal black hole, in which outer and inner horizons coincide. Again, the key equation (48) is characterized by one real singular point only, \(\xi =0\). First of all, we consider in more detail this intermediate case.

4.2 \(M^2=Q^2\): extremal black hole

4.2.1 General solution to the key equation

When \(M^2=Q^2\), the inner and outer horizons of the object coincide, and the equation (48) can be reduced to the Euler equation

$$\begin{aligned} \xi ^2 \Phi _{*}^{\prime \prime }(\xi ) + 2\xi \Phi _{*}^{\prime }(\xi ) - \frac{\mu ^2}{\Psi ^2_0 M^2} \left( \Phi _{*}- \frac{Q_{*}}{n_{*}\mu } \right) =0 . \end{aligned}$$
(58)

The general solution to (58) has the following form:

$$\begin{aligned} \Phi _{*}(\xi ) = \frac{Q_{*}}{n_{*}\mu } + C_1 \xi ^{\sigma _1} + C_2 \xi ^{\sigma _2} , \end{aligned}$$
(59)

where \(C_1\) and \(C_2\) are the constants of integration, and the power indices \(\sigma _1\) and \(\sigma _2\)

$$\begin{aligned} \sigma _1= & {} \frac{1}{2} \left( \sqrt{1+ \frac{4\mu ^2}{\Psi ^2_0 M^2}} -1 \right) \ge 0 ,\nonumber \\ \sigma _2= & {} -\frac{1}{2} \left( \sqrt{1+ \frac{4\mu ^2}{\Psi ^2_0 M^2}} +1 \right) <0 , \end{aligned}$$
(60)

are of the opposite sign. Since \(\sigma _1{+}\sigma _2={-}1\), we have only one parameter of modeling, which is associated, in fact, with the value of the ratio \(\frac{\mu ^2}{\Psi ^2_0 M^2}\).

4.2.2 The solution regular in the interval \(\xi _0<\xi <1\)

Here we assume that at \(\xi =\xi _{0}>0\) the function \(\Phi _{*}(\xi )\) takes the value \(\Phi _{*}(\xi _0)\). Such a boundary condition is typical for the case, when we deal with a magnetic star, and the parameter \(\xi _{0}\) relates to its radius. The spatial infinity \(r \rightarrow \infty \) relates to \(\xi \rightarrow 1\), and the corresponding value of the function \(\Phi _{*}\) is indicated as \(\Phi _{*}(1) \equiv \Phi _{\infty }\). Then the solution for the function \(\Phi _{*}(\xi )\) takes the analytical form

$$\begin{aligned} \Phi _{*}(\xi )= & {} \frac{Q_*}{n_{*}\mu }\left( 1-\xi ^{\sigma _2} \right) + \Phi _{\infty }\xi ^{\sigma _2} \nonumber \\&+ \left( \frac{\xi ^{\sigma _1} {-} \xi ^{\sigma _2}}{\xi ^{\sigma _1}_0 {-} \xi ^{\sigma _2}_0} \right) \left[ \Phi _*(\xi _0) {-} \frac{Q_*}{n_{*}\mu }\left( 1{-}\xi ^{\sigma _2}_0 \right) {-} \Phi _{\infty } \xi _0^{\sigma _2}\right] . \nonumber \\ \end{aligned}$$
(61)

There are two guiding parameters in this formula: the amplitude factor \(\frac{Q_*}{n_{*}\mu }\), which can be positive or negative with respect to the signs of the constants \(Q_*\) and \(\mu \), as well as, the positive factor \(\frac{\mu ^2}{\Psi ^2_0 M^2}\). Typical behavior of these profiles is illustrated in Fig. 1.

Fig. 1
figure 1

Typical profiles of the regular envelope function \(\Phi _*(\xi )\) for the case of extremal charged black hole (\(Q^2=M^2\)) given by the formula (61). The first number in the parentheses relates to the value of the parameter \(\frac{\mu ^{2}}{\Psi _{0}^{2}M^{2}}\), the second one corresponds to the parameter \(\frac{Q_{*}}{n\mu }\)

4.2.3 An example of the stepwise solution

Let us consider the illustration, which corresponds to the equilibrium function (50) with the delimiter \(\xi _{*}=\frac{1}{4}\), and the model parameters \(n_1=0\), \(n_2=1\), \(Q^{(1)}_{*}=0\)

$$\begin{aligned} \phi _{(\mathrm{eq})}= \Phi _{*}(\xi ) \times h \left( \xi - \frac{1}{4}\right) , \end{aligned}$$
(62)

where the envelope function in the interval \(\frac{1}{4}<\xi <1\) is of the form

$$\begin{aligned} \Phi _{*}(\xi ) = \frac{1}{81}\left( \frac{1}{\xi ^2} {+} 128 \xi {-} 48 \right) . \end{aligned}$$
(63)

This envelope function takes zero value at \(\xi =\frac{1}{4}\); its derivative

$$\begin{aligned} \Phi ^{\prime }_{*}(\xi ) = \frac{2}{81}\left( 64{-}\frac{1}{\xi ^3} \right) \end{aligned}$$
(64)

also vanishes at \(\xi =\frac{1}{4}\). This function is the exact solution to the key equation with the following values of the model parameters:

$$\begin{aligned} \frac{Q^{(2)}_*}{n_2 \mu }= {-} \frac{16}{27} , \quad \frac{\mu }{\Psi _0 M} =\sqrt{2} , \quad \sigma _1=1 , \sigma _2=-2 . \end{aligned}$$
(65)

and takes the value \(\Phi _{\infty }=1\) at the spatial infinity. Mention should be made that \(\Phi (\xi )=0\) is the exact solution to the key equation in the interval \(0<\xi <\frac{1}{4}\). On the delimiting sphere \(\xi =\frac{1}{4}\) the function (63) and the constant function \(\Phi (\xi )=0\) happen to be sewed. The profile of the function (62) with (63), as well as, two additional ones are depicted in Fig. 2.

Fig. 2
figure 2

Three examples of the stepwise envelope function (62); the profiles are distinguished by the parameters \(\sigma _{1}\) and \(\sigma _{2}\), indicated in the parentheses. Graphs of all the functions have the delimiter \(\xi _0 = \xi _* = \frac{1}{4}\), and tend to one at infinity

4.3 \(Q=0\): the Schwarzshild type black hole

4.3.1 Key equation and its general solution

When \(Q^2<<M^2\), the geometry of the background spacetime is close to the one of the Schwarzshild type; we mean that the magnetic charge \(\mu \) is non-vanishing and it plays an important role in the axionic halo formation, however, in the formation of the background gravitational field its contribution is negligible. Now we can obtain the reduced key equation from (48) as the limiting case \(Q \rightarrow 0\); we deal now with the Bessel equation

$$\begin{aligned} \xi ^2 \Phi _{*}^{\prime \prime }(\xi ) + \xi \Phi _{*}^{\prime }(\xi ) - \frac{\mu ^2}{\Psi ^2_0 M^2} \ \xi ^2 \left( \Phi _{*} - \frac{Q_{*}}{n_{*}\mu } \right) =0 . \end{aligned}$$
(66)

The general solution to (66) is

$$\begin{aligned} \Phi _{*}(\xi ) = \frac{Q_{*}}{n_{*}\mu } + C_1 I_0\left( \frac{\mu \xi }{\Psi _0 M}\right) + C_2 K_0\left( \frac{\mu \xi }{\Psi _0 M} \right) . \end{aligned}$$
(67)

Here \(I_0\) and \(K_0\) are the Bessel functions of the third and fourth type with indices equal to zero; they can be represented standardly as follows:

$$\begin{aligned} I_0(z)= & {} \sum _{m=0}^{\infty }\frac{\left( \frac{z}{2}\right) ^{2m}}{(m!)^2} , \quad I_0(z \rightarrow 0) \approx 1 + \frac{z^2}{4} , \end{aligned}$$
(68)
$$\begin{aligned} K_0(z)= & {} \int _0^{\infty } d \zeta e^{-z \cosh {\zeta }} , \quad K_0(z \rightarrow 0) \approx \log {\frac{z}{2}} . \end{aligned}$$
(69)

4.3.2 Regular solution in the interval \(0<\xi <1\)

If we consider the requirement \(\Phi _{*}(0) = 0\) we find immediately the regular profile

$$\begin{aligned} \Phi _{*}(\xi ) = \frac{Q_{*}}{n_{*} \mu } \left[ 1 - I_0\left( \frac{\mu \xi }{\Psi _0 M} \right) \right] . \end{aligned}$$
(70)

Clearly, \(\Phi ^{\prime }_{*}(0)=0\) and \(\Phi _{*}(1) = \frac{Q_{*}}{n_{*} \mu } \left[ 1 - I_0\left( \frac{\mu }{\Psi _0 M} \right) \right] \). In other words, the values of the function \(\Phi _{*}(\xi )\) and of its derivative are equal to zero at the Schwarzschild horizon, and the maximal value of the modulus of this function depends on the following two ratios: \(\left| \frac{Q_{*}}{n_{*} \mu }\right| \) and \(\Gamma \equiv \frac{\mu }{\Psi _0 M}\).

4.3.3 Stepwise equilibrium function

Again we assume that \(n_1=0\), \(n_2=1\) and \(Q^{(1)}_{*}=0\), so that the stepwise equilibrium function is of the form

$$\begin{aligned} \phi _{(\mathrm{eq})}= \Phi _{*}(\xi ) \times h \left( \xi - \xi _{*}\right) , \end{aligned}$$
(71)

with the envelope function

$$\begin{aligned} \Phi _*(\xi ) = \frac{Q_*}{n_* \mu } \left\{ 1 {-} \frac{K^{\prime }_0(\Gamma \xi _*) I_{0}(\Gamma \xi ) {-} K_0(\Gamma \xi ) I^{\prime }_{0}(\Gamma \xi _*)}{K^{\prime }_0(\Gamma \xi _*) I_{0}(\Gamma \xi _*) {-} K_0 (\Gamma \xi _*) I^{\prime }_0(\Gamma \xi )} \right\} ,\nonumber \\ \end{aligned}$$
(72)

which is the exact solution to the key equation. The prime denotes the derivative with respect to the argument of the function. We introduced here the auxiliary parameter \(\Gamma = \frac{\mu }{\Psi _0 M}\). Clearly, \(\Phi _*(\xi _*) =0\) and \(\Phi ^{\prime }_*(\xi _*) =0\). The delimiter value \(\xi _*\) satisfies the transcendent equation \(\Phi _*(1)= \Phi _{\infty }\). The typical behavior of the envelope function is presented in Fig. 3.

Fig. 3
figure 3

Illustrations of the stepwise equilibrium function (71) with the envelope function (72). For all profiles \(\frac{Q_{*}}{n_{*}\mu }=1\) and the delimiter is \(\xi _*=0.25\). The values of the parameter \(\Gamma = \frac{\mu }{\Psi _0 M}\) are written in parentheses

4.4 \(M=0\): the spacetime with naked singularity

When \(M=0\) the spacetime has no horizons, and the variable \(\xi \) takes values in the interval \((1,\infty )\). Now \(\nu ={-}1\), thus the Eq. (48) can be reduced to the Heun equation (56) with the parameters (57) by the replacement \(\xi \rightarrow x\). Now the key equation has no singular points in the interval \((1,\infty )\). Illustrations are presented in Fig. 4.

Fig. 4
figure 4

Typical profiles of the envelope function \(\Phi _*(\xi )\) for the magnetic naked singularity as the solutions to the Heun equation (48). The scalar \(\xi \) takes values in the interval \((1,\infty )\); the point \(\xi =1\) corresponds to the spatial infinity, and \(\xi \rightarrow \infty \), when \(r \rightarrow 0\). In the parentheses two parameters are presented: the first is \(\frac{\mu ^{2}}{\Psi _{0}^{2}Q^{2}}\), the second is \(\frac{Q_{*}}{n\mu }\). Asymptotically, at \(\xi \rightarrow \infty \), the solutions to the Heun equation tend to the constants \(\Phi _{0}\), which depend on the values of the chosen model parameters

4.5 Numerical analysis of the general case \(M^2>Q^2\)

For illustration of the results of the numerical analysis of the case \(\nu >0\), the most physically motivated, we studied systematically the Heun equation (48) for various values of the guiding parameters \(\frac{M^2{-}Q^2}{Q^2}>0\) and \(\frac{\mu ^2}{\Psi ^2_0 Q^2}\). Now, we deal with the interval \(0<\xi <1\), and for the outer zone of the object there are no real singular points in the key equation (48). Typical profiles of the envelope function \(\Phi _*(\xi )\) are presented in Fig. 5.

Fig. 5
figure 5

Typical profiles of the envelope function \(\Phi _*(\xi )\), the solutions to the key equation (48), for the case \(M^2>Q^2\) and \(\Phi (1)=0.5\); the numbers in parentheses correspond to the values of the parameters \(\frac{M^2}{Q^2}{-}1\) and \(\frac{\mu ^{2}}{\Psi _{0}^{2}M^{2}}\), respectively

5 Discussion and conclusions

The authors of the work [30] have formulated the idea that the dark matter axions form a Bose–Einstein condensate, and thus the behavior of the axion systems differs from the one of an ordinary cold dark matter, especially if the regime of interaction is non-linear and the external fields are strong. This idea emphasizes that the axion system is not a simple collisionless, pressureless cold gas; the axion system has to be characterized by some internal self-interaction. We also follow the idea that axionic systems are self-interacting, and we think that their internal structures are predetermined by the modified periodic potential (23). Minima of this potential predetermine the equilibrium states of the axions. Our modified periodic potential can be obtained from the standard one by introduction of the so-called envelope function \(\Phi _{*}\), which stands to emphasizes the fact that the equilibrium value of the axion field depends on the position in the spacetime.

Since the strict variation formalism requires the envelope function to depend on some scalar invariants, \(\Phi = \Phi _{*}(\xi ,\eta ,\zeta )\), we have linked these scalars with moduli of three Killing vectors, which are associated with static spherically symmetric spacetime under consideration.

This extension of the axionic potential, in its turn, led us to the necessity of the Lagrangian modification, in which we considered new terms associated with the Killing vector field. Based on analogy with the Einstein–Aether theory, we obtained the correspondingly extended master equations, thus providing the whole model to be self-consistent. The main conclusion of this first part of the work is the following: when the axion field is in the equilibrium state, for which the modified potential and its first derivative vanish, the presence of dynamically defined Killing vector fields does not violate the master equations for gravitational, electromagnetic and axion fields.

In Sect. 3 we considered the application of the developed formalism for the spacetime of the Reissner–Nordström type and have found exact solutions, which describe axionic halo profiles near the dyons and magnetic black holes. We have found the envelope functions for several values of guiding model parameters; the most interesting findings, from our point of view, are the solutions of the stepwise type, which describe two-level distributions of the axionic dark matter. In principle, such model distributions can be considered in the context of description of the magnetic black holes. Indeed, because of the gravitational attraction the axionic dark matter halo surrounding such a object should have empty zone from the outer horizon till to the first stable orbit of the rotating massive particle (see Figs. 2, 3).

The developed formalism also can be applied to the description of the dark matter filaments; for this purpose we can use the scalars \(\eta \) (38) and \(\zeta \) (41). We hope to consider these cosmological units in the nearest future.