Boundary layers and noise in magnetized plasmas line-tied at conducting surfaces

TL;DR

This paper investigates boundary layers and noise phenomena in magnetized plasmas with line-tied conducting boundaries, revealing that boundary layers are real and must be resolved for accurate simulations, especially as they scale with plasma resistivity.

Contribution

The study combines analytical and numerical methods to demonstrate the existence of boundary layers and noise in line-tied magnetized plasmas, emphasizing their importance for simulation accuracy.

Findings

Boundary layers exist near conducting boundaries in magnetized plasmas.

Short wavelength structures are present on long wavelength features.

Boundary layer widths scale with the square root of plasma resistivity.

Abstract

In magnetized plasma situations where magnetic fields intersect massive conducting boundaries, "line-tied" boundary conditions are often used, analytically and in numerical simulations. For ideal MHD plasmas, these conditions are arrived at given the relatively long time scales for magnetic fields penetrating resistively into good conductors. Under line-tied boundary conditions, numerical simulations often exhibit what could be construed as numerical "noise" emanating from the boundaries. We show here that this "noise" is real. By combining numerical and analytical methods, we highlight the existence of boundary layers near the conductors and confirm the appearance of short wavelength structures riding on long wavelength features. We conclude that for numerical fidelity the boundary layers need to be resolved. Boundary layer widths scale as the square root of the plasma \b{eta}.

Figures

Equations

$${\partial }_{t}n+\nabla \cdot n\mathbf{u}=0$$

$$nM{d}_t\mathbf{u}=-T\nabla n+j_y\hat{\mathbf{y}}\times \mathbf{B}$$

$$j_y=\hat{\mathbf{y}}\cdot \nabla \times \mathbf{B}$$

$${\partial }_t\mathbf{B}=\nabla \times (\mathbf{u}\times \mathbf{B})$$

$$\nabla \cdot \mathbf{B}=0$$

$${\partial }_{t}n+{\partial }_x{u}_x=-ik{u}_z$$

$${\partial }_t{u}_x=-T{\partial }_{x}n$$

$${\partial }_t{u}_z=-ikTn+(i/k){\nabla }^2{B}_x$$

$${\partial }_t{B}_x=-ik{u}_z,$$

$$({\gamma }^2-T{\partial }_x^2)n={\gamma }^2{B}_x$$

$$({\gamma }^2-{\nabla }^2)B_x=-k^{2}Tn$$

$${\nabla }^2{B}_x=k^{2}Tn.$$

$$({\gamma }^2-T{\partial }_x^2)n=0.$$

$$n=\sin (Kx),{\gamma }^2=-K^{2}T.$$

$$B_x=-T[\sin (Kx)-\sinh (kx)/\sinh (k)].$$

$${\gamma }^{2}n={\gamma }^2{B}_x.$$

$$({\gamma }^2-{\nabla }^2)B_x=0.$$

$$B_x=\cos (Kx),{\gamma }^2=-[K^2+k^2].$$

$$({\gamma }^2-T{\partial }_x^2)n={\gamma }^2{B}_x.$$

$$n(x)=\cos (Kx)-A\cos (\kappa x),\kappa \gg K.$$

$$n(x)=\cos (Kx)-[(K/\kappa )/\sin (\kappa )]\cos (\kappa x),\kappa \gg K.$$

$$n(x\to -1)=K[s-(1/\kappa )\sin (\kappa s)-(1/\kappa )(\cos \kappa /\sin \kappa )].$$

$$f=Ax+B{x}^2.$$

$$f_{-1}=-3f_1+f_2$$

$$f_{-2}=-8f_1+3f_2.$$

$$f=f_0+A{x}^2+B{x}^3.$$

$$f_{-1}=-(3/2)f_0+3f_1-(1/2)f_2$$

$$f_{-2}=-12f_0+16f_1-3f_2.$$

Full text

1 Introduction

There are many situations where strong magnetic fields in plasma intersect massive, inert, conducting boundaries. One example is emerging magnetic flux from the solar surface; here the solar surface and below is massive as well as conducting (Parker, 1972; Hood & Priest, 1979; Antiochos, 1987; Lionello et al., 1998; Priest, 2014; Huang et al., 2006). As another example are laboratory magnetized plasma experiments where magnetic lines meet vacuum vessel walls. In both these situations, the plasma MHD timescales are much shorter than resistive penetration times into the conductor. Thus, the field lines are effectively line-tied.

