© ESO 2022

1. Introduction

Ambipolar diffusion is an important, and at times crucial process at work in magnetised plasmas in which ionised and neutral populations coexist. For a long time, and increasingly in recent years, numerical models of the time-evolution of partially ionised plasmas that include ambipolar diffusion have been constructed. However, the ambipolar diffusion coefficient is proportional to B², and so, by its own nature, it can lead to the formation of sharp, mathematically singular current sheets in the neighbourhood of magnetic null points (Zweibel 1994; Brandenburg & Zweibel 1994). From the point of view of applied mathematics, singularities are a natural outcome of the evolution of nonlinear diffusion equations (see the book by Vázquez 2007). Interestingly, as detailed in the present paper, in the ambipolar diffusion case, the passage of magnetic flux across the singular current sheet occurs at a finite rate through the combination of the vanishing diffusion coefficient and the infinite slope of the magnetic profile. Now, if the time evolution of a system leads to the spontaneous formation of essential singularities that play a relevant role in the system, then the numerical solution of the problem may become quite complicated. Thus, a question arises with regard to the capabilities of state-of-the-art MHD numerical codes in coping with the expected singular current sheets and, in addition, whether it is possible to provide a suite of tests to probe their proficiency when dealing with such singularities.

The plasma in the low atmosphere of the Sun and other cool stars, in the interstellar medium, accretion disks, planetary ionospheres, and in various other cosmic environments is partially ionised – in fact, it is sometimes only weakly ionised (see, e.g. the reviews from Shu et al. 1987; McKee & Ostriker 2007; Zweibel et al. 2011; Crutcher 2012; Leake et al. 2014; Zweibel 2015; Ballester et al. 2018). For example, in crucial layers of the low solar atmosphere (e.g. the high photosphere and chromosphere), the standard VALC model (Vernazza et al. 1981) predicts the ionisation degree of Hydrogen to be below 10⁻². In dynamical models, such as those calculated with the radiation-MHD Bifrost code (e.g. Carlsson et al. 2016), the horizontal averages of the ratio between electron and Hydrogen number density, ⟨n_el/n_H⟩, are also typically below 10⁻² at those heights, in spite of the juxtaposition of cool pockets and hotter structures apparent in the 2D and 3D realistic numerical models that include those layers (see, e.g. Wedemeyer et al. 2004; Hansteen et al. 2007; Leenaarts et al. 2007, 2011; Nóbrega-Siverio et al. 2020a). In a magnetised, partially ionised plasma, it is only the charged species that are subject to the Lorentz force: this causes a drift between the neutrals and the electron-ion gas, which is countered by their mutual friction due to collisions, or by charge-exchange phenomena, or by a combination of both. When the neutral-charged coupling is strong enough, however, the plasma can still be treated as a single fluid, provided the standard Ohm’s law used in simple magnetohydrodynamics is augmented with extra terms (Cowling 1957; Braginskii 1965; Mitchner & Kruger 1973).

Following those basic references, in the simplest case such a generalised Ohm’s law can be written in the following form:

(1)

where u is the centre-of-mass velocity of the plasma, j_⊥ is the component of j perpendicular to B, σ_el is the electrical conductivity, e and n_e are the electric charge and number density of the electrons, respectively, and χ_a is a material coefficient, namely:

(2)

with ρ and ρ_n the total and neutral mass density, respectively, ν_in the ion-neutral collision frequency and m_in the reduced ion-neutral mass. The last two terms in Eq. (1) are the Hall term and the ambipolar term, respectively. To focus on the ambipolar diffusion problem, we can consider the case when the bulk plasma speed is zero (e.g. because the neutrals are at rest and the ionisation level is low) and the Hall and electric conductivity terms are negligible. In that case, one is left with the left-most and right-most terms in Eq. (1) alone. Using Faraday’s induction law, we obtain the corresponding equation for the evolution of the magnetic field:

(3)

with η_a, the ambipolar diffusion coefficient, defined by

(4)

Equations (3) and (4) constitute the pure ambipolar diffusion problem. The term “diffusion” in this context is standard in the literature, even though it is not clear that a general magnetic field configuration governed by Eqs. (3) and (4) should show what would intuitively be associated with simple diffusive behaviour. On the one hand, η_a depends on B²; thus, even the simplest diffusive behaviour we can expect from Eqs. (3) and (4) (e.g. when the field is parallel to itself everywhere) is intrinsically nonlinear. Furthermore, in a general case, the magnetic field has two or three nonzero components and j_⊥ is different from j, so Eq. (3) inextricably entangles the field components in a vector equation that is not a simple diffusion equation. There is extensive literature on this topic, starting with the seminal paper by Mestel & Spitzer (1956) and encompassing hundreds of papers, that studies various aspects of the ambipolar diffusion problem from a physical and astrophysical point of view; for more details, we refer to the reviews mentioned at the beginning of this section. Of particular interest with regard to the present article are those papers studying the sharp electric current sheets that form when the magnetic field profile approaches a null point, such as those of Zweibel (1994) and Brandenburg & Zweibel (1994) mentioned above. These authors have shown how the 1D B profiles governed by the ambipolar diffusion equation and containing a null tend to develop a singular shape in time of the form |B|∝|δ|^1/3, with δ the distance to the null, which could then remain as a stationary solution. They also showed how a differentially rotating magnetic slab tends to develop sharp current sheets between rotating layers – again because of the effect of ambipolar diffusion. Later papers showed how such singular profiles may affect the process of magnetic cancellation and reconnection in them (Brandenburg & Zweibel 1995; Heitsch & Zweibel 2003a,b).

In this paper we concentrate on a simple situation in which Eq. (3) leads to a nonlinear diffusion equation, namely, the case when B points everywhere in the same direction, such as, for instance, along the z-axis, B = B(x, y, t) e_z, for which the equation becomes

(5)

We aim to determine explicit self-similar solutions for this equation when χ_a is a constant; given the constraints, there are relevant cases in one and two dimensions, namely, when the system has either axial symmetry (the cylindrical case) or mirror symmetry (the plane-parallel or Cartesian case). Here, we are dealing with the cylindrical case, while the plane-parallel case will be treated in a follow-up study. Using coordinates (r, φ, z), the equation we solve here is therefore:

(6)

Equations (5) and (6) are of diffusion type with nonlinear diffusion coefficient χ_a B².

Diffusion problems in which the diffusion coefficient is proportional to a positive power (e.g., 2, or 5/2) of the diffusing quantity are important in many different fields of physics in addition to ambipolar diffusion, such as heat conduction in hot plasmas, flows in porous media, boundary layers, thin liquid film spreading, and many others (see the book by Vázquez 2007). In early studies (Zel’dovich & Kompaneets 1950; Barenblatt 1952; Pattle 1959), fundamental self-similar solutions were found in Cartesian, cylindrical, and spherical coordinates: we collectively refer to them as the Zeldovich-Kompaneets-Barenblatt-Pattle solutions (or ZKBP, for short¹). These solutions have compact support (i.e. the set of points where they are nonzero is contained within a finite distance of the origin) and are unsigned, that is, they do not have any internal nulls; the outer boundary has infinite slope and expands in time with a speed which is a natural combination of the parameters of the problem. Since those early works, the properties of the self-similar solutions of the nonlinear diffusion equation have been studied with great generality using phase-plane techniques for the spatial part (see, e.g. Grundy 1979; Hulshof 1991 and the encompassing review by Vázquez 2007). In such works, similarity solutions were considered for the so-called ‘porous medium equation’, both for the 1D and for the radially-symmetric multidimensional cases. A general way of writing that equation is the following:

