text:
# A practical spectral method for hyperbolic conservation laws

## I Introduction

It is well known that due to their accuracy and efficiency, spectral methods have great advantages over local methods in solving applicable scientific and engineering problems . Given a hyperbolic system of nonlinear conservation laws

$`\text{u}_t+\text{f}(\text{u})_x=0,`$ (1)
with initial condition

$`\text{u}(x,0)=\text{u}_0(x),`$ (2)
its solution may not exist in the classical sense because of possible discontinuities. The direct use of spectral methods to this problem will encounter the Gibbs oscillations , which lead to blow-up in the time integration. Therefore, it has been of great interest in making spectral methods applicable to the hyperbolic conservation law systems in the past two decades. The expectation, as point out by Gottlieb , is that the use of spectral methods will enable us not merely to capture the shock, but also to capture the delicate features and structures of the flow. Obviously, spectral approaches will be extremely valuable to a class of aerodynamical problems that involve the interaction of both turbulence and shock . Such an interaction is some of the most challenging problems in computation fluid dynamics.

The main objective of previous studies is to recover smooth solutions from those that are contaminated by Gibbs oscillations and, meanwhile, to improve the rate of convergence. For shock-capturing, many up-to-date local methods, such as weighted essentially non-oscillatory (WENO) scheme , central schemes , arbitrary-order non-oscillatory advection scheme , Gas kinetic , and image processing based schemes , perform well. The success of these local shock-capturing schemes lies in their appropriate amount of intrinsic numerical dissipation, which is introduced either by explicit artificial viscosity terms, or by upwinding, or by appropriate local average (non-oscillatory central schemes), or by relaxation . The characteristic decomposition based on Roe’s mean matrix can also be considered as a local averaging on the Jacobian matrix $`\text{A(u)}=\text{f(u)}/\text{u}`$. In his recent work, LeVeque  outlined the relation between some approximate Riemann solvers and relaxation schemes. The relative amount of numerical dissipations may explain why the local characteristic decomposition is not necessary in low-order methods while it seems to be indispensable in high-order methods . However, when these local methods are used in the cases where flow structures with fine details are needed to be resolved together with shocks, numerical dissipations are usually found to be so large that they smear the fine details . The spectral methods, on the contrary, contribute very little numerical dissipation and dispersion in principle when used to approximate spatial derivatives. Nevertheless, when they are applied to the approximation of spatial derivatives on a domain containing discontinuities, the Gibbs oscillations again have to be suppressed by appropriate means, for which there are two general routines: (1) explicit artificial viscosity, e.g. spectral viscosity (SV) method proposed by Tadmor  and (2) filtering, which is the central issue of the present study.

At present, most filters are constructed in the spectral domain and they can be called spectral filters. Hussaini et al.  suggested a list of typical spectral filers, i.e. Lanzos filter, raised cosine filter, sharpened raised cosine filer and exponential cutoff filter. More sophisticate and effective filters in spectral domain are Vandeven’s $`p`$-th order filter , Cai et al.’s sawtooth function  and Gottlieb and Tadmor’s regularized Dirichlet function .

Apart from methods that are implemented in the spectral domain, filters in the physical domain are also developed. To design an appropriate filter in the physical domain is generally more difficult than in the Fourier domain. A straightforward procedure is to make use of numerical dissipation partially contained in high-order shock-capturing schemes , such as the ENO scheme, where, actually, numerical dissipations were introduced both in the Fourier domain (via an exponential filter) and in the physical domain (via ENO polynomial interpolation). Such a strategy was generalized by Yee et al.  to the so-called characteristic filter method (with finite difference or finite volume) which has been successfully applied to shock-capturing and turbulence simulations. Another procedure was initiated by Gottlieb  by using the Gegenbauer polynomial to resolve the oscillatory partial Fourier summation. Some promising numerical results from the filter approach can be found in Refs. .

Nevertheless, some prominent issues may hinder the success of filtered spectral methods for practical numerical computations involving shocks. In our view, there are two very important issues that are vita to the success of a filtered spectral scheme. The first issue is how to select optimal filters for shock-capturing. The issue is very complex and cannot have a unique answer at this point because there are so many properties to mind, such as flatness, ripple, filter length, effective frequency range and length of transition band, to name only a few. Loosely speaking, it is desirable to have filters that are free of dispersion errors, flat while having very small transition band, short in length while having high resolution. Moreover, what adds to the complexity is that solutions to different hyperbolic conservation law systems may have different Fourier spectral distributions. Consequently, one faces the difficulty that an optimal filter for a problem may not even work for another one. Therefore, it is desirable to have a filter which accounts for this change by an adjustable parameter.

The second issue also concerns the change of solution, but is due to the time evolution for a given hyperbolic conservation law system. As the Fourier spectral distribution of a given problem changes with time, the question is how to implement the filter? It will be too dissipative for many systems if a filter is applied at each time step of the integration. Therefore, an adaptive implementation, which is controlled by a sensor during the time integration, is appropriate. Given the complexity, it is unlikely that there will be an ideal solution to these issues in the immediate future. Obviously, these challenging and interesting problems call for the further study of spectral filter approaches.

The objective of this work is to construct, analyze and implement a class of high-order lowpass filters which can be used either in the Fourier domain or in the physical domain. These filters are constructed within the framework of the discrete singular convolution (DSC) algorithm  and they are used to suppress Gibbs oscillations arising in the Fourier pseudospectral method (FPM). Such lowpass filters have a wide effective wavenumber range that makes it possible to capture the fine flow structures, which is a desirable objective of spectral methods for hyperbolic conservation law systems. It is this effective wavenumber range that controls the resolution of the overall method. In our design, this effective range can be varied according to the resolution requirement by a parameter. Different DSC kernels have different magnitude responses and adjustability in the Fourier domain, which in turn influence the accuracy and resolution. The DSC filters used in this work were originally designed in a conjugate filter oscillation reduction (CFOR)  scheme. By ‘conjugate filters’ it means that the effective wavenumber range of the lowpass filter is largely overlapped with that of the highpass filter which is used for the approximation of spatial derivatives. The resulting Fourier domain algorithm is essentially a windowed Fourier pseudospectral method for shock-capturing.

The rest of this paper is organized as follows. In Section II, a brief introduction will be devoted to the construction of the DSC lowpass filters followed by a discussion of two types of implementations of these filters to spectral methods. Extensive numerical experiments are carried out in Section III, including Burgers’ equation, shock-tube problems, shock-entropy wave interaction, shock-vortex interaction. The shock-capturing ability and the high resolution feature of the present method are addressed. Concluding remarks end this paper.

## II Theory and algorithm

### A DSC lowpass filter

In the context of distribution theory, a singular convolution is defined by

$`f(x)=(T\eta )(x)={\displaystyle _{\mathrm{}}^{\mathrm{}}}T(xt)\eta (t)𝑑t,`$ (3)
where $`T(x)`$ is a singular kernel and $`\eta (x)`$ is an element of the space of test functions. Interesting singular kernels include that of Hilbert type, Abel type and delta type. The former two play important roles in the theory of analytical functions, processing of analytical signals, theory of linear responses and Radon transform. Since delta type kernels are the key element in the theory of approximation and the numerical solution of differential equations, we focus on singular kernels of the delta type

$`T(x)=\delta ^{(q)}(x),q=0,1,2,\mathrm{},`$ (4)
where superscript $`(q)`$ denotes the $`q`$th-order ‘derivative’ of the delta distribution $`\delta (x)`$, with respect to $`x`$, which should be understood as generalized derivatives of distributions. When $`q=0`$, the kernel $`T(x)=\delta (x)`$ is important for the interpolation of functions. In this work, only the case of $`q=0`$ will be involved. One has to find appropriate approximations to the above singular kernel, which can not be directly realized in computers. To this end, we consider a sequence of approximations

$`\underset{\alpha \alpha _0}{lim}\delta _\alpha ^{(q)}(x)=\delta ^{(q)}(x),q=0,1,2,\mathrm{}`$ (5)
where $`\alpha `$ is a parameter which characterizes the approximation with the $`\alpha _0`$ being a generalized limit. Among various candidates of approximation kernels , regularized Shannon kernel (RSK)

$`\delta _{\sigma ,\mathrm{\Delta }}(xx_k)={\displaystyle \frac{\mathrm{sin}\frac{\pi }{\mathrm{\Delta }}(xx_k)}{\frac{\pi }{\mathrm{\Delta }}(xx_k)}}\mathrm{exp}\left({\displaystyle \frac{(xx_k)^2}{2\sigma ^2}}\right)`$ (6)
is frequently used due to its close link to the Shannon sampling theorem. In this formula, $`\mathrm{\Delta }`$ is the grid spacing and $`\sigma `$ determines the width of the Gaussian envelop. For a given $`\sigma 0`$, the limit of $`\mathrm{\Delta }0`$ reproduces the delta kernel (distribution). With the RSK, a function $`u`$ can be approximated by a discrete convolution