In numerical solutions of magnetized plasma for such situations, so-called “line-tied” boundary conditions are implemented by assuming that both the tangential electric field and the perturbed normal magnetic field be zero at the boundary. For MHD, since the flows are dominantly $E\times B$, this translates to zero tangential flows (where strong fields cut into the boundary). In addition, if the boundary is assumed effectively impervious, the normal flow is also set to zero, and the normal derivative of the density is zero (“hard-wall” boundaries). Finally, the normal derivative of the tangential perturbed magnetic field is assumed to be zero. [(Hood, 1986) has considered line-tied boundary conditions in detail. He distinguishes between impervious boundaries and boundaries that allow normal flow through. In this paper, we assume Hood’s impervious choice.]

As an example, consider a 1D shear Alfvén wave set up between two conducting plates at $x=-1$ and $x=+1$, as shown in Fig. 1. There is an equilibrium uniform magnetic field, $\boldsymbol{B}_{0}=B_{0x}\boldsymbol{\hat{x}}$.

An initial condition on the $z$-directed plasma flow, $u_{z}=\sin(Kx)$, $K=\pi/2$, is given. In a code we describe later (Guzdar et al., 1993), which uses a fourth order finite difference scheme, two “ghost points” outside of the domain need to be specified. The simplest prescription, consistent with Line-tied/Hard-wall (LTHW) boundary conditions, is to symmetrize $B_{z}$ and antisymmetrize $u_{z}$ across the boundaries. Further, $u_{z}(x=\pm 1)$ are set to zero, while $B_{z}(x=\pm 1)$ are stepped. Corresponding boundary conditions for other possibly involved variables would be antisymmetric $u_{x}$ and $B_{x}$, and symmetric density, $n$. The 1D Alfvén wave, in a numerical solution as initialized above, works very well, and, in particular, is noise-free.

The situation, however, is more complex in the case of a 2D or 3D situation. A system set up such as in Fig. 1, if subjected to 2D or 3D perturbations, can be shown to generate small scale structures. While the structures may remain small amplitude, they sometimes devolve into apparent numerical instability. The latter consequence is especially true for nonlinear problems. It is well recognized by workers in the field that this type of "noise" is often observed. A remedy sometimes suggested is to place a strong but narrow diffusion layer very close to the boundary. While this is apparently sometimes workable, there is a risk of slippage in line-tying (to the extent that line-tying is an important part of the matter at hand).

In this paper, we explore a system, such as in Fig. 1, both analytically and numerically, to understand the origin of the "noise". In particular, we study the linear ideal MHD normal modes of the system of Fig. 1, with LTHW boundary conditions. From analytic solutions, we find two characteristics of the normal modes:

(1) the eigenfunctions show the presence of sharp boundary layers near the conducting walls;

(2) the eigenfunctions show a two-space scale structure, specifically a high wave number "ripple" riding on low-wavenumber envelope solutions.

This kind of multiscale structure could easily be misconstrued as "noise". In fact, if the boundary layers are not sufficiently numerically resolved, this could result in real numerical noise. Our conclusion from these findings is that high resolution numerical analysis may be generally necessary for these types of problems.

In the next four sections, we analytically calculate the eigenmodes of the system. This calculation, we will show, requires boundary layer and multiscale methods. In Sec. 6, we run several numerical simulations to bolster and confirm the analytic solutions. We compare different numerical solution methods to quantify the grid resolution necessary to address line-tied problems.

2 Equations