(7)

with u(x, t) a scalar function of the radial coordinate x = |x|; x the radial vector in n ≥ 1 dimensions; m a fixed parameter generally fulfilling m ≥ 1; and Δ the n-dimensional Laplacian operator. Except for a constant, the particular case for m = 3 and two dimensions coincides with Eq. (6). The results of Hulshof (1991) are of particular interest: that author considered solutions with sign changes, that is, internal nulls, within their domain, but otherwise of the standard self-similar form u(x, t) = t^αU(η), with η = x t^−γ; α and γ two constants; and U(η) the (undetermined) spatial part for the solution. For the case of a solution that expands in space and decays in time (α < 0 and γ > 0), he proved the existence of a countable series of self-similar solutions with compact support, all with a finite number of internal nulls; the first solution in the series is the ZKBP solution, which has no nulls, and all the others have one internal null more than the preceding one in the series. The exponents α and γ in the series are ordered such that the solutions decay faster but expand more slowly the greater the number of nulls in their interior. In his paper, Hulshof (1991) provides this classification, proves the mathematical theorems supporting it, and studies the phase-plane properties of the solutions; however, their specific spatial shapes U(η) are not shown.

In this paper, we are interested not only in the mathematical properties of the ambipolar diffusion as a nonlinear diffusion process, but also in the inclusion of ambipolar diffusion terms in MHD codes. In astrophysics, over the past few decades, multidimensional MHD computer codes have been developed that model a variety of physical processes including ambipolar diffusion. Representative examples of such codes and simulations outside solar physics can be found in Basu & Mouschovias (1994), Mac Low et al. (1995), Padoan et al. (2000), Basu & Ciolek (2004), Kudoh & Basu (2008), Choi et al. (2009), Gressel et al. (2015), Tomida et al. (2015), O’Sullivan & Downes (2007), Masson et al. (2012), Viganò et al. (2019), and Grassi et al. (2019). The consideration of ambipolar diffusion processes in solar physics started many decades ago (e.g. Parker 1963), but it has undergone a true explosion in terms of its use in large numerical models (e.g. Leake et al. 2005; Leake & Arber 2006; Arber et al. 2007; Cheung & Cameron 2012; Leake & Linton 2013; Martínez-Sykora et al. 2012, 2017a,b, 2020a,b; Ni et al. 2015, 2016, 2021; Khomenko et al. 2017, 2018, 2021; González-Morales et al. 2018, 2020; Nóbrega-Siverio et al. 2020a,b; Popescu Braileanu & Keppens 2021). Such numerical calculations often encounter a problem: given the comparatively high values of χ_a in different cosmic environments, the advance in time may grind to a halt in magnetised regions when a standard Courant-Friedrichs-Lewy condition is adopted for the timestep based on η_a. The recent papers by González-Morales et al. (2018) and Nóbrega-Siverio et al. (2020b) describe the construction of so-called super-time-stepping (STS) modules for the Mancha code and for the Bifrost code, respectively, designed with the aim to overcome that stiffness problem. In at least three of the papers cited above (Masson et al. 2012; Viganò et al. 2019; Nóbrega-Siverio et al. 2020b), the basic ZKBP solution in cylindrical coordinates was used to test the ambipolar diffusion module, given its simplicity and the analytical expression available for it. However, as already mentioned, the ambipolar diffusion problem tends to give rise to sharp current sheets and singularities. The ZKBP solution is comparatively smooth in that sense and it would be good to have some other canonical solutions on hand that include current sheets having higher degrees of singularity, which naturally occur in the ambipolar diffusion problem.

The objective of this paper is to calculate signed self-similar solutions of Eq. (6) with compact support and propose that they be used as simple, but nonetheless demanding tests for the ambipolar diffusion modules in MHD codes that include the generalised Ohm’s law. First, we obtain, through separation of variables and numerical calculation, a set of such solutions. As expected from the general considerations above, they have a finite number of nulls within their domain which turn out to be locations of unavoidable singularities; hence, we devised a method that permits a smooth calculation of the solutions across the singularities by numerical means. Secondly, we used some of those solutions to test the capabilities of the new ambipolar diffusion module for the Bifrost code implemented by Nóbrega-Siverio et al. (2020b), thereby showing that these solutions constitute a suitable test for benchmarking MHD codes.

The layout of the paper is: in Sect. 2, we describe the basic equation, the boundary conditions, and the method to calculate the solutions numerically for the cylindrical case. The corresponding solutions are presented and briefly discussed in Sect. 3. Then, in Sect. 4, the use of all those solutions as tests for numerical codes is explained on the basis of the specific example in the Bifrost code. Sect. 5 presents a discussion and conclusions.

2. Finding self-similar solutions

2.1. The self-similar ansatz

To find self-similar solutions of the nonlinear diffusion equation for the case with cylindrical symmetry, Eq. (6), we first note that, just as in the case of heat flux in the nonlinear heat conduction equation, we can here speak of a diffusive flux given by

(8)

which tends to carry magnetic field from more strongly magnetised regions to weaker-field ones. Now, for the separation of variables, we follow standard applied mathematical procedures (e.g. Vázquez 2007), by seeking solutions of the form

(9)

and with L(t) and B₀(t) being global scales for the spatial coordinate and for the magnetic field, respectively. For simplicity, we choose f(ξ = 0) = 1, so B₀(t) corresponds to the value of the solution at the cylindrical axis. If a dot and a prime indicate total derivatives with respect to t and ξ, respectively, Eq. (6) becomes:

(10)

while the diffusive flux defined by Eq. (8) takes the form:

(11)

To proceed, all variables are normalised by choosing values for B₀ and L at a given initial time, t₀₀, say B₀₀ and L₀₀. Also, a natural time unit is defined to be , which is equivalent to setting χ_a = 1 in the foregoing equations. Then we request that all terms in Eq. (10) be time-independent by setting

(12)

with K and H two constants which are positive for decreasing B₀ and growing L, respectively. From Eq. (12), we have:

(13)

with

(14)

and τ_d(t) is an evolving timescale for the diffusion problem:

(15)

Using K and H, Eq. (10) becomes an equation for the spatial part of B, namely,

(16)

To find self-similar solutions we must therefore solve Eq. (16) subject to appropriate boundary conditions. This equation is akin to those found for the general problem in the mathematical literature mentioned in the introduction. In fact, introducing τ_d as time variable instead of t would bring our solutions to the standard self-similar form B(r, τ) = τ^α f(r τ^−γ) (e.g. Zel’dovich & Raizer 1967; Vázquez 2007).

2.2. Constraints at the boundaries

At the axis, the boundary conditions to be imposed are:

(17)

The second condition in Eq. (17) is the usual prescription to prevent the gradient of the solution (i.e. the electric current in this case) from being undefined at the axis. The conditions of Eq. (17) suffice to specify the problem, and, for arbitrary K and H, the solutions of Eq. (16) extend to infinity. Here, we want to find solutions with compact support. Fixing r = L (i.e. ξ = 1) as the outer edge of the solution, an extra boundary condition must be specified there. For continuity, we first impose

(18)

Then, to find the value of the derivative at the outer boundary, we use the condition that the integral of the physical solution remains constant in time; this makes sense at least in the ambipolar diffusion case, in which the integral of the axial magnetic field, namely, the total magnetic flux Φ, is conserved. So:

(19)