$`u(x){\displaystyle \underset{k=W}{\overset{W}{}}}\delta _{\sigma ,\mathrm{\Delta }}(xx_k)u(x_k),q=0,1,2,\mathrm{},`$ (7)
where $`\{x_k\}_{k=W}^W`$ are a set of discrete grid points which are centered at $`x`$, and $`2W+1`$ is the computational bandwidth, or effective kernel support, which is usually smaller than the computational bandwidth of the spectral method - the entire domain span. Equation (7) is often referred as a DSC and the RSK is used as the prototype of a DSC lowpass filter.

The feasibility of the DSC-RSK as lowpass filters in numerical computation is based on an important observation: by adjusting the ratio $`r`$ ($`r=\frac{\sigma }{\mathrm{\Delta }}`$) of the regularizer, one can control the effective wavenumber region and thus, control the resolution of the lowpass filter. We plot the frequency response profile of the DSC-RSK lowpass filter in FIG. 1. A distinct property of the RSK lowpass filter is that there is no ripples in its passband; thus the application of this filter will not impart the low frequency. In contrast, some lowpass filters constructed with certain polynomials may have ripples in their passband, which would rule out their application to numerical integrations, especially for long time evolution, where the solution may be distorted after repeated filtering.

The maximum effective wavenumber (frequency) range of the DSC-RSK is determined by the ratio $`r`$, which in turn is determined by computational bandwidth $`W`$. The philosophy is simple since with the increase of $`W`$, the Gaussian window should also be enlarged to fully utilize the stencil. The optimal values $`r`$ at different $`W`$ can be estimated and they characterize the effective wavenumber range. The larger the optimal $`r`$ is, the wider the effective wavenumber range will be. When $`r`$ value is chosen larger than the optimal value for a given $`W`$, the accuracy of the resulting filter will be degraded by oscillations. Hence we do not recommend the use of $`r`$ beyond the optimal value unless in some circumstances where the filer is utilized to stabilize a long-time integration.

### B Fourier domain algorithm

We denote $`\zeta (x_j)=\delta _{\sigma ,\mathrm{\Delta }}(x_j)`$ the proposed DSC-RSK lowpass filter in the physical domain $`[0,L]`$ and $`\widehat{\zeta }(\omega _n)`$ its image in the Fourier domain, where $`\mathrm{\Delta }=L/N`$, $`x_j=j\mathrm{\Delta }`$ and $`\omega _n=\frac{2\pi n}{L}`$. The obvious relation between $`\zeta (x_j)`$ and $`\widehat{\zeta }(\omega _n)`$ is

$`\widehat{\zeta }(\omega _n)=\mathrm{\Delta }{\displaystyle \underset{j=0}{\overset{N1}{}}}\zeta (x_j)e^{i\omega _nx_j}.`$ (8)
As in the previous work , filtering in the Fourier domain is usually accomplished by multiplying the given Fourier coefficients $`\widehat{u}(\omega _n)=\mathrm{\Delta }_{j=0}^{N1}u(x_j)e^{i\omega _nx_j}`$ by $`\widehat{\zeta }(\omega _n)`$ directly, i.e., a windowed Fourier transform,

$`u^\zeta (x_j)={\displaystyle \frac{1}{L}}{\displaystyle \underset{n=\frac{N}{2}}{\overset{\frac{N}{2}1}{}}}\widehat{\zeta }(\omega _n)\widehat{u}(\omega _n)e^{i\omega _nx_j}.`$ (9)
Therefore, the modified Fourier coefficients $`\widehat{\zeta }(\omega _n)\widehat{u}(\omega _n)`$ in Eq. (9) are the Fourier coefficients of a localized version of the original function $`u(x_j)`$: $`u^\zeta (x_j)`$. Hence, errors arising from the discontinuity are also localized and the accuracy away from the discontinuity can be ensured. In practical applications, one only needs to carry out the filtering and differentiation $`D_xu=u_x`$ in one step

$`u_x^\zeta (x_j)={\displaystyle \frac{1}{L}}{\displaystyle \underset{n=\frac{N}{2}}{\overset{\frac{N}{2}1}{}}}i\omega _n\widehat{\zeta }(\omega _n)\widehat{u}(\omega _n)e^{i\omega _nx_j}.`$ (10)
A standard forth-order Runge-Kutta (RK-4) scheme is employed for the time advance scheme in this work. We use a TVD switch to adaptively activate the application of the DSC-RSK lowpass filter. Whenever the total variation of the approximate solution in the two consecutive time steps exceeds a prescribed criteria, the lowpass filter is activated. The DSC-RSK lowpass filter is only applied at the end of each full RK-4 cycle. The reader is referred to our earlier work  for the detail of this adaptive filtering.

### C Physical domain algorithm

It is more difficult to apply the filter in the physical domain than in the Fourier domain. On one hand, it is not easy to quantitatively analyze the influence of numerical dissipation in the physical domain. On the other hand, the convolution operation

$`u^\zeta (x)={\displaystyle u(x\xi )\zeta (\xi )𝑑\xi },`$ (11)
complicates the physical domain application of the filter. However, in some practical applications, it turns out to be quite necessary to implement the filter in the physical domain, especially for problems whose boundaries require some special treatments. The reflective boundary is one of such examples that need to be treated by using grid points outside the computational domain. In this work, we propose the method to apply our DSC-RSK lowpass filter in the physical domain. Due to the fact that the RSK lowpass filter is an interpolatory kernel, we consider a two-step procedure. Step one is the prediction of mid-mesh values from the nodal values

$`u_{i+\frac{1}{2}}^\zeta ={\displaystyle \underset{j=W}{\overset{1}{}}}u_{i+j+1}\zeta ((j+{\displaystyle \frac{1}{2}})\mathrm{\Delta })+{\displaystyle \underset{j=1}{\overset{W}{}}}u_{i+j}\zeta ((j{\displaystyle \frac{1}{2}})\mathrm{\Delta }).`$ (12)
The other is the reconstruction of the nodal values from the mid-mesh values

$`u_i^{\zeta \zeta }={\displaystyle \underset{j=W}{\overset{1}{}}}u_{i+j+\frac{1}{2}}^\zeta \zeta ((j+{\displaystyle \frac{1}{2}})\mathrm{\Delta })+{\displaystyle \underset{j=1}{\overset{W}{}}}u_{i+j\frac{1}{2}}^\zeta \zeta ((j{\displaystyle \frac{1}{2}})\mathrm{\Delta }),`$ (13)
where, $`u_i`$ denotes $`u(x_i)`$ and $`\mathrm{\Delta }`$ is the grid spacing. The lowpass filter with a $`r`$ value smaller than the optimal value is used in one of these two steps. The activation of the RSK lowpass filter is controlled in the same manner as that of the frequency domain algorithm.

As described, the DSC-RSK lowpass filter can be applied either in the physical domain or in the Fourier domain. In general, filtering in the Fourier domain is much easier and more convenient than in the physical domain because no additional summation is required to implement the filter in the Fourier spectral method. In contrast, the physical domain implementation requires two additional summations. Therefore, we have used the filter in the Fourier domain in most of our numerical experiments except for the cases in which the reflective boundary condition is involved. We have checked that whenever applicable, results obtained from the Fourier domain approach and the physical domain approach are identical.

### D Post processing filter

Once the time integration has been completed, a strong filter is usually required for post processing the solution so as to make the solution more presentable. In principle, we only need to reduce the $`r`$ value of the DSC-RSK lowpass filter at the last time step of the computation. As an alternative, we provide a series of lowpass filters based on the Lagrange delta kernel which is given by

$`\delta _W(x,x_j)={\displaystyle \underset{k=W,kj}{\overset{W}{}}}{\displaystyle \frac{xx_k}{x_jx_k}},j=W,\mathrm{},W.`$ (14)
FIG. 2 shows the magnitude response of the Lagrange delta kernel. In fact, the highpass versions of these filters are graphically identical to centered finite difference (FD) schemes. In our numerical experiments, we adopt Lagrange-2 or Lagrange-4 as post processing filters. Their efficiency is demonstrated by our numerical results.

## III Numerical experiments

In this section, we examine the validity and demonstrate the performance of proposed windowed (filtered) Fourier spectral scheme for conservation law systems. A large collection of standard linear and nonlinear benchmark problems , including the linear advection equation, 1D inviscid Burgers’ equation, shock tube (Sod and Lax) problems, 1D shock-entropy interaction, 2D shock-entropy interaction, 2D shock-vortex interaction, 2D advection of an isentropic vortex and flow past a cylinder, are considered in the present study. Some of these problems have periodic boundary conditions therefore the windowed Fourier pseudospectral method will be applied directly. For problems whose boundary conditions are not periodic, we symmetrically double the computational domain in both $`x`$ and $`y`$ directions to compute the flux derivative. It is noted that the extension is necessary for constructing periodic problem feasible for the Fourier spectral method. Our filtering procedures, as described earlier, can be applied with general boundary conditions, such as Dirichlet or Neumann. In the appendix, we list the filter parameter $`r`$ used for each numerical example.

### A Scalar conservation law systems

We begin our numerical experiments with the scalar conservation law system, which is given by

