1 Introduction

The Standard Model of accretion in astrophysics [1, 2] is based on the idea that the angular momentum transport across the plasma disk is ensured by effective viscosity, due to the turbulence phenomena triggered by the magneto-rotational instability (MRI) [3], see also Refs. [4,5,6]. However, the necessity to account for a magnetic field, due to the central object, together with the long-wavelength response of the plasma disk, implies that the azimuthal component of the generalized Ohm law is consistent only postulating an anomalous plasma resistivity, see Refs. [7, 8]. This effective resistivity makes the disk magnetic field limited in amplitude, because of its diffusive nature [7]. As a result, the possibility to deal with the high magnetic field strength required for the jet generation leads to investigate specific scenarios [9].

In Refs. [10, 11] (see also Ref. [12]), it has been investigated the possibility for a small scale backreaction of the plasma disk to the magnetic field of the central object, outlining the emergence of a crystalline morphology of the magnetic field microscales, i.e., a local radial oscillation of the magnetic flux function. In particular, in the limit of a linear backreaction, in Ref. [10] it has been written down a Master Equation for the magnetic flux function, which is the starting point of the present investigation. The relevance of this new paradigm relies on the possibility to use the O-points of the resulting configuration as the sites where the jet formation could start, as discussed in Refs. [13, 14].

Here, we develop a detailed general analysis of the linear Master Equation for the crystalline profile of the plasma disk, outlining two main significant features: i) performing a general Fourier analysis, we clarify how the decreasing behavior of the plasma backreaction with the vertical quote (discussed in Refs. [10, 11], see also Ref. [15]) is actually a peculiar property of the fundamental harmonics in the radial oscillation; ii) by analyzing the Master Equation in the framework of a separable solution, we clarify the existence of regions of the parameter space, where a relevant sequence of O-points of the magnetic configuration can take place.

Both the results mentioned above have a precise physical implication. The point i) shows how the higher order harmonics of the radial oscillation of the flux function are actually associated to an inverse behavior of the vertical gradient between the background and the linear perturbation. This issue calls attention for its possible implications on the morphology of the nonlinear backreaction which is associated with the fragmentation of the disk in a ring-like morphology, see Ref. [11] and also Refs. [16,17,18]. The point ii) demonstrates that the topology of the induced radial component of the magnetic field can be characterized by a sequence of O-points, where such a radial component vanishes (noting that the dipole-like magnetic field of the central object is essentially vertical on the equatorial plane). This feature, absent in the analysis in Refs. [10, 11], suggests that the magnetic micro-structure, emerging when the perturbation scale is sufficiently small (i.e., for sufficiently high values of the plasma \(\beta \) parameter), could be an interesting mechanism in triggering sites for the jet formation. It is worth noting that we also outline the existence of a parameter region for the model in which both the radial and the vertical dependence of the backreaction magnetic surfaces oscillate.

2 Standard model for accretion

The basic idea for accretion onto a compact object has been formulated in Ref. [1] and Ref. [2], where essentially thin disk configurations were considered. Actually, thin accretion disks are commonly present around stellar and black hole systems and they have the advantage to elucidate the accretion features in a one-dimensional model, i.e., essentially the radial dependence on the r coordinate is described, being the vertical one (z) averaged out of the problem and the azimuthal dependence on \(\phi \) is forbidden by the axial symmetry of the steady disk configuration [7].

For such a steady axially symmetric one-dimensional representation of a thin disk configuration, accreting material onto a compact object, having a mass M and an almost vertical magnetic field \(\mathbf{B} = B\hat{\mathbf{e }}_z\), is determined via the implementation of the mass conservation equation and of the momentum conservation ones.

If we define the disk accretion rate as

$$\begin{aligned} {\dot{M}} \equiv - 2\pi r\varSigma u_r \; , \end{aligned}$$
(1)

where \(\varSigma \) and \(u_r\) denote the superficial density and the radial infalling velocity of the fluid, respectively, then the mass conservation equation reduces to the following condition

$$\begin{aligned} d{\dot{M}}/dr = 0\; . \end{aligned}$$
(2)

The radial component of the momentum conservation equation splits into two components, one stating that the disk angular velocity \(\omega \) is essentially equal to the Keplerian one \(\omega _K\) and the advection and the pressure gradient balance each other separately, i.e.,

$$\begin{aligned} \omega (r) \simeq \omega _K= & {} \sqrt{GM/r^3}, \end{aligned}$$
(3)
$$\begin{aligned} \rho _0 u_r\frac{\hbox {d}u_r}{\hbox {d}r} + \frac{\hbox {d}p_0}{\hbox {d}r}= & {} 0 \; , \end{aligned}$$
(4)

\(\rho _0(r)\) and \(p_0(r)\) denoting the values of the mass density and of the pressure, taken on the equatorial plane of the disk.

The separation of the gravo-static equilibrium into Eqs. (3) and (4) is a natural assumption when the thinness of the disk is taken into account, see Ref. [7], and it relies on the idea that the advection and pressure gradient terms balance each other, on a different scale with respect to the Keplerian motion of the disk. In Ref. [19], it has been systematically investigated the equilibrium configuration of a disk, demonstrating that, only when the disk is thick enough, the differential angular velocity significantly deviates from the Keplerian value (a configuration very far from the applicability presented in Eqs. (3) and (4) is recovered in the so-called ADAF [20] and see also Refs. [21, 22] for the hydrostatic equilibrium in strong gravitational field).

When, in the next section, we will analyze the magnetic micro-structures generated at small spatial scales, we will see that a deviation from the Keplerian rotation of the disk can be also induced by the plasma backreaction and in that case, it is compensated by the Lorentz force. However, there, the advection term is not present, but, as it has been discussed in Ref. [8], such contributions, at high \(\beta \) values of the plasma, live on a macroscopical spatial scale, decoupling from the microscale balance of the Lorentz force with the correction to the centripetal Keplerian force (see also the discussion in the following Sect. 3.1).

The vertical component of the momentum equation determines the decay law of the mass density far away the equatorial plane. The details of such decay depend on the particular equation of state we adopt for the fluid, but a reliable expression, valid both for the general polytropic and isothermal case, reads as

