2 Equations

As mentioned, we investigate the system shown in figure 1. Ideal MHD equations in two dimensions are considered. In usual notation, the governing equations are:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}n+\unicode[STIX]{x1D735}\boldsymbol{\cdot }n\boldsymbol{u}=0 & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle nMd_{t}\boldsymbol{u}=-T\unicode[STIX]{x1D735}n+j_{y}\hat{\boldsymbol{y}}\times \boldsymbol{B} & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle j_{y}=\hat{\boldsymbol{y}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{B} & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{B}=\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B}) & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0, & \displaystyle\end{eqnarray}$$

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 figure 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}=1=M$ . Thus, energies are normalized to the magnetic energy, and, in normalized units, $T$ is of order the plasma $\unicode[STIX]{x1D6FD}$ .) 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 (2.1)–(2.5), we obtain a linearized system in the variables $n$ , $u_{x}$ , $u_{z}$ and $B_{x}$ , given by

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}n+\unicode[STIX]{x2202}_{x}u_{x}=-\text{i}ku_{z} & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}u_{x}=-T\unicode[STIX]{x2202}_{x}n & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}u_{z}=-\text{i}kTn+(\text{i}/k)\unicode[STIX]{x1D6FB}^{2}B_{x} & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}B_{x}=-\text{i}ku_{z}, & \displaystyle\end{eqnarray}$$

where, here, $k$ is assumed non-zero, and $\unicode[STIX]{x2202}_{x}B_{x}=-\text{i}kB_{z}$ . The Laplacian $\unicode[STIX]{x1D6FB}^{2}$ is defined as $\unicode[STIX]{x2202}_{x}^{2}-k^{2}$ . Here, we have assumed solutions of the form $\exp (\text{i}kz)$ . (For $k=0$ , equations (2.8) and (2.9) are replaced by $\unicode[STIX]{x2202}_{t}u_{z}=\unicode[STIX]{x2202}_{x}B_{z}$ , and $\unicode[STIX]{x2202}_{t}B_{z}=\unicode[STIX]{x2202}_{x}u_{z}$ .) The quantities $n$ , $u_{x}$ , $u_{z}$ , $B_{x}$ , $B_{z}$ , are all small perturbations, compared to $n_{0}$ and $B_{0x}$ . These equations are fourth order in $\text{d}/\text{d}t$ 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.

We note here that, clearly, the above set of equations can be solved as an eigenvalue problem with solutions of the form $\exp (\unicode[STIX]{x1D6FE}t)$ , $\unicode[STIX]{x1D6FE}$ being constant. Homogeneous Dirichlet boundary conditions on $u_{x}$ and $u_{z}$ at $x=\pm 1$ ensure a well-posed problem. In this paper, we will, rather, solve these equations asymptotically by using two small parameters: we will generally assume that $\unicode[STIX]{x1D6FD}\ll 1$ , $T\sim O(\unicode[STIX]{x1D6FD})$ , leading to ‘reduced equations’, and we will assume that we are considering long wavelength modes in $x$ , i.e. the $x$ -scale of the modes considered is of $O(1)$ . The parameter $k/K$ , then, could be taken to be large, in subsidiary expansion, when needed. These approximations help to separate the slow and fast modes and to highlight the physics. We reiterate that 2-D ( $x$ and $z$ ) displacements of the system in figure 1 lead to slow waves and fast waves. We have not considered here perturbations into the page, $u_{y}$ and $B_{y}$ ; the latter would result in pure Alfvén waves.

3 Mismatch

In this section, we explore some physics of the two types of modes in the system. While the physics is well known, we describe the modes in the context of how the spatial solutions of the various variables interact to possibly cause ‘mismatches’, precipitating sharp structures.

We begin by noting that if $k$ is set to zero in (2.6)–(2.9), the system is one-dimensional. In that case, the modes clearly separate into a pure sound wave and a pure Alfvén 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 behaviour in $u_{x}$ , i.e.  $\cos (Kx)$ . The latter is consistent with the hard-wall anti-symmetric boundary conditions for $u_{x}$ .