$`u_t+f(u)_x=0,`$ (15)
where $`f(u)`$ is a function of $`u`$. It is generally believed that spectral methods are not very suitable for simple shock profiles and should not be used if one is interested primarily in the shock profiles rather than other detailed structures of the flow. In this study, we choose scalar conservation law systems to illustrate that the proposed global spectral scheme performs well for a class of localized shock problems.

Example 1. We first solve a linear advection equation, which is of the form

$`\begin{array}{cc}\hfill u_t+u_x=& 01<x<1,\hfill \\ \hfill u(x,0)=& u_0(x)\text{periodic,}\hfill \end{array}`$ (18)
where, $`u_0`$ is an initial value

$`u_0(x)=\{\begin{array}{cc}\frac{1}{6}(G(x,\beta ,z\delta )+G(x,\beta ,z+\delta )+4G(x,\beta ,z))\hfill & 0.8x0.6\hfill \\ 1\hfill & 0.4x0.2\hfill \\ 1|10(x0.1)|\hfill & 0x0.2\hfill \\ \frac{1}{6}(F(x,\alpha ,a\delta )+F(x,\alpha ,a+\delta )+4F(x,\alpha ,a))\hfill & 0.4x0.6\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}`$ (24)
Here, functions $`G`$ and $`F`$ are defined as

$`\begin{array}{cc}\hfill G(x,\beta ,z)=& e^{\beta (xz)^2}\hfill \\ \hfill F(x,\alpha ,a)=& \sqrt{\text{max}(1\alpha ^2(xa)^2,0)},\hfill \end{array}`$ (27)
where $`a=0.5`$, $`z=0.7`$, $`\delta =0.005`$, $`\alpha =10`$ and $`\beta =\frac{\text{log}2}{36\delta ^2}.`$

The solution of this problem contains a combination of a Gaussian, a square wave, a sharp triangle wave and a half ellipse. The exact solution at any time can be easily obtained by a translation of the initial solution at speed 1. In FIG. 3 and FIG. 4, we display the numerical results for this problem at $`t`$=8 with 128 and 256 grid points, respectively. We observe that our method works well in both cases.

Example 2. Secondly, we test the proposed method by considering a moving W-shape wave, i.e., a piecewise continuous initial value, for the linear advection equation with an initial value

$`u_0(x)=\{\begin{array}{cc}1\hfill & 0x0.2\hfill \\ 4x\frac{3}{5}\hfill & 0.2x0.4\hfill \\ 4x+\frac{13}{5}\hfill & 0.4x0.6\hfill \\ 1\hfill & 0.6x0.8\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}`$ (33)
This case is very similar to Example 1. It contains the so-called contact discontinuity and is quite difficult to solve in hyperbolic conservation laws. In particular, the interaction of two lines at $`x=0.2`$ and $`x=0.6`$ is difficult to resolve. Hence, it is a good test for the shock-capturing ability of the present method. We compute the solution up to $`t`$=8. The results are shown in FIG. 5 and FIG. 6.

Example 3. Having solved the linear equation with different initial values, we consider the solution of the most popular inviscid Burgers’ equation, in which $`f(u)=u^2/2`$ and the Riemann type initial value is given by

$`u(x,0)=\{\begin{array}{cc}1\hfill & x0\hfill \\ 0\hfill & x>0.\hfill \end{array}`$ (36)
This is a standard benchmark problem in hyperbolic conservation laws and has been considered by numerous researchers. The exact solution is a shock wave with a constant velocity

$`u(x,t)=\{\begin{array}{cc}1\hfill & xSt<0\hfill \\ 0\hfill & xSt>0,\hfill \end{array}`$ (39)
where the speed of the shock front is

$`S={\displaystyle \frac{1}{2}}.`$ (40)
We should note that this case is a non-periodic boundary problem. We must symmetrically double the computational domain in the $`x`$-direction when we calculate the derivative of the flux. This operation is necessary because it generates the periodic boundary condition so that the Fourier spectral method can be used. The numerical results are plotted in FIG. 7 at time $`t`$=2.

Example 4. Another Riemann type initial value for the inviscid Burgers’ equation with a flux of $`u^2/2`$ is given by

$`u(x,0)=\{\begin{array}{cc}0\hfill & x<0\hfill \\ 1\hfill & x0.\hfill \end{array}`$ (43)
The exact solution of this problem is a rarefaction wave

$`u(x,t)=\{\begin{array}{cc}0\hfill & \frac{x}{t}<0\hfill \\ \frac{x}{t}\hfill & 0<\frac{x}{t}<1\hfill \\ 1\hfill & \frac{x}{t}>1.\hfill \end{array}`$ (47)
Numerical results are plotted in FIG. 8 at $`t`$=2 with different numbers of nodal points. Particularly, one can see that the end of the rarefaction fan is well resolved.

Example 5. Finally, we use non-convex flux to test the convergence to the physically correct solution. The flux is a non-convex function

$`f(u)={\displaystyle \frac{1}{4}}(u^21)(u^24)`$ (48)
with a Riemann type initial value

$`u(x,0)=\{\begin{array}{cc}3\hfill & x<0\hfill \\ 3\hfill & x0.\hfill \end{array}`$ (51)
The reader is referred to  for more detailed information about this problem. The results are displayed in FIG. 9 at $`t`$=0.04. Again, present spectral solver yields a satisfactory resolution.

### B Euler systems in one dimension

In this subsection, we perform numerical experiments by using the proposed scheme for the 1D Euler equation of gas dynamics. In one dimension, the Euler equation takes the form

$`U_t+F(U)_x=0`$ (52)
with

$`U=\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right);F(U)=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ u(E+p)\end{array}\right),`$ (59)
where, $`\rho ,u,p`$ and $`E`$ denote the density, velocity, pressure and total energy per unit mass $`E=\rho (e+(u^2)/2)`$, respectively. Here, $`e`$ is the specific internal energy. For an ideal gas with the constant specific heat ratio ($`\gamma =1.4`$) considered here, one has $`e=p/(\gamma 1)\rho `$. We consider the following two well-known Riemann problems.

Example 6. We solve for the solution of two shock tube problems with Sod’s and Lax’s initial conditions , i.e.,

$`(\rho ,u,p)_{t=0}=\{\begin{array}{cc}(1,0,1),\hfill & x<0\hfill \\ (0.125,0,0.1),\hfill & \text{otherwise}\hfill \end{array}`$ (62)
for the Sod problem and

$`(\rho ,u,p)_{t=0}=\{\begin{array}{ccc}(0.445,0.698,3.528),\hfill & x<0\hfill & \\ (0.5,0,0.571),\hfill & \text{otherwise}\hfill & \end{array}`$ (65)
for the Lax problem. The numerical results for these two problems are shown in FIG. 10 and FIG. 11, respectively. We can see that the present method gives a good resolution in both cases. Similar to all the methods depending on the numerical dissipation, the resolution for the linear degenerate contact discontinuity is slightly lower than that for the generic nonlinear shock wave. If combined with Harten’s artificial compression method, present spectral solver should have better resolution for the contact discontinuity.

Example 7. In this example, the interaction of an entropy wave of small amplitude with a Mach 3 right-moving shock in a one-dimensional flow is investigated. The computational domain is taken as $`[0,9]`$ and the flow field is initialized with

$`(\rho ,u,p)_{t=0}=\{\begin{array}{cccc}(3.85714,\hfill & 2.629369,& 10.33333)& x0.5\hfill \\ (e^{ϵ\text{sin}(\kappa x)},\hfill & 0,& 1.0)& x>0.5\hfill \end{array}`$ (68)
where $`ϵ`$ and $`\kappa `$ are the amplitude and the wave number of the entropy wave before the shock. This problem is significant due to its relevance to the interaction of shock-turbulence. Spectral methods, by nature of their high accuracy and negligible dispersive and dissipative errors, are ideally suited to the numerical study of turbulence and have already been widely applied. Here we would like to argue the feasibility of the present spectral method to the shock-turbulence interaction. In our test, we vary the wavenumber of the pre-shock wave while keeping its amplitude unchanged ($`ϵ=0.01`$). As the wavenumber increases, the problem becomes more and more challenging because the amplified high-frequency entropy waves after the shock are always mixed with spurious oscillations. As the wavenumber increases, it is very difficult to distinguish the high-frequency waves from the spurious oscillations. A low-order scheme may dramatically damp the transmitted high frequency waves. Even some popular high-order schemes also encounter the difficulty in preserving the amplitudes of the entropy waves due to their excessive dissipation with a given mesh size. Therefore, a successful shock-capturing method should be able to eliminate Gibbs’ oscillation while capturing the shock and preserving the high-frequency entropy waves.

In our test, we let the shock move from $`x=0.5`$ to $`x=8.5`$. For the purpose of comparison with the previous results, we only give the results at the interval $`[4.0,9.0]`$. Also, in order to discharge the transient waves due to the non-numerical initial shock profile, we plot the length of the amplified entropy waves in the same manner as that in Ref. . Furthermore, we must point out two nontrivial remarks for the following cases with different $`\kappa `$. One is that we have replotted the final results in a denser grid, otherwise the bounded strip will not be fully spanned due to the fact that no enough grid values locate near all the peaks /valleys of the wave. Another is that we have used a local filter in order to eliminate the oscillations near the shock when we present final results after the computation.