In problems for which B does not change sign, this constant differs from zero. This implies B₀ L² = const, so that, from Eq. (13), K/H = 2, which would directly lead to the ZKBP solution. However, new self-similar solutions arise from imposing an alternative possibility, namely,

(20)

for which f(ξ) must pass through one or more nulls. A derivative condition at ξ = 1 can then be obtained by integrating Eq. (16) and applying Eq. (17):

(21)

from which

(22)

This condition is non-trivial, since the spatial self-similar profile can have an infinite slope at the outer boundary. Physically speaking, this is equivalent to saying that the diffusive flux vanishes at the outer edge; that is to say that the magnetic flux remains confined in the (expanding) domain of the solution. To summarise the boundary conditions at the outer edge of the solution, we have:

(23)

2.3. Integration of the equations: usage of the phase plane

Equation (16) is an ordinary differential equation for which standard numerical techniques (such as a second-order extended Euler method or a fourth-order Runge-Kutta scheme) work well almost everywhere. However, our solutions have infinite slope wherever f = 0; the numerical integration can therefore be inaccurate even for small steps of ξ near such nulls. We have, however, devised a change of variable that renders the equation non-singular at the nulls, so that the numerical procedure becomes smooth everywhere: when nearing the nulls one can change the independent variable from ξ to f (admissible wherever f′≠0), use the diffusive flux D as auxiliary variable, and split Eq. (16) into two first-order equations:

(24)

From Eq. (24), at an arbitrary null located at ξ_n, the slope of the D(f) function reaches a finite value, ξ_n H, which proves the adequacy of (24) to carry the numerical integration through the nulls. In turn, however, Eq. (24) becomes singular at the extrema of f, where D = 0. Hence, when calculating solutions, the procedure to follow is to switch from Eq. (16) to Eq. (24), or viceversa, at intermediate locations between each extremum and node (chosen not too close to the singular points or nodes). Using this strategy, the numerical integration becomes smooth.

3. Solutions

3.1. The eigenvalue problem

The solutions of Eq. (16) subject to the inner boundary conditions of Eq. (17) for arbitrary positive values (K, H) in general extend to infinity, asymptotically behaving like f ∝ ξ^−K/H (since the f² term decreases with ξ much more rapidly than the others). To understand how to extract solutions that also fulfil the outer boundary condition of Eqs. (18) and (22) at a finite radius ξ = 1, we change the independent variable to , so that Eq. (16) becomes

(25)

the prime now understood as the derivative with respect to . Figure 1 shows solutions for Eq. (25) for five values of in the interval [4.0, 7.2] that are suitably chosen so as to illustrate the behaviour of the solutions. Starting from , all curves have an oscillatory range, with one or two crossings of the horizontal axis, together with an asymptotic tail extending to infinity. The change from one regime to the other is quite abrupt, apparent almost as a corner when it occurs for small f, that is, for near the value for which the asymptotic tail goes from below to above the axis.

Fig. 1. Selected solutions of the fundamental equation in an infinite domain subject to the inner boundary conditions of Eq. (17) with asymptotic shape toward infinity that illustrate the identification of the eigenvalues and eigenfunctions for the fundamental equation subject to the full boundary conditions of Eqs. (17), (18), and (22).

In fact, for (red curve), the asymptotic tail is exactly zero and the end of the oscillatory range occurs at a precise point (, dashed vertical line) with an infinite slope. At the junction, the condition (f³)′ = 0 applies, since (f³)′ changes sign when varying across the critical . The solution for the critical therefore obeys the equation and both the inner and outer boundary conditions, thus constituting an eigenfunction with just one internal crossing of the axis (at ). Repeating this procedure for solutions with more crossings, we find an infinite, but countable collection of eigenfunctions which solve the full problem (as predicted by Hulshof 1991). To transform the eigenfunctions back to unstretched coordinates, we set , , , so that these solutions fulfil the original Eq. (16) and all boundary conditions of Eqs. (17), (18), and (22).

3.2. The eigenfunctions

Figure 2 (top panel) presents the solutions of the eigenvalue problem for the fundamental and the first four harmonics; the fundamental eigenfunction is the ZKBP solution, which has the following simple profile:

(26)

Fig. 2. First five eigenfunctions as a function of the stretched coordinate . The actual solutions are shown in the upper panel, while the corresponding dimensionless diffusive flux D is displayed in the lower panel. The locations of the null crossings of the solutions are indicated with vertical dashed lines. In both panels, the eigenvalues are indicated.

In the bottom panel, the behaviour of −f² f′ (and so, of F_d) is presented for the five solutions shown in the upper panel, with the null crossings marked by vertical dashed lines. Detailed runs of f and −f² f′ as functions of for the first five harmonics is provided in tabular form as attached files at a public repository² Table 1 provides the eigenvalues and the locations of the end point, , for those harmonics and the fundamental. For clarity in the diagram, the fifth harmonic is not shown in Fig. 2.

The functions resemble standard Sturm-Liouville solutions, except that the crossings of the horizontal axis and the external boundary point have an infinite slope. These singularities are essential, as can be seen by considering the integral of the magnetic field in Eq. (19) restricted to the interval between two successive internal nodes, (ξ_n, ξ_n + 1), or between the origin and the first node. Since the positions of the nodes are fixed, the spatial part of the integral, , is constant. However, the product in front of the integral, B₀(t) L²(t), is a decreasing function of time, as easily proved from Eq. (13) and noting that increases upward of 2 as the harmonic number is increased. Now, when the diffusing physical variable is a conserved quantity, the only possibility for its integral between two nodes (or between the origin and the first node) to decrease in absolute value is through diffusive flux to or from the neighbouring lobe: , so that positive and negative values from neighbouring lobes mutually cancel. The need for a singularity at the nodes follows from that condition.

We can understand the nature of the singularities at internal nulls through an asymptotic approximation. Choosing any of the nodes, say the one located at ξ_n, and calling δ the distance from it, namely, , we can show that close to the node (more precisely: for |δ|/ξ_n ≪ 1, |fK|≪|ξ_n Hf′| and ), the following approximate solution is valid:

(27)

with . This implies that |f|∝|δ|^1/3 near the node. At the outer edge, the solution also has infinite slope, but the power law is different, since there the diffusive flux vanishes. Indeed, it is found that on the inside of (and near) the edge: |f|≈|2 H δ|^1/2.

Some of the properties of the diffusive flux at the null crossings are clearly reflected in the curves shown in the bottom panel of Fig. 2. For instance, the diffusive flux at those locations is different from zero, as expected. Secondly, the singularities at the nulls are apparent through the vertical slopes of the curves. Finally, the diffusive flux vanishes at the outer edge, as it must to fulfil the boundary condition.

4. The time-dependent solutions as tests of multidimensional MHD codes

4.1. Scope of the tests

The fundamental solutions of the self-similar problem presented in the previous section can be used as basic tests to check the correctness and capabilities of modules designed to solve the ambipolar diffusion problem in general magnetohydrodynamics codes. Given the presence of essential singularities within their domain, and the important physical role played by the latter in the ambipolar diffusion equation, these solutions can provide relevant, and demanding, tests that go beyond the case used in the past, namely the ZKBP solution, which has no internal nulls. To show the potential of these new tests, in the following, we present examples of application to a specific numerical code, namely the Bifrost code. First, the code is briefly introduced (Sect. 4.2). Then, a number of tests are carried out (Sect. 4.3). Finally, the asymptotic properties of the solutions for large times are discussed (Sect. 4.4).

4.2. The Bifrost code and the ambipolar diffusion module