Similarly, the $k=0$ counterparts of (2.8) and (2.9) set up an Alfvén wave, coupling $u_{z}$ to $B_{z}$ . Here, $u_{z}$ is co-sinusoidal, i.e.  $\cos (Kx)$ , consistent with the anti-symmetric line-tied boundary condition, and $B_{z}$ is sinusoidal (and so symmetric at the boundaries).

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’ in the coupling between slow and fast mode spatial sinusoidal parities. To facilitate discussion, let us assume that $T\ll 1$ . In this case, we will show, the slow wave is almost sound like, and the fast wave is almost Alfvén like. Now, consider the slow wave. The sinusoidal parities of this sound-like mode were described above. However, if one examines (2.6), within the assumption of the parities of $n$ and $u_{x}$ just mentioned, the parity that would ‘fit’ the variable $u_{z}$ , from (2.6), would be sinusoidal, i.e.  $\sin (Kx)$ . But this parity violates the assumption of anti-symmetric boundary conditions for $u_{z}$ . From (2.9), a similar mismatch also occurs for $B_{x}$ . This violation is what we refer to as a ‘mismatch’. We will show that this type of mismatch gives rise to a sharply varying solution for $B_{x}$ , since $B_{x}=0$ is enforced at the plates.

A similar mismatch also arises in the fast mode for $k$ non-zero. The fast mode, for low $T$ , is largely 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 anti-symmetric at the plates; while $B_{x}\sim \cos (Kx)$ , consistent with $B_{x}$ also being anti-symmetric at the plates. However, from (2.8), the parities as deduced for $u_{z}$ and $B_{x}$ would force $n$ to be $\cos (Kx)$ , which is anti-symmetric at the boundaries, and so violates the symmetric boundary condition requirement.

We will now show below that these mismatches play a role in formations leading to multiscale structure.

4 Analytical solution – slow mode

We proceed to obtain analytically the slow and fast eigenmodes for the system in figure 1. We begin by rewriting the four equations (2.6)–(2.9) as two coupled equations. We use (2.7) in (2.6) and (2.9) in (2.8) to get the following 2-field equations for $n$ and $B_{x}$

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6FE}^{2}-T\unicode[STIX]{x2202}_{x}^{2})n=\unicode[STIX]{x1D6FE}^{2}B_{x} & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6FE}^{2}-\unicode[STIX]{x1D6FB}^{2})B_{x}=-k^{2}Tn, & \displaystyle\end{eqnarray}$$

where solutions of the form $\exp (\unicode[STIX]{x1D6FE}t)\exp (\text{i}kz)$ are assumed, and $\unicode[STIX]{x1D6FE}$ is the eigenvalue.

We now make the low- $\unicode[STIX]{x1D6FD}$ assumption $T\sim O(\unicode[STIX]{x1D6FD})\ll 1$ . 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 condition $\unicode[STIX]{x1D6FD}\ll 1$ ensures that the sound speed is subdominant to the Alfvén speed.

We begin with the slow mode, wherein we anticipate that $\unicode[STIX]{x1D6FE}^{2}\ll 1$ (the slow mode is sub-Alfvénic). In this case, from (4.2), we can discard the $\unicode[STIX]{x1D6FE}^{2}$ term; but we assume the mode is long wavelength in $x$ , i.e.  $\text{d}/\text{d}x\sim O(1)$ . Further, in ‘optimal ordering’, we assume $k$ , for the moment, is also of order unity, and, so, the Laplacian is ordered to be $O(1)$ . Since the right-hand side is proportional to $T\sim \unicode[STIX]{x1D6FD}\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 (4.2)

(4.3) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}B_{x}=k^{2}Tn.\end{eqnarray}$$

From (4.1), comparing the two $\unicode[STIX]{x1D6FE}^{2}$ terms, we see that the term on the right-hand side must be neglected since $B_{x}$ is small, making the right-hand side doubly small. The remaining two terms on the left-hand side are of the same order, and so (4.1) reduces as

(4.4) $$\begin{eqnarray}(\unicode[STIX]{x1D6FE}^{2}-T\unicode[STIX]{x2202}_{x}^{2})n=0.\end{eqnarray}$$

Equation (4.4) is now an eigenvalue equation for $n(x)$ with eigenvalue $\unicode[STIX]{x1D6FE}^{2}$ . The $n$ eigenfunction, and the eigenvalue, invoking the appropriate hard-wall boundary condition, are

(4.5a,b ) $$\begin{eqnarray}n=\sin (Kx),\quad \unicode[STIX]{x1D6FE}^{2}=-K^{2}T.\end{eqnarray}$$

This is a long wavelength sound wave. (The boundary condition, $\unicode[STIX]{x2202}_{x}n=0$ , is inferred from (2.7), given that $u_{x}=0$ at the boundary.)

Turning to (4.3), 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 Dirichlet boundary condition, is

(4.6) $$\begin{eqnarray}B_{x}=-T[\sin (Kx)-\sinh (kx)/\sinh (k)].\end{eqnarray}$$

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