As mentioned, we investigate the system shown in Fig. 1. Ideal MHD equations in 2D are considered. In usual notation, the governing equations are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the temperature, $T$, is assumed isothermal. The vector $\boldsymbol{B}$ is assumed to be in the $x-z$ plane. (For the purposes of this paper, the system in Fig. 1 will be normalized as follows: we set the distance between the plates to 2, where the domain of $x$ is [-1,1]. Also, we set the Alfvén speed to unity. In particular, $B_{0x}=1$, and $n_{0}M=1$. Thus energies are normalized to the magnetic energy, and, in normalized units, $T$ is the plasma $\beta$.) While these equations are nonlinear, and what we present in this paper pertains also to nonlinear conditions, it suffices to demonstrate our findings in the linear regime. Accordingly, from Eqns (1)–(5), we obtain a linearized system in the variables $n$, $u_{x}$, $u_{z}$, and $B_{x}$, given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\partial_{x}B_{x}=-ikB_{z}$. The Laplacian $\nabla^{2}$ is defined as $\partial_{x}^{2}-k^{2}$. Here, we have assumed solutions of the form $\exp(ikz)$, where $k=2\pi/L$. The quantities $n$, $u_{x}$, $u_{y}$, $u_{z}$, $B_{x}$, $B_{y}$, $B_{z}$, are all small perturbations, compared to $n_{0}$ and $B_{0x}$. These equations are fourth order in $d/dt$ and thus yield four ideal MHD eigenmodes, two fast modes and two slow modes. Generally speaking, with the LTHW boundary conditions as above, this system will exhibit the type of noise discussed above.

3 Mismatch

We begin by noting that if $k$ is set to zero in Eqs (6)–(9), the system is 1-D. In that case, the modes clearly separate into the sound wave and the Alfven wave. In particular, the first two equations can yield a standing sound wave, parallel to the B field, coupling $n$ and the parallel flow $u_{x}$. The boundary conditions are "purely" hard-wall: in particular, a sinusoidal, i.e., $\sin(Kx)$, density perturbation, (given the symmetric boundary conditions at $x=-1$ and $x=+1$) “forces” a co-sinusoidal behavior in $u_{x}$, i.e., $\cos(Kx)$. The latter is consistent with the hard-wall antisymmetric boundary conditions for $u_{x}$. Similarly, the last two equations set up an Alfven wave coupling $uz$ to $B_{x}$. Both these variables are co-sinusoidal, ie $\cos(Kx)$, consistent with the antisymmetric line-tied boundary conditions for each.

If, however, $k$ is non-zero, the eigenfunctions are more intricate. We will show this in the next two sections. Here, we first point out that there is a “mismatch”, between slow and fast modes, in their spatial sinusoidal parities. Consider the slow wave. This is a sound like mode, propagating parallel to the B field. The sinusoidal parities were described above. However, if one examines Eq (6), within the assumption of the parities of $n$ and $u_{x}$ just mentioned, the parity that would “fit” the variable $u_{z}$, from Eq (6), would be sinusoidal, ie, $\sin(Kx)$. But this parity violates the assumption of antisymmetric boundary conditions for $u_{z}$. This violation is what we refer to as a “mismatch” between the slow and fast modes. We will show that this type of mismatch gives rise to a boundary layer in $u_{z}$, since $u_{z}=0$ is enforced at the plates.

A similar mismatch also arises in the fast mode for $k$ nonzero. The fast mode is generally an Alfvénic mode, involving largely $B_{x}$ (or $B_{z}$) and $u_{z}$. As already discussed, this has $u_{z}\sim\cos(Kx)$, being antisymmetric at the plates; while $B_{x}\sim\cos(Kx)$, consistent with $B_{x}$ also being antisymmetric at the plates. However, from Eq (8), the parities as deduced for $u_{z}$ and $B_{x}$ would force $n$ to be $\cos(Kx)$, which is antisymmetric at the plates, and so violates the symmetric boundary condition requirement for $n$ at the plates.

We will now show below that these mismatches play a role in boundary layer formation in these modes, leading to short scale structure.

4 Analytical Solution — slow mode

We proceed to obtain analytically the slow and fast eigenmodes for the system in Fig. 1. We begin by rewriting the four Eqs (6)–(9) as two coupled equations. We use (7) in (6) and (9) in (8) to get the following 2-field equations for $n$ and $B_{x}$

[TABLE]

[TABLE]

where solutions of the form $\exp(\gamma t)\exp(ikz)$ are assumed, and $\gamma$ is the eigenvalue.

We now make the low-beta assumption $T\ll 1$. We obtain a reduced description for the slow mode in the $T\ll 1$ regime. Such reduction has the advantage of practically decoupling the slow and fast modes, allowing for easier equations that need to be solved. (Our conclusions in this paper, we will show, are independent of this reduction.) Note, the limit $T\ll 1$ ensures that the sound speed is subdominant to the Alfvén speed.