The numerical experiments are performed in double precision with the 3D radiation-MHD Bifrost code (Gudiksen et al. 2011), whose core uses a staggered mesh in cartesian coordinates and a sixth-order differential operator. The code was used in a minimal configuration in which only the induction equation with ambipolar diffusion term is solved (i.e. no mass, motion, or energy equations are solved) in two dimensions (x, y), imposing cylindrical symmetry in the initial condition only. For the ambipolar diffusion term, we used the module developed by Nóbrega-Siverio et al. (2020b), which employs the super time-stepping (STS) technique (Alexiades et al. 1996). This module also includes an automatic selection of the free parameters of the STS method to obtain the best performance within a range that does not compromise the precision nor the STS stability. In addition, hyperdiffusion terms are included to guarantee the stability of the code in regions with strong gradients. The coefficients for the hyperdiffusion terms are chosen to be as small as possible, so that the stability is preserved while minimizing the diffusion impact in the numerical solution. The details of the automatic selection of the STS parameters and of the hyperdiffusion coefficients can be found in Sects. 5 and 6 of the paper by Nóbrega-Siverio et al. (2020b). In that paper, the validation of the module was achieved by reproducing the ZKBP solution in a situation with cylindrical symmetry. The decay in time of the maximum field and the spatial expansion of the outer front were found to follow an approximate power law, with the exponents matching with 10⁻³ accuracy those of the self-similar solution in Eq. (13). Thanks to the STS implementation, a speed-up factor of 8 compared with the simple implementation based on the Courant-Friedrichs-Lewy criterion for the ZKBP test could be achieved.

4.3. Tests for different harmonics

4.3.1. First harmonic

To test Bifrost’s STS module in the simplest configuration, we start by focusing on the first harmonic discussed in Sect. 3.2, which is shown as the blue curve of Fig. 2³.

For the test, we used the Bifrost code in two spatial dimensions (x, y), with the magnetic field pointing in the z-direction; the integration domain is a square of side 5.12 units. As an initial condition, we used the theoretical solution for that harmonic centred at the centre of the square; the outer front of the solution is located at a distance of unity from the centre, which implies L(t = 0) = 1, so as to fulfil the normalisation condition in a simple way. Beyond the outer front, the domain is padded with zeros. Three levels of spatial resolution are used, namely n_x = n_y = 512, n_x = n_y = 1024, and n_x = n_y = 2048, with grid sizes of 0.01, 0.005, and 0.0025, respectively. The calculation is carried out until the maximum of the solution is 10⁻² times the initial value; this occurs at t = 4.8 × 10³, which, using the values of H and K for the first harmonic in Table 1, corresponds to the advanced diffusion time τ_d = 5.3 × 10⁴. Before any comparison is made with the theoretical solutions of the previous section, we show in Fig. 3 the B(x, y, t) distribution at the initial (upper panel) and final (lower panel) times of the run for the case with the lowest resolution (n_x = n_y = 512). As shown in the figure, the symmetry of the initial condition is maintained to a high degree until the advanced time: a movie containing the time evolution would just show equal expansion in all directions. Simultaneously, from the field values of the grey-scale maps as given in the colour bars, we see that the magnetic field is decreasing by a large factor. The question, however, is whether the radial expansion and the field decrease accurately follow the power laws seen in the past sections.

Fig. 3. 2D maps B(x, y) showing the result of the Bifrost runs for the first harmonic at the initial time (t = 0; τd = 1; upper panel) and at the end of the run (t = 4.8 × 103; τd = 5.3 × 104; lower panel). In the figure, the results for the case with the lowest resolution, nx = ny = 512, are presented. The decrease of the field strength with time can be clearly appreciated through the colour bar scales.

For the comparison of the Bifrost solutions with the theoretical ones, we first determine, at each timestep, the location of the maximum field, [x_max(t),y_max(t)]⁴. Then, the radial coordinate is defined as

(28)

For the presentation of the results in the paper, a zonal averaging is carried out from the full 2D solution B(x, y, t) at any given time, namely, a 1D radial distribution is obtained by defining a collection of bins in the r-variable and averaging the value of B(x, y, t) for those (x, y) whose radial distance to (x_max, y_max) falls within each given bin. Concerning the theoretical solutions, the time-dependence (13)–(15) and the scalings and must be used with H and K coefficients given in the row corresponding to the first harmonic in Table 1 (second row after the header).

In Fig. 4, the initial function is shown in the left panel; the final state at diffusive time τ_d = 5.3 × 10⁴ for the lowest-resolution case n_x = n_y = 512 can be seen in the middle panel. The full time evolution for that case is presented in an accompanying movie (Animation 1); shown are the theoretical solution described in Sect. 3 (red) and the Bifrost solution (black) – or, more precisely, the solution averaged over radial bins of width equal to the grid size in the coordinate directions (i.e. Δr = 0.01 for this case). For all animations in the paper, we used a snapshot cadence corresponding to equal jumps in the logarithm of the diffusive time variable, τ_d. The (numerically determined) position of the outer front is marked with a dashed blue vertical line; that of the singular current sheet with a red dashed line. By construction, the theoretical and numerical curves lie on top of each other at time t = 0 (left panel). In the middle panel, we notice a small mutual deviation at the final time. To provide a quantitative measure for the deviation which is independent of the binning used to draw the figures, the absolute value of the relative deviation between the maxima of the theoretical and the Bifrost solutions for each of the resolution levels is shown in the right panel as a function of time. The relative deviation first grows and then stays at a roughly constant level in the advanced stages: this is due to the improvement of spatial resolution that naturally comes about through the expansion of the spatial support of the function as time progresses.

Fig. 4. Illustration of the result of the Bifrost runs for the first harmonic. Left panel: theoretical (red) and numerical (black) initial condition for the lowest-resolution case, nx = 512. The positions of the outer front and of the singular current sheet for the numerical solution are marked with dashed blue and red lines, respectively. Two horizontal lines are also superimposed at B(r) = 0 and at the maximum of the field of the numerical curve. Middle panel: configuration at diffusive time τd = 5.3 × 104 (the full time evolution is presented in the accompanying movie, Animation 1). Right panel: relative mutual deviation of the maxima of the theoretical and Bifrost solutions for the three resolution levels used in the test.

The general conclusion based on Fig. 4 and the accompanying animation is that the match is excellent. Clearly, the most sensitive location where one can check the mutual approximation between the theoretical and numerical solutions is the internal current sheet. Figure 5 shows the situation in the immediate neighbourhood of the sheet for the lowest- and highest-resolution cases (black for n_x = n_y = 512 and blue for n_x = n_y = 2048; the theoretical solution is drawn in red) and for the most advanced time of the experiment. The locations of the centres of the radial bins have been marked with symbols along the curves, for better comparison. The three solutions are expanding outwards, but we can see that the Bifrost solutions have a slight delay compared with the theoretical one. Also, the case with highest resolution does a very good job at matching the infinitely steep profile of the theoretical current sheet, whereas the n_x = n_y = 512 case has smoothed the profile to a certain extent.

Fig. 5. Solutions for the cylindrical problem in the immediate neighbourhood of the internal current sheet for the lowest- and highest-resolution cases (nx = ny = 512, black, and nx = ny = 2048, blue, respectively), in either case for the final time of the experiment. The red solid line is the theoretical curve.