Figure 2. Value of $B_{x}(x)$ plotted from (4.6);  $k/K$ chosen to be 20 to highlight boundary layers at the $x=\pm 1$ plates.

We note here that we have not used any asymptotic analysis, i.e. separating equation sets between ‘outer equations’ and ‘inner layer equations’, to arrive at the solution for $B_{x}$ (4.6). However, we clearly see a boundary layer type formation in figure 2, provided $k/K\gg 1$ . This is because for our calculation we had assumed an ‘optimal ordering’, $k\sim K$ , and we were fortunate enough to be able to solve the full equation. Alternatively, if we had assumed $k/K$ large, at the outset, 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 traditional boundary layer analysis, leading to an ‘inner solution’ matching asymptotically with the ‘outer solution’ (Bender & Orszag Reference Bender and Orszag1999). The ‘inner solution’ would have exhibited the sharp variation seen in figure 2. (In comment, we add that we have indeed verified this behaviour. 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 traditional multiscale asymptotic methods.

It is worth noting here that the emergence of a rapid variation in the $k/K$ large case has an important consequence. In particular, the ‘reduced equations’ of MHD (Strauss Reference Strauss1976) are derived under the assumption of long parallel wavelengths compared with perpendicular wavelengths, generally stated as $k_{\Vert }\ll k_{\bot }$ . But this is exactly our condition above. Thus, the application of the reduced equations at conducting surfaces should require a multiscale layer analysis near the conductors. This is a situation indeed encountered. The appearance and necessity of sharply varying solutions for reduced equations has been reported previously (Scheper & Hassam Reference Scheper and Hassam1999). In the latter paper, it was shown that low frequency motions near the solar conducting surface would require a boundary layer type 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 ( $\unicode[STIX]{x2202}_{t}=0$ ), with density having $k$ and $K$ variations, can be readily found. This too exhibits sharp variations at the boundaries. Emergent magnetic flux from the solar surface indeed is thought to exhibit sharply varying structures (Parker Reference Parker1979, p. 269, figure 10.12).

5 Analytical solution – fast mode

We now turn to the fast mode. We begin as in § 4 with (4.1) and (4.2). This time, we order $\unicode[STIX]{x1D6FE}^{2}$ as $O(1)$ , as we anticipate an Alfvén mode. Since $T\sim \unicode[STIX]{x1D6FD}\ll 1$ , and we assume $\unicode[STIX]{x2202}_{x}^{2}$ is $O(1)$ , for long wavelengths in the outer equations, equation (4.1) reduces to

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}^{2}n=\unicode[STIX]{x1D6FE}^{2}B_{x}.\end{eqnarray}$$

This yields $n=B_{x}$ . Using the latter in (4.2), the right-hand side of (4.2) must be discarded, being small, leading to

(5.2) $$\begin{eqnarray}(\unicode[STIX]{x1D6FE}^{2}-\unicode[STIX]{x1D6FB}^{2})B_{x}=0.\end{eqnarray}$$

The latter equation is an eigenvalue equation for $B_{x}$ , again for long wavelengths and optimal ordering for $k/K$ . With line-tied boundaries, the eigenfunction and eigenvalue are

(5.3) $$\begin{eqnarray}B_{x}=\cos (Kx),\quad \unicode[STIX]{x1D6FE}^{2}=-[K^{2}+k^{2}].\end{eqnarray}$$

From (5.1), 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 layer type behaviour as $x\rightarrow \pm 1$ . In what follows, we show that indeed boundary layer type analysis is needed. A rapid variation layer in solution generally indicates rapid variation in the $x$ -direction, i.e.  $\text{d}/\text{d}x\gg 1$ . As a consequence, there may be sharp variations in the eigenfunction, or, alternatively, a multiple scale character to the solution (Bender & Orszag Reference Bender and Orszag1999). For this paper, we have investigated both possibilities and we find that the latter character, multiscale behaviour, is the apt description, at least in the case of the fast mode. We proceed below with a multiscale analysis.

We go back to (4.1). We note that $\unicode[STIX]{x1D6FE}^{2}$ is already known to be $O(1)$ . But this time we also allow $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x$ large. In that case, the term $T\unicode[STIX]{x2202}_{x}^{2}$ must be taken to be optimally $O(1)$ , and, so, retained in the equation. Thus, the multiple scale equation from (4.1) remains

(5.4) $$\begin{eqnarray}(\unicode[STIX]{x1D6FE}^{2}-T\unicode[STIX]{x2202}_{x}^{2})n=\unicode[STIX]{x1D6FE}^{2}B_{x}.\end{eqnarray}$$

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