First, we consider $`\kappa =13`$. In this case, the computational domain is deployed with 513 grids points. Such a mesh is suitable for the Fourier spectral method. The numerical results of the amplitude of the entropy waves are shown in FIG. 12(a). It can be seen that the shock and the generated entropy waves are ideally captured. The entropy waves fully span the strip bounded by two solid lines at $`\pm `$ 0.08690716, which show the amplitude of the amplified entropy waves predicted by the linear analysis. Since post-shock waves are monochromic, a simple calculation indicates the point per wavelength (PPW) value of these waves is 7.125. Such value is larger than the spectral resolution (PPW=2.0) but it is much smaller than that obtained by most shock-capturing schemes. The transmitted waves suffer the similar smearing as the contact waves. The shock wave, on the contrary, is self-repairable due to its compression nature. This observation verified that the numerical dissipation needed for shock-capturing will degrade the resolution of the spectral method. It also verified that this degradation can be isolated and alleviated by the DSC-RSK filter.

We next double the wave number $`\kappa `$ up to 26. A mesh of 513 grid points is not enough for simultaneous shock capturing and high-frequency wave resolving. Therefore, we increase the mesh size to $`N=1025`$ for this computation. The results for this case are presented in FIG. 12(b). It is seen that no obvious deterioration happens to the results and the resolution in the post-shock waves is still satisfactory. The fullness of the wave is comparable with that of $`\kappa =13`$. But the first two post-shock waves have been slightly polluted.

When the wave number $`\kappa `$ is increased to 39 and 52, we use 2049 grid points to resolve such high-frequency waves. FIG. 12(c) and FIG. 12(d) show the results of $`\kappa =39`$ and $`\kappa =52`$, respectively. The same PPW value $`7.125`$ is maintained in this high-frequency case. The good resolution of the present method in simulating this 1D shock entropy interaction indicates that our scheme is capable of and efficient in distinguishing the high-frequency entropy waves from spurious oscillations.

Example 8. Results are now shown for the problem of Shu and Osher which also describes the interaction of an entropy sine wave with a Mach 3 right-moving shock. The computational domain is taken as $`[1,1]`$ and the flow field is initialized with

$`(\rho ,u,p)_{t=0}=\{\begin{array}{cccc}(3.85714,\hfill & 2.629369,& \hfill 10.33333)& \hfill x0.8\\ (1.0+ϵ\text{sin}(\kappa \pi x),\hfill & 0,& \hfill 1.0)& \hfill x>0.8,\end{array}`$ (71)
where $`ϵ`$=0.2 and $`\kappa `$=5. FIG. 13(a) and FIG. 13(b) show the solution of the density with 128 and 256 cells, respectively. The ‘Exact’ solution here is the solution computed by the third-order ENO with 1200 cells. As evident in the solution, the complicated flow field behind the shock is well resolved and the shock remains sharp. However, small oscillations near the shock still exist due to the fact that we do not attempt to locate this shock intentionally. If we predict the position of the shock and apply a local filter near the shock at the final time step, it is believed that more satisfactory results can be obtained.

### C Euler systems in two dimensions

In two dimensions, the Euler equation for gas dynamics in a vector notation takes the conservation form

$`U_t+F(U)_x+G(U)_y=0`$ (72)
with

$`U=\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ E\end{array}\right);F(U)=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ u(E+p)\end{array}\right);G(U)=\left(\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ v(E+p)\end{array}\right),`$ (85)
$`p=(\gamma 1)(E{\displaystyle \frac{1}{2}}\rho (u^2+v^2)),`$ (86)
where $`(u,v)`$ is the 2D fluid velocity and $`\gamma =1.4`$ is used in all computations.

Example 9. As the first example, we use our spectral scheme to study the interaction between a normal shock and a weak entropy wave which makes an angle $`\theta (0,\pi /2)`$ against the $`x`$-axis. If $`\theta =0`$, we have essentially the 1D problem (see example 7). The initial conditions are defined as follows: the right state of the shock is given as $`(\rho _r,u_r,v_r,p_r)`$=(1, 0, 0, 1) with a given shock Mach number $`M_s=3`$. We add a small entropy wave to the flow on the right of the shock which is equivalent to changing only the density of the flow on the right of the shock. In our problem we change the density $`\rho _r`$ in the right state of the shock by multiplying $`\rho _ϵ`$

$`\begin{array}{cc}\rho =\hfill & \rho _r\rho _ϵ\hfill \\ \rho _ϵ=\hfill & e^{(ϵ/p_r)\mathrm{sin}(\kappa (x\mathrm{cos}\theta +y\mathrm{sin}\theta ))}\hfill \end{array}`$ (89)
with $`ϵ`$ and $`\kappa `$ being the amplitude and wavenumber, respectively. In order to carry out a long time integration and to enforce periodic boundary condition in $`y`$-direction, the computational domain is taken to be $`[0,9]\times [0,\frac{2\pi }{\kappa \mathrm{sin}\theta }]`$. We initially position the normal shock at $`x=0.5`$ and allow it to move up to $`x=8.5`$. In our simulation, we adopt $`ϵ=0.1`$, $`\kappa =15`$, and $`\theta =30^{}`$. The performance is measured by drawing the maximum amplitudes of the amplified entropy waves in the $`y`$-direction for all fixed grid values $`x[7.4,8.4]`$, and comparing them with the amplitude predicted by the linear analysis, which is 0.08744786.

We use 513 grid points in the $`x`$-direction and 32 grid points in the $`y`$-direction in our test. In order to avoid peak/valley losses of the entropy waves due to insufficient grid points near these positions, we interpolate our final results to a denser grid when plotting them. FIG. 14 displays the amplitudes of the amplified entropy waves. Obviously, the compressed entropy waves are well resolved in the post-shock regime even though small oscillations remain near the location of the shock. It is believed that such a trivial flaw is acceptable.

This example has further demonstrated the capability of our spectral scheme for capturing small scale waves in the presence of a shock.

Example 10. This model problem describes the interaction between a stationary shock and a vortex. It has many potential applications and has been treated by many researchers using various techniques. Early numerical studies in this area using spectral methods focused primarily on shock-fitting techniques. However, for the cases in which a vortex causes a shock to deform at the point of bifurcation, it is generally difficult for a shock-fitting method to apply. This problem was successfully simulated in Refs.  by using the Chebyshev spectral method. In this paper, we also adopt this model to access the capability of the proposed spectral shock-capturing method. This problem is defined as follows. The initialization is imposed on the domain $`[0,2]\times [0,1]`$ with a stationary normal shock at $`x=0.5`$ and a Mach 1.1 flow at the inlet. The right state of the shock is given as $`(\rho _r,u_r,v_r,p_r)=(1,1.1\sqrt{\gamma },0,1)`$. A vortex is generated by a perturbation to the velocity $`(u,v)`$, temperature $`T`$ and entropy $`S`$ which is centered at $`(x_c,y_c)=(0.25,0.5)`$ and in the forms

$`u^{}`$ $`=`$ $`ϵ\tau e^{\alpha (1\tau ^2)}\mathrm{sin}\theta `$ (90)
$`v^{}`$ $`=`$ $`ϵ\tau e^{\alpha (1\tau ^2)}\mathrm{cos}\theta `$ (91)
$`T^{}`$ $`=`$ $`{\displaystyle \frac{(\gamma 1)ϵ^2e^{2\alpha (1\tau ^2)}}{4\alpha \gamma }}`$ (92)
$`S^{}`$ $`=`$ $`0`$ (93)
where, $`\tau =\frac{r}{r_c}`$, $`r=\sqrt{(xx_c)^2+(yy_c)^2}`$ and $`\theta =\text{tan}^1((yy_c)/(xx_c))`$. Here $`ϵ`$ denotes the strength of the vortex, $`\alpha `$ is the decay rate of the vortex. Note that we can obtain perturbations in $`\rho `$ and $`p`$ according to the relations of $`T=p/\rho `$ and $`S=\mathrm{ln}\frac{p}{\rho ^\gamma }`$. In our numerical experiment, we choose $`ϵ=0.3`$, $`r_c=0.05`$ and $`\alpha =0.204`$.

A reflective boundary is imposed on both the upper and lower boundaries. We solve the Euler equations on a 257 $`\times `$ 129 Cartesian grid, which is uniform in $`y`$-direction and refined in $`x`$-direction around the shock using the Roberts transformation . Once again, we emphasize that the DSC-RSK lowpass filter is used in the physical domain directly for this problem due to the presence of the reflective boundaries. FIGs. 15(a-e) show a time evolution of shock-vortex interaction in terms of the pressure contours. We can see that even at $`t`$=0.8 the solution (FIG. 15(e)) is almost free of oscillations and the fine scale features are captured very well. Particularly, the shock bifurcation reflected on the top boundary is shown clearly.

In order to explore the potential of our algorithm further, we increase the flow Mach number up to 3.0 and the strength of the vortex $`ϵ`$ up to 0.6. The results at $`t`$=0.2 are shown in FIG. 16. Deformation and bifurcation of the shock are observed clearly. These results indicate that the applicability of our spectral shock-capturing scheme is not limited to low Mach number and low vortex strength. Flows of high Mach number and high vortex strength can also be successfully simulated.