The quality of the Bifrost solution can also be tested by checking how accurately it follows the theoretical power laws of decay and expansion in Eqs. (13)–(15). For the comparison we use the values of the exponents of those laws that correspond to the first harmonic, namely α = −0.422 and γ = 7.78 × 10⁻² (Table 1). The comparison is shown in Fig. 6 for the worst-resolution case n_x = n_y = 512 (black solid: Bifrost solution; red dashed: theoretical laws). The left panel is for the maximum B₀(t); more demanding tests are provided by the positions of the outer front, L(t), and of the internal current sheet, L_cs(t), both shown in the central panel. We note that the determination of the singular points for the numerical curves in the central panel was done using the binned solution explained in the previous paragraph. The match is excellent: a minimum-square fit to the Bifrost curves yields power-law exponents which match those of the theoretical law of Eqs. (13)–(15) with a precision of at least a few significant digits; in relative terms, the approximation is 2 × 10⁻³ for B₀(t) (left panel), 2 × 10⁻² (middle panel, L), and 1 × 10⁻² (middle panel, L_cs). The corresponding precision for the highest resolution case n_x = n_y = 2048 is 7 × 10⁻⁴ for B₀(t), 5 × 10⁻³ for L(t), and 4 × 10⁻³ for L_cs(t), respectively. Of the three quantities, the location of the outer front is the most difficult to determine, given the sharp corner between the solution (which falls to zero as | 2 H δ |^1/2, see Sect. 3.2), and the zeros beyond the front. In contrast, the location of the internal current sheet is easier to ascertain when one has the solution at hand, since it cuts the horizontal axis with a steep slope. As apparent in the figure and in the match of the power law exponents, both are well determined even in the worst-resolution case. We also note that the precision increases roughly linearly with the number of grid cells, which reflects the first-order accuracy of the STS algorithm.

Fig. 6. Time evolution of various significant variables in the calculation for the first harmonic illustrated in Fig. 4. Left panel: maximum of the field. Middle panel: location of the singularities, with the upper curve corresponding to L(t), and the lower one to Lcs(t). Right panel: flux integral in the central lobe of the solution. In each panel, the result of the Bifrost calculation is shown as a black solid line and the power laws derived from Eqs. (13)–(15) as red dashed lines.

Throughout the paper we have pointed out that the different harmonics that solve the eigenvalue problem associated with Eq. (16) provide a more stringent test for MHD numerical codes than the basic ZKPB solution: the reason is that the internal current sheets are non-avoidable singularities but with finite diffusive flux across them. This fact is essential for the time evolution of the solution; the numerical calculation of the diffusive flux is challenging precisely because of the infinite slope of the theoretical solution at those points. A further test of the capability of a code to deal with this difficulty is obtained, therefore, by calculating the integral of B(x, y, t) in any of the lobes of the solution (i.e. calculating an integral such as Eq. (19), but limited to a single cylindrical shell of constant magnetic field sign) and studying its time dependence. The integral of B(x, y, t) in the ring between the nulls or in the circle inside the internal current sheet must evolve in time exactly as B₀(t) L²(t), and so, using the analytical laws Eqs. (13)–(15), as a power-law of τ_d(t) with exponent . This implies (see Table 1) that the integral is conserved in time for the ZKBP solution () and decays increasingly rapidly for successively higher harmonics. For the first harmonic, in particular, the power law is . A comparison between numerical solution and analytical power law appears in the right panel of Fig. 6, again for the worst-resolution case (n_x = 512). The power-law exponent of the minimum-square fit to the numerical curve matches the analytical value with a accuracy of 8 × 10⁻³. For the highest-resolution case, n_x = n_y = 2048, the fit is four times more accurate, 2 × 10⁻³.

4.3.2. Higher harmonics

To complete the section on the tests for various harmonics, we here study one with a higher number of null crossings. As a relevant example, we consider the third harmonic (light-green curve in Fig. 2), which has a total of three internal null crossings in addition to the outer front. Checking with the theoretical solution, we expect (see Table 1) a faster decay of the maximum in this case (α = −0.474) and a much slower expansion of the spatial support (γ = 2.58 × 10⁻²). The result of the calculation for the lower-resolution case (n_x = n_y = 512) is presented in Fig. 7 (left panel: initial condition; middle panel: final snapshot in the series). The calculation was run until a similar time t as in the previous section; given the higher values of the sum K + H, this implies a more advanced diffusion time τ_d(t), namely τ_d = 1 × 10⁵, with the maximum of the solution decreasing to 4 × 10⁻³ of the initial value. We see that the numerical solution deviates from the theoretical one more markedly than in the corresponding case for the first harmonic (Fig. 4): here, the effective resolution for each of the lobes is roughly twice as bad as for the run for the first harmonic. This is also reflected in the corresponding curves for the relative deviation of the maxima of the solution between the Bifrost and theoretical solutions (right panel of the figure): for the highest-resolution case (n_x = n_y = 2048), the maximum relative deviation, while still quite good (on the order of 10⁻²), is a factor 2 greater in the experiment for the third harmonic than in that for the first. For the low-resolution case, the relative deviation is in the range of a few percent.

Fig. 7. Illustration of the result of the Bifrost runs for the third harmonic. Left and middle panels: initial condition and the final configuration in the calculation for nx = 512, respectively. The black thick curve is the numerical solution, whereas the red one is the theoretical solution. The position of the outer front and of the singular current sheet for the numerical solution are marked with dashed blue and red lines, respectively. The full time evolution is presented in the accompanying movie, Animation 2. Right panel: relative mutual deviation of the maxima of the theoretical and Bifrost solutions for the three resolution levels used in the test.

The determination of the location of the singular features and of the corresponding power-law exponents is a challenging test for the lower-resolution case (n_x = n_y = 512) of this harmonic or higher (Fig. 8). The location of the outer front L(t) is calculated by the numerical solution with a relative maximum deviation of 7 × 10⁻³ from the theoretical value; the fit of L(t) to a power law yields an exponent that matches the theoretical value, γ, within 5 × 10⁻². For the innermost current sheet, the situation is more demanding for the numerical solution: L_cs(t) is determined with accuracy better than 4 × 10⁻³ throughout the calculation; however, the fit of L_cs(t) to a power law yields an exponent that matches the theoretical value only within 4 × 10⁻². The fitted power laws for the time decay of the maximum of the solution and of the integrated magnetic flux within the innermost lobe have exponents that match the theoretical ones with an acceptable accuracy of 3 × 10⁻³ in the first case and 9 × 10⁻³ in the second. For the higher-resolution cases the Bifrost calculation does a much better job. For the highest resolution case (n_x = n_y = 2048), the determination of the power law exponents matches the theoretical values with accuracy 1 × 10⁻³ (decay of the maximum), 2 × 10⁻² (location of the outer front), 1 × 10⁻² (location of the innermost current sheet) and 3 × 10⁻³ (flux integral).

Fig. 8. Time evolution of various significant variables in the calculation for the third harmonic illustrated in Fig. 7. Left panel: maximum of the field. Middle panel: location of the singular points, where the upper curve is for L(t), and the lower one for Lcs(t). Right panel: value of the flux integral in the central lobe of the solution. In each panel, the result of the Bifrost calculation is shown as a black solid line; the power laws derived from Eqs. (13)–(15) are shown as red dashed lines.

4.4. Tests for more general initial conditions: the asymptotic evolution in time

In the absence of perturbations, the pure eigenmodes must evolve exactly fulfilling the self-similar shape and the power laws of Eqs. (13)–(15). However, given the errors associated with the numerical solution of the equation (both those associated with the discretisation of the initial condition and with the accumulated error along the update in time), the numerical solution will always deviate (if only by a small amount) from the pure theoretical eigenmode. However, in the previous section we show that the match of the numerical solution to the theoretical eigenmodes is excellent all the way to very advanced diffusion times, such as τ_d = O(5 × 10⁴); this is a qualitative indication both of the validity of the numerical solution and of the fact that the eigenmodes of the problem themselves, when subjected to small perturbations, remain in the close neighbourhood of the exact eigenmode up to advanced diffusion times.

