{'fraction_non_alphanumeric': 0.048700498101297035, 'fraction_numerical': 0.014351235389850569, 'mean_word_length': 4.121874210659257, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fokas's Uniform Transform Method for linear systems

We demonstrate the use of the Unified Transform Method or Method of Fokas for boundary value problems for systems of constant-coefficient linear partial differential equations. We discuss how the apparent branch singularities typically appearing in the global relation are removable, allowing the method to proceed, in essence, as for scalar problems. We illustrate the use of the method with boundary value problems for the Klein-Gordon equation and the linearized Fitzhugh-Nagumo system. The case of wave equations is treated separately in an appendix.


Introduction
The Unified Transform Method (UTM) or Method of Fokas presents a new approach to the solution of boundary value problems (BVPs) for integrable nonlinear equations [6,8], or even more successfully for linear, constant coefficients partial differential equations (PDEs) [5,8]. In this latter context, the method has resulted in new results for one-dimensional problems involving more than two spatial derivatives [7,8], for elliptic problems [3,8], and for interface problems [4,11], to name but a few. A more complete picture can be obtained from the online repository at [1].
Fewer studies exist of how the method applies to linear systems of equations or to higherorder scalar evolution equations that may be rewritten as such. In fact, to our knowledge, the application of the UTM to the standard wave equation presented in the appendix here is new, albeit not surprising. Systems of equations are examined in [9,13,15,14] and [16]. In [16], a problem from thermoelastic deformation is examined, while in [9] two water waverelated systems are treated. In [15], the UTM is applied to the wave equation in a moving domain, while in [13,14] the free Klein-Gordon equation is examined. The treatment in these papers is non-generic in the sense that the different branches of the dispersion relations of the systems considered are either polynomial or rational in the wave number (as in [9,15]) or are conveniently parameterized to avoid radicals, as in [13,14,16]. In addition, all five papers start from a Lax pair and construct a Riemann-Hilbert problem almost exclusively, with almost no hints as to how to generalize the more accessible UTM for scalar systems based on the use of Green's Theorem. Generically, the different branches of the dispersion relations for a first-order evolution system of dimension N depend on radicals as they are the roots of an N -th order polynomial. Below we investigate the general case of an N -dimensional first-order linear system, using two examples.
The first example is the (free: no potential) Klein-Gordon (KG) equation [10,13,15] u tt − u xx + αu = 0, where α is a constant parameter and indices denote partial derivatives. For our purposes, the KG equation is rewritten as a two-dimensional first-order system: where q = u, p = u t . Our approach to the KG equation differs from [13,15] in that all our considerations are based on so-called local relations, to which Green's Theorem is applied. No parameterization of the dispersion relation branches is used. Rather we use the branches of the dispersion relation is their original form. Lastly, in [13,14] only α > 0 is considered. The second example is the linearized Fitzhugh-Nagumo (FN) system of partial differential equations [12] v where β is a constant parameter. In general, we consider systems of the form where Q is an N -dimensional vector and Λ is a linear-operator matrix of size N × N and of order n.
Example 1. For the KG equation, we have so that n = 2, and Example 2. Similarly, for the FN system, allowing us to read off Λ(−i∂ x ) and n = 2 as well.
Throughout, we contrast the general systems case with the general scalar case, where the equation can be written in the form for a scalar-valued function u(x, t).
Although different nonlocal systems can be considered by allowing Λ to depend rationally on its argument (see e.g. [17] for the scalar setting), in this paper we restrict to the case where Λ (and λ) depends polynomially on its argument. As a consequence, the dispersion relation for the scalar case, easily found by equating u = exp(ikx−ωt), is given by ω = λ(k), and ω depends polynomially on k. It should be remarked that in the above calculation of the dispersion relation, we have followed the convention in the literature on the UTM for the dispersion relation. Thus the dispersion relation as used here differs by a factor of i from the standard use. For instance, in the UTM dispersive equations are characterized by a purely imaginary dispersion relation for k ∈ R.
Similarly, for the systems case (3), we let so that ω satisfies where I is the N × N identity matrix. Thus the different branches Ω 1 , . . . , Ω N of the dispersion relation are roots of an N -th order polynomial. Generically, they depend on radicals of order up to N (N -th roots) whose arguments are polynomials in k. It follows that these dispersion branches are sheets of an N -valued function and they have branch point singularities in the complex k plane. This is not always the case, as is illustrated in the appendix for the wave equation and in the examples treated in [9]. Throughout this paper, we assume that all branches Ω j (k) are distinct, except at isolated values of k ∈ C.
Example 1. For the KG equation, which has branch points at the square roots of −α. We choose the branch cut that connects these branch points.
Example 2. For the FN system, which has 4 branch points.
One of the main advantages of the UTM is its ability to characterize the number and the type of boundary conditions required to ensure wellposedness of a given initial-boundary value problem [8]. We wish to see to what extent the same can be done for systems of PDEs. To this end, we do not specify boundary conditions for our examples at this point. Rather, we will see how the application of Fokas's UTM determines the information that should be provided on the boundary of our domain. We limit ourselves to problems posed on the half line x > 0. It is anticipated that extending our results to problems on the finite interval x ∈ (0, L) is comparable to extending the UTM for scalar problems on the half line to scalar problems on the finite interval. Of course, in all cases, we specify initial conditions Q(x, 0) = Q 0 (x).
In the next sections, we go through the extension of Fokas's UTM, as applied to systems of linear, constant-coefficient evolution PDEs. Each section deals with a different step of the method, so as to present the method in an algorithmic way.

The local relation
The first step for the application of the UTM is to rewrite the system of equations in divergence form. We refer to this form as the local relation. Following the appendix of [9], we rewrite (3) as where Ω(k) = diag(Ω 1 , . . . , Ω N ), the diagonal matrix with the different branches of the dispersion relation as diagonal elements. The matrix A(k) diagonalizes the matrix Λ(k): Lastly, the vector X(x, t, k) is a differential matrix operator of degree at most n − 1, polynomial in k, defined by and the last equality defines the matrix-valued polynomials c j (k). Equation (12) is verified by working out the product rules of both terms. For the scalar case, both equations are identical. For the case of systems, an important difference between (12) and (13) is that (13) contains Ω(k) and A(k), which are typically branched in the complex k plane. On the other hand, the left-hand side of (12) depends only on Λ(k) and X(x, t, k), which are not branched as functions of k. In practice, (13) is more useful, as (12) requires the calculation of a matrix exponential with non-diagonal exponent. Even in specific examples where the matrix exponential is easily calculated, it is useful to have the local relation in terms of the branches of the dispersion relation.
Example 1. For the KG equation, exp(Λ(k)t) is easily computed directly: which is not branched since all functions above with square root arguments are even. With the local relations (13) are Example 2. For the FN equation, and the local relations (13) are after multiplication by Ω j (k), and using that Ω 1 (k)Ω 2 (k) = β. Remarks.
• For the scalar case, there is no difference between the two forms (12) and (13). The local relation is easily obtained by multiplying the PDE (7) by exp(−ikx + ωt) and using integration by parts to get to the divergence form. The dispersion relation is found during this process [5] as well. The same process can be used for systems of PDEs, but additional linear algebra is required to get to the divergence form, as only specific linear combinations of the equations allow for this form.
• In practice, it may not always be possible to write the set of local relations as compactly as above, using index notation: the first local relation may depend on Ω 2 (k) and so on. In fact, this is the case in the second example, but the use of the dispersion relation for the FN system allows the further simplification.
• In [16] and to a lesser extent in [9], the Lax pair formalism for the UTM is used. For constant-coefficient systems linear PDEs this is not necessary, and we stay within the framework of using the local relation and Green's Theorem (in the next step).

The global relation
For the next step, we integrate each local relation over an infinite strip in the (x, t) plane, cornered at the origin, see Fig. 1. Using Green's Theorem, we obtain using (12), and using (13). Here Note that we have replaced the upper limit T (see Fig. 1) of the temporal variable by t. The relations (21) and (22) are referred to as the Global Relations. The same comments can be made as for the local relations: (21) is written entirely in terms of quantities without branch points, while (22) contains branched quantities throughout, but is easier to write down in practice. For specific examples, it is more convenient to start directly from the derived local relations.
Example 1. For the KG equation, wherê for f = q, p, κ = 0, 1 and the index κx denotes κ derivatives with respect to x.
Example 2. For the FN system, where the hat denotes the half-line Fourier transform, as above, and e Ω j s v κx (0, s)ds, κ = 0, 1.
As in the case of constant-coefficient scalar evolution equations posed on the half line, the global relation is valid in the lower half of the complex plane, Im k ≤ 0. Indeed, the halfline Fourier transforms in (23a) are defined only if the exponential does not grow, assuming integrable or square-integrable initial conditions and solutions. The other quantities are defined by proper integrals and do not impose any further restriction.

A "solution" formula
The global relations contain the Fourier transform of the quantities we wish to solve for. At this point, we solve for these transforms and invert using a regular Fourier inversion, interpreting the half-line transforms as whole-line transforms of functions that are zero for x < 0. Thus using (21), and, clearly equivalent, These equations do not provide expressions for the solutions of (3). Although the first term contains only known quantities, the second term depends on m j − 1 derivatives of Q j evaluated at the boundary x = 0. Here Q j is the j-th component of Q, j = 1, . . . , N and m j is the maximal degree as a function of k occurring in the j-th column of Λ(k). If m j = 0, the function Q j nor any of its x derivatives appear in the second term and no boundary data involving Q j should be prescribed. We know from experience with scalar problems [5,8] that it is unlikely that all these unknown boundary functions need to be specified. In fact, specifying them freely is likely to result in inconsistencies: e.g., it is not possible to prescribe the Dirichlet and Neumann data for the heat equation on the half line independently.
using the same notation as above.
In the following sections, we aim to turn (28) and (29) into genuine solution formulas, depending on the correct number of boundary functions.

Deformation of the integration path
Following the approach to solve constant-coefficient scalar equations, we wish to deform the integration path in the expressions (28) and (29) as far away from the real line as possible as k → ∞ [5,8]. In some cases all or part of the real line will remain in place, and other curves in the upper-half plane resulting in zero contributions might be added to the integration path, as in the scalar case. As in that case, we have no intent to deform the path of integration of the first term: the integrand of this term is known explicitly and it requires no further manipulation.
Thus we focus on the second term. Specifically, we wish to elucidate the role of the branch points in the systems case. Since (28) is written in terms of quantities that are not branched, it is clear that the branch points of A(k) and Ω(k) that are apparent in (29) are all removable. This is also immediately clear from the expressions (30a-b) and (31a-b). Encircling any of the branch points results in an interchange of the two sheets, effectively switching the indices on Ω 1 and Ω 2 . This permutation leaves the integrands invariant, confirming the removability of their branch points. As stated above, this is true in general, due to the non-branched nature of the integrands in (28). For any given example, the "solution formula" is conveniently written using (29), resulting in integrands that are symmetric functions of the different branches of the dispersion relation. Encircling a branch point amounts to a permutation of the indices of these branches, producing no change due to the symmetry of the integrands. It follows that the process of deforming the path of integration proceeds very much as in the scalar case, with minor modifications due to the vector structure as discussed below.
We examine the k dependence of the integrand in the second term of (28). Since x and t are independent variables, both defined on their respective half lines, different restrictions on k ∈ C are imposed by the need to control the x and t dependence of the integrand separately. Since the x dependence is confined to the exponential, it follows that necessarily k ∈ C + = {C : Im k ≥ 0} for the integral to be defined. Similarly, using (23b) and following the reasoning from the scalar case [5,8], we need Re Ω(k) ≥ 0, i.e., the matrix Ω(k) is positive definite. Thus, we define the inaccessible region D + in the upper-half plane: Further, for k ∈ C + /D + , the integrand decays exponentially as k → ∞. Thus the contribution to the second term of any path that tends to infinity in C + /D + vanishes, by Jordan's Lemma [2]. It follows from Cauchy's Theorem [2] that we may write using (21), and, clearly equivalent, The integrand of the second term is a linear combination of different exponentials exp (ikx − Ω j t). One may consider splitting the integral into different parts, each depending on one of these exponentials only. If this is done, the different integrands are not symmetric functions of the branched quantities, and the branch points are not removable. As a consequence, much greater care is required for the deformation of the integration path away from the real line. To facilitate this, it is convenient to deform the integration path above all finite and branch points is trivially allowed by Cauchy's Theorem. Deforming around branch points at infinity should be considered on a case-by-case basis, as convergence issues come into play. We will not consider this case further.
After the initial finite deformation, further deformations into the region not containing the branch points can be made without these points playing any role at all. Once the branch points are below the path of integration, it may be possible and desirable that different integration paths ∂D + j are used for different components Q j , if these components do not depend on all Ω j .
Example 1. For the KG equation q and p depend on both Ω 1 and Ω 2 . In general, D + is conveniently found by first determining its boundary, followed by examining the different regions it defines. For the KG equation, (9) reads Splitting this expression in its real and imaginary parts and imposing that Re ω = 0 determines the boundary of D + . This boundary is shown in red in Figure 3. This boundary does not separate the upper-half plane in different regions. Since Ω 1 = −Ω 2 , there is no subset of the upper-half plane where the real part of both Ω 1 and Ω 2 is positive. Thus D + consists of the whole upper-half plane and no deformation from the integration along the real line to the upper-half plane in (30a-b) is done.
Remark. Alternatively, one may deform the integration to go above the branchpoint at i √ α (α > 0) or those at ± √ −α (α < 0). After doing so, the integral on the second term may be distributed, resulting in one integral where the whole upper-half plane is accessible. By Cauchy's Theorem, the contribution from this integral vanishes. In other words, the net contribution from this branch of the dispersion relation is hiding in the other integral as the contribution from the branch point(s) as the other integral is deformed back to the real line.
Example 2. For the FN equation the situation is more straightforward. Figure 3: The region D+, with boundary shown in red, covering the whole upper-half plane, for α > 0 (left) and α < 0 (right). Figure 4: The region D + for the FN equation and the deformed path of integration along its boundary in red.
which is illustrated in Fig. 4. Depending on the value of β, some of the different branch points of Ω 1 and Ω 2 may lie between the real line and ∂D + . As long as the integrals are not distributed over their integrands, this is not a concern. Thus v(

Symmetries of the dispersion relation
In order to eliminate the superfluous dependence on unknown boundary conditions, we need the symmetry group of the dispersion relation. This is the collection of transformations k → ν(k) that leave the dispersion relation invariant. Applying these transformations allows us to transform one branch of the dispersion relation into another one. The symmetries are found by solving det(Λ(ν) − Ω j (k)I) = 0, j = 1, . . . , N.
which determines ν as a function of k, since Ω j depends on k through (9). It is clear that the identity transformation ν 1 = k is a solution of (38). Below we show that for both the KG and FN system, the only nontrivial symmetry is ν 2 : k → −k. The appendix contains an example of a different symmetry.

The elimination of unknown boundary functions
The deformed solution formula (34) depends on the quantities where Q j,lx denotes l derivatives with respect to x of Q j (x, t). Through these quantities, the solution formula exhibits dependence on a number of boundary functions Q j,lx . In order to eliminate such unnecessary boundary functions from the deformed solution formulae (34), new global relations (22) are obtained by using the symmetries from the previous section.
Since the matrix exp(Ωs) is diagonal, the m-th component of (22) contains only g (m) j,k , for varying j and k. In other words, Ω m is the only branch of the dispersion relation appearing in the m-th component. If Ω m (ν l (k)) = Ω n (k) ∈ {Ω j , j = 1, . . . , N }, the new global relation below is considered: withG m (ν l (k)) = t 0 e Ωn(k)s (A(ν l (k))X(0, s, ν l (k))Q(0, s)) m ds, where the subindex m indicates the m-th component is taken. It follows that these new global relations depend on the same quantities g (m) j,l as does the deformed solution formula (34). The new global relation is valid for all k for which ν l (k) ∈ {k ∈ C : Imk ≤ 0}. In other words, k is in the closed lower-half plane transformed under ν −1 l . This process results in a number of new global relations. If these global relations are valid in regions of C where all or a part ∂D + lies, they can be used to eliminate unwanted boundary functions by solving the new global relations for these unwanted boundary conditions, as in the scalar case [5,8]. The above procedure determines the exact number of boundary conditions that needs to be prescribed in order for the boundary-value problem to be well posed: this number is N j=1 (m j + 1) minus the number of boundary functions that can be eliminated.
Example 1. For the KG equation, there is only the nontrivial symmetry ν 2 (k) = −k. Both Ω 1 (k) and Ω 2 (k) are invariant under k → −k, thus we obtain two additional global relations, both valid in the closed upper-half plane. Importantly, both are valid on the real line, thus both can be used in (30a). The two new global relations are The solution formula (30a) depends on g 0 , g 1 , and g 1 . Using the two new global relations above, we expect to eliminate two of these. Assuming that Dirichlet boundary conditions are specified, we wish to eliminate the dependence in (30a) on g (1) 1 and g (2) 1 , which encode the Neumann data. Note that no boundary-value dependence on p(x, t) arises. Using the expressions above, (30a) becomes Here the contributions in the right-hand side of (45a) that exhibit dependence on the left-hand side.
We have not used that Ω 1 = −Ω 2 or that Ω 2 1 = Ω 2 2 = −(α + k 2 ), so as to exhibit the symmetry of the solution. It is a straightforward check from the above formulae that q t = p. Example 2. As above, the FN equation only has the nontrivial symmetry ν 2 (k) = −k.
Both Ω 1 (k) and Ω 2 (k) are invariant under k → −k, thus two additional global relations are obtained, both valid in the closed upper-half plane, which is where ∂D + is. Thus both can be used in (37a). The two new global relations, solved for g The solution formula (31a) depends on g (1) 0 , g 0 , g 1 , and g (2) 1 . Using the two new global relations above, we expect to eliminate two of these. Assuming that Neumann boundary conditions are specified, we wish to eliminate the dependence in (30a) on g (1) 0 and g (2) 0 , which encode the Dirichlet data. No boundary-value dependence on w(x, t) arises. Using the expressions above, 1 − Ω 1 e −Ω 1 t g (1) 1 where In (48a), the path of integration along ∂D + for the terms involving initial conditions has been deformed back to the real line, to allow their combination with the initial-condition terms already present in (37a).
In the above examples R 1 , R 2 , S 1 and S 2 represent the right-hand terms depending on the solution we wish to obtain. In the next section, we show these terms are zero.

A solution formula
The new global relations (42) depend on the solution throughQ j (ν l (k)). Thus solving these relations for boundary functions and substitution in the solution formula (34) introducesQ j (ν l (k)) in the right-hand side of (34). It appears no effective solution formula has been obtained yet. At this point, we single out the right-hand terms depending on the solution, and we wish to show that their contribution is zero, as in the examples below. In general, in these culprit terms the time-dependent part of the exponential cancels, allowing the deformation of the contour into the previously inaccessible region. Further, for all examples we have examined, the right-hand side terms involvingQ j (ν l (k)) are symmetric under a permutation of the indices of the different branches Ω j of the dispersion relation, see Example 1 and 2, above. As a consequence, the integral over ∂D + may be distributed to separate off these terms without introducing branching. This is to be expected: the solution of the original boundary-value problem should be independent of the choice of indices on {Ω j , j = 1, . . . , N }: switching our labels on Ω 1 and Ω 2 , for instance, should not affect the solution. If the solution is to be symmetric under these permutations, we should expect the same for the terms containingQ j (ν l (k)).
Example 1. We show that R 1 = 0 = R 2 in (45a). Recall that the original solution formulae (30a) could not be deformed in the upper-half plane due to the presence of either e −Ω 1 t or e −Ω 2 t . These exponentials canceled, and they are absent in (46). On the other hand, the factor e ikx decays in the upper half plane and it follows from Jordan's Lemma that both R 1 and R 2 are zero. Thus (45a) with R 1 = 0 = R 2 represent the final form of the solution of the KG equation posed on the positive half line with Dirichlet boundary data. The solution formula shows that no boundary information on p(x, t) is required.
Example 2. The fact that S 1 = 0 = S 2 follows in exactly the same way: the region D + was previously inaccessible, but because the time-dependent exponential is absent, we may apply Cauchy's Theorem around D + , with the integral contribution from the path at infinity vanishing, due to Jordan's Lemma. Thus S 1 = 0 = S 2 since no singularities are present in D + .

Remarks.
• For these particular example, one can observe that R 1 = 0 = R 2 , S 1 = 0 = S 2 in another way. Consider R 1 . Replacing k → −k in the integral, we get which is the inverse Fourier transform of the Fourier transform of q(x, t), but evaluated at −x. Since x is positive, it follows that R 1 = 0. A similar argument works for R 2 , demonstrating that R 1 = 0 = R 2 without contour deformation. For S 1 and S 2 , we deform back to the real line, after which the above argument can be repeated. This method works for these examples because the lone nontrivial symmetry ν 2 (k) = −k indicates the presence of a mirror symmetry in (1) and (3a-b): the substitution x → −x leaves both equations invariant and one expects that the half-line boundary-value problems for these equations may be solved using the method of images.
• In the first example, we eliminated g (1 1 and g (2) 1 , the time transforms of the Neumann data. Similary, in the second example we eliminated g (1 0 and g (2) 0 . One may solve the first new global relations for g (1) 0 and the second one for g (2) 1 , for instance. When one does so, the terms depending on the Fourier transformsQ j (ν l (k)) of the solution evaluated at −k are not symmetric under permutation of the indices, as we expect. Nonetheless, as the reader easily verifies for the KG equation with g (1) 0 and g (2) 0 eliminated (i.e., g (2) 0 and g (1) 0 are assumed known), the contributions from the term containingq(−k, t) still vanish, due to Jordan's Lemma. Thus the resulting solution formula for q(x, t) and p(x, t) is not symmetric under a permutation of the dispersion branch indices. This may be explained by noting that specifying boundary conditions by supplying g (2) 0 and g (1) 0 breaks the symmetry at the level of the problem statement already: when one specifies g (2) 0 = g 0 (Ω 2 , t) = t 0 e Ω 2 s q(0, s)ds, one needs to include the information whether this function of t originates from integration involving Ω 1 or Ω 2 .

A three-dimensional example due to Fokas and Treharne
In [16], Treharne and Fokas study a three-dimensional linear system on the half line, originating from a problem in elasticity. We discuss how their systems fits within the algorithm outlined in this paper. Although [16] is a seminal paper, as the first work where the UTM is considered for a system, the arguments presented there are specific to that system (the same can be said for the considerations in [9]) and do not highlight generic features of the UTM as applied to systems. In particular, we have emphasized the role played by expressions involving symmetric functions of the roots of the dispersion relations, allowing us to bypass the need for any branch cuts. Doing so, all expressions involved in the formal manipulations are analytic and much of the flavor of the scalar case is retained. Indeed the solution itself should be invariant under a relabeling of the branches of the dispersion relation, assuming the boundary conditions preserve this indifference.
The Treharne-Fokas system has the dispersion relation Here α, β and λ are real parameters. This dispersion relation is cubic in ω and quartic in k with a branch point at infinity in the k variable. Other than by using Cardano's formulae, it is not possible to obtain the branches Ω j (k) explicitly. The use of Cardano's formulae results, as expected, in horrendous expressions whose manipulation is not helpful. However, by a straightforward asymptotic analysis of the dispersion relation, one can identify the regions in the complex k plane for ω with simple asymptotic behavior as k → ∞, allowing us to label the different branches Ω 1 , Ω 2 , Ω 3 . From the knowledge of the branches one readily obtains the associated eigenvectors (in terms of Ω j , j = 1, 2, 3) that diagonalize the relevant operator. Much of the analysis presented above follows immediately and is not presented. Note that knowledge of the asymptotic behavior of Ω j , j = 1, 2, 3 is sufficient to deform the paths of integration. At the next step one requires the symmetries that transmute one branch of the dispersion relation into another. Perhaps surprisingly, these symmetries are easily obtained in terms of ω from the dispersion relation: since the dispersion relation is a bi-quadratic in k, there are four symmetries: . With the symmetries in hand, the arguments presented earlier follows and a solution expression is obtained. Of course, one does not require every combination of ω j and ν j and a judicious choice should be made to eliminate as many boundary values as possible. It is here that the peculiarities of the PDE problem at hand come into play.
The solution formula obtained as described above depends on the expressions for Ω 1 , Ω 2 and Ω 3 : in effect we have reduced the solution of a PDE boundary-value problem to that of solving a third-order polynomial. For a general system, the solution expression depends on the roots Ω j (k) and on the symmetries ν j (k). Neither one can be written down explicitly for general dispersion relations. The complexity of solving boundary-value problems for linear systems of constant-coefficient PDEs is thus reduced to that of solving one additional polynomial problem compared to the scalar case. For scalar problems, Ω(k) is given explicitly, but the determination of the symmetries requires the solution of a polynomial of degree one less than the order of the problem (after the trivial symmetry ν 1 (k) = k has been divided out). For systems, one faces the additional task of obtaining the branches of the dispersion relation Ω j (k).
The analysis in [16] differs from ours in that the authors of [16] parameterize the dispersion relation, thereby diagonalizing the operator and eliminating any need to analyze expressions which may be branched. In our work, by employing symmetric functions of the branches of the dispersion relation, one yet again avoids the detailed analysis of branched expressions. The benefit of parametrization is that expressions are obtained that are more explicit. Of course, the fact that this problem may be suitably parametrized is not generic. Only dispersion relations given by an algebraic curve of genus zero can be parameterized. The genus of the dispersion relation and its corresponding parametrization may be readily obtained using the symbolic algebra package Maple and its package algcurves.

Conclusions
We have demonstrated how Fokas's Unified Transform Method can be generalized from scalar evolution PDEs with constant coefficients to systems of such equations. Our main goal has been to show that the method continues its applicability much as in the scalar case, with added complexity relative to the system under consideration. It is difficult to make general statements for such a large class of problems, but some trends are clear. Although the different branches of the dispersion relation typically contain radicals, the different stages of the solution process can be executed without the need for branch cuts to be introduced, as the functions that arise are symmetric functions of the dispersion relation branches. Further, as in the scalar case, the symmetries of the dispersion relation give rise to new global relations, allowing for the elimination of unknown boundary data. Alternatively, this process allows one to determine the amount of boundary information that needs to be supplied to have a well-posed problem.
Here a ∈ R. Note that the wave equation is recovered if a = 0. We follow the algorithm for the UTM for systems and we rewrite this equation as This leads to and Note that α 1 > 0, α 2 < 0. For the wave equation, α 1,2 = ±1. Since this implies that one of the characteristic speeds is positive (information is carried away from the boundary), and one is negative (information is carried towards the boundary), we expect that one boundary condition needs to be prescribed. With the local relations are This leads to the global relations for j = 1, 2, valid in Im k ≤ 0. Herê e Ω j s u(0, s)ds, g e Ω j s u t (0, s)ds.
It should be noted that specifying the Dirichlet boundary condition u(0, t) results in both g with 1 + ikg We do not write down a solution formula for v(x, t) = u t (x, t). The first term in (58) is known and no manipulation of it is necessary. The second term brings in dependence from both the Dirichlet and Neumann boundary data and we expect to be able to eliminate one of these. We wish to deform the path of integration used for I 2 as far as possible from the real line in the k plane. As for the KG equation, D + = {k ∈ C : Im k > 0}, thus no deformation away from the real line is allowed.
Next, we examine the symmetries of the dispersion relation. We solve where ω = Ω 1 or ω = Ω 2 . We obtain These symmetries are the reason we proceed with Example (50) instead of the regular wave equation: these symmetries are reminiscent of those used to solve scalar interface problems, see [4], for instance. More importantly, they go beyond k → −k, which is used in the main sections of the paper.
Since Ω 1 (ν 3 (k)) = Ω 2 (k) and Ω 2 (ν 2 (k)) = Ω 1 (k), we obtain two new global relations, both valid on the real line: the first one is obtained by substituting k → ν 3 (k) in (56) with j = 1, while the second is obtained by letting k → ν 2 (k) in (56) with j = 2. These new global relations are and both are valid for Im k ≥ 0, since α 1 /α 2 < 0. Next, we wish to eliminate either the Dirichlet or Neumann data from (59b). If Dirichlet data is given, we wish to use (62a-b) to eliminate the Neumann data functions g (j) 1 , j = 1, 2. On the other hand, if Neumann boundary conditions are provided, we use (62a-b) to eliminate g (j) 0 , j = 1, 2. As in the scalar case [5,8], if a Robin boundary condition is given, all of g (j) m , j = 1, 2, m = 0, 1 are eliminated using (62a-b) and the two time transforms of the Robin condition. For exposition sake, let us assume Dirichlet data is given. We substitute g (j) 1 , j = 1, 2 obtained from (62a-b) in (59b). This results in where Note that T 1 (x, t) is determined by the initial data. It can be combined with I 1 (x, t). The function T 2 (x, t) depends only on the Dirichlet data, assumed known. Lastly, T 3 (x, t) exhibits dependence on the solution, but we show below that T 3 (x, t) ≡ 0. Thus is a solution of (50), posed on x > 0 with Dirichlet data specified at x = 0. This solution is valid for all x and t in the quarter plane x > 0, t > 0. To show that T 3 (x, t) ≡ 0, it suffices to note that due to the absence of the exponential, deformation of the path of integration into D + (the upper-half plane) is allowed, and the conclusion follows from Jordan's Lemma and Cauchy's Theorem. Observe that the singularity at k = 0 is removable.
Lastly, we restrict to the case of the wave equation (a = 0, Ω 1,2 = ±ik) and we illustrate how the solution obtained here can be reduced to d'Alembert form. We examine the Dirichlet, Neumann, and Robin problem. The characteristics for the wave equation are drawn in Fig. 5. One family of characteristics moves initial data away from the line t = 0, to the left. The second family takes information to the right, either initial data coming from t = 0 or else boundary data from x = 0. It follows that points in the blue region in Fig. 5 depend on initial data only. Their solution expression should be the same, independent of which boundary conditions are imposed. Points above the blue region depend on both initial data and boundary data. In this region we expect a different solution expression, depending on which boundary conditions are specified.
Examining I 1 (x, t) with a = 0, we find x t 0 Figure 5: In the blue region, the solution depends only on initial data. Above this region, the solution depends on both initial and boundary data.
using well-known calculational properties of the Fourier transform. This is the d'Alembert form of the solution, as expected in the blue region x > t. This implies that we should find I 2 (x, t) = 0 in the blue region. Indeed, with a = 0 and x > t, 1 + ikg 1 + ikg Both x − t + s > 0 and x + t + s > 0 if x > t. Thus the integrands are decaying exponentially in the upper-half k plane and by Cauchy's Theorem and Jordan's Lemma, I 2 (x, t) ≡ 0. Note that the singularity at k = 0 in the above integrals is removable. Thus in x > t, This formula may be extended to include x = t, provided the initial and boundary data are compatible at (x, t) = (0, 0). Next, we examine x < t.
For the Dirichlet problem, I 2 (x, t) = T 1 (x, t) + T 2 (x, t). With a = 0, Next, again with a = 0, where we have used the definition of the delta function and the fact that the Dirichlet boundary condition u(0, t) = 0 for t < 0. Using that u 0 (x) = 0, v 0 (x) = 0 for x < 0 in (66), we finally obtain It is a straightforward check that this solution satisfies the wave equation, the initial conditions, and the Dirichlet boundary condition. Further, if initial and boundary data are compatible at (x, t) = (0, 0) then the solution is continuous at x = t.
For the Neumann problem, I 2 (x, t) = T 1 (x, t) + T 2 (x, t). With a = 0 and x < t, the new global relations (62a-b) are ikg (1) 1 , which are substituted in (59b), resulting in We have omitted the terms containingû(−k, t) andv(−k, t), which are easily shown to have zero contributions. Here the Heaviside function. As before, we have used that u 0 (x) = 0, v 0 (x) = 0 for x < 0. Combining our results, we obtain As for the Dirichlet problem, it is straightforward to check that this solution satisfies the wave equation, the initial conditions, and the Neumann boundary condition. If initial and boundary data are compatible at (x, t) = (0, 0) then the solution is continuous at x = t. Lastly, we consider Robin boundary conditions: with given real constants a and b and a time-dependent function f (t). We evaluate this boundary condition at t = s, multiply by exp Ω j s and integrate, to obtain with e Ω j s f (s)ds, j = 1, 2.
The two equations (77) are valid for all k ∈ C. We solve these equations combined with the new global relations (42a-b) for g (j) m , j = 1, 2, m = 0, 1, obtaining g (j) 0 = f (j) + (Ω jû0 (−k) −v 0 (−k)) − e Ω j t (Ω jû (−k, t) −v(−k, t)) γ + ik , (79a) g (j) 1 = ikf (j) − γ(Ω jû0 (−k) −v 0 (−k)) + γe Ω j t (Ω jû (−k, t) −v(−k, t)) γ + ik . (79b) These expressions are substituted in (59b) with a = 0. This results in with which depends on the boundary conditions only, the dependence on the initial condition u 0 (x), and the contribution from the initial condition v 0 (x). In these definitions, the contour C accommodates the presence of the singularity at k = iγ: in the integral obtained after substituting (79a-b) in (59b), the singularity at k = iγ is removable, since the integrand of (59b) is analytic in all of C. If we distribute the integral so as to isolate the dependence ofû(−k, t) (the term withv(−k, t) cancels), both integrals have a pole singularity. If γ > 0, this singularity is in the upper-half plane, which is where we wish to close the contour to show that the contribution fromû(−k, t) vanishes. To avoid this issue, we deform the path of integration to pass above k = iγ along path C, prior to distributing the integral. Thus C is the real line if γ < 0, and C is asymptotic to the real line but passes above iγ if γ > 0. Using this contour, it immediately follows that the terms containingû(−k, t) vanish, as desired. Thus we have omitted these terms in (80).
In what follows, we need 1 2π which are obtained using the definition of the delta function and standard contour integration. As elsewhere in this paper, all integrals over k are principal-value integrals. Using the definition of f (j) , switching the order of integration, and using (84), we find