Example 11. In this case, we investigate the advection of an isentropic vortex in a free stream. As the exact solution of this problem is available, it is a good test for examining the accuracy and stability of shock capturing schemes.

Consider a mean flow of $`(\rho _{\mathrm{}},u_{\mathrm{}},v_{\mathrm{}},p_{\mathrm{}},T_{\mathrm{}})=(1,1,1,1,1)`$ with a periodic boundary condition in both directions. At $`t_0`$, the flow is perturbed by an isentropic vortex $`(u^{},v^{},T^{})`$ centered at $`(x_0,y_0)`$, having the forms of

$`u^{}`$ $`=`$ $`{\displaystyle \frac{\lambda }{2\pi }}(yy_0)e^{\eta (1s^2)},`$ (94)
$`v^{}`$ $`=`$ $`{\displaystyle \frac{\lambda }{2\pi }}(xx_0)e^{\eta (1s^2)},`$ (95)
$`T^{}`$ $`=`$ $`{\displaystyle \frac{(\gamma 1)\lambda ^2}{16\eta \gamma \pi ^2}}e^{2\eta (1s^2)}.`$ (96)
Here, $`s=\sqrt{(xx_0)^2+(yy_0)^2}`$ is the distance to the vortex center; $`\lambda `$ is the strength of the vortex and $`\eta `$ is a parameter determining the gradient of the solution. In our test, we choose them as $`\lambda =5`$ and $`\eta =1`$ unless we further specify. Note that for an isentropic flow, relations $`p=\rho ^\gamma `$ and $`T=p/\rho `$ are valid. Therefore, the perturbed $`\rho `$ is required to be

$`\rho =(T_{\mathrm{}}+T^{})^{1/(\gamma 1)}=\left[1{\displaystyle \frac{(\gamma 1)\lambda ^2}{16\eta \gamma \pi ^2}}e^{2\eta (1s^2)}\right]^{1/(\gamma 1)}.`$ (97)
The computational domain is taken to be $`[0,10]\times [0,10]`$ with the center of the vortex being initially located at $`(x_0,y_0)=(5,5)`$, the geometrical center of the computational domain.

Since there is no presence of shock in this problem, spectral method on its own can already provide excellent results if the time integration period is small enough. As such, the lowpass filter is not needed in an initial time period. However, as time progresses, errors would accumulate rapidly and the computation could become unstable if the lowpass filter is not activated to efficiently control the dramatically nonlinear growth of the errors. In our numerical experiment which is not shown here, the computation blows up at around $`t`$=13 if no filter is used.

Two experiments are conducted in this study. One is to examine the accuracy of our spectral code and the other is to investigate the stability of our scheme for a long time integration. For the first experiment, we compute the solution of $`\rho `$ at the time $`t`$=2 using three sets of meshes ($`N=32,64,128`$). To make the spatial error dominant, we optimize the CFL number to be 0.01.

Two error measures, $`L_1`$ and $`L_2`$, have been used. To be consistent with the literature , two errors used in this problem are defined as

$`L_1`$ $`=`$ $`{\displaystyle \frac{1}{(N+1)^2}}{\displaystyle \underset{i=0}{\overset{N}{}}}{\displaystyle \underset{j=0}{\overset{N}{}}}|f_{i,j}\overline{f}_{i,j}|`$ (98)
$`L_2`$ $`=`$ $`{\displaystyle \frac{1}{(N+1)}}\sqrt{{\displaystyle \underset{i=0}{\overset{N}{}}}{\displaystyle \underset{j=0}{\overset{N}{}}}|f_{i,j}\overline{f}_{i,j}|^2},`$ (99)
where $`f`$ is the numerical result and $`\overline{f}`$ the exact solution. Note that they are not the standard definitions. The errors for the density with respect to the exact solution are listed in Table I in which we also list the results of some other schemes for a comparison. As the table highlights, the accuracy of the Fourier spectral method is comparative to that of the DSC method , while both Fourier spectral and DSC methods are much more accurate than the fourth-order accurate, conservative central scheme (C4).

Our next numerical experiment concerns the performance of our spectral code for a long time integration. Here, we compute the solution of $`\rho `$ at the time $`t`$=100, 200, 400, 600, 800 and 1000, with the grid $`N=64`$. Note that we choose the gradient of the solution $`\eta =0.5`$ in this case. We list the $`L_1`$ and $`L_2`$ errors for the density in Table II. It can be seen that our spectral code is highly accurate even after such a long time evolution. The results indicate that our spectral scheme together with the DSC lowpass filter is reliable for the numerical approximation of conservation laws.

Example 12. Finally, we consider the problem of a supersonic flow past a cylinder . This example is used to mainly examine the ability of our DSC filter spectral scheme in handling a non-rectangular domain.

The physical domain is on the $`xy`$ plane, and the computational domain is chosen to be  $`\times `$  on the $`\xi \eta `$ plane. The mapping between the computational domain and the physical domain is

$`\begin{array}{cc}x=\hfill & \hfill (R_x(R_x1)\xi )\text{cos}(\theta (2\eta 1))\\ y=\hfill & \hfill (R_y(R_y1)\xi )\text{sin}(\theta (2\eta 1))\end{array}`$ (102)
where we take $`R_x=3`$, $`R_y=6`$ and $`\theta =\frac{5\pi }{12}`$.

The flow field is initialized with a Mach 3 shock moving from the inlet, i.e., $`\xi =0`$, toward the cylinder. The state in front of the shock is given as $`(\rho _r,u_r,v_r,p_r)=(1.4,1.0,0,1.0)`$. We impose a reflective boundary condition on the surface of the cylinder, i.e., $`\xi =1`$, and the outflow boundary condition at $`\eta =0`$ and $`\eta =1`$. In our simulation, a uniform mesh of $`65\times 129`$ in the computation domain is used. FIG. 17(a) depicts a diagram of the mesh (drawing every other grid line) in the physical domain. Like the case on a rectangular domain, we double the computational domain in both $`\xi `$ and $`\eta `$ directions when we approximate the first order derivative of the flux using the Fourier spectral method. The numerical result given in terms of pressure contours (obtained at $`t`$=4.5) is shown in FIG. 17(b). In summary, our spectral scheme can be easily applied to the nonrectangular domain, as long as the domain can be smoothly transformed to the rectangular one.

## IV Conclusion Remarks

In this work, we introduce the use of discrete singular convolution (DSC) filters in the Fourier pseudospectral method (FPM) for the simulation of hyperbolic conservation law systems. The fast Fourier transform (FFT) is used as the basic scheme, while DSC low-pass filters are used to suppress unphysical oscillations. The DSC filters are implemented in either the Fourier domain or the physical domain, depending on the physical problem under consideration. The Fourier domain implementation is easy and cost efficient, whereas the physical domain implementation allows more flexibility to handle some special boundary conditions. A sensor technique developed in the conjugate filter oscillation reduction (CFOR) scheme  is adopted in the present work to appropriately activate the DSC filter. The fourth-order Runge-Kutta scheme is utilized for the time integration.

Extensive numerical experiments are considered to validate the proposed approach and to demonstrate its usefulness. Excellent numerical results are obtained for all problems examined. The DSC filters turn out to be very effective in removing the Gibbs oscillations produced by the spectral approximation while maintaining a high accuracy. Moreover, it is also highlighted that these filters can dramatically improve the stability of the spectral method for the long time integration. The main features of proposed method are the follows:
(1) The basic scheme, the FPM is of spectral accuracy, which is important for problems that require not only capturing shock waves, but also resolving fine flow structures.
(2) The use of FFT endows the method high efficience, which is desirable for large scale problems in science and engineering.
(3) Numerical dissipation of the DSC filters is adjustable for a given problem, which endows the proposed method the ability to handle problems that are sensitive to numerical dissipation. The interaction of turbulence and shock is arguably a problem of this kind.
(4) The Fourier domain algorithm can be regarded as a windowed Fourier pseudospectral method for shock-capturing. The alternative use of the DSC filter in the physical domain makes the proposed method feasible to problems requires special boundary conditions, e.g., reflective boundary condition.
(5) The proposed scheme by-passes finding the characteristics which is required in many other shock-capturing schemes.
(6) The use of the method in higher dimensions is straitforward.

It is note that the state of art shock-capturing schemes, such as WENO and central schemes, do not depend on any adjustable parameter. To make the proposed scheme robust and user friendly, future work will focus on developing sophisticated sensors so that the filter parameter will be chosen automatically for each given problem. We believe that this can be accomplished by analyzing the Fourier power spectrum distribution of the approximate solution at selected time steps.

Acknowledgments This work was supported in part by Michigan State University.

APPENDIX
The optimal values ($`r`$) for DSC-RSK lowpass filters in our numerical experiments.


format: text

text:
# Uncovering Low-Dimensional Topological Structure in the QCD Vacuum

## Abstract