(5.5) $$\begin{eqnarray}n(x)=\cos (Kx)-A\cos (\unicode[STIX]{x1D705}x),\quad \unicode[STIX]{x1D705}\gg K,\end{eqnarray}$$

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

(5.6) $$\begin{eqnarray}n(x)=\cos (Kx)-[(K/\unicode[STIX]{x1D705})/\sin (\unicode[STIX]{x1D705})]\cos (\unicode[STIX]{x1D705}x),\quad \unicode[STIX]{x1D705}\gg K.\end{eqnarray}$$

This solution is plotted in figure 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 $\unicode[STIX]{x1D705}/K$ , but also that the amplitude of the ripple is small, by the same factor. The factor $\unicode[STIX]{x1D705}/K$ , using (5.3) for $\unicode[STIX]{x1D6FE}$ , can be seen to be of $O(1/\unicode[STIX]{x1D6FD}^{1/2})$ , and thus the scale of the ripple and its amplitude both go as $\unicode[STIX]{x1D6FD}^{1/2}$ . It can also be noted that the $A\cos (\unicode[STIX]{x1D705}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 $(\unicode[STIX]{x1D705}/K)$ as above. The structure of the boundary layer type behaviour can be seen upon expanding (5.6) near the boundaries. At the left boundary, we let $x=-1+s$ , where $s\ll 1$ and $\unicode[STIX]{x1D705}s\sim 1$ . Equation (5.6) then becomes

(5.7) $$\begin{eqnarray}n(x\rightarrow -1)=K[s-(1/\unicode[STIX]{x1D705})\sin (\unicode[STIX]{x1D705}s)-(1/\unicode[STIX]{x1D705})(\cos \unicode[STIX]{x1D705}/\sin \unicode[STIX]{x1D705})].\end{eqnarray}$$

From this solution, we see that the scale of the sharp variation is given by $\unicode[STIX]{x1D705}s\sim 1$ .

Figure 3. Value of $n(x)$ plotted from (5.6), at $T=0.2$ and $k=20$ . Sharp structured ripple riding on long wavelength envelope of $n(x)$ evident. Note $\text{d}n/\text{d}x$ is forced to be zero at the $x=\pm 1$ plates.

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 Reference Trefethen2000). We subjected the linear equations, (2.6)–(2.9), to an initial condition given by $n=((1+\cos (\unicode[STIX]{x03C0}x))/2)^{4}\exp (\text{i}kz)$ . 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 discretization; hence, we kept the time step $\unicode[STIX]{x1D6E5}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 figure 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.

Figure 4. Solutions of $n$ and $u_{x}$ at $t=0,5,9,12$ from the $N=200$ Chebyshev run.

Figure 5. Finite-difference solutions at $t=9$ with various numbers of grid points $N$ and boundary conditions. Black dashed lines are the $N=200$ Chebyshev solution.

The initial condition for $n$ in the above simulation was chosen to demonstrate the fact that the early time solution can be relatively noise free, but that ‘noise’ can evolve for later times. This is consistent with the fact that the initial condition is more slow wave type in early times but, being a general initial state, couples to and excites fast wave modes at later times.

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 §§ 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

(6.1) $$\begin{eqnarray}f=Ax+Bx^{2}.\end{eqnarray}$$

This gave the following conditions for the ghost points, viz.

(6.2) $$\begin{eqnarray}\displaystyle & \displaystyle f_{-1}=-3f_{1}+f_{2} & \displaystyle\end{eqnarray}$$

(6.3) $$\begin{eqnarray}\displaystyle & \displaystyle f_{-2}=-8f_{1}+3f_{2}. & \displaystyle\end{eqnarray}$$

Likewise, for $(\text{d}f/\text{d}x)(b)=0$ boundary conditions, we fitted to

(6.4) $$\begin{eqnarray}f=f_{0}+Ax^{2}+Bx^{3}.\end{eqnarray}$$

This gave the following conditions for the ghost points, viz.

(6.5) $$\begin{eqnarray}\displaystyle & \displaystyle f_{-1}=-(3/2)f_{0}+3f_{1}-(1/2)f_{2} & \displaystyle\end{eqnarray}$$

(6.6) $$\begin{eqnarray}\displaystyle & \displaystyle f_{-2}=-12f_{0}+16f_{1}-3f_{2}. & \displaystyle\end{eqnarray}$$

We then re-ran our test example with the two sets of boundary conditions. Figure 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 § 3. In contrast, the solution with polynomial-fit boundary conditions is noticeably smoother.

Figure 6. Scalings of errors with respect to $N$ for different implementations of spatial discretization.

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. Figure 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.