$$\begin{aligned} \rho (r,z) \simeq \rho _0 (r) \Big ( 1 - \frac{z^2}{H^2}\Big ) \; , \end{aligned}$$
(5)

where H(r) denotes the half-width of the disk and if it is thin we must have \(H/r \ll 1\). The spatial scale of the disk vertical extension is fixed in the vertical gravostatic equilibrium via the ratio between the sound velocity to the Keplerian angular velocity, i.e., \(H\sim v_s/\omega _K\). The sound velocity can be reliably estimated via the relation \(v_s\sim \sqrt{K_\mathrm{B}T/m}\), being T the disk temperature, \(K_\mathrm{B}\) the Boltzmann constant and m the proton mass.

Finally, the azimuthal component of the momentum conservation equation, which regulates the angular momentum transport across the disk, provides the following expression of the accretion rate:

$$\begin{aligned} {\dot{M}} = 6\pi \eta _v H\; , \end{aligned}$$
(6)

here \(\eta _v\) denotes the shear viscosity coefficient, associated to the friction of the disk layers in differential rotation with \(\omega _K\). This viscosity effect can not be due to the kinetic conditions of the plasma, which is actually quasi-ideal (see Ref. [23]), and in Ref. [1] \(\eta _v\) has been justified via a turbulent behavior emerging in the disk spatial microscales. Such a coefficient can be estimated by comparing the exact expression of the \((r,\,\phi )\) component of the viscous stress tensor \(\tau _{ij}\) with its interpretation in terms of the correlation function of the turbulent radial and azimuthal velocity component \(v_r\) and \(v_{\phi }\), respectively, i.e., we have

$$\begin{aligned} \tau _{r\phi } = \eta _vr\frac{\hbox {d}\omega _K}{\hbox {d}r} = - \frac{3}{2}\eta _v\omega _K= - \langle \rho _0 v_rv_\phi \rangle = - \alpha \rho _0 v_s^2,\nonumber \\ \end{aligned}$$
(7)

from which we get the basic Shakura expression

$$\begin{aligned} \eta _v = \frac{2}{3}\alpha \rho _0 v_sH\; . \end{aligned}$$
(8)

The parameter \(\alpha \le 1\) has been introduced to phenomenologically account that all the supersonic velocity fluctuations are unavoidably damped in the turbulent regime.

By combining this expression for \(\eta _v\) with the relation (6), we arrive to the final expression for the accretion rate of the disk in terms of basic quantities (particle density, temperature, central body mass), namely

$$\begin{aligned} {\dot{M}} = 4\pi \alpha \rho _0 v_sH^2 = 4\pi \alpha \rho _0 \omega _K H^3. \end{aligned}$$
(9)

The possibility to deal with a turbulent behavior of the plasma requires that a sufficiently significant instability can affect the equilibrium configuration. Since in the case of a Keplerian disk the steady profile turns out to be stable [3], after some years of investigation, the source of turbulence has been identified in the MRI [24, 25]. Firstly derived in Ref. [26] and Ref. [27], MRI is an Alfvénic instability taking place when the magnetic field intensity within the plasma is below a given threshold. In particular, for a Keplerian disk, the simplest condition to deal with MRI reads as follows

$$\begin{aligned} \omega _A \equiv kv_A < \sqrt{3}\;\omega _K\; , \end{aligned}$$
(10)

\(v_A\) being the Alfvén speed and k the perturbation wave-number. In correspondence to an assigned value of the magnetic field, a minimal spatial scale exists \(\lambda = 2\pi /k\) for which MRI holds. Such a scale, say \(\lambda _{\min }\) is easily identified as

$$\begin{aligned} \lambda _{\min } \simeq \frac{2\pi }{\sqrt{3}} \frac{v_A}{v_s} H \sim \frac{H}{\sqrt{\beta }}\; , \end{aligned}$$
(11)

\(\beta \) denoting the standard plasma parameter, which in astrophysics can take also rather large values.

In Ref. [10], it has been shown that just on the small scale \(\lambda _{\min }\), there stated as \(\sim v_A/\omega _K\), the plasma backreaction to the central body can also take a very different nature with respect turbulence. In fact, a steady perturbation to the equilibrium is allowed for which the correction of the centripetal force associated to the differential rotation, is directly linked to the perturbed Lorentz force raised in the plasma by the emergence of current densities on very small scales. We observe that the background magnetic field is current-free since it is due to the central body magnetosphere and, within the disk, it is a vacuum field. In Ref. [11], this idea of a microscopic crystalline morphology of the perturbed magnetic field, i.e., its radial oscillating behavior on the scale \(\lambda _{\min }\), is implemented also in the limit of a nonlinear beackreaction. As a result, it is possible to show that the disk morphology is decomposed in a series of microscopic ring-like structures, i.e., the mass density of the disk acquires periodic nodes, see also the global model developed in Ref. [12].

In the simplest case of a linear small scale backreaction within the disk, the magnetic flux function is currugated and this perturbed profile is associated to a linear two-dimensional partial differential equation, dubbed in what follows as Master Equation for the disk crystalline structure. We will discuss in detail the morphology of the solutions associated to this equation in order to extract information about the small scale properties of the disk steady profile on which MRI can develop. We conclude by observing that in Ref. [8], it was argued how the presence of a small scale crystalline disk backreaction and the associated relevant micro-current density can have significant implications on the problem of the “anomalous” disk resistivity, required to account for the observed accretion rates in systems like X-ray binary stars.

3 Physical model

We consider a steady and axialsymmetric thin disk configuration, characterized by a magnetic field \(\mathbf{B} \) having the following poloidal form

$$\begin{aligned} \mathbf{B} = -r^{-1}\partial _z \psi \hat{\mathbf{e }}_r + r^{-1}\partial _r\psi \hat{\mathbf{e }}_z \, , \end{aligned}$$
(12)

