{'fraction_non_alphanumeric': 0.02714447663433682, 'fraction_numerical': 0.009635974304068522, 'mean_word_length': 4.615738284703802, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Toward a Consistent Framework for High Order Mesh Refinement Schemes in Numerical Relativity

It has now become customary in the field of numerical relativity to couple high order finite difference schemes to mesh refinement algorithms. To this end, different modifications to the standard Berger-Oliger adaptive mesh refinement algorithm have been proposed. In this work we present a fourth order stable mesh refinement scheme with sub-cycling in time for numerical relativity. We do not use buffer zones to deal with refinement boundaries but explicitly specify boundary data for refined grids. We argue that the incompatibility of the standard mesh refinement algorithm with higher order Runge Kutta methods is a manifestation of order reduction phenomena, caused by inconsistent application of boundary data in the refined grids. Our scheme also addresses the problem of spurious reflections that are generated when propagating waves cross mesh refinement boundaries. We introduce a transition zone on refined levels within which the phase velocity of propagating modes is allowed to decelerate in order to smoothly match the phase velocity of coarser grids. We apply the method to test problems involving propagating waves and show a significant reduction in spurious reflections.


Introduction
Long term stable evolution of non linear hyperbolic partial differential equations often require techniques to efficiently deal with vast length scales while simultaneously resolving fine scale features. Indeed, many numerical simulations in computational astrophysics, cosmology, numerical relativity and fluid dynamics are confronted with processes that span a wide range of time and length scales. In the context of numerical relativity, these simulations are often performed in full three spatial dimensions without symmetry assumptions. For these cases, running fine unigrid integrations is computationally expensive and often impractical. Recent establishments in the numerical solution of partial differential equations by finite difference techniques has seen an increasing use of nested grids and mesh refinement techniques in order to minimize the truncation error incurred with minimal computational and memory requirements [5,4,28].
The principle behind mesh refinement schemes is to recursively refine areas of the computational grid that are likely to induce higher discretization errors. This approach efficiently focuses computational effort and resources in places where it is needed compared to refining the entire grid. Extensive theory has since been developed for the method in different contexts, and mesh refinement algorithms have been widely adopted in the literature [20,23,29,12,34]. One of the key aspects in the implementation is the way in which the decision to add or remove levels of refinement is made. In this work we focus on the concept of fixed mesh refinement (FMR), where the grid hierarchy is created once and remains fixed for the duration of the computation [6,29,10]. This differs from adaptive mesh refinement (AMR), where the algorithm is endowed with error estimation routines that dynamically determines which areas need refinement. For most applications in numerical relativity, one can know before hand which areas within the computational domain require more refinement. As a result the FMR configuration is widely adopted in the field. In the context of numerical relativity, mesh refinement has seen application ranging from studies in critical phenomena, single black hole space times, binary black hole collisions, black hole lattice universes, neutron stars and core-collapse supernovae, within the non linear framework of general relativty, see for example [19,29,7] and references therein.
Another component of the scheme is the inter-level coupling among nested grids. While the coarse grid has to supply boundary conditions to finer grids during evolution, one can choose not to update the coarse grid solution with the fine grid solution. In this case, the coarse grid solution is independent of that of the nested grid, a configuration that is refered to as one way (parasitic) nesting. In this work we employ two way (interactive) nesting, where we, in addition, update coarser levels with finer levels once the finer levels have been integrated to the same time level as the coarser ones. See [17] for a comparison between parasitic and interactive coupling within the context of the shallow water equations.
Traditionally mesh refinement techniques were coupled to second order convergent methods. On the other hand, recent trends in the numerical simulation community has seen the coupling of higher order finite difference methods to the mesh refinement framework [30,19,11,36,25]. This combines the efficiency of local mesh refinement with the robustness and accuracy of higher order methods. However, there is an inherent incompatibility between high order time discretization schemes with the standard mesh refinement algorithm that may result in loss of convergence or even instabilities [19]. This issue is related to how the computation of boundary data for the refined grids is handled. A search for a stable high order mesh refinement implementation has resulted in several modifications to the standard method in an effort to address this subject. For example, [29,7] employs the idea of buffer zones, where additional points are added closer to the boundary of the refined grids 1 . In this setup, boundary conditions are not prescribed explicitly in the refined levels, the integration is only applied to a progressively smaller domain in the refined grids and the buffer zone is ultimately discarded. Another approach is the tapered boundary approach [19,11]. Here, one performs the integration at level l using the past domain of dependence of the child grid only. Other approaches have been to refine only in space and using the same time step for all levels [10]. In this work we use a framework where we refine both in space and time and the treatment of interface boundaries is dictated by the time marching algorithm, fourth order accurate Runge Kutta algorithm in this case.
In addition to issues of convergence and stability, one has to address the problem of spurious reflections off refinement boundaries that arise when propagating waves cross refinement boundaries. This is essential for gravitational wave source simulations as the waves are normally extracted at a large radius. As a result, propagating waves will have crossed several refinement boundaries, before reaching the radius of extraction. In [10], the idea of derivative matching was proposed in order to minimize spurious reflections for second order convergent schemes. In addition, the concept of mesh adapted stencils (MAD) was introduced in [2]. However these implementations do not involve refinement in time. Other methods that have been applied in the Advanced Weather Research and Forecasting Model, Euler equations and Maxwell equations involve the use of sponge layers in the refined levels to ensure that the solution in the refined levels will be nudged towards that of the coarser grids at refinement boundaries. This may involve the addition of artificial damping and dissipation terms in the system under consideration [17,32]. See also the treatment of [27] in the case of first order convergent schemes. In this work we propose a simple scheme that is adopted from the animation and image processing community [14], to deal with transitions from fine to coarse grid solutions at refinement boundaries.
Finally, we note that high order mesh refinement implementations are often endowed with several lower order approximations in order to reduce the overall computational demand [7,25,36]. In the case where there is refinement in time, this often includes the use of lower order interpolations in time. This results in some savings in memory usage as the storage of past time levels is lower for lower order interpolations. For example, one typically needs n + 1 past points in order to have a n-th order polynomial interpolant. In addition, one can reduce the number of buffer zones in the refinement boundaries by successively lowering the order of finite differencing as one approaches the refinement boundary. Other approaches in the case where artificial dissipation is needed, include the use of Kreiss-Oliger dissipation operators of order lower than that required by the finite differencing operator. This lowers the number of required ghost points as well as being cheaper computationally. Ultimately, these approximations have a bearing on the convergence order, the amount of spuri-ous reflections induced at interface boundaries and the overall performance of a given implementation. In order to keep the overall scheme consistent, we make no such simplifying assumptions in this work. This paper is organized as follows: in Section 2 we convey the framework of our FMR approach. We review our boundary application method in §2.4.1 and introduce the transition zone in §2.4.2. We present our results in Section 4 and finally concluding in Section 5.

Grid layout
The grid hierarchy is arranged by first discretizing the spatial domain into a relatively coarse uniform mesh that covers the entire computational domain. This constitutes the base or root grid H 0 0 with mesh size h 0 . Finer grid patches of mesh sizes h l /r are then overlaid as required to the base grid with each grid at level l having mesh size h l−1 /r. More than one grid patch can be added in a given level. This forms a tree or a hierarchy of grids H l p , where the indices l and p represent the level and patch number respectively. This configuration is depicted in Figure 1. It is at this point that we emphasize that each grid in a given level l has its own solution vector and is evolved independently of all other grids. Of course it has to depend on the parent grid, within which it is nested, for boundary data.
Each grid patch added at a given level must satisfy certain conditions. Among them, the idea of proper nesting: A fine grid at level l must start and end at the corner of a cell belonging to level l − 1. Moreover, grids at higher levels cannot 'float'. This means that if there is a grid at level l + 2, it must be contained in a grid at level l + 1 that is itself properly nested on a grid at level l. The refinement factor r must be an integer, and is the same for all levels. This results in a constant CFL for all added levels, thus the same integration routine of the base level is stable on all levels, if it is stable on level 0. It also implies that grids at higher levels require r l time steps to catch up with a single time step of the base grid.

Inter-level communication
Each grid in a given level l can be indexed independently by its own (i l , j l , k l ) coordinate system. However, for reasons of inter-level communications, there is a mapping from the H l p coordinates (i l , j l , k l ) to the H l−1 p coordinates (i l−1 , j l−1 , k l−1 ) and viceversa. This relation is expressed as, for a staggered hierarchy. Such communications are necessary for the computation of initial conditions, boundary conditions and for updating the coarse grid with the fine Fig. 1 A grid hierarchy demonstrating proper nesting. A single grid H 0 0 covering the entire domain is marked by 'a'. There is one refinement grid H 1 0 at level one marked by 'b'. Two disjoint refinement grids are H 2 0 and H 2 1 are marked by 'c' and 'd' respectively. Note that the ghost zones are not included in the grids. This figure is used to emphasize that the grid hierarchy is not some complex data structure, but that the overlaid grids are independently stored in memory. grid solution. Initial data can be generated by spatial interpolation from the previous grid level (prolongation) or by calling the same initialization routine that was used to initialize the base grid. All levels are added and initialized at the same initial time t 0 . Once all levels have been integrated to the same time, data in the finer meshes is used to update data in the coarse levels through the use of interpolation operators, a process called restriction.
It is important to note that with the mapping of indices 1, all fine grid points are staggered about coarse grid points at lower levels. This has the advantage that if the base grid is discretized strategically to 'avoid' certain points (x, y, z) ∈ R 3 , (e.g. by using a cell centered grid to avoid dealing explicitly with the points on the edges) such points remain avoided in all refined levels. However, from a computational standpoint, this may be expensive because communicating data across levels requires three dimensional interpolation all the time. This is different from cases where some points in the fine grids are allowed to coincide with the ones from coarser levels. For such cases, inter-level communication can occur via injection, where data is simply copied to corresponding points at a given level.

Numerical Integration
Spatial differentiation is handled through fourth order convergent finite differencing. The first derivative is given by the operator, while the second derivative is given by, The cross derivative is given by sequential application of 2. In the case of the xy derivative, this results in, We use these centered stencils for all derivatives, except advection terms for which we use the lop sided formulas, Time integration is carried out using a fourth order accurate Runge Kutta scheme through the method of lines framework. A Runge Kutta step from y n to y n+1 is accomplished with, where the quantities Y i defined by, serve to store boundary data. The control algorithm that is responsible for the evolution of the entire grid hierarchy is orchestrated by a recursive procedure, shown in Algorithm 1 Algorithm 1: A simple illustration of the AMR integration algorithm with refinement factor r.

Boundary Treatment
The base grid consist of a ghost zone that serve to store boundary data. For refined grids, in addition to ghost zones, there is also a transition zone that ensures a smooth transition from the well resolved solution of the finer grids to the less resolved solution of the coarser grids. We discuss each in turn.

Ghost zone
For a fourth order accurate computation of centered derivatives, we use two ghost zones on each side of a grid 2 . Outer boundary conditions on the base grid are supplied by the user, for example, the user may specify periodic or sommerfeld type boundary conditions. We distinguish between two types of boundary conditions for the finer levels, those that coincide with the outer boundary and those that simply border a cell from the underlying coarse grid, which we refer to as refinement boundaries. If the ghost zones of refined grids coincide with those of the base grid, the boundary is filled using the prescribed procedure for outer boundaries. Otherwise, we use coarse grid data to fill the fine grid boundaries. Imposing boundary data on refinement boundaries can be reduced to the problem of imposing time dependent boundary conditions within a method of lines framework. An interesting peculiarity of higher order Runge Kutta schemes is that under certain conditions, they have a tendency to behave like lower order schemes, a phenomenon which is refered to as order reduction. This phenomenon is caused by the stiffness of the system under consideration in the case of implicit schemes [26,15], or inconsistent application of boundary data for the intermediate Runge Kutta stages in the case of explicit schemes [8,24]. Of course, in both cases the problem is exacerbated by the non linerity of the problem. Although not identified as such, order reduction effects in the context of mesh refinement were reported in [29,19] when the conventional method of imposing boundary conditions was used, i.e, simply applying boundaries corresponding to the solution evaluated at the intermediate times of the Runge Kutta stages. In this work, we instead use the Runge Kutta method itself to fill the ghost zones. We use Equations 9, where the k i are given by, [22] k 1 = hy , In the equations above, y , y and y are time derivatives of the quantities under evolution while f y is the Jacobian matrix of the PDE system. See [22,21] for details on the derivation. The time derivatives can easily be obtained by polynomial interpolation methods since the coarse grid points at the advanced time will already have been computed before advancing the refined levels. However, a subtle issue arises in this case. To evolve the finer grids, at least four past points of the coarse grid solution are needed in order to obtain third order interpolants. This means each finer grid can only be initialed after the coarser grid has evolved at least four time steps. This is undesirable in the context of FMR. We opt to use the fact that the classical Runge Kutta method has a built-in interpolant, termed dense output. This interpolant is given by, [16] y(t n + θh) = y n + with 0 ≤ θ ≤ 1 and the b i are polynomials in θ, with 1 ≤ m ≤ 3. One does not need to compute the Jacobian matrix f y explicitly since we are only interested in the product f y y which can be computed from the system 10 as, See also, [21,9,30]. The implication is that we do not store the solution history of coarser grids as is currently done in some implementations. Instead, we store the four intermediate k i values corresponding to the current time, for all levels. This configuration allows for time interpolation through (11). Of course this is followed by spatial interpolation, for which we employ fourth order barycentric Lagrange interpolation; higher than fourth order was found to be unreliable in some of the runs.

Transition zone
To complete the specifications on treatment of the boundary, we examine what happens close to the refinement boundary. Consider a grid hierarchy with two levels l 0 and l 1 . Parametrize the solution F (x) on such a hierarchy as, where f (x, l 0 ) and f (x, l 1 ) are the solutions on the base and refined grids respectively, and w is a binary weight function which takes the value w = 1 if x is within the refined region and w = 0 everywhere else. A plot of the weight function is depicted in Figure 2.
Note the discontinuity at the transition points x = 10 where w transitions from w = 0 to w = 1 indicating a switch from the solution F (x) = f (x, l 0 ) to F (x) = f (x, l 1 ). Also at x = 90, w transitions from w = 1 to w = 0, indicating the switch from F (x) = f (x, l 1 ) back to F (x) = f (x, l 0 ). Because of the dispersion relation for propagating waves ( §7.1), there is a difference in phase speeds of propagating modes in the coarse and fine grid levels. As a result of the discontinuous transition in the weight function w, waves propagating from refined regions abruptly change their phase velocities when crossing refinement boundaries, creating a glitch that will seed spurious reflections. To circumvent this problem, we introduce a transition zone on the refined levels, within which the weight function w = w(x) is allowed to vary smoothly from w = 0 to w = 1 across the refinement boundary. This can be accomplished by Hermite interpolation. For a transition beginning at x = a and ending at x = b, one can derive the following profiles, where the variable t is defined as, These profiles are shown in Figure 3. In this case the weight function varies continuously from w = 0 to w = 1, allowing the solution F (x) to vary smoothly across the transition zone. This can also be interpreted as a smooth acceleration and deceleration of the associated phase speeds of propagating modes in the solution F (x). See [14] for a discussion on transition profiles.
The procedure we follow to ensure a smooth transition from the fine grid solution to the coarse grid solution is as follows.
-Fill ghost zone points using the method outlined in §2.4.1.
-Fill the transition zone by blending values from the refined grid, with that of the coarse grid according to the weight function w(x).
The transition zone is evolved along with the fine grid solution to ensure a smooth coupling with the refinement boundary and thus the coarse grid solution. However, for reasons of stability, we do not use transition zone values when updating the coarse grid solution with the fine grid solution. Unless otherwise specified, we take the width of the transition zone to be three through out this work. For this size, the smooth profiles given above are equivalent.

Evolution system
We adopt the BSSN formulation of the Einstein field equations [31,3]. The evolution equations are given in terms of the variables, subject to the constraints,γ = Det γ ij = 1,Ã i j = 0. Additional variablesΓ i = −γ ij ,j are also introduced. The evolution equations for these variables are derived from the ADM equations and are given by, The superscript [· · · ] T F denotes the trace-free part with respect to the physical metric γ ij , and The Ricci tensor R ij is now written as a sum of two pieces where R φ ij is given by

Numerical results
For ease of exposition, we restrict to the case β i = 0. Where refinement is used, we restrict to a refinement factor r = 2. All runs employ three ghost points and, where appropriate, three transition points.

Wave equation: Periodic boundaries
In this section, we carry out evolutions of the wave equation with mesh refinement. The wave equation in flat Cartesian coordinates is given by, We instead cast it in an alternative form, by introducing a new auxiliary variable Π = ∂ t φ, As initial data, we choose a sinusoidal profile with periodic boundary conditions on the domain x ∈ [−0.5, 0.5]. The region x ∈ [−0.25, 0.25] is refined by a factor of r = 2. Although the wave propagates essentially in one dimension, we evolve it using the full 3D grid with periodic boundary conditions for the outer boundaries. In Figure 4 we plot the solution errors for evolutions with resolution dx = 1/25ρ for ρ = 1, 2, 3. The errors show fourth order convergence as desired.

Gauge Wave
We now evolve a non linear gauge wave using the BSSN system. The Gauge wave test is characterized by a line element which results from a non linear gauge transformation of the flat Minkowski space time in Cartesian coordinates, resulting in where the function H is given as, for some constant A, and d is the wavelength. We evolve the above metric using the BSSN formulation with the Harmonic gauge condition, We choose the amplitude A = 0.1 and the wavelength d = 1. As in the last case, the simulation domain covers the range x ∈ [−0.5 : 0.5] with refinement boundaries in the region x ∈ [−0.25 : 0.25]. Evolving the gauge wave initial data with the BSSN formulation requires the addition of artificial dissipation to achieve stable evolutions. This is true even for unigrid runs, see for example [1]. Our dissipation operator takes the form for rth order accurate finite difference stencils. In Figure 5 we plot the solution errors for evolutions with resolution dx = 1/25ρ for ρ = 1, 2, 3, 4. The errors scale according to fourth order convergence. The non linearity of this problem makes it ideal to demonstrate the effects of order reduction. Figure 6 shows a similar case, but with boundary data of the refined levels imposed in the conventional way. Notice the loss of convergence at high resolution. The loss of convergence at high resolutions was also identified in [29].

Wave equation: Gaussian pulse
In the following test we show how our proposed algorithm handles artificial reflections that often arise when a propagating wave crosses a mesh refinement boundary. We are interested in waves propagating outward from the fine grid across mesh refinement boundaries into the coarser grid. We evolve the wave equation, with initial data given by a Gaussian pulse centered at the origin, with σ = 0.25 and A = 1. Although the wave propagates essentially in one dimension, we evolve it using the full 3D grid with periodic boundary conditions for the outer boundaries. The simulation domain covers x ∈ [−4, 4]. The region x ∈ [−1, 1] is further refined by a factor of r = 2. The solution is shown in Figure 7(a). The pulse starts initially at x = 0 with amplitude one and produces two pulses each with amplitude 0.5 traveling in opposite directions. In Figure 7(b) we show the result after the pulses have crossed refinement boundaries at x = ±1. When each pulse crosses a refinement boundary, spurious reflections are generated. These travel in a direction opposite that of the inducing pulse. When no transition zone is used, the spurious reflections reinforce at x = 0 and can exceed the discretization error in amplitude. Employing a transition zone significantly reduces these artificial reflections.

Teukolsky Wave
In this section we evolve the Einstein field equations in three space dimensions. As initial data, we use the Teukolsky solution for a quadrupole l = 2, even parity m = 0 waves [33]. The metric for the quadrupole modes is given by, The coefficients A, B and C are constructed via a generating function F (x) and are given by, We take F (x) to be a superposition of ingoing (x = t + r) and outgoing (x = t − r) waves, and, where we have chosen the particular case F 1 (x) = −F 2 (x) = Ae −x 2 . The angular functions f uv for even parity m = 0 modes are given by We note that this represents time symmetric data and so K ij = 0 and K = γ ij K ij = 0 at the initial slice t = 0. We use an amplitude of A = 10 −6 to complete the specification of initial data.
Because of the symmetries of the problem, we impose mirror symmetry boundary conditions along the planes x = 0, y = 0 and z = 0. We thus evolve the Octant propagates radially outward crossing mesh refinement boundaries along the x = 4, y = 4 and z = 4 planes, eventually reaching the radiation boundary and leaving flat Minkowski spacetime. Although the Teukolsky wave is a routine problem for testing numerical relativistic codes, it is especially challenging for a mesh refinement code. This is mainly because the refined region is Cartesian, while the wave propagates spherically outward. As a result, the wavefront will not encounter the refinement boundaries at the same time. In Figures 8 we show the result at t = 8 for a run without a transition zone. Spurious ripples are generated when the wave initially hits the refinement boundary; the reflections continue to be generated until the wave has fully crossed the refinement boundary. These ripples are reflected toward the origin as expected. In Figure 9, we show a similar run using a transition zone. In this case, spurious reflections are significantly minimized.

Concluding Remarks
In this work, we have presented a fourth order convergent mesh refinement scheme that also significantly minimizes spurious reflections off refinement boundaries that are caused by differing levels of accuracy between two successive refinement levels. This is an important issue for the field of numerical relativity where the use of higher order finite differencing is becoming increasingly common [25,13,35]. For these higher order methods, the truncation error can become so small that the dominant error comes from spurious reflections. Our method is not restricted to any formulation of the Einstein field equations. Indeed one can apply it to any hyperbolic system of partial differential equations. Because we are not using the buffer zone approach as detailed in [29,7], our method requires a total of six points on each boundary, irrespective of the time integration method used. This differs from employing buffer zones in that, for a fourth order accurate Runge Kutta algorithm, along with fourth order finite differencing involving lop-sided advection stencils, twelve points are needed along each boundary (the situation could be worse for higher order finite differencing) [18,7]. This is a significant saving in memory usage, especially in three space dimensions where the buffer zone can be a significant part of the grid. In addition, the blending operation we employ to fill the transition zone is cheaper than having to repopulate the entire buffer zone after every time step. We also note that the use of a transition zone is computationally cheaper than the sponge boundary method since one has to populate the sponge boundary at every intermediate Runge Kutta step, while the transition zone is only populated at the end of the time step. Moreover, there is often a level of experimentation required to determine how large a sponge zone one should use. Although one can have the transition zone as large as desired, we have found satisfactory results with a size that spans only three fine grid points.
The implementation described here uses Runge Kutta dense output formulas to interpolate in time, which avoids any potential issues with polynomial interpolation. In particular, because one does not have to couple fine grid solutions to the solution history of coarser grids, fine grids can be immediately initialized along with the base grid. This is especially attractive since fine grids can be initialized to the same accuracy as the base grid, by using the same initialization routine as the base grid. In addition, the use of high order polynomial interpolation in time requires the use of asymmetric stencils, which results in oscillations in the refinement boundaries and ultimately, unstable modes at higher resolutions [7]. This dense output method of time interpolation, which is routinely used in the solution of ordinary differential equations [16], was also used in the context of mesh refinement [21] for conservation laws and [9] for solving Maxwell's equations.
It would be interesting to investigate the efficiency of the transition zone implementation in an adaptive context or whether it would minimize spurious reflections that are a result of shock waves crossing refinement boundaries. There is also the question of how it would affect conservation along interface boundaries in the context of conservation laws. These and other issues involving black hole spacetimes are a subject of further study.

Acknowledgments
I thank Denis Pollney for discussions and comments. BM acknowledges financial support from the National Research Foundation (South Africa).

Dispersion relation
In this section we calculate the dispersion relation and phase velocity resulting in discretizing the wave equation with a fourth order stencil along with Runge Kutta time marching. We follow the approach (and notation) of [10] where a similar calculation was given using second order finite differences and iterative Crank-Nicholson time marching. The wave equation (31) can be written in matrix form as where we have defined the vector V as With this identification, we denote by V n j the solution at time step n and grid point j. The second derivative operator ∂ xx appearing in (53) is given by the stencil, For our analysis, we consider plane wave solutions for V , of the form, for some constant vector W . Using the classical fourth order Runge Kutta scheme to advance (53) in time, results in the update rule, Plugging in (56) to the above rule results in the relation, where we have defined Λ as, The system (58) represents an eigenvalue problem. In particular, W is an eigenvector corresponding to the eigenvalue e iωdt for the matrix in (58). Further analysis shows the eigenvalues to be the pair, This expression represents the dispersion relation, relating the frequency ω with the wave number k. For completeness we calculate the phase velocity, v p (λ) = ξ/k for ξ = Re (ω). From (60) we get, Therefore, v p (λ) = λ 2πdt arcsin 2Λ(3 − 2Λ 2 ) 9 − 4Λ 6 (2 − Λ 2 ) . (62)