We begin with the slow mode, wherein we anticipate that $\gamma^{2}\ll 1$ (the slow mode is sub-Alfvénic). In this case, from Eq (11), we can discard the $\gamma^{2}$ term; but we order the Laplacian to be $O(1)$, i.e., long wavelengths. Since the RHS is proportional to $T\ll 1$, we conclude that $B_{x}$ must be small to allow an optimal balance in that equation. Thus, for $T\ll 1$, we have from Eq (11)

[TABLE]

From Eq (10), comparing the two $\gamma^{2}$ terms, we see that the term on the RHS must be neglected since $B_{x}$ is small, making the RHS quadratically small. The remaining two terms on the LHS are of the same order, and so Eq (10) reduces as

[TABLE]

Eq (13) is now an eigenvalue equation for $n(x)$ with eigenvalue $\gamma^{2}$. The $n$ eigenfunction and the eigenvalue, invoking the appropriate hard-wall boundary condition, are

[TABLE]

This is a long wavelength sound wave.

Turning to Eq (12), this is now an inhomogeneous equation for $B_{x}$. The solution, constructed from a homogeneous solution and an added particular solution, and satisfying the appropriate line-tied boundary condition, is

[TABLE]

Note in the solution the appearance of hyperbolic behaviour, which brings into doubt the use of simple symmetric/antisymmetric boundary conditions. In fact, if we assume $k\gg K$, the appearance of a boundary layer is clearly evident, as shown in Fig. 2. The Figure shows $B_{x}(x)$ plotted from Eq. (15). $k/K$ is chosen to be 20 to highlight the boundary layer at the $x=\pm 1$ plates. Note that for $k/K\gg 1$, $B_{x}$ approaches the boundaries at $x=\pm 1$ as symmetric but then there is a boundary layer near the plate enforced by the zero boundary condition at $x=\pm 1$. As is evident, the "mismatch" discussed earlier results in the boundary layer.

We note here that we have not used any boundary layer analysis to arrive at the solution for $B_{x}$ (Eq. 15). However, we clearly see a boundary layer in Fig. 2, provided $k/K\gg 1$. This is because for our calculation we assumed an "optimal ordering", $k\sim K$, and we were fortunate enough to be able to solve the equation. Alternatively, if we had assumed $k/K$ large, right from the start, we would have found the $\sin(Kx)$ solution but not the $\sinh(kx)$ part. We would have then had a "mismatched" boundary condition. This would have precipitated a boundary layer analysis, leading to an "inner solution" matching asymptotically with the "outer solution" (Bender & Orszag, 1999). We have indeed verified this behavior. We do not show this calculation here, except to mention it, as we already have a uniformly valid solution everywhere. In the next section (the fast mode) we will indeed be forced to resort to boundary layer methods.

It is worth noting here that the emergence of a boundary layer in the $k/K$ large case has an important consequence. In particular, the "reduced equations" of MHD (Strauss, 1976) are derived under the assumption of long parallel wavelengths compared with perpendicular wavelengths, generally stated as $k_{\parallel}\ll k_{\perp}$. But this is exactly our condition above. Thus, the application of the reduced equations at conducting surfaces should require boundary layer analysis near the conductors. This is a situation indeed encountered. The appearance and necessity of a boundary layer for reduced equations has indeed been reported previously (Scheper & Hassam, 1999). In the latter paper, it was shown that low frequency motions near the solar conducting surface would require a boundary layer analysis (thus necessitating the use of non-reduced equations in the layer.)

As a final comment on the slow mode, we note that an equilibrium solution ($\partial_{t}=0$), with density having $k$ and $K$ variations, can be readily found. This too has boundary layers at the plates.

5 Analytical Solution — fast mode

We now turn to the fast mode. We begin as in Sec 4 with Eqs (10) and (11). This time, we order $\gamma^{2}$ as $O(1)$, as we anticipate an Alfvén mode. Since $T\ll 1$, and we assume $\partial_{x}^{2}$ is $O(1)$ for long wavelengths, Eq (10) reduces to

[TABLE]

This yields $n=B_{x}$. Using the latter in Eq (11), the RHS of (11) must be discarded, being small, leading to

[TABLE]