A further question concerns the time evolution of the (theoretical and numerical) solutions when the initial condition deviates from a pure eigenmode by a finite perturbation: we consider whether the exact eigenmodes act as attractors for neighbouring functions – at least for those with the same number of zero-crossings of the solution. This question is of interest from the point of view of applied mathematics, of the physics of ambipolar diffusion, and of the construction of AD modules in numerical codes. In this section, we first obtain hints from the Bifrost solutions of the problem for the first (Sect. 4.4.1) and higher (Sect. 4.4.2) harmonics, which represent two different patterns of behaviour. Then we briefly discuss a few results from the literature (Sect. 4.4.3).

4.4.1. The first harmonic

Figure 9, which shows a selection of four snapshots of the accompanying Animation 3, presents (black curves) the results of a test run with Bifrost for the first harmonic in which the initial condition deviates by a significant amount from the pure eigenmode. For that example, the initial condition is given by:

(29)

Fig. 9. Four snapshots corresponding to the time evolution for initial condition with zero net flux, compact support, and a single internal null point, but otherwise of arbitrary shape. These snapshots are taken from the accompanying Animation 3, where the full time evolution can be seen. In both the figure and animation, the shape of the first harmonic corresponding to the unsigned flux of the solution toward the end of the evolution is superimposed as a green curve.

and zero outside of that interval. The parameters r₁ and W₁ were chosen as W₁ = 0.06 and r₁ = 0.48432, so that the function (top left panel) has zero net magnetic flux, only one internal null and essentially zero derivative at the origin; intermediate resolution (n_x = n_y = 1024) was chosen for this example. By letting the code run, we observe that the shape of the solution quickly approaches that of the first eigenmode: we added to the panels (green curves) a first harmonic calculated with the parameters (unsigned flux and maximum field) which hold toward the end of the evolution: we see that the initial solution quickly tends toward that shape; for instance, at the comparatively early diffusive time τ_d = 2.8 (lower left panel in the figure; see also the accompanying Animation 3) the solution is already quite near the first harmonic toward which it is advancing asymptotically. The run is carried out until a very advanced diffusive time, τ_d = 6.9 × 10⁴ (lower right panel); from about τ_d = O(10³) onward, the two curves are almost superimposed on each other, the largest discrepancy being found near the outer front.

We ran further cases of finite perturbation to the first harmonic. In all cases, the solution quickly approaches the shape of the exact first eigenfunction. This is highly suggestive that the first harmonic is not just stable against small perturbations (as shown in Sect. 4.3.1) but also in the sense that it is an attractor at least for initial functions that have zero net flux and a single zero-crossing in the middle of the domain.

Of special importance is the condition of zero net flux: when running those cases one has to make sure that the initial condition is in flux balance, at least up to a few significant digits; in other words, the integrated net flux must be orders of magnitude below the unsigned flux. From the theory (see Sect. 4.4.3), we know that an initial condition which is not in flux balance will evolve into the simple ZKBP solution corresponding to the same net flux. However, this is not going to modify the solution to any major extent until it has evolved so much that the unsigned flux of its lobes is comparable to the net signed flux of the initial condition. With the integrated flux evolving as B₀(t) L(t)², that is, as , the evolution will need to be calculated up to large values of τ_d for the ZKBP form to begin to become apparent in the solution; for instance, if the initial flux imbalance is 10⁻⁴ in relative terms, τ_d ought to be enormously large, possibly on the order of 10¹³, for the ZKBP function to begin to dominate the shape of the solution.

4.4.2. Higher harmonics

To test the behaviour of the higher harmonics when supplemented with a finite perturbation, we reran an experiment for the third harmonic, but this time adding a zero-flux perturbation such that the resulting function has the same number of zero-crossings along the radial direction as the original eigenmode, namely, three internal ones and the outer front. The evolution in this case is markedly different to all previous cases (see Fig. 10 and the accompanying Animation 4): the time evolution leads to a transformation of the third harmonic through merging of pairs of internal nullpoints⁵. In the figure and animation, the numerical solution is the black solid line. For the initial condition (top left panel in the figure) we additionally show its constituents, namely, a pure third eigenmode (red line) and the zero-flux, finite perturbation added to it, namely, 0.2 sin(2 π r²) (red dashed). As time proceeds, the shape of the solution increasingly deviates from that of a third harmonic: in the top-right panel, the two innermost nulls are about to merge, a process which is already complete at the time of the bottom-left panel. Finally, as time advances, the solution tends to the shape of the first harmonic with remarkable precision. For proper comparison, we have determined an exact first harmonic with the parameters (unsigned flux and maximum magnetic field) corresponding to the final state of the numerical solution and calculated its profile backward in time using the standard power-laws for that harmonic. That profile is superimposed (as a solid green curve) to all snapshots with τ_d > 10.

Fig. 10. Selected snapshots from the time evolution of the third harmonic when a non-small perturbation is added initially (see also Animation 4 for a high-cadence display of the evolution). Top left panel: initial condition (black solid), and its components, namely, an exact third harmonic (red solid) and a non-small perturbation (red dashed). The other three panels highlight the gradual conversion of the solution toward the shape of the first harmonic. The green curve corresponds to the time evolution of the first harmonic with unsigned flux and maximum field matching the solution in the final snapshot. The value of τd in the images is calculated with the H and K of the first harmonic (second row in Table 1 after the header). The corresponding values with the H and K for the third harmonic would be: τd = 1; 157; 619 and 1.26 × 105, respectively.

Finally, we also carried out a number of further experiments with different finite perturbations to the third harmonic, and also to the fifth harmonic. The results are similar in terms of null-point merging and transformation to the first harmonic after a finite time as in the experiment in this subsection.

The empirical results obtained here are compatible with the possibility that initial conditions consisting of a zero net-flux finite perturbation to harmonics higher than the first (and such that the superposition has compact support) end up with an exact first-harmonic shape. Also, judging by our results, it is not necessary to go to very large values of the diffusive time, τ_d, for this process to be well advanced, thus constituting a case of what is known in the mathematical literature as ‘intermediate asymptotics’. On the other hand, when using really poor spatial resolution to solve the problem of the evolution of the higher harmonics with no explicit finite perturbation, the numerical errors are likely to lead to a transformation of those harmonics into the first harmonic in a comparatively short integration time. Therefore, testing the higher harmonics with an MHD code adds a layer of difficulty to the test: whenever the initial discretisation or the numerical error made during the update in time is not small, the higher harmonics will evolve toward the first one and the test will fail.

4.4.3. The mathematical theory