Recently, we have pointed out that sign-coherent 4-dimensional structures can not dominate topological charge fluctuations in QCD vacuum at all scales. Here we show that an enhanced lower-dimensional coherence is possible. In pure SU(3) lattice gauge theory we find that in a typical equilibrium configuration about 80% of space-time points are covered by two oppositely-charged connected structures built of elementary 3-dimensional coherent hypercubes. The hypercubes within the structure are connected through 2-dimensional common faces. We suggest that this coherence is a manifestation of a low-dimensional order present in the QCD vacuum. The use of a topological charge density associated with Ginsparg-Wilson fermions (“chiral smoothing”) is crucial for observing this structure.

An important aspect of the QCD vacuum structure relates to the local patterns in topological charge fluctuations. Phenomena such as a large $`\eta ^{}`$ mass, the $`\theta `$-dependence, and possibly also spontaneous chiral symmetry breaking (S$`\chi `$SB) are directly related to these fluctuations. Contrary to the situation in the past, the latest developments related to lattice chiral symmetry allow us to study these issues without a bias related to a particular model picture. This approach was initiated in Ref.  and is based on the use of the topological charge density $`q_x=\text{tr}\gamma _5(1\frac{1}{2}D_{x,x})`$ associated with Ginsparg-Wilson fermions . This density (a) satisfies the index theorem on the lattice , (b) its renormalization properties are analogous to those in the continuum , (c) exhibits smoother behavior than the underlying gauge fields due to non-ultralocality (“chiral smoothing”, and (d) can be naturally eigenmode-expanded .

Due to (d) it is possible to define an effective density  encoding the topological charge fluctuations up to low-energy scale $`\mathrm{\Lambda }`$, i.e. $`q_x^{(\mathrm{\Lambda })}_{|\lambda |\mathrm{\Lambda }a}(1\frac{\lambda }{2})c_x^\lambda ,`$ where $`c_x^\lambda =\psi _x^{\lambda +}\gamma _5\psi _x^\lambda `$ is the local chirality of the mode with eigenvalue $`\lambda `$. Using this framework it was demonstrated unambiguously that the topological charge fluctuations at low energy are not dominated by coherent unit lumps (e.g. instantons), and some interesting quantitative information about the typical amounts of charge in the local “hot spots” was obtained .

Motivated by (c), we now concentrate on the possibility that ordered behavior exists in the full $`q_x`$ containing all fluctuations up to the lattice cutoff. Indeed, if one is interested in the fundamental structure in the QCD vacuum relevant at all scales, then the full $`q_x`$ has to be examined. On the other hand, this strategy might appear rather futile since in the continuum  $`q(x)q(0)0,|x|>0`$, and thus a 4-d coherence associated with objects of finite physical size should not be prevalent in $`q_x`$ . However, such arguments do not exclude the dominance of an ordered coherent structure carried by lower-dimensional manifolds embedded in the 4-d Euclidean space. We demonstrate that a long-range coherence of this type exists.

If physically relevant 4-d coherent structures existed, they could be identified by finding all elementary sign-coherent 4-d hypercubes on the lattice and determining all connected regions formed by such cubes. Consistency of building a 4-d manifold requires that different hypercubes be connected through common 3-d faces. We have done such analysis using the overlap operator  and Wilson gauge action at two lattice spacings ($`a=0.110`$ fm and $`a=0.082`$ fm). The isolated islands of such coherence were found shrinking to mere points in the continuum limit as expected . Proceeding to the lower-dimensional manifolds, we have applied the same procedure using 3-d elementary cubes with structures defined by maximal connected regions (cubes connected through 2-d faces). For a given configuration, we have found all such structures and ordered them by the number of space-time points they contain. The typical (very stable) situation is shown in Fig.1 (left) where we plot the fraction of total number of space-time points taken by a given structure. We also show (right) what happens for the same configuration after performing a random permutation of sites.

We learn that: (i) $`q_x`$ is dominated by two oppositely-charged connected coherent structures, covering together about 80% of space-time points (even after the continuum limit extrapolation). (ii) The random permutation of sites destroys this long-distance coherence. We thus conclude that the observed behavior is a manifestation of an inherent low-dimensional order present in a typical snapshot of the QCD vacuum (equilibrium configuration) and absent in a completely disordered medium.

Finally, we wish to emphasize the following points. (A) From the lattice point of view, the two dominating structures represent 3-d hypersurfaces that are intertwined, embedded in the 4-d lattice, and almost filling the space-time. However, this is not sufficient to conclude that inside such sheets there are 3-d coherent regions of finite physical size surviving the continuum limit. The nature and dimensionality of such low-dimensional regions will be considered elsewhere . (B) The use of $`q_x`$ associated with GW fermions (chiral smoothing) is crucial for uncovering this coherent structure. Indeed, simple ultralocal (even improved) gauge operators do not appear to exhibit such behavior. (C) We speculate that the long-range nature of this low-dimensional topological structure might have crucial consequences for understanding the mechanism of S$`\chi `$SB in QCD .


text:
# Discovery of a High-Redshift (z=0.96) Cluster of Galaxies Using a FIRST Wide-Angle Tailed Radio SourceBased in part on observations obtained at the W. M. Keck Observatory.

## 1 Introduction

Clusters of galaxies trace the largest gravitationally-bound mass concentrations in the Universe; their epoch of formation, and thus their abundance as a function of redshift, is a strong constraint on $`\mathrm{\Omega }_M`$. High-redshift clusters are very difficult to find – the faint optical magnitudes of the associated galaxies and the X-ray emission from the intracluster medium (ICM) will often be missed by current surveys at redshifts of $`z1`$. Any new clusters revealed at these distances add significantly to the few that are known, and are useful for studies of cosmology and galaxy evolution.

Using a unique method based on radio galaxy morphology, we have confirmed the existence of dozens of clusters with redshifts up to $`z=0.85`$ (Blanton et al. 2000, 2001). Bent double-lobed radio sources (particularly FR I-type sources \[Fanaroff & Riley 1974\], including wide-angle tails \[WATs\] and narrow-angle tails \[NATs\]) are thought to achieve their distinct “swept-back”, C-shaped morphologies as a result of interaction with dense gas in which they are embedded. The canonical model describes the radio source’s host galaxy as having a significant peculiar velocity relative to the cluster with which it is associated – the ram pressure of the dense ICM then pushes back the lobes of the radio source giving it a “bent” appearance (Owen & Rudnick 1976). Often, the host galaxies of WAT radio sources are cD ellipticals, thought to reside at the bottom of cluster potential wells, and are thus not expected to have significant peculiar velocities. An alternate explanation (Burns et al. 1994), posits bulk motion of the ICM as a result of a recent or ongoing cluster-cluster (or sub-cluster) merger. The ram pressure of the merging gas then bends the lobes. Either way, bent-double radio sources are often pointers to clusters, and they may be flags in particular for clusters that have recently formed or merged.

We have selected 384 bent-double (mostly WAT) radio sources from the first 3000 deg2 of the VLA FIRST (Faint Images of the Radio Sky at Twenty-cm) survey (Becker, White, & Helfand 1995). Imaging and multi-slit mask spectroscopic follow-up on a moderate-to-high redshift subset of ten of these revealed eight clusters with redshifts as high as $`z=0.85`$ (Blanton et al. 2000). A complete (area- and magnitude-limited) low-redshift subset including 40 of the 384 objects was observed optically, and $`50\%`$ of them were found to be in richness class Abell 0 or greater clusters (Blanton et al. 2001). The majority of the low-redshift bent doubles that are found in clusters, based on our optical richness measurements, are detected in the ROSAT All Sky Survey (five out of six detected for $`z<0.2`$). Those that are not found in clusters, based on our richness measurements, are not detected in the RASS (only one out of nine detected for $`z<0.2`$). In other words, the combination of radio and optical data gives us a very clear indication for the presence, or lack of, a cluster, as confirmed with X-ray observations.

In order to identify high-$`z`$ candidates, we obtained $`R`$-band images of fields surrounding bent doubles that were blank to the limit of the Digitized Sky Survey (DSS). Some of these had no optical identifications even down to the limit of our images ($`m_R2223`$). We then performed near-infrared (NIR) observations at the positions of these sources in hopes of uncovering very red, high-$`z`$ cluster galaxies.

NIR observations have been successful in the past in identifying high-$`z`$ clusters. Stanford et al. (1997, 2001) found one of the most distant clusters known at $`z=1.27`$ using NIR and optical imaging, while Rosati et al. (1999) used a combination of X-ray, NIR, and optical observations to confirm the existence of a cluster at $`z=1.26`$ that may be part of a superstructure including the Stanford et al. (1997) cluster. Using the NIR alone, or in conjunction with radio and/or X-ray data, gives the potential for locating clusters that are at such high-$`z`$ they would be missed by all but the deepest optical observations. Recently, Gladders & Yee (2000) have performed a high-redshift cluster search by using optical and NIR photometry and identifying clusters by their “red sequence,” or clustering on a color-magnitude diagram. Chapman et al. (2000) have used a combination of radio observations with optical and NIR photometry to identify likely clusters at redshifts greater than 0.8.

Here we present optical and NIR observations taken at the position of the bent-double radio source, BD1137+3000. A combination of imaging and spectroscopy confirms that this source is located in a distant cluster of galaxies with a redshift $`z=0.96`$.

We assume $`H_{}=70`$ km s-1 Mpc-1, $`\mathrm{\Omega }_M=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ throughout.

## 2 BD1137+3000: Radio Data

The radio source is detected in the FIRST survey as a three component object, including a core and two lobes (FIRST J113733.6+300009, FIRST J113734.0+300002, and FIRST J113732.8+300003). The position of the core is $`\alpha `$ = 11h 37m 33.62s, $`\delta =+30^{}00\mathrm{}9\stackrel{}{\mathrm{.}}5`$. The source is clearly bent with an opening angle between the two lobes of $`80^{}`$. The size of the source, measured from the outer edge of one lobe to the core and out to the outer edge of the other lobe is $`30\mathrm{}`$ (240 kpc). The total flux density, as measured by FIRST is 9.1 mJy. The southeast lobe is the brightest of the three components, with a flux density of 4.1 mJy. The core has a comparable flux density (3.8 mJy), while the southwest lobe is 1.2 mJy. It is difficult to classify the source as an FR I or II based on its morphology. The southeast lobe is bright at the end of the lobe (typical of an FR II source). The southwest lobe brightens approximately halfway between the core and the edge of the lobe, bends at this bright area, and then fades toward the edge of the lobe (typical of an FR I source). We expect that there is additional extended emission related to the source that is not seen because (1) there is surface brightness dimming $`(1+z)^4`$ and (2) it is resolved out by the high resolution of the FIRST images. The flux density measured from the low-resolution NRAO VLA Sky Survey (NVSS; Condon et al. 1997), which resolves out less extended emission, 9.8 mJy, is slightly higher than the FIRST value. Using the NVSS flux density value, and assuming $`P\nu ^\alpha `$ with $`\alpha =0.8`$, the radio power of the source at 1440 MHz is $`P_{1440}=4.0\times 10^{25}`$ W Hz-1. This power is slightly below the FR I/II break for radio galaxies (Fanaroff & Riley 1974), and is typical of the class of bent double radio sources known as wide-angle tails (WATs; O’Donoghue et al. 1993) which are very often associated with cD galaxies in the centers of clusters. An overlay of the FIRST radio contours onto a Keck near-infrared camera (NIRC) image in the K filter of the cluster center is shown in Figure 1. The figure is 280 kpc on a side for our assumed cosmology. The radio contours are logarithmically spaced in the range 0.45 to 4.0 mJy. The radio core is coincident with an elliptical galaxy, as seen in the near-infrared image. This NIR image, which covers only the inner $`35\mathrm{}`$ square region of the cluster, is used for display purposes only in this paper. We include analysis of NIR data from the KPNO 2.1m telescope covering a much wider field of view in the following section.

## 3 Optical/NIR Imaging Data and Reductions

Optical data taken in the R-band at the KPNO 2.1m under partly cloudy conditions on UT 6 April 1997 revealed no identification for the radio galaxy to a limit of $`R22`$. This motivated us to observe the object in the NIR, and later to reobserve it to deeper limits in the optical under photometric conditions.

### 3.1  NIR

BD1137+30 was observed with IRIM at the KPNO 2.1m on UT 25 and 26 May 1997; the image scale was $`1\stackrel{}{\mathrm{.}}09`$/pixel. Observations were taken in a 16-point grid pattern. The total exposure in the $`K`$-band was 64 minutes (64 frames of 4$`\times `$15 s coadds), and that in the $`J`$-band was 32 minutes (16 frames of 4$`\times `$30 s coadds). The seeing was approximately $`1\stackrel{}{\mathrm{.}}6`$ for the $`J`$-band image and $`1\stackrel{}{\mathrm{.}}7`$ for the $`K`$-band image. A non-linearity correction was applied to the data to account for the changing response of the instrument when high fluxes are encountered. The correction was kindly provided by F. R. Chromey, derived from a non-linearity sequence taken the night immediately following our run, and amounted to a 3.3% correction at 30,000 counts. The data were corrected for the effects of dark current and flat fields using IRAF. The flat fields were constructed for each band by median combining the dithered target images and rejecting the objects contained therein. These types of flats were superior to dome flats (which we also constructed). The images were then processed with the IRAF contributed package DIMSUM111Deep Infrared Mosaicking Software, um… Created by P. Eisenhardt, M. Dickinson, S. A. Stanford, & J. Ward. We used DIMSUM to sky-subtract the data by using the sky derived from eight images taken immediately before and after each frame. The frames were then aligned and combined to produce a ’first pass’ combined image. During this first pass, object masks were constructed for the frames to help eliminate object residuals in the final combined image. A ’mask pass’ run incorporated these masks and produced the final combined image. We modified DIMSUM to attain better alignment of the images before combining them. The final mosaicked image reduced the image scale by a factor of four (to $`0\stackrel{}{\mathrm{.}}27`$/pixel – see Stanford, Eisenhardt, & Dickinson (1995)). During the reduction process, a few bad frames were identified and rejected, leaving us with a total of 30 min in the $`J`$-band and 56 min in the $`K`$-band. The final ’mask pass’ combined images were normalized so that the counts are counts per second. UKIRT faint standards were used to derive the photometry.

### 3.2 Optical

A total exposure in the $`R`$-band of 7200 s ($`6\times `$ 1200s) at the MDM 1.3m telescope with the Templeton CCD obtained on UT 21 May 1998 was kindly provided by J. Kemp. This setup produced an image scale of $`0\stackrel{}{\mathrm{.}}5`$/pixel. The night was photometric and the seeing was approximately $`1\stackrel{}{\mathrm{.}}21\stackrel{}{\mathrm{.}}3`$. The data were processed following standard IRAF procedures to remove the effects of the bias and flat field. The images were aligned and median-combined with cosmic rays rejected. The Landolt standard field PG1323 was used to calibrate the data.

## 4 Photometry / Image Analysis

FOCAS was used to measure magnitudes for the standard stars and objects in the target fields for the NIR and optical data. FOCAS ’total’ magnitudes were measured for the standards. The error for the optical photometry, estimated from the scatter in the offset between the measured and Landolt-catalogued magnitudes for five stars in PG1323, is 0.04 mag. The error in the NIR magnitudes, estimated from the difference in the zero point derived for the two nights, is 0.1 mag.

Aperture magnitudes were measured for objects in the target field. To assure that we were measuring the same physical aperture in each band, as well as to make the analysis more straightforward, adjustments were made to the final $`R`$, $`J`$, and $`K`$ images. They were all scaled to have the same pixel scale ($`0\stackrel{}{\mathrm{.}}5`$/pix) and the same orientation (N up, E to the left). The $`R`$ and $`J`$ images were smoothed with a Gaussian to match the seeing in the $`K`$-band image (the worst of the three, with 1$`\stackrel{}{\mathrm{.}}7`$). We then used FOCAS to measure magnitudes within an aperture with a diameter equal to 2$`\times `$FWHM (3$`\stackrel{}{\mathrm{.}}4`$ = 6.8 pixels). This is close to a total magnitude and will measure the same part of the galaxy in each of the bands. At $`z=0.96`$ (see §6), this aperture corresponds to 27 kpc for our assumed cosmology. The detection limit was 2.5$`\sigma `$ above sky and the minimum area was 12 pixels2 = 2.9 arcsec2 (the square of the FWHM). The sky was calculated by taking the average pixel value in a square annulus (default width 3 pixels) some buffer distance (default 2 pixels) away from the edge of the objects’ last isophotes. A small difference in the sky value had a very significant effect on faint objects. To assure that we weren’t including galaxy light in the sky subtraction, since many of the galaxies are very close together, we increased the buffer distance until the change in distance no longer affected the measured aperture magnitude (i.e., the sky value was no longer changing). We used a buffer of 15 pixels in $`R`$ and 20 pixels in $`J`$ and $`K`$.

The aperture magnitudes, measured as described above, for the radio source counterpart are 22.90$`\pm 0.04`$, 19.46$`\pm 0.1`$, and 17.78$`\pm 0.1`$, in $`R`$, $`J`$, and $`K`$, respectively. This gives colors of ($`RK`$) = 5.12$`\pm 0.11`$, ($`JK`$) = 1.68$`\pm 0.14`$, and ($`RJ`$) = 3.44$`\pm 0.11`$. The absolute magnitude of the host galaxy is $`M_V=22.29`$, assuming no evolution, a redshift $`z=0.96`$ (see §6), and a $`K`$-correction from Coleman, Wu, & Weedman (1980). This is typical of a bright elliptical galaxy, but fainter than a cD.

## 5 Color-Magnitude Diagrams

Color magnitude diagrams (CMDs) were constructed for objects in the entire 15.6 arcmin2 field as well as for only the objects found within a $`1\mathrm{}`$ radius aperture (480 kpc at $`z=0.96`$, see §6) from the radio source. No separation was made between stars and galaxies (all were included). At magnitudes close to that of the radio source, star/galaxy separation is very difficult, and not very important since number counts at these faint magnitudes are dominated by galaxies. The CMDs are presented in Figures 2 and 3. The top panels include the objects within $`1\mathrm{}`$ from the radio source, while the bottom panels include the objects from the entire field. The objects that are confirmed to be at the cluster redshift (see §6) are circled. Three obvious, foreground stars are seen in the diagrams – the three blue objects apparent in the CMDs constructed from the objects within one arcmin of the radio source. There is a much larger range in the colors of objects when looking at the CMDs for the whole field. The objects in the $`1\mathrm{}`$ radius aperture have less scatter (they are much more concentrated on the diagram) than those in the whole field, and are redder, consistent with membership in a distant cluster. Most of the objects within $`1\mathrm{}`$ from the radio source have colors in the range $`4<(RK)<6`$ and $`1.5<(JK)<2.5`$. Comparing the CMDs to those presented in Stanford, Eisenhardt, & Dickinson (1998), which includes diagrams for 19 clusters from $`z=0.30.9`$, and Stanford et al. (1997), which displays a CMD for a $`z=1.27`$ cluster, suggests that this cluster has a redshift between $`z=0.9`$ and 1.27, and closer to 0.9.

The $`R`$-band and $`K`$-band images are presented in Figure 4 and 5. Objects that were detected with at least 2.5$`\sigma `$ significance in both $`R`$ and $`K`$, with a minimum area of 12 pixels2 (2.9 arcsec2), and that have $`4<(\mathrm{R}\mathrm{K})<6`$ are circled. Spectroscopically confirmed cluster members (see §6) are labeled with letters and numbers; details for these objects are given in Table 1. The radio source host galaxy is marked ’A’. Source ’C’ was detected at only $`2\sigma `$ in the $`R`$-band, and was not detected in the $`J`$-band, but is circled because it has a measured redshift.

## 6 Spectroscopy

### 6.1 Long-slit Observation

A long-slit observation (2 $`\times `$ 30 min.) of the cluster that included the radio counterpart and two other galaxies was made at the Keck II 10m telescope with the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) on 13 April 1998. The spectra were corrected for bias and flat field effects using standard IRAF procedures. They were wavelength calibrated using comparison Ne + Ar spectra and flux calibrated using a standard star (Feige 34). Redshifts were derived from cross-correlating the spectra with an elliptical galaxy template (M32 shifted to the rest frame).

### 6.2 Multi-slit Mask Observation

A multi-slit mask observation of the field surrounding the radio source was made at the Keck II 10m telescope with LRIS on 29 May 2001. The field was observed for a total exposure of 60 min. (3 $`\times `$ 20 min.). We reduced the data from the LRIS-R, which was used with a 150 lines mm-1 grating. The central wavelength was $`8495\mathrm{\AA }`$, the dispersion was $`4.8\mathrm{\AA }`$/pix, and the pixel scale was 0.2109 arcsec/pix. Objects were chosen for the slit mask based on their $`(RK)`$ colors. Objects that had $`4<(RK)<6`$ were candidates for inclusion on the mask, since they provide us the best opportunity of confirming galaxies at redshifts $`z1`$. The mask included a total of twelve program object slits.

The data were reduced following standard IRAF procedures. The frames were corrected for bias using the overscan region. An internal halogen flat taken through the mask following the data frames was used to correct for flat field effects. The three object frames were then aligned and combined with cosmic rays rejected.

Spectra were successfully extracted for seven of the program objects; the remaining five were too faint for successful extractions. The spectra were wavelength-calibrated using comparison spectra of Hg + Ar lamps. Redshifts were initially estimated by identifying a few obvious lines, including the Ca II H+K break, and the \[O II\] $`\lambda 3727`$ and \[O III\] $`\lambda \lambda 4959,5007`$ emission lines. The spectra were then Fourier cross-correlated with an elliptical galaxy template (M32 shifted to the rest frame) using FXCOR in IRAF. The errors in the redshifts are computed in FXCOR and are based on the fitted peak height and the antisymmetric noise, or “r-value” (Tonry & Davis 1979). The average (mean) error for the redshifts from both the long-slit and multi-slit mask observations is 150 km s-1.

### 6.3 Spectroscopic Results

All ten spectra that were successfully extracted had redshifts in the range $`0.950<z<0.965`$, giving us a $`100\%`$ success rate of identifying cluster galaxies based on their NIR/optical colors. The objects and their redshifts are listed in Table 1. The spectrum of the radio host galaxy is shown in Figure 6, and is typical of an elliptical galaxy. The spectrum has been smoothed with an 11 pixel boxcar. We calculated the line-of-sight velocity dispersion using $`\sigma _{}=\sqrt{(N1)^1_{i=1}^N\mathrm{\Delta }v_i^2}`$, where $`\mathrm{\Delta }v_i=c(z_i\overline{z})/(1+\overline{z})`$. The 68% confidence uncertainty is calculated following Danese, De Zotti, and di Tullio (1980), and includes errors due to the small number of member galaxies in the field and the measurement errors. We find $`\sigma _{}=530_{90}^{+190}`$ km s-1. This is a typical value for an Abell richness class 0 cluster (Lubin & Bahcall 1993).

## 7 Limit on the X-ray Luminosity

An observation in the ROSAT  public data archive serendipitously includes the position of BD1137+3000. The ROSAT  PSPC observation was centered at $`\alpha `$ = 11h36m33.6s, $`\delta =+29^{}48\mathrm{}00\stackrel{}{\mathrm{.}}0`$ leaving the position of the bent-double radio source 17$`\stackrel{}{\mathrm{.}}`$8 from the center, very close to the inner support structure of the PSPC. The observation was performed on 1991 May 30, for a total duration of 33447 seconds.

We processed the data using the Snowden routines222“Cookbook for Analysis Procedures for ROSAT XRT/PSPC Observations of Extended Objects and the Diffuse Background”, S. L. Snowden (1994), with support from the USRSDC (Snowden et al. 1994) which correct for exposure, vignetting, variations in the detector quantum efficiency and non-cosmic background components. After applying these routines, 17884 seconds worth of good data remained. We extracted the X-ray count rate within a $`1\mathrm{}`$ (480 kpc) radius aperture centered on the position of the radio source, and took background from two large circular apertures with similar offsets from the center of the PSPC field as the source region. The source region contained counts above the background at the $`2.1\sigma `$ level, with a $`90\%`$ confidence upper limit on the count rate of $`1.1\times 10^3`$ ct s-1 in the $`0.42.0`$ keV energy band. We converted this to a bolometric X-ray luminosity with PIMMS, using a model including Galactic absorption of $`N_H=1.93\times 10^{20}`$ erg cm-2 (Dickey & Lockman 1990), a Raymond-Smith plasma with a temperature of 2.4 keV (estimated from the $`\sigma T`$ relation \[Lubin & Bahcall (1993)\]), a chemical abundance of 0.4 times the solar value, and a redshift of $`z=0.96`$. The upper limit on the bolometric X-ray luminosity is then $`L_{X,bol}3\times 10^{44}`$ erg s-1, which is consistent with the values expected from the $`L_x\sigma `$ relation ($`3\times 10^{44}`$ erg s-1 \[Edge & Stewart 1991\]; $`5\times 10^{43}`$ erg s-1 \[Borgani et al. 1999\]). Chandra and XMM-Newton observations are planned to confirm the presence of an X-ray emitting ICM in this cluster, and to study it in more detail.

## 8 Conclusions

We have provided evidence at optical and near-infrared wavelengths for a high-redshift cluster of galaxies associated with a wide-angle tailed radio galaxy at $`z=0.96`$. Color-magnitude diagrams reveal an over-density of red objects near the radio source with colors consistent with elliptical galaxies at this redshift. Objects are clustered on the color-magnitude diagrams in the range $`4<(RK)<6`$ and $`1.5<(JK)<2.5`$.

Using spectroscopy from the Keck II LRIS, we have confirmed a total of ten galaxies at the cluster redshift. The line-of-sight velocity dispersion is $`\sigma _{}=530_{90}^{+190}`$ km s-1, typical of an Abell richness class 0 cluster. Using data from the ROSAT public archive, we limited the cluster’s bolometric X-ray luminosity to $`L_{X,bol}3\times 10^{44}`$ erg s-1, which is consistent with the value expected from the $`L_X\sigma `$ relation.

The majority of high-redshift clusters found in traditional X-ray and optical surveys select the most X-ray luminous and optically rich clusters. The radio selection technique employed here is an important complement to those methods, in part because it can locate high-z clusters with a wide range of X-ray luminosities and optical richnesses, thus aiding in the study of galaxy evolution and cluster formation.

E. L. B. thanks Adam Stanford and Craig Sarazin for helpful discussions regarding the NIR and X-ray data reduction, respectively. Support for E. L. B. was provided by NASA through the Chandra Fellowship Program, grant award number PF1-20017, under NASA contract number NAS8-39073. Part of the work reported here was done at the Institute of Geophysics and Planetary Physics, under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under contract W-7405-Eng-48. Support for M. D. G. was provided by a grant from the National Science Foundation (AST-99-70884). The FIRST project is supported by grants from the National Geographic Society, the National Science Foundation (AST 00-98259 and AST 00-98355), NASA (NAG 5-6035), NATO, IGPP, Columbia University, and Sun Microsystems.