The latter Eq is an eigenvalue equation for $B_{x}$. With line tied boundaries, the eigenfunction and eigenvalue are

[TABLE]

From Eq (16), as above, $n$ has the same parity as $B_{x}$. However, this parity violates the zero derivative boundary condition required for $n$. (This is a manifestation of the "mismatch" discussed earlier). This violation indicates that solutions near the boundaries must be reconsidered: that is, there may be boundary layers as $x\rightarrow\pm 1$. In what follows, we show that indeed there is a boundary layer. A boundary layer in solution generally indicates rapid variation in the x-direction, i.e., $d/dx\gg 1$. As a consequence, there may be sharp variations in the eigenfunction, or, alternatively, a multiple scale character to the solution (Bender & Orszag, 1999). For this paper, we have investigated both possibilities and we find that the latter character, multiscale behaviour, is the apt description. We proceed below with a multiscale analysis.

We go back to Eq (10). We note that $\gamma^{2}$ is already known to be $O(1)$. But this time we order $\partial/\partial x$ large. In that case, the term $T\partial_{x}^{2}$ must be taken to be optimally $O(1)$, and, so, retained in the equation. Thus, the multiple scale equation from Eq (10) remains

[TABLE]

Assuming even $n(x)$, the general solution to this inhomogeneous equation can be written

[TABLE]

(Note that the particular solution is approximate: a term of $O(K^{2}/\kappa^{2})$ has been discarded.) The constant $A$ is fixed upon demanding zero derivative of $n$ at $x=-1$. Thus,

[TABLE]

This solution is plotted in Fig. 3, showing clearly the short wavelength ripple riding on the long wavelength solution for $n$. Note also the zero derivative at the boundaries. Note further that the ripple is of short wavelength, by a factor $\kappa/K$, but also that the amplitude of the ripple is small, by the same factor. It can also be noted that the $A$cos$(\kappa x)$ part of the solution is sub-dominant in the region away from $x\rightarrow\pm 1$, but that it becomes comparable to the particular solution, cos$(Kx)$, in the boundary layers. This gives a narrow layer width, narrower by the same factor $(\kappa/K)$ as above. The structure of the boundary layer can be seen upon expanding (21) near the boundaries. At the left boundary, we let $x=-1+s$, where $s\ll 1$ and $\kappa s\sim 1$. Equation (21) then becomes

[TABLE]

We conclude that the ripple is part of the eigenmode, and the “noise” must be concluded to be real. In particular, in a numerical solution, it must be resolved.

6 Numerical solutions

As a more in-depth examination of the noise, we conducted a series of numerical experiments with high resolution, using various boundary conditions in a finite-difference code, and increasing the number of grid points. We also compared the finite-difference code with a much more accurate solution using a Chebyshev spectral method, with which the discretization error decreases exponentially when the number of collocation points $N$ increases (Trefethen, 2000). We subjected the linear equations, (6)-(9), to an initial condition given by $n=((1+\cos(\pi x))/2)^{4}\exp(ikz)$. The temperature $T=0.2$ and the wavenumber $k=15$. Our primary objective was to quantify the numerical errors due to different implementations of spatial disretization; hence, we kept the timestep $\Delta t=5\times 10^{-5}$ for all calculations, with a trapezoidal leapfrog scheme for time stepping. Here we first report the solutions from a simulation using the Chebyshev code with $N=200$. By varying $N$ for a convergence test, we confirmed that the discretization error was reduced to approximately the level of floating point round-off error at this resolution. We show in Fig. 4 traces of $n$ and $u_{x}$ at $t=0$, $5$, $9$, and $12$. We note the early time traces are fairly smooth, and of long wavelength, but the late time solutions can be construed as “noisy”. Of particular note is the highly “noisy” trace of $u_{x}$, and the short-scale “ripple” riding on top of the long wavelength envelope of the density profile.

Now we turn to the finite-difference implementations of the linear problem. Here we use a fourth order finite-difference scheme that requires two "ghost" points outside the boundary. In one set of calculations, we implemented the symmetric ($n$)/anti-symmetric ($u_{x}$, $u_{z}$, $B_{x}$) boundary conditions for the ghost points. However, we recall, our analysis of Secs. 4 and 5 above had revealed boundary layers, thus rapidly varying functions, bereft of any simple symmetries across the boundaries. To address this behaviour, we tried also, in another set of calculations, to numerically fit the two ghost points to a solution using polynomials. For example, for the boundary condition $f(b)=0$, we used