A number of fundamental results have been rigorously proved in the mathematical literature concerning the asymptotic behaviour in time of some of the solutions of the porous medium equation and related equations (e.g. Kamin & Vázquez 1991; Bernis et al. 1993; Hulshof et al. 2001). What is of interest for us here is, primarily, the results that can be applied to the cylindrically symmetric case with diffusion coefficient which is proportional to the square of the dependent variable (n = 2, m = 3 in the notation of Eq. (7)). The most basic result, already mentioned in Sect. 4.4.1, is that initial conditions which have a finite nonzero flux integral (called ‘the mass’ in the mathematical literature for the PME) converge toward the ZKBP solution with the same flux integral (‘mass’) asymptotically in time (Vázquez 2007, Theorem 18.2); here, allowance is made for either a positive or negative flux integral by globally changing the sign of the ZKBP solution; also, ‘convergence’ is meant in the sense that the L^p norm of the difference between the actual solution and the ZKBP function tends to zero as t → ∞ faster than a negative power of the time with an exponent which is a function of n, m, and p (e.g. −1/3 for n = 2 and m = 3 in the L² norm; see details in the book by Vázquez 2007). A complementary result is the following: when the initial condition has positive net flux and its negative part has compact support, then the whole solution evolves into a positive function after a finite time (Vázquez 2007, Theorem 18.29). Since we are dealing with signed functions which have zero flux integral, these results are of interest mainly because they impose a strict condition on the possible flux imbalance caused by numerical errors (as discussed in Sect. 4.4.1, final paragraph): if it is not small, the numerical solutions will approach the ZKBP solution in a comparatively short time. However, the flux imbalance in all the Bifrost experiments discussed in the present paper is small enough that they have not shown this behaviour even though they have been run until a very long diffusive time.

Of special interest in the present paper would be mathematical results concerning the first harmonic in the eigenfunction series, meaning the harmonic with a single zero crossing within its spatial support. Bernis et al. (1993) deal with the plane-parallel 1D problem; the paper is devoted to the study of the so-called ‘dual porous medium equation’ (DPME), which is obtained through a double integration in space of the standard PME. These authors prove the existence of an asymptotic trend of the solutions of that equation: upper and lower bounds are given in the form of power-laws in time with a gamma exponent corresponding to the first zero-mass symmetric harmonic. The bounds bracket the global size of the solution (as measured by the supremum) indicating that it decreases in time in the same way as that harmonic. However, there does not seem to be any extension of those results to prove the actual convergence, for instance, in some L^p metric, of the difference between the solution of the PME equation and the first harmonic.

5. Summary and discussion

In the present paper, we first deal with the problem of ambipolar diffusion of an axial magnetic field, imposing the condition of cylindrical symmetry around the central axis of the domain, B = B(r) e_z. By solving an eigenvalue problem, we calculate explicit solutions that are self-similar, have compact support, and pass through a finite number of nulls. The resulting eigenfunctions have an intrinsic singularity at the nulls, so phase-plane techniques have been used to calculate the passage through them. The singularities constitute a finite collection of sharp current sheets through which magnetic flux is transported in spite of the vanishing magnetic field thanks to the infinite field gradient. The existence and various properties of this type of eigenfunction had been proved in the applied mathematics literature devoted to the porous medium equation (Hulshof 1991); to our knowledge, their explicit shapes have not been published yet. We can compare the form of the solutions in the neighbourhood of the singularities calculated in the present paper (Eq. (27)) with the properties of the solutions of the PME: both Grundy (1979) and Vázquez (2007, Chap. 16, and references therein) study in detail the nature of the critical points for the spatial part of the self-similar solutions of the PME in 1D, 2D, or 3D. These authors consider a phase-plane version of the equation; Grundy (1979), in particular, provides a list of all possible critical points in the problem and a table with a classification of the solution segments linking those critical points. We find that the individual lobes of our eigenfunctions can be matched to three different categories in his table (namely, to those of the third, fourth and fifth rows on page 276), corresponding, respectively, to (a) the link between our central point and the first internal null; (b) the link between the outer front and the nearest internal null; and (c) the link between two internal nulls (for the second or higher harmonics). It can be seen that those critical points possess singularities of the same order as those found in our text.

In the second part of the paper, we propose the set of self-similar solutions as tests for MHD numerical codes with ambipolar diffusion capabilities. To show their usefulness and validity, a battery of tests was carried out for the Bifrost code in two spatial dimensions starting from initial conditions with cylindrical symmetry (Sect. 4). We showed that the ambipolar diffusion module in Bifrost can cope with the passage of the solutions through the current sheets, with the level of accuracy increasing the higher the spatial resolution and in spite of the intrinsic singularity in them. Vice versa, the tests show that these functions can probe the capabilities of ambipolar diffusion modules to a larger extent than the simple ZKBP solution that has been used thus far (e.g. Masson et al. 2012; Viganò et al. 2019; Nóbrega-Siverio et al. 2020b). As test functions, the various harmonics proposed in our paper have the comparative advantage that they combine the B ∝ | δ |^1/3 singularity at the internal nulls (with δ the distance to the null) with the finiteness of the nonlinear diffusive flux, and this combination must be sufficiently well reproduced by the code if it is to pass the test. The ZKBP solution, instead, has a null just at the outer front, and the singularity there is of a lower order (∝|δ|^1/2), with zero diffusive flux across it. On the other hand, the scarcity of tests for the ambipolar diffusion term until now is in contrast with, for instance, the case of HD shocks, for which a whole category of exact solutions is available (the solutions of the Riemann problem) that have been used to develop sophisticated numerical schemes and tests (see Laney 1998; Toro 2009). The contrast to shocks, in fact, is interesting because of the differences in their mathematical and physical nature: in shocks, it is the (magneto)hydrodynamic evolution of the hyperbolic components of the PDE that leads to the formation of the singularities, which is then smoothed through simple diffusive phenomena (typically viscosity). In the ambipolar diffusion problem, wherever there is a null, it is the diffusive phenomenon itself that creates and maintains the singularity.

The shape of the solutions near the nulls will be modified when the problem also includes standard Ohmic diffusion (or any other linear diffusion term for the magnetic field or artificial diffusion designed for the smooth operation of the code). Let us consider a problem that includes Ohmic diffusivity in addition to the ambipolar term, but in which the latter is dominant in a given region. Under such circumstances, the approach to the singularity is expected to have the standard nonlinear diffusion profile f ∝ | δ |^1/3, except in the immediate neighbourhood of the null where even a small amount of linear (e.g. Ohmic) diffusivity will prevent an infinite slope. Parker (1963) considered that problem in the case of a stationary situation (see also Cheung & Cameron 2012). If the slope of the field is not infinite at the null, then the diffusive flux −χ_a B² B′ will be zero there. However, the total diffusive flux (Ohmic plus ambipolar), −B′/(μ₀ σ)−χ_a B² B′, will be nearly uniform around the null and fixed by the conditions imposed by the dominant agent outside of the null, that is, the ambipolar diffusion. For the time-dependent, self-similar solutions of the present paper, we expect a similar situation: the slope of the function will be finite in the neighbourhood of the null (as forced by the simple Ohmic diffusion), but the diffusive flux at the null will be as given by the pure ambipolar problem. The time evolution of the global solution is therefore expected to be quite close to that of the pure ambipolar problem. A behaviour of the numerical solutions around the null as just described also fits the results of the numerical tests in Sect. 4, even though the hyperdiffusion terms used by the Bifrost code are not of the simple Ohmic diffusion kind (see Nóbrega-Siverio et al. 2020b). We have seen that for a sufficient spatial resolution, the accuracy of the exponent in the power laws is high and, more to the point, the integral of the diffusing variable (the magnetic flux) in the central lobe follows the analytical law of time evolution quite precisely (Fig. 6). That time dependence is physically and mathematically determined by the rate of exchange of magnetic flux across the nulls, which is given by the diffusive flux. If the presence of hyperdiffusivities (or Ohmic diffusivities) changed the diffusive flux by any important amount, the power-law behaviour would not be fulfilled.