where we use the aforementioned standard cylindrical coordinates \((r\, ,\phi \, ,z)\) (\(\hat{\mathbf{e }}_r\), \(\hat{\mathbf{e }}_{\phi }\) and \(\hat{\mathbf{e }}_z\) being their versors) and \(\psi (r,z)\) is the magnetic flux function. The disk also possesses a purely azimuthal velocity field \(\mathbf{v} \), i.e.,

$$\begin{aligned} \mathbf{v} = \omega (\psi ) r \hat{\mathbf{e }}_{\phi } \, , \end{aligned}$$
(13)

where the angular velocity \(\omega \) is a function of \(\psi \) because we are in the range of validity of the so-called co-rotation theorem [28], i.e., Eqs. (12) and (13) holds together with the stationary and axial symmetry hypotheses.

Since, we are assuming the plasma disk is quasi-ideal (actually it is true in many range of observed mass density and temperature), we neglect the poloidal velocities, especially the radial component, which are due to effective dissipation, according to the Shakura idea of accretion discussed in the previous section. Furthermore, we observe that the co-rotation condition (13) prevents the emergence of an azimuthal magnetic field component via the dynamo effect.

We now split the magnetic flux function around a fiduciary radius \(r = r_0\), as follows

$$\begin{aligned} \psi = \psi _0(r_0) + \psi _1 (r_0, r - r_0, z) \, , \end{aligned}$$
(14)

where \(\psi _0\) is the vacuum contribution of the central object around which the disk develops (essentially a vertical magnetic field comes out from the dipole-like nature of the field and from the thinness of the disk), while \(\psi _1\) is a small (still steady) correction, here considered of very small scale with respect to the background quantities. By other words, we are studying a small backreaction of the plasma which is embedded in the central object magnetic field, whose spatial (radial and vertical) scales are sufficiently small to produce non-negligible local currents in the disk.

According to Eqs. (14) and (12), also the magnetic field is expressed as \(\mathbf{B} = \mathbf{B} _0 + \mathbf{B} _1\). According to the validity of the co-rotation theorem, at any order of perturbation of the steady configuration, we expand the angular velocity as follows

$$\begin{aligned} \omega (\psi ) \simeq \omega _0(\psi _0) + \left( d\omega /d\psi \right) _{\psi = \psi _0}\psi _1 \, . \end{aligned}$$
(15)

In Ref. [10] (see also Refs. [11, 12]), it was shown that, in the linear regime, i.e., \(|\mathbf{B} _1|\ll |\mathbf{B} _0|\), the equilibrium configuration near \(r_0\) reduces to the radial equilibrium only, which, at the zeroth and first order in \(\psi \), gives the following two equations

$$\begin{aligned}&\omega _0(\psi _0) = \omega _K \, , \end{aligned}$$
(16)
$$\begin{aligned}&\partial ^2_r \psi _1 + \partial ^2_z\psi _1 = -k_0^2( 1 - z^2/H^2)\psi _1\, , \end{aligned}$$
(17)

respectively. We recall that H denotes the half-depth of the disk and \(k_0\sim 1/\lambda _{\min }\) is the typical wave-number of the small scale backreaction, taking the explicit form

$$\begin{aligned} k_0^2 \equiv 3\omega _K^2/v_A^2 \, , \end{aligned}$$
(18)

and, for a thin isothermal disk, it results \(k_0H = \sqrt{3}\beta \equiv 1/\epsilon \). As postulated above, in order to deal with small scale perturbations, we have to require that the value of \(\beta \) is sufficiently large.

In order to study the solutions of Eq. (17), it is convenient to introduce dimensionless quantities, as follows

$$\begin{aligned} Y \equiv \frac{k_0\psi _1}{\partial _{r_0}\psi _0} \, ,\;\; x \equiv k_0 \left( r - r_0\right) \, ,\;\; u = z/\delta \, , \end{aligned}$$
(19)

where \(\delta ^2 = H/k_0\). Hence, we get

$$\begin{aligned} \partial ^2_xY + \epsilon \partial ^2_uY = -\left( 1 - \epsilon u^2\right) Y \, , \end{aligned}$$
(20)

which, in what follows we dub the Master Equation for the crystalline structure of the plasma disk. We recall that, according to Eq. (5), the u variable takes a finite range across the thin disk configuration. Moreover, it is also easy to check the validity of the relations

$$\begin{aligned}&B_z = B_{0z}\left( 1 + \partial _xY\right) \, , \end{aligned}$$
(21)
$$\begin{aligned}&B_r\equiv B_{1r} = -B_{0z} \sqrt{\epsilon }\partial _uY \, , \end{aligned}$$
(22)

where the existence of linear perturbation regime requires \(|Y|\ll 1\).

It is relevant to investigate the solutions of Eq. (20) in view of determining the physical conditions under which the crystalline structure, discussed in Refs. [10, 29], can actually take place.

3.1 Reconciling accretion with magnetic micro-structures

Let us now clarify how the magnetic micro-structures can be relevant for the accretion picture of plasma disk. Actually, in Refs. [10,11,12] the profile of the backreaction on small scale has been analyzed in the presence of differential angular velocity only, i.e., without any poloidal matter flux and therefore accretion features are not addressed.

In the original idea, proposed in Refs. [30, 31], the formation of magnetic micro-structures was intended to be an intriguing paradigm in which accretion could be sustained also in absence of dissipative effects, only on the basis of the ideal plasma morphology. The crucial point of this reformulation consists of the possibility that, in the case of nonlinear plasma backreaction, a large number of X-points can form. In fact, when the vertical magnetic field is dominated by the backreaction, it structure is intrinsically characterized by a radial oscillation and therefore X-points appear with the same periodicity of the magnetic surface oscillation. Clearly, near such points, the plasma disk manifests its porosity, simply because the prescription of the ideal electron force balance \(v_rB_z = 0\) allow now for nonzero radial velocity also in the absence of turbulent viscosity. Of course, in order such a plasma porosity becomes an efficient tool for the disk accretion, a mechanism able to pump the plasma into the X-point is required. In Ref. [30], it was postulated that the pumping of plasma could be supported by intermittent ballooning instabilities, i.e., modulating a naturally mechanism observed in laboratory plasma physics to the astrophysical context scenario.