[TABLE]

This gave the following conditions for the ghost points, viz

[TABLE]

[TABLE]

Likewise, for $(df/dx)(b)=0$ boundary conditions, we fitted to

[TABLE]

This gave the following conditions for the ghost points, viz

[TABLE]

[TABLE]

We then re-ran our test example with the two sets of boundary conditions. Fig. 5 shows the solutions at $t=9$ with various numbers of grid points $N$. Here, the black dashed line in each panel shows the Chebyshev solution with $N=200$, the blue line shows the solution with symmetric/anti-symmetric boundary conditions, and the orange line shows the solution with polynomial-fit boundary conditions. Panel (a) shows the cases with $N=50$. Here, the polynomial-fit solution exhibits high-amplitude fluctuations near the boundaries, because polynomial extrapolation can be very unreliable when the solution is under-resolved. Panel (b) shows the cases with $N=75$. Now the two sets of boundary conditions yield very similar solutions, and the solutions appear to be "noisy". Note that the fluctuations in the finite-difference solutions are out of phase with respect to the Chebyshev solution, due to the finite-difference discretization error. Panel (c) shows the cases with $N=150$, where the finite-difference solutions now approach the Chebyshev solution. While the two finite-difference solutions are rather similar, the close-up view near the boundary at $x=-1$, shown in panel (d), reveals a notable difference. Here we note that the solution with symmetric/anti-symmetric boundary conditions exhibits low-amplitude fluctuations at the grid scale, as a consequence of the "mismatch" discussed in Sec. 3. In contrast, the solution with polynomial-fit boundary conditions is noticeably smoother.

To quantify the error of finite-difference solutions, we measured the root-mean-square error of the density profile $n$ relative to the Chebyshev solution at $N=200$, which was taken as a proxy for the accurate solution. Fig. 6 shows the scaling of errors as a function of $N$, for both sets of finite-difference solutions, as well as the errors of Chebyshev solutions. As a consequence of the exponential convergence of spectral methods, the error of Chebyshev solution decreases rapidly as $N$ increases. On the other hand, we expect the error of a fourth-order finite-difference method to decrease as $N^{-4}$. The $N^{-4}$ scaling is observed for both sets of boundary conditions at low $N$. However, the "mismatch" problem degrades the convergence rate to approximately $N^{-2}$ at high $N$ for the symmetric/anti-symmetric boundary conditions. In comparison, the $N^{-4}$ scaling is retained in the implementation of polynomial-fit boundary conditions. We conclude that the polynomial-fit boundary conditions are superior to the symmetric/anti-symmetric boundary conditions, provided that the fluctuations of the solution are well-resolved.

7 Conclusion

We have shown the appearance of rapidly varying functions in MHD situations which involve line-tying at boundaries. Such short scale variations appear either as boundary layers, near the conducting boundaries, or as small amplitude ripples riding on the solution envelope. The width of the boundary layers scales as the aspect ratio of the disturbance, ie, the ratio of layer width to parallel wavelength scales as $k_{\parallel}/k_{\perp}$, where $k_{\parallel}$ and $k_{\perp}$ are parallel and perpendicular wavenumbers, and $k_{\parallel}/k_{\perp}\ll 1$. As far as the ripples, the ratio of ripple wavelength to parallel wavelength scales as $\sqrt{1/\beta}$. Clearly, these smaller scales would have to resolved. Nonlinear simulations we have performed, preliminarily, show that numerical instability can be triggered from the "noise" if there is insufficient resolution. We note also that the ripple amplitude is small compared to the envelope, scaling as $\sqrt{\beta}$.

Our findings exemplify the caution of Gresho & Lee (1981); to quote these authors: "Don’t suppress the wiggles – they’re telling you something!". Since rapidly varying solutions are a real consequence of line-tied boundary conditions, it is imperative that sufficient grid resolution be ensured and appropriate boundary conditions implemented in numerical calculations. In particular, diffusion layers near the boundary, if implemented, should be used with caution, depending on the importance of line-tying physics to the problem at hand.