The final point in this discussion concerns the asymptotic behavior in time (i.e. for large values of τ_d) of the eigenfunctions presented in this paper (Sect. 4.4). By using the Bifrost code, we checked that a small enough perturbation (e.g. due to the discretisation error of the initial condition and the accumulated error through the update in time) does not bring the eigenfunctions out of their initial shape at least until a very advanced diffusive time, τ_d, except for a small amount that depends on the spatial resolution used for the numerical calculation. On the other hand, through experiments in which the initial condition, while still of zero net flux and with compact support, deviates from the exact eigenfunction shape by a non-small amount, we have seen that harmonics higher than the first one decay into the latter through merging of the internal nulls (or, possibly, via the ejection of the outermost ones through the outer front). Instead, if that kind of finite perturbation was given to the first harmonic, the solution ends up adopting an exact first-harmonic shape. This can be of interest from a twofold perspective: on the one hand, from the mathematical point of view, our empirical results (obtained numerically) suggest the possibility that the first harmonic is an attractor for initial conditions with cylindrical symmetry which have zero net magnetic flux. On the other hand, from the point of view of testing MHD codes, the higher harmonics pose a double challenge when given as initial conditions with no finite perturbation to them: the code must be able to maintain the shape of the harmonic until a sufficiently long diffusive time is reached.

To conclude, as a general recipe, when testing an MHD code with ambipolar diffusion capabilities, the following steps are recommended

1. First to test for the ZKBP solution, which has been the standard practice in the past;

2. Then to check for the first harmonic and see that it keeps its theoretical shape for large diffusive times and following the standard power laws with sufficient precision.

3. And, finally, to carry out tests with the higher harmonics to check (a) whether small numerical perturbations do not bring them out of their initial shape and (b) whether explicit, finite initial perturbations lead them to transform ino the first harmonic.

In this paper, we have presented a new family of exact self-similar solutions to the problem of ambipolar diffusion in a cylindrical geometry and tested them with the Bifrost code. The extension of these results to problems with different geometries will be provided in a future paper.

Movies

Movie 1 associated with Fig. 4 (animation_1) Access here

Movie 2 associated with Fig. 7 (animation_2_512) Access here

Movie 3 associated with Fig. 9 (animation_3_1st_harm_arbitr_2) Access here

Movie 4 associated with Fig. 10 (animation_4_3rd_harm_arbitr) Access here

Acknowledgments

The authors are grateful to an anonymous referee for the careful reading of the manuscript and useful suggestions. F.M.I. thanks Prof. J. Hulshof for his guidance concerning the mathematical literature on the Porous Medium Equation and for comments to the manuscript; bibliography suggestions by Prof. J.L. Vázquez are also appreciated. This research has been supported by the Spanish Ministry of Science, Innovation and Universities through projects AYA2014-55078-P and PGC2018-095832-B-I00. The authors are also grateful to the European Research Council for support through the Synergy Grant number 810218 (ERC-2018-SyG). D.N.S. acknowledges support by the Research Council of Norway through its Centres of Excellence scheme, project number 262622, and through grants of computing time from the Programme for Supercomputing. A.W.H. gratefully acknowledges the financial support of STFC through the Consolidated grant, ST/S000402/1, to the University of St Andrews.

References

All Tables

All Figures

	Fig. 1. Selected solutions of the fundamental equation in an infinite domain subject to the inner boundary conditions of Eq. (17) with asymptotic shape toward infinity that illustrate the identification of the eigenvalues and eigenfunctions for the fundamental equation subject to the full boundary conditions of Eqs. (17), (18), and (22).
In the text

	Fig. 2. First five eigenfunctions as a function of the stretched coordinate . The actual solutions are shown in the upper panel, while the corresponding dimensionless diffusive flux D is displayed in the lower panel. The locations of the null crossings of the solutions are indicated with vertical dashed lines. In both panels, the eigenvalues are indicated.
In the text

Fig. 3. 2D maps B(x, y) showing the result of the Bifrost runs for the first harmonic at the initial time (t = 0; τd = 1; upper panel) and at the end of the run (t = 4.8 × 103; τd = 5.3 × 104; lower panel). In the figure, the results for the case with the lowest resolution, nx = ny = 512, are presented. The decrease of the field strength with time can be clearly appreciated through the colour bar scales.

In the text

Fig. 4. Illustration of the result of the Bifrost runs for the first harmonic. Left panel: theoretical (red) and numerical (black) initial condition for the lowest-resolution case, nx = 512. The positions of the outer front and of the singular current sheet for the numerical solution are marked with dashed blue and red lines, respectively. Two horizontal lines are also superimposed at B(r) = 0 and at the maximum of the field of the numerical curve. Middle panel: configuration at diffusive time τd = 5.3 × 104 (the full time evolution is presented in the accompanying movie, Animation 1). Right panel: relative mutual deviation of the maxima of the theoretical and Bifrost solutions for the three resolution levels used in the test.

In the text

	Fig. 5. Solutions for the cylindrical problem in the immediate neighbourhood of the internal current sheet for the lowest- and highest-resolution cases (nx = ny = 512, black, and nx = ny = 2048, blue, respectively), in either case for the final time of the experiment. The red solid line is the theoretical curve.
In the text

Fig. 6. Time evolution of various significant variables in the calculation for the first harmonic illustrated in Fig. 4. Left panel: maximum of the field. Middle panel: location of the singularities, with the upper curve corresponding to L(t), and the lower one to Lcs(t). Right panel: flux integral in the central lobe of the solution. In each panel, the result of the Bifrost calculation is shown as a black solid line and the power laws derived from Eqs. (13)–(15) as red dashed lines.

In the text

Fig. 7. Illustration of the result of the Bifrost runs for the third harmonic. Left and middle panels: initial condition and the final configuration in the calculation for nx = 512, respectively. The black thick curve is the numerical solution, whereas the red one is the theoretical solution. The position of the outer front and of the singular current sheet for the numerical solution are marked with dashed blue and red lines, respectively. The full time evolution is presented in the accompanying movie, Animation 2. Right panel: relative mutual deviation of the maxima of the theoretical and Bifrost solutions for the three resolution levels used in the test.

In the text

Fig. 8. Time evolution of various significant variables in the calculation for the third harmonic illustrated in Fig. 7. Left panel: maximum of the field. Middle panel: location of the singular points, where the upper curve is for L(t), and the lower one for Lcs(t). Right panel: value of the flux integral in the central lobe of the solution. In each panel, the result of the Bifrost calculation is shown as a black solid line; the power laws derived from Eqs. (13)–(15) are shown as red dashed lines.

In the text

Fig. 9. Four snapshots corresponding to the time evolution for initial condition with zero net flux, compact support, and a single internal null point, but otherwise of arbitrary shape. These snapshots are taken from the accompanying Animation 3, where the full time evolution can be seen. In both the figure and animation, the shape of the first harmonic corresponding to the unsigned flux of the solution toward the end of the evolution is superimposed as a green curve.

In the text

Fig. 10. Selected snapshots from the time evolution of the third harmonic when a non-small perturbation is added initially (see also Animation 4 for a high-cadence display of the evolution). Top left panel: initial condition (black solid), and its components, namely, an exact third harmonic (red solid) and a non-small perturbation (red dashed). The other three panels highlight the gradual conversion of the solution toward the shape of the first harmonic. The green curve corresponds to the time evolution of the first harmonic with unsigned flux and maximum field matching the solution in the final snapshot. The value of τd in the images is calculated with the H and K of the first harmonic (second row in Table 1 after the header). The corresponding values with the H and K for the third harmonic would be: τd = 1; 157; 619 and 1.26 × 105, respectively.

In the text