Furthermore, in Ref. [8], it was emphasized how, limiting attention to a one-dimensional model, a natural inconsistency appears between the Shakura idea of accretion and the emergence of magnetic micro-structures. The nature of this incoherent formulation of the small scale plasma backreaction into a visco-resistive accretion scenario can be easily recognized by observing that, in the Ohm law, the largest contribution to the radial infalling velocity comes from the induced oscillating current densities and therefore, it has, in turn, an oscillating character too (which contradicts the smooth Shakura radial infalling). In a more recent two-dimensional formulation [32] (see also Ref. [33]), it has been reconciled the disk magnetic micro-structure with the Shakura scenario of accretion, by making use of the smallness of the poloidal velocity field with respect to the toroidal differential rotation of the disk. Addressing the generalized Grad-Shafranov equation for an accretion disk (see Ref. [19]), we have overcome the inconsistency observing that the contribution of the microscale phenomenon must be averaged out on the disk macroscales. As a result, the effect of the current induced by the backreaction allows a proper balance of the Ohm law, even in the absence of resistivity high values, but the induced poloidal velocity is averaged to zero on a macroscopical scale.

This scenario suggests that the emergence of magnetic micro-structures across the accreting disk could provide a viable alternative paradigm to the necessity of the so-called disk anomalous resistivity, without rejecting the solid scenario of angular momentum transport as driven by the effective viscosity associated to the turbulence that MRI is able to trigger. However, the presence of the small scale backreaction requires very high values of the plasma \(\beta \) parameter and therefore there is a disk temperature interval in which MRI can still survive, but the disk is too cold to manifest magnetic micro-structures. In such a region of the disk parameters, the necessity of an effective large resistivity to account for accretion onto compact objects, like X-rays binary stars, can probably no way avoided.

4 A viable solution to the Master Equation

We investigate the solution of Eq. (20) for the crystalline morphology where we consider the regime \(\epsilon \ll 1\), corresponding to large values of the \(\beta \) plasma parameter. Since the crystalline structure is associated to a periodicity in the x-dependence, we naturally search a solution for the dimensionless first-order magnetic flux function Y in the form

$$\begin{aligned} Y(x,u) = \sum _{n=1}^{n=\infty } F_n(u) \sin (nx + \varphi _n(u)) \, , \end{aligned}$$
(23)

where \(F_n\) and \(\varphi _n\) are amplitude and phase, respectively. Above, we included a u-dependent phase \(\phi _n(u)\) in order to properly weight both \(\sin \) and \(\cos \) functions in the Fourier expansion. Once set to zero the coefficients of the different trigonometric functions, this expression provides the two following equations

$$\begin{aligned}&2\frac{\hbox {d}F_n}{\hbox {d}u}\frac{\hbox {d}\varphi _n}{\hbox {d}u} + F_n\frac{d^2\varphi _n}{\hbox {d}u^2} = 0 , \end{aligned}$$
(24)
$$\begin{aligned}&\epsilon \frac{\hbox {d}^2F_n}{\hbox {d}u^2} -\epsilon F_n \Big (\frac{\hbox {d}\varphi _n}{\hbox {d}u}\Big ) ^2 = -\left[ (1-n^2) - \epsilon u^2\right] F_n . \end{aligned}$$
(25)

Equation (24) admits the solution

$$\begin{aligned} \hbox {d}\varphi _n/\hbox {d}u = \varTheta /F_n^2, \end{aligned}$$
(26)

with \(\varTheta =const.\), which reduces Eq. (25) to the following closed form in \(F_n\)

$$\begin{aligned} \epsilon \frac{\hbox {d}^2F_n}{\hbox {d}u^2} -\epsilon \frac{\varTheta ^2}{F_n^3} = -\left[ (1-n^2) - \epsilon u^2\right] F_n . \end{aligned}$$
(27)

We observe that, if the periodicity is not associated to the fundamental wave-number \(k_0\) but to a close one \(k = \chi k_0\) (where \(\chi \simeq 1\)), the equation above is simply mapped by replacing the integer n by \(\chi n\). For \(\epsilon \ll 1\), it is clear that we must have \(n = 1\) and then we assume \(1-\chi ^2n^2 = \epsilon \). By this choice, Eq. (27) becomes independent of \(\epsilon \), reading as (noting \(F\equiv F_1\))

$$\begin{aligned} \frac{\hbox {d}^2F}{\hbox {d}u^2} - \frac{\varTheta ^2}{F^3} = - \left( 1 - u^2\right) F . \end{aligned}$$
(28)

We underline how we set \(1-\chi ^2 n^2 = \epsilon \), thus removing a free parameter from the model, in order to properly recover, for \(\varTheta = 0\), the same behavior discussed in Ref. [10], and it corresponds, from a physical point of view, to deal with a radial wave-number near \(k_0\). In the next section, we will study the Master Equation by using the separation variable method and limiting therefore attention to linear behaviors only. The present case is recovered when the separation constant \(\gamma \) obeys the relation \(\gamma = 1 - \epsilon \), again ensuring that the radial wave-number, fixed by the \(\gamma \) value, remains close to \(k_0\). There we will provide a systematic analysis of the mathematical problem for different positive values of \(\gamma \) in order to span all the parameter space, without focusing attention to physical constraints. However, it remains clear that, for sufficiently low values of \(\gamma \), the short wave-length approximation could break down and, on the contrary, for too large \(\gamma \) values, the backreaction scale would become smaller than the Debye length, invalidating the MHD scheme. Such peculiar situations take place in rather extreme conditions and they are fixed by the details of the disk configuration, so that we will not account for them in the next section.

Let us consider the problem of solving Eq. (28) with boundary condition \(F(0)={\bar{F}}\), being \({\bar{F}}\) a generic constant (we consider \(0<{\bar{F}}\ll 1\) since such an equation is invariant for \(F\rightarrow -F\) and because the linear regime must be preserved) and \(F'(0)=0\) (here and in the following the prime denotes u derivatives). This last condition has been selected in order to guarantee the symmetry of the solution with respect to the equatorial plane. The plots of the numerical solution, in correspondence of different values of \(\varTheta \) at fixed \({\bar{F}}\), are in Fig. 1. Such a solution is a Gaussian for \(\varTheta =0\) as derived in Ref. [10], while for \(\varTheta \ne 0\) it starts to increase for increasing |u| values. It is worth recalling that, when F is greater than unity, the linear approximation fails and the Master Equation is no longer predictive.

Fig. 1
figure 1

Solution of Eq. (28) for the amplitude F of the perturbed flux function as a function of the normalized vertical height of the disk u. We fix \({\bar{F}}=0.01\) varying the parameter \(\varTheta \) as indicated in the plot

Full size image

We note that, only for \(\varTheta = 0\), an exponentially decaying behavior of the linear perturbation above the equatorial plane is obtained, just like the background profile [7]. Such a natural similarity in the vertical behavior of the background and of the perturbations was also at the ground of the analyses in Refs. [10, 11, 29], but the solution derived above demonstrates how this feature is a peculiar property of the simplest case only, i.e., when the phases \(\varphi _n\) are constant quantities.

We observe that in Eq. (26), in the linear regime for \(F_n\ll 1\), we must require that the constant \(\varTheta \) be correspondingly small in order to avoid an exceedingly large value of the u-derivative of the phase \(\phi _n\). Despite the effective \(k_z\) value is rather large in the considered scenario (since it is estimated as \(k_z \sim \delta ^{-1}\sim 1/H\sqrt{\epsilon }\)), however, we can explore only region of u-values \(u\sim k_z\) limited by the request that the linearity of the Master Equation is not broken. In this respect, the increasing behavior of F with large absolute values of u is physically acceptable. Such a behavior can be considered as predictive in the present configurational scheme, only up to the validity of the linearity \(Y\ll 1 \Rightarrow F\ll 1\), so that such an increasing behavior is never associated to a diverging non-physical feature: our model is a local one, we are exploring a small region around a fiducial \(z_0\) value and only in that neighborhood, our linear equilibrium makes sense. Similar considerations hold also for the radial dependence, when the oscillating behavior were mapped into an increasing one with the radial coordinate x, namely the trigonometric function could become hyperbolic ones (see the next section).

We also that, from a physical point of view, the increasing behavior is not surprising. In fact, it can be explained because, while the background profile is fixed by a gravostatic equilibrium, the perturbation morphology has nothing to do with gravity and it is mainly produced by the balance of the centripetal (say centrifugal) and the Lorentz forces.

However, it is important to stress that for sufficiently small values of \(\varTheta \) a peak of the function F around the equatorial plane is still present. This feature takes place because for small values of \(\varTheta \) the effect of the nonlinear term \(\varTheta /F^3\) is minimized.

From a physical point of view, a process able to induce a crystalline structure in the plasma disk (for instance a sound or a gravitational wave), i.e., a boundary condition able to select the solution associated to a radial corrugation of the magnetic flux function, should be vertically coherent in phase to produce the same vertical gradient of the background profile in the perturbation too. On the contrary, a vertical shear also in the phase of the underlying process is responsible for an inversion of the vertical gradient profile between the background and the perturbation.

The most important physical conclusion of the present analysis is to demonstrate that the crystalline structure is a very general feature of the linearized perturbed equilibrium, but, being associated to small values of \(\epsilon \) (i.e., large values of the plasma \(\beta \) parameter), the radial oscillation can be associated to a basic wave-number only, very close to the natural one \(k_0\).

5 Separable solution of the Master Equation

Since the approach in the previous section gave us the solution \(Y = F(u)\sin (x)\) as the only viable case for \(\epsilon \ll 1\), it is natural to investigate solutions relying on the separation of variable method. In particular, we will be able to explore the situation in which the wave-number of the radial oscillation (x-dependence) is a generic constant times the fundamental wave-number \(k_0\).

In order to investigate solutions relying on the separation of variable method, the Master Equation (Eq. 20) can be rewritten using \(Y(x,u)=Y_1(x) Y_2(u)\) in the following decoupled form:

$$\begin{aligned}&(\hbox {d}^2Y_1/\hbox {d}x^2)/Y_1=-\gamma , \end{aligned}$$
(29)
$$\begin{aligned}&\epsilon (\hbox {d}^2Y_2/\hbox {d}u^2)/Y_2+(1-\epsilon u^2)=\gamma , \end{aligned}$$
(30)

where we have introduced the arbitrary constant \(\gamma \). The first equation has two possible kind of solutions, namely

$$\begin{aligned}&Y_1(x)= C_1 \sin (\sqrt{|\gamma |} x)+ D_1 \cos (\sqrt{|\gamma |} x),\; \gamma >0, \end{aligned}$$
(31)
$$\begin{aligned}&Y_1(x)= C_1 \sinh (\sqrt{|\gamma |} x) + D_1 \cosh (\sqrt{|\gamma |} x) ,\; \gamma <0. \end{aligned}$$
(32)

Equation (30) can be rewritten in the form of the Kummer equation of parabolic cylinder functions using the scaling \(u\rightarrow {u}/\sqrt{2}\). In particular, we get

$$\begin{aligned} \frac{\hbox {d}^2Y_2}{\hbox {d}u^2} - \Big ( \frac{1}{4} {u}^2 + A \Big ) Y_2=0, \end{aligned}$$
(33)

where \(A=(\gamma -1)/(2\epsilon )\).

Fig. 2
figure 2

Solution of Eq. (33) for the vertical (u) dependent part \(Y_2\) of the perturbed flux function. We use \(Y_2(0)=0.01\) and \(Y_2'(0)=0\), by varying the parameter A in \([1,\,0,\,-0.4,\,-6.5,\,-25]\) as indicated in the plot

Full size image
Fig. 3
figure 3

Plot of the perturbed magnetic flux function \(Y(x,u)=Y_1(x)Y_2(u)\) solution of Eq. (31) (with \(C_1=1\) and \(D_1=0\)) and Eq. (33) for \(Y_2(0)=0.01\), \(Y_2'(0)=0\), \(A=-11\) (i.e., \(\gamma =0.78\))

Full size image
Fig. 4
figure 4

Cross sections of the function \(Y(x,u)=Y_1(x)Y_2(u)\) form Fig. 3 as function of the vertical coordinate u for fixed values of the radial coordinate x (upper panel) and as function of x at given u (lower panel), as indicated in the plots

Full size image
Fig. 5
figure 5

Cross section of the current density distribution \(J=\partial _x^2 Y+\epsilon \partial _u^2Y\) provided by Eq. (20) with the solution of Fig. 3. Panel structure as Fig. 4

Full size image

Let us now analyze the properties of the solution by numerically integrating Eq. (33) using initial conditions \(Y_2(0)=0.01\) and \(Y_2'(0)=0\). We also fix the parameter \(\epsilon =0.01\) and run the \(\gamma \) parameter in order investigate the morphology of the solution for different values of the constant A. From the numerical analysis, a range of the parameter A emerges (dependent on the choice of \(\epsilon \)) for which the behavior of the function \(Y_2\) becomes to oscillate, as well sketched in Fig. 2. Since these values of A still corresponds to positive values of the constant \(\gamma \) (for our setup, \(\gamma \lesssim 0.9\)), we obtain a range of the solution of the Master Equation where the radial and vertical dependence of the perturbed magnetic flux function both have an oscillating profile as represented in Fig. 3 and in the detailed cross sections in Figs. 4 and 5.

Clearly, when \(\gamma \) takes negative values, the x dependence of the function Y is no longer oscillating (while the u-dependence still remains oscillating) and the crystalline morphology of the disk is suppressed. As we shall discuss below, the oscillating nature of the vertical dependence of the perturbed magnetic flux function has the relevant physical implication that, according to the expression for the magnetic field in Eqs. (21) and (22), a series of O-points appears in the disk configuration, in which the radial component \(B_r\) vanishes.

6 Whittaker equation and O-points

The many different behaviors of the solution of Eq. (33) can be better understood by investigating peculiar properties of the corresponding analytic solutions. In particular, using the following change of variables: \(y={u}^2/2\) and \(Y_2=(2y)^{-1/4} W(y)\), Eq. (33) rewrites now as

$$\begin{aligned} \frac{\hbox {d}^2 W}{\hbox {d} y^2}+\Big (-\frac{1}{4}- \frac{A}{2 y}+\frac{3}{16 y^2}\Big )W=0, \end{aligned}$$
(34)

which is a particular case of Whittaker’s equation, i.e.,

$$\begin{aligned} \frac{\hbox {d}^2 W}{\hbox {d} y^2}+\Big (-\frac{1}{4}+ \frac{K}{y}+\frac{1/4 -\mu ^2}{y^2}\Big )W=0, \end{aligned}$$
(35)

with \(K=-A/2\) and \(\mu =\pm 1/4\). The linearly independent solutions of this equation are given by

$$\begin{aligned} W^{(1)}_{K,\mu }&= e^{-y/2} y^{\mu +1/2} M(1/2+\mu -K,1+2\mu ,y), \end{aligned}$$
(36)
$$\begin{aligned} W^{(2)}_{K,\mu }&= e^{-y/2} y^{\mu +1/2} U(1/2+\mu -K,1+2\mu ,y), \end{aligned}$$
(37)

where M and U are the Kummer and Tricomi functions [34], also known as confluent hypergeometric functions of the first and second kind, respectively. We will use only the function M because the function U is multivalued, and we study the solutions only for finite y since the model holds for a limited range of u values. Moreover, for \(\gamma \le 1\), the function \(W_{K,\mu }\) is bounded since \(M(1/2+\mu -K,1+2 \mu ,y)\) behaves like

$$\begin{aligned} \varGamma (1+ 2 \mu ) \frac{e^{y/2} y^{-1/2-\mu -K}}{\varGamma (1/2+\mu -K)}+\frac{(-y)^{-1/2-\mu +K}}{\varGamma (1/2+\mu +K)}, \end{aligned}$$

where \(\varGamma \) denotes the standard gamma function. Specifically, the Kummer function M(a, b, y) is formally defined as

$$\begin{aligned} M(a,b,y)=1+ \frac{a_{(1)}}{b_{(1)}}y+ \dots + \frac{a_{(n)}}{b_{(n)}}\frac{y^n}{n!}, \end{aligned}$$

where the rising factorial denoted with the index (n) is defined as

$$\begin{aligned}&X_{(0)}=1,&\quad X_{(n)}= X (X+1)(X+2)\dots (X+n-1). \end{aligned}$$

The sum is converging everywhere, while the derivative of M with respect to y is given by the following recurrent equation: \(d M(a,b,y)/dy=(a/b) M(a+1,b+1)\).

Using this formalism, we get the following solutions of Eq. (33):

$$\begin{aligned} Y_2^+(u)&=c_1^{+} u \,e^{-u^2/4} M(A/2+3/4, 3/2; u^2/2), \end{aligned}$$
(38)
$$\begin{aligned} Y_2^-(u)&=c_1^{-} e^{-u^2/4}\; M(A/2+1/4,1/2;u^2/2), \end{aligned}$$
(39)

where \(c^{\pm }_{1}\) are integration constants and we denoted with the choice \(\mu =\pm 1/4\).

The zeros of the derivative of the solutions above correspond to that one of the radial component of the magnetic field from Eq. (22). Using recursive expressions, we get

$$\begin{aligned}&\hbox {d}Y_2^+/\hbox {d}u\nonumber \\&\quad =c_1^{+} e^{-u^2/4}[(1-u^2/2) M(A/2+3/4,3/2; u^2/2)+\nonumber \\&\qquad +u e^{-u^2/4}(A/3+1/2) u M(A/2+7/4, 5/2; u^2/2)], \end{aligned}$$
(40)
$$\begin{aligned}&\hbox {d}Y_2^-/\hbox {d}u\nonumber \\&\quad =c_1^- e^{-u^2/4}(-u/2) M(A/2+1/4,1/2; u^2/2)+\nonumber \\&\qquad +c_1^{-} u e^{-u^2/4}(A+1/2) M(A/2+5/4, 3/2; u^2/2). \end{aligned}$$
(41)
Fig. 6
figure 6

Plot of the magnetic flux function component \(Y_2^{-}(u)\) in Eq. (39), corresponding to \(\mu =-1/4\), for different values of the parameter A as in the numerical analysis of Fig. 2. The integration constant \(c_1^{-}\) is set to match the boundary condition in \(u=0\)

Full size image
Fig. 7
figure 7

Plot of the function \(Y_2^{+}(u)\) in Eq. (38), corresponding to \(\mu =1/4\) (parameter choice as in Fig. 6)

Full size image
Fig. 8
figure 8

Behavior of \(\hbox {d}Y_2^-/\hbox {d}u\) in Eq. (41) corresponding to parameter choice of Fig. 6

Full size image

Since in the general solution both \(\mu =-1/4\) and \(\mu =1/4\) have to appear in a linear combination, we plot in Figs. 6 and 7 the function \(Y_2(u)\) for both these values, respectively. Comparing these two pictures with the numerical integration in Fig. 2, we see that the choice \(\mu =-1/4\) better addresses the profile associated to the considered physical boundary conditions of the numerical integration. Moreover, it is worth underline the following behavior depending on the values of the A parameter (we focus only on the half-plane \(u\geqslant 0\)): for \(-25<A<-16\), \(Y_2^-\) oscillates; for \(-16<A<0\), the oscillations disappear and the solution has only some maxima or minima; for \(A\geqslant 0\) it is an increasing function. The transitions of the analytic solution are in agreement with those of the numerical solution of the previous section. In Fig. 8 we instead plot the graph of \(\hbox {d}Y_2^-/\hbox {d}u\). The O-points of the radial component of the magnetic field on the vertical axis clearly emerge, placed in symmetric positions with respect the equatorial plane. For the specific case \(A=-25\), the positive zeros of \(B_r\) can be found at \(u=\)[0, 0.6, 1.3, 1.9, 2.55, 3.2, 3.9, \(\ldots \)]. In this respect, we note that algorithms are present for systematically finding the roots of Kummer functions [34].

6.1 Asymptotic solutions

Taking into account the solutions for \(Y_2(u)\) in Eqs. (38) and (39), let us now define the following functions

$$\begin{aligned} Z(A,u)=&\varPhi _1(u)\cos (\pi (1/4+A/2))+\nonumber \\&-\varPhi _2(u)\sin (\pi (1/4+A/2), \end{aligned}$$
(42)
$$\begin{aligned} V(A,u)=&\frac{1}{\varGamma (1/2-A)}(\varPhi _1(u)\sin (\pi (1/4+A/2))+\nonumber \\&+ \varPhi _2(u)\cos (\pi (1/4+a/2)), \end{aligned}$$
(43)

where

$$\begin{aligned}&\varPhi _1(u)=\frac{\varGamma (1/4-A)}{\sqrt{\pi }\;2^{A/2+1/2}}Y_2^+, \end{aligned}$$
(44)
$$\begin{aligned}&\varPhi _2(u)=\frac{\varGamma (3/4-A/2)}{\sqrt{\pi }\;2^{A/2-1/2}}Y_2^-. \end{aligned}$$
(45)

In the asymptotic limit \(-A\gg u^2\) (for \(A<0\)) and finite u (we recall that the thin-disk hypothesis limit the range of u values), we obtain

$$\begin{aligned}&Z(A,u)+i \varGamma (1/2-A) V(A,u)\nonumber \\&=\frac{e^{i \pi (1/2+1/2 A)}}{2^{(1/2+1/2 A)}\sqrt{\pi }}\varGamma (1/4- A/2)e^{\theta _r+i \theta _i}e^{i p u}, \end{aligned}$$
(46)

where \(p=\sqrt{-A}\), while \(\theta _r(u)\) and \(\theta _i(u)\) are assigned polynomials. Thus, the function Z(A, u), the solution \(Y_2(u)\) and its derivative with respect to u oscillate giving rise to a sequence of O-points across the vertical axis.

For positive values of A, we get instead the following asymptotic behavior for \(A\gg u^2\):

$$\begin{aligned} Z(A,u)=\frac{\sqrt{\pi }}{2^{(1/2+1/2 A)}\varGamma (3/4+ A/2)}e^{-\sqrt{A} u+ \theta _1}, \end{aligned}$$
(47)

where \(\theta _1\) is again an assigned polynomial. This case for \(A>0\) is not relevant for the present physical discussion.

6.2 Laguerre polynomials

We have shown how, in the case \(-A\gg u^2\), there is a countable set of O-points, but it is also possible to find other zeros series in the correspondence of smaller |A|. It is easy to verify that the following relation takes place between the Kummer function M(a, b, y) and the Laguerre polynomials \(L_n^{(\alpha )}\):

$$\begin{aligned} M(-n,\alpha +1,y)=L_n^{(\alpha )}(y)/\left( {\begin{array}{c}n+\alpha \\ n\end{array}}\right) . \end{aligned}$$
(48)

For \(\mu =-1/4\), Eq. (39) contains the function \(M(A/2+1/4,1/2;u^2/2)\), and we easily get the specific values: \(n=-A/2-1/4\) and \(\alpha =-1/2\). We thus obtain the following general expression for the derivative of \(Y_2^-\):

$$\begin{aligned}&\hbox {d}Y_{2}^-/\hbox {d}u=-c_1^- e^{-u^2/4}\nonumber \\&\quad \times \Big (u L_{n-1}^{(1+\alpha )}(u^2/4)+(u/2)L_{n}^{(\alpha )}(u^2/4)\Big )/\left( {\begin{array}{c}n+\alpha \\ n\end{array}}\right) . \end{aligned}$$
(49)

The equation above is a particular form of Eq. (41) and, since it is well-known that for each n there is a finite number of zeros of the Laguerre polynomials which increases with n, we have shown how a countable set of zeros on the z axis emerges. In particular, the zeros of Laguerre polynomial \(L_n^{(\alpha )}\) belong to the interval \((0, n+\alpha +(n-1) \sqrt{n+\alpha })\), and in Fig. 9 we plot, as an example, \(\hbox {d}Y_{2}^-/\hbox {d}u\) from Eq. (49) for the first 10 integer n values.

Fig. 9
figure 9

Plot of \(\hbox {d}Y_{2}^-/\hbox {d}u\) from Eq. (49) (zoom of the half-plane \(u>0\)) by setting integers n in [1, 10]. Standard color scheme from red (\(n=1\)) to blue (\(n=10\)). The number of oscillations increases with increasing n

Full size image

7 Discussion

We have shown how the z-dependence of the magnetic field is associated to O-points of the corresponding configuration. In such points, the radial component of the magnetic field, due to the plasma backreaction only, vanishes. Here, we outline a new feature induced by the crystalline morphology of the plasma response to the dipole field of the central object.

Let us now discuss the relevance of the O-points derived above, in view the physical properties of the plasma disk thought as an accreting structure. Our analysis has been performed for an ideal plasma disk in the absence of poloidal components of the velocity field, i.e., \(v_r = v_z \equiv 0\). However, in the Standard Model of accretion [1], the plasma must have a finite electric conductivity and the azimuthal component of the generalized Ohm law reads

$$\begin{aligned} v_zB_r - v_rB_z = J_{\phi }/\sigma , \end{aligned}$$
(50)

where \(J_{\phi }\) is the corresponding azimuthal component of the current density and \(\sigma \) is the constant coefficient of electric conductivity. In the standard mechanism for accretion onto a compact object, both the quantities \(B_r\) and \(v_z\) are essentially negligible and, where \(\sigma \) takes very large values, the accretion is suppressed, i.e., since \(B_z\ne 0\), we must get \(v_r\simeq 0\).

In the spirit of the present analysis, in those regions where the plasma is quasi-ideal, we find the relation

$$\begin{aligned} v_z = v_r B_z/B_r. \end{aligned}$$
(51)

Clearly, nearby the numerous O-points of the magnetic profile, where \(B_r\) is nearly vanishing, we get very high values of the vertical velocity, which suggests the existence of privileged plasma sites for the jet formation. This perspective has been investigated in some detail in Refs. [13, 14].

It is worth observing that, when an anomalous resistivity is considered, the term \(v_z B_r\) can be actually safely neglected since the balance of the Ohm law is guaranteed by the relation \(v_rB_z = - J_{\phi }/\sigma \). Of course, if we can state that \(v_zB_r\sim v_rB_z\), we naturally have large values of \(v_z\) (due to the smallness of \(B_r\) in a thin disk). However, this situation can not be regarded as the most general one, unless some other information on the velocity field is assigned by extra physics in the disk. On the contrary, when the values of sigma are very large, the balance between the two terms is mandatory and \(v_z\) is very large near O-points. Clearly, also the solution \(v_r\sim v_z = 0\) is viable, but this situation is just that one discussed in Sect. 3, where differential disk angular velocity is included only. However, in order to get the equality \(v_rB_z = v_z B_r\) and the O-points, simultaneously, \(J_{\phi }/\sigma \) must be a negligible contribution, i.e., \(k_0\) can not exceed a critical magnitude fixed by the value of the electric conductivity \(\sigma \). If not, we can still have a large value of \(v_z\), but again in specific cases only.

Actually, the study of the crystalline micro-structures in the presence of poloidal velocity has not yet been fully investigated, see Ref. [32] for recent developments (the most important difficulty relies on the nonlinear advection terms in the poloidal velocity components). However, the picture associated to the Master Equation is reliable in view of the smallness characterizing the poloidal velocity components with respect to the toroidal rotation velocity, as expected in a thin accretion disk.

The main merit of the present analysis consists of having outlined the general character that the O-points takes in the linear profile of the perturbed magnetic field, as soon as the short wavelength backreaction of the plasma is excited in disk with high \(\beta \) values.

8 Conclusions

We analyzed general features of the small scale morphology of the backreaction that is generated in a thin plasma, embedded in the magnetic field of the gravitationally confining central object. In particular, we studied the Master Equation of the so-called crystalline structure [10], searching for a satisfactory characterization of the admissible profile in the linear backreaction limit.

We have followed two different, but complementary, approaches. The first one is based on a general Fourier expansion, while the second one relies on a separated variable procedure. This analysis offered a complete spectrum of the available solutions, outlining new behaviors in the plasma equilibrium, absent in the basic treatment in Refs. [10, 11].

The main merit of this systematic study of the linear backreaction of the plasma disk consists of the relevant physical implications of the obtained solutions, as discussed above, i.e., the inverse behavior of the vertical gradient with respect to the background one (as the highest harmonics are considered) and the possible emergence of a consistent sequence of O-points, where the radial magnetic field component vanishes.

The first of these issues is of interest in view of the nonlinear scenario of the plasma backreaction discussed in Refs. [11, 12], where, see also Ref. [29], the possibility for the radial fragmentation of the disk into a microscopic array of rings is demonstrated. In this respect, the linear inversion of the vertical gradients suggests that, in the nonlinear regime, a separation of the disk into two symmetric components, up and down the equatorial plane, could take place, with significant implications on the transport features across the disk.

The second result, concerning the appearance of O-point series of the magnetic configuration, is relevant in view of the realization of conditions for jet emission. In fact, nearby the O-points, when the plasma resistivity is suppressed, the vertical component of the plasma velocity can take very large values, see the discussion in Refs. [13, 14] about the seeds of a vertical matter flux from the disk. Using the solution given in terms of hypergeometric confluent functions, we have also analytically shown the existence of the O-point series and such a series critically depends on the parameter of the separation of variable technique.

The systematic analysis here pursued of the Master Equation for the crystalline structure of a stellar accretion disk, constitutes a well-grounded starting point for the investigation of nonlinear features of the magnetic confinement of the plasma disk. In fact, as discussed in Ref. [11], the extreme nonlinear regime is characterized by a dominant influence of the MHD-force in confining the plasma with respect to the radial gravitational field. Furthermore, in Ref. [35] it is outlined how the request that the Master Equation holds in the nonlinear case too turns out to be a proper choice to close the equilibrium system, preserving the crystalline morphology.