text:
# Interaction cost of non-local gates

## Abstract

We introduce the interaction cost of a non-local gate as the minimal time of interaction required to perform the gate when assisting the process with fast local unitaries. This cost, of interest both in the areas of quantum control and quantum information, depends on the specific interaction, and allows to compare in an operationally meaningful manner any two non-local gates. In the case of a two-qubit system, an analytical expression for the interaction cost of any unitary operation given any coupling Hamiltonian is obtained. One gate may be more time-consuming than another for any possible interaction. This defines a partial order structure in the set of non-local gates, that compares their degree of non-locality. We analytically characterize this partial order in a region of the set of two-qubit gates.

An elementary concern in quantum information theory is to stablish the trade-off between different physical resources that are relevant for information processing. A controlled Hamiltonian interaction between quantum systems is one instance of useful resource. It can be employed, for example, to simulate the dynamics of another multipartite quantum system. On the other hand multi-particle unitary gates are a requirement for universal quantum computation. In particular, two-qubit gates —together with one-qubit gates— can be taken as the building block of quantum computers.

A detailed study of the connections existing between non-local Hamiltonians and non-local gates is thus of interest from a quantum information perspective, but this issue is also relevant in other areas. For instance, the synthesis of multipartite gates from Hamiltonian interactions —and, in particular, time–minimizing schemes— has been recently analyzed in the context of quantum control theory . Whereas the requirements for arbitrary manipulation of single qubits are presently met in a number of experimental schemes, the engineering of two-qubit gates can be only (partially) achieved with very few systems . In real experiments not only an interaction Hamiltonian between the qubits, but also considerable command on them in order to process the interaction, are required. For instance, mechanisms to switch on and off the interaction, as well as to accurately drive the systems towards the desired joint evolution, are needed. But even from a theoretical perspective, a description of two-qubit gates in terms of interactions able to prescribe optimal protocols for gate synthesis was so far missing. Here we shall provide such a description.

More generally, we consider a set of subsystems with a given Hamiltonian $`H`$, and assume that arbitrarily fast local unitaries (LU) can be performed to properly tailor the evolution that $`H`$ induces. The aim is to perform some joint unitary transformation $`U`$ on the systems. This is the setting considered in  and corresponds to the so-called gate simulation under LU of . Two definitions are needed to specify the problems that we shall address.

Definition 1: The interaction cost $`𝒞_H(U)`$ of a non-local gate $`U`$ given a Hamiltonian $`H`$ denotes the minimal time needed in order to perform $`U`$ using the interaction $`H`$ and fast LU.

Definition 2: We say gate $`U`$ is more non-local than gate $`V`$, and write $`VU`$, when for all interactions $`H`$ the interaction cost of U is never smaller than that of $`V`$,

$$VU𝒞_H(V)𝒞_H(U)H.$$
(1)
First we shall show how the interaction cost $`𝒞_H`$ can be explicitly computed for any gate and any interaction of a two-qubit system. This is possible by considering results recently developed in the areas of quantum control  and quantum information . In  considerable progress towards the solution was made, and only a final optimization was left unsolved. The results of  provide the tools needed to perform such an optimization and thereby complete the results of .

Definition 2 introduces a partial order structure in the set of non-local gates. This structure captures the intuition, in terms of the resources needed to perform a gate, that one gate may be “more non-local” than another. Our second result is an analytical characterization of this structure in a region of the set of two-qubit gates.

We start by describing known facts concerning the simulation of non-local Hamiltonians and the synthesis of non-local gates.

($`i`$) Optimal simulation of two-qubit Hamiltonians under LU. Any Hamiltonian acting on two qubits is uniquely represented, for the purposes of simulation under LU, by its canonical form

$$H=\underset{i}{}h_i\sigma _i\sigma _i,h_1h_2|h_3|,$$
(2)
where $`\sigma _i`$, $`i=1,2,3`$, stand for the Pauli matrices. In the rest of the paper $`H`$ denotes a Hamiltonian written in its canonical form and $`\stackrel{}{h}`$ denotes the vector $`(h_1,h_2,h_3)`$ with its properly ordered coefficients. The special majorization relation $`\stackrel{}{x}_s\stackrel{}{y}`$ between three dimensional real vectors $`\stackrel{}{x}`$ and $`\stackrel{}{y}`$ is relevant in this context. It is given by the set of inequalities

$`x_1`$ $``$ $`y_1,`$ (3)
$`x_1+x_2x_3`$ $``$ $`y_1+y_2y_3,`$ (4)
$`x_1+x_2+x_3`$ $``$ $`y_1+y_2+y_3.`$ (5)
where the components $`x_i`$ and $`y_i`$ are assumed to fulfill $`x_1x_2|x_3|`$ and $`y_1y_2|y_3|`$ .

Fact 1 (Theorem of ): The minimal time overhead $`t_{H^{}|H}`$ (i.e., the inverse of efficiency $`s_{H^{}|H}`$) required to simulate Hamiltonian $`H^{}`$ by Hamiltonian $`H`$ and fast LU is the minimal value of $`c0`$ such that the vectors $`\stackrel{}{h}^{}`$ and $`\stackrel{}{h}`$ satisfy $`\stackrel{}{h}^{}_sc\stackrel{}{h}`$. Protocols for optimal simulation are known.

($`ii`$) Optimal synthesis of two-qubit gates under LU. Any two-qubit gate $`U`$ can be written in terms of local unitaries $`u_i`$ and $`v_i`$ and Pauli matrices $`\sigma _k`$ as

$$U=(u_1v_1)e^{i_k\lambda _k\sigma _k\sigma _k}(u_2v_2).$$
(6)
In  it is shown how to obtain this decomposition. Notice that gate $`U`$ is equivalent, up to local unitaries performed on the qubits before and after $`U`$, to

$$U_\stackrel{}{\lambda }e^{i_k\lambda _k\sigma _k\sigma _k}.$$
(7)
Since we assume the ability to perform instantaneous (i.e., sufficiently fast) LU operations, the synthesis of $`U`$ is as time-consuming as that of $`U_\stackrel{}{\lambda }`$, and we need only focus on the later. In addition , to each $`U`$ there corresponds a unique $`U_{\stackrel{}{\lambda }^0}`$ with $`\stackrel{}{\lambda }^0(\lambda _1^0,\lambda _2^0,\lambda _3^0)`$ such that $`\lambda _1^0\lambda _2^0|\lambda _3^0|`$, $`\lambda _1^0,\lambda _2^0[0,\pi /4]`$ and $`\lambda _3^0(\pi /4,\pi /4]`$, that we shall call its canonical form. In what follows we will often represent any two-qubit gate $`U`$ by its canonical form $`U_{\stackrel{}{\lambda }^0}`$ or by its corresponding (unique) vector $`\stackrel{}{\lambda }^0=(\lambda _1^0,\lambda _2^0,\lambda _3^0)`$. Recall that all commutators $`[\sigma _j\sigma _j,\sigma _k\sigma _k]`$ vanish, and that $`\mathrm{exp}(i\pi /2\sigma _j\sigma _j)=i\sigma _j\sigma _j`$ is a local gate. This implies that for any vector $`\stackrel{}{n}=(n_1,n_2,n_3)`$ with integer components $`n_j`$,

$$U_{\stackrel{}{\lambda }^0}=_{LU}U_{\stackrel{}{\lambda }^0}U_{\frac{\pi }{2}\stackrel{}{n}}=U_{\stackrel{}{\lambda }^0+\frac{\pi }{2}\stackrel{}{n}},$$
(8)
with $`\stackrel{}{\lambda }^0+\pi /2\stackrel{}{n}`$ essentially exhausting all vectors compatible with the gate $`U_{\stackrel{}{\lambda }^0}`$ . \[$`=_{LU}`$ is used to denote equivalence under LU.\].

In Theorem 10 of  the problem of time-optimally producing a two-qubit gate $`U`$ using interaction $`H`$ is shown to reduce to a specific minimization over all possible decompositions of $`U`$ of the form (6). Here we rephrase the theorem in terms of the notion of Hamiltonian simulation and the concepts introduced before. Without loss of generality, we refer only to unitary operations that can be written as in (7), and associate a self-adjoint operator $`H_\stackrel{}{\lambda }_i\lambda _i\sigma _i\sigma _i`$ to each possible decomposition.

Fact 2 (Theorem 10 of , readapted): The time-optimal way to synthesize gate $`U`$ with interaction $`H`$ and fast $`LU`$ consists of simulating, among all Hamiltonians $`H_\stackrel{}{\lambda }`$ such that $`U=\mathrm{exp}(iH_\stackrel{}{\lambda })`$, the one with smallest time overhead $`t_{H_\stackrel{}{\lambda }|H}`$. The minimal interaction time (i.e., the interaction cost $`𝒞_H(U)`$) is given by the smallest time overhead $`t_{H_\stackrel{}{\lambda }|H}`$.

Our first aim is to perform the optimization described in Fact 2. This is feasible because we have an analytical characterization both of all (infinitely many) decompositions of $`U`$ (cf. Eq. (8)) and of the time overhead $`s_{H_\stackrel{}{\lambda }|H}`$ for any decomposition (cf. Fact 1), as expressed in the following lemma.

Lemma: The interaction cost $`𝒞_H(U)`$ is the minimal value $`c0`$ such that a vector $`\stackrel{}{n}`$ of integers exists satisfying

$$\stackrel{}{\lambda }^0+\frac{\pi }{2}\stackrel{}{n}_sc\stackrel{}{h}.$$
(9)
It is useful to introduce, for each $`\stackrel{}{n}`$, the pre-cost $`c_\stackrel{}{n}`$ as the minimal value $`c0`$ such that $`\stackrel{}{\lambda }^0+\pi /2\stackrel{}{n}_sc\stackrel{}{h}`$. Pre-cost $`c_\stackrel{}{n}`$ is the overhead needed to simulate $`H_{\stackrel{}{\lambda }^0+\pi /2\stackrel{}{n}}`$ by $`H`$ or, equivalently, the minimal time $`t`$ needed to travel, in the set of non-local gates, from the identity operator to $`U`$ along the path defined by $`\stackrel{}{\lambda }^0+\pi /2\stackrel{}{n}`$. Intuitively, a large $`\stackrel{}{n}`$ corresponds to a “long” —and therefore non-optimal— path. Following this intuition we arrive at our first result.

Theorem 1: The interaction cost $`𝒞_H(U)`$ or minimal time needed to create gate $`U`$ by using Hamiltonian $`H`$ and fast LU is given by

$$𝒞_H(U)=\mathrm{min}\{c_{(0,0,0)},c_{(1,0,0)}\},$$
(10)
that is, the minimal of two pre-costs, one corresponding to the canonical vector $`\stackrel{}{\lambda }^0=(\lambda _1^0,\lambda _2^0,\lambda _3^0)`$ of $`U`$ and the other to the vector $`(\frac{\pi }{2}\lambda _1^0,\lambda _2^0,\lambda _3^0)`$ . The time-optimal protocol consists in simulating the corresponding Hamiltonian (either $`\stackrel{}{h}_1(\lambda _1^0,\lambda _2^0,\lambda _3^0)`$ or $`\stackrel{}{h}_2=(\frac{\pi }{2}\lambda _1^0,\lambda _2^0,\lambda _3^0)`$) by $`H`$ for time $`t=𝒞_H(U)`$.

Remark. Thus, in order to time-optimally perform gate $`U`$ with Hamiltonian $`H`$, we can proceed as follows. Using Ref. , we compute $`\stackrel{}{\lambda }^0`$ from $`U`$, and using Refs.  we compute $`\stackrel{}{h}`$ from $`H`$. Theorem 1 gives the minimal time of simulation and the Hamiltonian (either $`\stackrel{}{h}_1`$ or $`\stackrel{}{h}_2`$) to be simulated, and finally Ref.  describes an optimal protocol for simulating the convenient Hamiltonian by $`H`$ and LU.

Proof: We need to see that $`𝒞_H(U)`$ as given by Eq. (10) is the minimal pre-cost, i.e. $`𝒞_H(U)c_\stackrel{}{n}`$ for all $`\stackrel{}{n}`$. It is straightforward to check from Eq. (5) that (T1.$`i`$) for any two vectors $`\stackrel{}{x}`$ and $`\stackrel{}{y}`$, with components $`x_1x_2|x_3|`$, $`y_1y_2|y_3|`$, the minimal $`c0`$ such that $`\stackrel{}{x}_sc\stackrel{}{y}`$ satisfies $`c3x_1/y_1`$; (T1.$`ii`$) if $`\stackrel{}{x}_s\stackrel{}{x}^{}`$, then $`\stackrel{}{x}^{}_s\stackrel{}{y}\stackrel{}{x}_s\stackrel{}{y}`$ ($`_s`$ is a partial order!). In particular, let $`c^{}0`$ be the minimal value such that $`\stackrel{}{x}^{}_sc^{}\stackrel{}{y}`$. Then $`\stackrel{}{x}_s\stackrel{}{x}^{}\stackrel{}{x}_sc^{}\stackrel{}{y}`$, so that the minimal $`c0`$ such that $`\stackrel{}{x}_sc\stackrel{}{y}`$ always satisfies $`cc^{}`$. Now, recall that by definition $`\pi /4\lambda _1^00`$, and notice that if some component $`n_j`$ of $`\stackrel{}{n}`$ fulfills $`|n_j|>1`$, then the maximal component of the reordered version  of $`\stackrel{}{\lambda }^0+\pi /2\stackrel{}{n}`$ is at least $`3\pi /4`$. Thus, because of (T1.$`i`$), $`\stackrel{}{\lambda }^0_s\stackrel{}{\lambda }^0+\pi /2\stackrel{}{n}`$. Then (T1.$`ii`$) implies that $`c_{(0,0,0)}c_\stackrel{}{n}`$. Therefore we can restrict our attention to vectors $`\stackrel{}{n}`$ with $`|n_j|1`$. A case by case check shows that the pre-costs $`c_\stackrel{}{n}`$ with $`\stackrel{}{n}\{(1,1,1),(0,1,0),(0,0,1),(0,0,1)\}`$ fulfill $`c_{(1,0,0)}c_\stackrel{}{n}`$, since (cf. point (T1.$`ii`$) above) $`\stackrel{}{\lambda }^0+\pi /2(1,0,0)_s\stackrel{}{\lambda }^0+\pi /2\stackrel{}{n}`$ . Similarly, we obtain that for the remaining vectors $`\stackrel{}{n}`$ with $`|n_j|1`$ the pre-costs satisfy $`c_{(0,0,0)}c_\stackrel{}{n}`$, because $`\stackrel{}{\lambda }^0+\pi /2(0,0,0)_s\stackrel{}{\lambda }^0+\pi /2\stackrel{}{n}`$. The only remaining configurations, with vectors $`\stackrel{}{n}\{(1,0,0),(0,0,0)\}`$, are incomparable according to the $`_s`$ relation —unless $`\lambda _1^0+|\lambda _3^0|\pi /4`$, in which case we always obtain $`c_{(0,0,0)}c_{(1,0,0)}`$—, and this is why the optimization of Eq. (10) has to be performed. $`\mathrm{}`$

Corollary: ($`a`$) When U is such that $`\lambda _1^0+|\lambda _3^0|\pi /4`$, then the interaction cost is always given by

$$C_H(U)=c_{(0,0,0)}.$$
(11)
($`b`$) If, instead, $`\lambda _1^0+|\lambda _3^0|<\pi /4`$, then Hamiltonians $`H`$ and $`H^{}`$ always exist such that $`C_H(U)=c_{(0,0,0)}^\stackrel{}{h}<c_{(1,0,0)}^\stackrel{}{h}`$ and $`C_H^{}(U)=c_{(1,0,0)}^\stackrel{}{h}^{}<c_{(0,0,0)}^\stackrel{}{h}^{}`$.

Proof: ($`a`$) follows from the fact that $`\lambda _1^0+|\lambda _3^0|\pi /4\stackrel{}{\lambda }^0_s\stackrel{}{\lambda }^0+\pi /2(1,0,0)`$. (b) can be checked by considering $`H`$ and $`H^{}`$ given by $`\stackrel{}{h}=\stackrel{}{\lambda }^0`$ and $`\stackrel{}{h}^{}=(\pi /2\lambda _1^0,\lambda _2^0,\lambda _3^0)`$. $`\mathrm{}`$

In order to analyze Eq. (10) we first consider some examples. For the Ising interaction $`H=h\sigma _3\sigma _3`$ (equivalently $`h\sigma _1\sigma _1`$) and an arbitrary gate $`U`$, Eq. (10) reads (cf. Theorem 2 of ),

$$C_{h\sigma _1\sigma _1}(U)=\frac{\lambda _1^0+\lambda _2^0+|\lambda _3^0|}{h}.$$
(12)
Let us now instead focus on three specific gates and arbitrary interactions. By $`|m|n`$ ($`m,n=0,1`$) we denote the computational basis of two-qubits. The CNOT gate is defined as

$$|m|n|m|nm,$$
(13)
where $``$ is sum modulo 2. Using the method described in Ref.  we obtain its canonical vector, $`\stackrel{}{\lambda }^0=(\pi /4,0,0)`$. Similarly, the SWAP gate,

$$|m|n|n|m,$$
(14)
has vector $`\stackrel{}{\lambda }^0=(\pi /4,\pi /4,\pi /4)`$. We also consider a third, intermediate gate $`U_{XY}`$ with $`\stackrel{}{\lambda }^0=(\pi /4,\pi /4,0)`$, that corresponds to

$$|m|ni^{|mn|}|n|m.$$
(15)
For these three gates we find

$`𝒞_H(\text{CNOT})`$ $`=`$ $`{\displaystyle \frac{\pi }{4}}{\displaystyle \frac{1}{h_1}},`$ (16)
$`𝒞_H(U_{XY})`$ $`=`$ $`{\displaystyle \frac{\pi }{4}}{\displaystyle \frac{2}{h_1+h_2|h_3|}},`$ (17)
$`𝒞_H(\text{SWAP})`$ $`=`$ $`{\displaystyle \frac{\pi }{4}}{\displaystyle \frac{3}{h_1+h_2+|h_3|}}.`$ (18)
With these examples at hand we make the following two observations. First, to any gate $`U`$ there corresponds a natural interaction $`H_U`$, with vector either $`\stackrel{}{h}_1`$ or $`\stackrel{}{h}_2`$ as defined in Theorem 1. This natural interaction allows to perform gate $`U`$ optimally without need to intermediately simulate another Hamiltonian and therefore the time inefficiencies inherent in the process of simulation are avoided. In this sense the natural interactions for the CNOT gate, gate $`U_{XY}`$ and the SWAP gate are, respectively, the Ising interaction $`\sigma _1\sigma _1`$, the XY–model interaction $`\sigma _1\sigma _1+\sigma _2\sigma _2`$ and the Heisenberg or exchange interaction $`\sigma _1\sigma _1+\sigma _2\sigma _2+\sigma _3\sigma _3`$.

The second observation is that for any fixed interaction $`H`$, e.g. $`H=h\sigma _1\sigma _1`$ as in Eq. (12), the interaction cost induces an order in the set of gates. For instance, according to Eq. (12), a SWAP is the most time–consuming gate when the Ising interaction is available. Eqs. (16-18) also show, however, that such an order depends on the available interaction. Using the exchange interaction $`H=\sigma _1\sigma _1+\sigma _2\sigma _2+\sigma _3\sigma _3`$, $`U_{XY}`$ is twice as time-consuming as a SWAP gate.

Let us move to Definition 2. It endows the set of non-local gates with a partial order structure based on the notion of interaction cost, but which is independent of any particular interaction. By comparing the resources required to perform two gates, such a partial order captures the intuition that some gates are more non-local than others.

We have already argued that no gate more non-local (i.e., more time-consuming for all interactions) than all the others exists. It is also easy to see that a gate $`\alpha \stackrel{}{\lambda }^0`$ is always less non-local than $`\stackrel{}{\lambda }^0`$ for any $`\alpha [0,1]`$, since the pre-costs are linear in $`\alpha `$. Next we present an analytical characterization of the partial order relation $`VU`$ in a region of the set of two-qubit gates .

Theorem 2: Let $`U`$ and $`V`$ be two two-qubit gates with corresponding ordered vectors $`\stackrel{}{\lambda }_U^0`$ and $`\stackrel{}{\lambda }_V^0`$ such that in both cases the restriction $`\lambda _1^0+|\lambda _3^0|\pi /4`$ holds. Then gate $`U`$ is more non-local than gate $`V`$ if and only if $`\stackrel{}{\lambda }_V^0_s\stackrel{}{\lambda }_U^0`$,

$$VU\stackrel{}{\lambda }_V^0_s\stackrel{}{\lambda }_U^0.$$
(19)
Proof: Recall that the restrictions on $`\stackrel{}{\lambda }_U^0`$ and $`\stackrel{}{\lambda }_V^0`$ imply, because of Corollary 1, that the interaction costs $`𝒞_H(U)`$ and $`𝒞_H(V)`$ are given, respectively, by the smallest $`c_U,c_V0`$ such that

$`\stackrel{}{\lambda }_U^0`$ $`_s`$ $`c_U\stackrel{}{h},`$ (20)
$`\stackrel{}{\lambda }_V^0`$ $`_s`$ $`c_V\stackrel{}{h}.`$ (21)
Suppose first that $`VU`$, that is, that for any Hamiltonian $`H`$ we have $`C_H(V)C_H(U)`$. Then we also have $`\stackrel{}{\lambda }_V^0_s𝒞_H(U)\stackrel{}{h}`$. In particular, if we choose the interaction $`H`$ to have vector $`\stackrel{}{h}=\stackrel{}{\lambda }_U^0`$, we have $`𝒞_H(V)𝒞_H(U)=1`$ and $`\stackrel{}{\lambda }_V^0_s𝒞_H(U)\stackrel{}{h}=\stackrel{}{\lambda }_U^0`$, which proves the direct implication. The inverse implication follows from (T1.$`ii`$) of the proof of Theorem 1, which shows $`\stackrel{}{\lambda }_V^0_s\stackrel{}{\lambda }_U^0`$ implies that $`𝒞_H(V)𝒞_H(U)`$ for all $`H`$. $`\mathrm{}`$

As an example of this result, we see that the $`U_{XY}`$ gate is more non-local than the CNOT gate, and that, as it was to be expected, gates with sufficiently small components $`|\lambda _i^0|`$ are less non-local than those with large $`|\lambda _i^0|`$ .

Finally, recall that entanglement can be used as a catalyzer for Hamiltonian simulation . In particular, if to each of the two interacting qubits $`A`$ and $`B`$ we attach an extra qubit $`A^{}`$ and $`B^{}`$, where the pair $`A^{}B^{}`$ is in a maximally entangled state, Hamiltonian $`H`$ between $`A`$ and $`B`$ can be used to perform, without consuming the entanglement of $`A^{}B^{}`$, more powerful simulations than before, provided that fast LU are allowed in $`AA^{}`$ and $`BB^{}`$. Consequently, the interaction cost of gates is modified when not only LU, but also entangled ancillas are available.

In this work we have characterized the time-optimal synthesis of two-qubit unitary transformations using an arbitrary two-qubit Hamiltonian. In particular, the interaction cost $`𝒞_H(H)`$ has been computed and optimal protocols have been described. We have also characterized, in a region of the space of two-qubit gates, a partial order structure related to their degree of non-locality. These results can be applied to the study of the interaction cost for particular processes, such as the creation of a maximally entangled state  or the transmition of a classical or quantum bit of information from one qubit to another . All these discussions involve only two interacting qubits. It would be desirable to obtain a generalization to higher–dimensional systems. The lack of an analog to decomposition (6) in these cases is a serious drawback. Another interesting generalization consists in considering the asymptotic scenario, where the aim is to perform a large number of copies of the same gate.

This work was supported by the European Community under project EQUIP (contract IST-1999-11053), the ESF and the Institute for Quantum Information GmbH. G.V. is supported by grant HPMF-CT-1999-00200 (Marie Curie fellowship) of the European Community.


format: text

text:
# Untitled Document

Roman O. Popovych and Irina A. Yehorchenko

Institute of Mathematics, National Academy of Science of Ukraine,

3 Tereshchenkivs’ka Street, 01601, Kyiv-4, Ukraine

$``$ E-mail: rop@imath.kiev.ua

$``$ E-mail: iyegorch@imath.kiev.ua

Group Classification of Generalised Eikonal Equations

## Abstract

A new approach to the problem of group classification is applied to the class of first-order non-linear equations of the form $`u_au_a=F(t,u,u_t)`$. It allowed complete solution of the group classification problem for a class of equations for functions depending on multiple independent variables, where highest derivatives enter nonlinearly. Equivalence groups of the class under consideration and algebraic properties of the symmetry algebra are studied.

The class of equations considered presents generalisation of the eikonal and Hamilton Jacobi equations. The paper contains the list of all non-equivalent equations from this class with symmetry extensions, and proofs of such non- equivalence.

New first order non-linear equations possessing wide symmetry groups were constructed.

1. Introduction. The problem of group classification of differential equations is determination of all particular equations from the class possessing wider symmetry when the class in general. Mathematical models, and in particular, differential equations of mathematical physics, often contain parameters (numerical or functional), that are determined experimentally. For this reason the parameters are not fixed exactly. Theoretical selection of mathematical models can be done on the basis of the symmetry principle, that is to select among given values of parameters ones that lead to a mathematical model with given properties, or with the widest possible symmetry group .

We considered a class of mathematical models given by a rather general set of non-linear first-order equations that are invariant under the Euclid group – the group of space rotations and translations. We did not fix any symmetry with respect to time and dependent variables. The set of mathematical models considered incorporates models with relativistic and Galilei symmetries that are well-known and widely used in physics, geometric optics and other applications.

The basics for the group classification were constructed by S. Lie himself  in his classification of integrable ordinary differential equations. The start of systematic study of similar problems for partial differential equations is linked to papers by L.Ovsyannikov (see  and references wherein) and by his collaborators. In particular, Ovsyannikov introduced such important concepts as arbitrary element, kernel of principal groups, equivalence group and a number of others.

At present there are quite a few papers on group classification of important classes of models in mechanics, physics, biology and other sciences (see e.g. , and these present only a small part of existing papers!).

Until recently the principal method for solving of group classification problems for partial differential equations was direct integration of determining equations for a Lie symmetry operator with sorting out all possible cases of integration. Application only of the above method substantially narrowed the scope of solvable group classification problems, as determining equations may be extremely difficult to solve, if it is possible at all. Most of the group classification problems solved in the literature are problems for classes of equations with arbitrary functions on only one variable each. Arbitrary functions with multiple variables result in determining equations being partial differential equations, rather than ordinary. It is obvious that this complicates integration considerably.

Recently new approaches appeared having allowed solution of group classification problems that could not be solved by means of the methods known before. Let us mention in more detail a few papers whose ideas are the closest to this paper.

I.S. Akhatov, R.K. Gazizov and N.K. Ibragimov , developing the approach by L.Ovsyannikov, systematically employed the equivalence group for simplification of classifying conditions.

R. Zhdanov and V. Lahno  suggested a new approach to group classification of PDEs that presents a generalisation of S. Lie’s approach to group classification of ODEs  (see also ). This approach is based on classification of low-dimensional abstract algebrae and works best for equations with two independent variables. It allows avoiding direct solving of cumbersome determining equations and finding complete solutions for non-linear group classification problems for parameter functions depending on many variables. Determining equations to be solved are partial differential equations with respect to the parameter functions. Further development for this approach presented in this paper is applied for equations with arbitrary number of independent variables. R. Popovych and R. Cherniha  employed for classification of systems of coupled non-linear Laplace equations a combined method in which the number of cases that are considered during integration of determining equations, is essentially decreased at the expense of the study of possible structure of the symmetry algebra. There are a number of other papers suggesting new approaches to group classification problems.

In the present paper the combined method is used for group classification of generalised eikonal and Hamilton-Jacobi equations. The results of the paper show that this method that includes investigation of determining equations, of equivalence transformation groups and of the symmetry algebra structures is an efficient tool for finding equations from a certain class with symmetry extensions. However, the method is not strictly algorithmic, and the steps to be used would be different for different classes of equations.

We succeeded to give a full description of subclasses for the class considered where symmetry extension arises compared to the general symmetry group for the class in total. We expect that the method suggested would allow solving other group classification problems, especially for higher order equations and for systems of PDEs, and extending the variety of non-linear equations with wide symmetry that could be used for construction of mathematical models in physics and other natural sciences.

2. Formulation of the Problem. Let us consider the equation of the form

$$u_au_a=F(t,u,u_t)$$
(1)
for a real function $`u=u(t,x)`$ of $`n+1`$ independent variables $`t=x_0`$ and $`x=(x_1,x_2,\mathrm{},x_n)`$, $`n2`$. Here and below a lower index of a function will designate differentiation with respect to the corresponding variable. The indices $`a`$ and $`b`$ take values from 1 to $`n`$. The summation over repeated indices is implied unless stated otherwise. The class of equations (1) is the general class of non-linear first-order equations invariant under the Euclid group – the group of space rotations and translations. It is remarkable also as it includes all the following well-known equations:

$$\begin{array}{cc}F=2mu_t\hfill & \text{ is the Hamilton-Jacobi equation for a free particle;}\hfill \\ F=u_t^21\hfill & \text{ is the relativistic Hamilton equation;}\hfill \\ F=u_t^2\hfill & \text{ is the eikonal equation.}\hfill \end{array}$$

The symmetry properties of these equations and relations between them were investigated in . Equations of the form $`u_t^2=f\left(u\right)u_au_a`$ that are equivalent to the eikonal equation, were used in  for studying of solutions of non-linear wave equations.

We perform group classification for equations of the form (1) by an arbitrary element — a smooth function $`F=F(t,u,u_t)0.`$ (Here and below smoothness means continuous differentiability.)

3. Remarks on Formulation of Classification Problems and Schemes for their Solving. We would like to explain our choice of the class we consider for group classification. This choice was governed by the following rules:

1. A class has to generalise physically interesting equations.

2. A problem has to be symmetric by itself: a class under consideration should have a wide equivalence group.

3. A class has to preserve at least some of important symmetry properties of equations being generalized.

The second rule was the basis of the approach suggested in . Both second and third rules were the reason why we included dependence on the variable $`t`$ for the arbitrary function $`F`$. That preserved the possibility of symmetries interchanging $`t`$ and $`u`$ that the listed well-known equations have. Introduction of $`t`$ substantially increases the equivalence group, and despite most resulting classes with symmetry extensions could be represented by $`F=F(u,u_t)`$, or $`F`$ does not depend on $`t`$, the equivalence transformations employed in the classification process belonged to the wide equivalence group corresponding to the class with $`F=F(t,u,u_t)`$.

The general outline of our classification could be presented as follows.

1. Find determining equations and the kernel of main groups.

2. Find the equivalence group.

3. Integrate a part of the system of determining equations and obtain a possible form for invariance algebra operators.

4. Classify possible symmetry algebrae that are lowest-order extensions of the main algebra and that possess a certain structure identified as a separate case at the step 3. Distinguish the corresponding non-equivalent cases.

5. Impose conditions for the invariance algebrae having larger dimension or different structure than those considered at the step 4, find new cases. The requirement of a different structure than was considered before allowed to study remaining cases by means of the standard approach.

4. Kernel of Main Groups and Equivalence Group. Let an infinitesimal operator

$$Q=\xi ^0(t,x,u)_t+\xi ^a(t,x,u)_a+\eta (t,x,u)_u$$

generate a one-parameter symmetry group of local transformation for the equation (1). Then from the infinitesimal invariance criterion  after transition to the manifold given by the equation (1) in the prolonged space, and splitting by non-related variables we derive the following determining equations for the coefficients of the operator $`Q`$:

$$\xi _b^a+\xi _a^b=0,\xi _a^a=\xi _b^b,ab$$
(2)
(here there is no summation by $`a`$ and $`b`$),

$$\left(\xi _t^a+\xi _u^au_t\right)F_{u_t}2\xi _u^aF+2\eta _a2u_t\xi _a^0=0,$$
(3)
$$\xi ^0F_t+\eta F_u+\left(\eta _t+\left(\eta _u\xi _t^0\right)u_t\xi _u^0u_t^2\right)F_{u_t}=2\left(\eta _u\xi _1^1\xi _u^0u_t\right)F.$$
(4)
If the function $`F`$ is not fixed, then splitting in (3) and (4) by the “variables” $`F`$, $`F_t`$, $`F_u`$, $`F_{u_t}`$, $`u_t`$, we obtain that $`\eta =\xi ^0=0`$, $`\xi _1^1=0`$, $`\xi _t^a=\xi _u^a=0`$, whence taking into account (2) we come to the following statement.

Statement. The kernel of the main groups for the equation (1) is the Euclid group $`E\left(n\right)`$, the algebra of which is $`A^{\mathrm{ker}}=e\left(n\right)=_a,J_{ab}=x_a_bx_b_a.`$

The equivalence group for the equation (1) coincides with the group generated by the set of one-parameter groups of local symmetries of the system

$$u_au_a=F,w=u_t,F_a=0,$$
(5)
whose infinitesimal operators have the form

$$\begin{array}{c}\widehat{Q}=\widehat{\xi }^0(t,x,u)_t+\widehat{\xi }^a(t,x,u)_a+\widehat{\eta }(t,x,u)_u\hfill \\ +\widehat{\theta }(t,x,u,w)_w+\widehat{\chi }(t,x,u,w,F)_F.\hfill \end{array}$$

From the infinitesimal invariance criterion for the system (5) after splitting with respect to non-related variables we obtain determining equations for coefficients of the operator $`\widehat{Q}`$, from which it follows that an infinitesimal operator of any one-parameter equivalence group for the equation (1) is a linear combination of operators

$$_a,J_{ab},x_a_a2F_F,\overline{\xi }_t+\overline{\eta }_u+2\left(\overline{\eta }_u\overline{\xi }_uu_t\right)F_F,$$
(6)
where $`\overline{\xi }`$ and $`\overline{\eta }`$ are arbitrary smooth functions of variables $`t`$ and $`u`$. Thus, equivalence transformations, that act non-trivially on the parameter $`F`$, have the form

$$\begin{array}{c}\stackrel{~}{t}=\zeta (t,u),\stackrel{~}{u}=\phi (t,u),\stackrel{~}{x}=\delta x,\hfill \\ \stackrel{~}{F}=\delta ^2\left(\phi _u\zeta _u\stackrel{~}{u}_{\stackrel{~}{t}}\right)^2F=\delta ^2\left(\frac{\zeta _t\phi _u\zeta _u\phi _t}{\zeta _t+\zeta _uu_t}\right)^2F,\hfill \end{array}$$
(7)
where $`\delta `$ is a non-zero constant, $`\zeta `$ and $`\phi `$ are arbitrary smooth functions of variables $`t`$ and $`u`$, for which $`\zeta _t\phi _u\zeta _u\phi _t0`$.

If we limit the class of equations (1), having imposed an additional condition $`F_t=0`$ (this subclass is separated in the process of classification), then for calculation of the corresponding equivalence group it is necessary to add this condition to the system (5). As a result, the equivalence group is narrowed: an infinitesimal operator of any one-parameter equivalence group for the equation (1) with $`F=F(u,u_t)`$ is a linear combination of the operator $`t_t`$ and of operators (6), where the functions $`\overline{\xi }`$ and $`\overline{\eta }`$ now depend only on the variable $`u`$, and therefore equivalence transformations that act non-trivially on the parameter $`F`$ have the form

$$\begin{array}{c}\stackrel{~}{t}=\widehat{\delta }t+\zeta \left(u\right),\stackrel{~}{u}=\phi \left(u\right),\stackrel{~}{x}=\delta x,\hfill \\ \stackrel{~}{F}=\delta ^2\left(\phi _u\zeta _u\stackrel{~}{u}_{\stackrel{~}{t}}\right)^2F=\delta ^2\left(\frac{\widehat{\delta }\phi _u}{\widehat{\delta }+\zeta _uu_t}\right)^2F,\hfill \end{array}$$
(8)
where $`\delta ,`$ $`\widehat{\delta }`$ are non-zero constants, $`\zeta `$ and $`\phi `$ are arbitrary smooth functions of the variable $`u`$, and $`\phi _u0`$.

Further limitation of the class of equations (1) may lead to appearance of equivalence transformations of the form (7), that are different from (8) (see the proof).

5. Result of Classification. All possible cases of extension of the maximal invariance algebra in the Lie sense for the equation (1) up to equivalence transformations (7) are exhausted by the cases listed in the Table 1.

In the Table 1 $`F=F(t,u,u_t)`$, $`f=f(u,u_t)`$, $`h=h(u_t)`$ are arbitrary smooth functions of their arguments, $`\delta `$ is constant, $`\epsilon _1,\epsilon _2=\pm 1`$, $`(\epsilon _1,\epsilon _2)(1,1)`$, $`J_{ab}=x_a_bx_b_a`$, $`D=t_t+u_u+x_a_a`$. In the case 7 $`g_{\mu \nu }`$ is metric tensor of Minkowsky space $`^{1,n}`$, i.e. $`g_{00}=g_{11}=g_{22}=g_{33}=1`$, $`g_{\mu \nu }=0`$, $`\mu \nu `$; $`c^\mu `$, $`b^{\mu \varkappa }`$, $`d`$, $`a^\varkappa `$, $`\eta `$ are arbitrary smooth functions of the variable $`u`$ (for the set of such operators to form an algebra, it is necessary to require that these functions were infinitely differentiable or real analytical); the indices $`\mu `$, $`\nu `$ and $`\varkappa `$ take values from 0 to 3.

Let us note that the known Hamilton-Jacobi and eikonal equations, and the relativistic Hamilton equation are selected in the family of equations (1) as representatives of the classes of equivalent equations having the widest symmetry.

Table 1

|  | $`F`$ | Basis Symmetry Operators |
| --- | --- | --- |
| $`0`$ | $`F(t,u,u_t)`$ | $`_a`$, $`J_{ab}`$ |
| $`1`$ | $`e^{\delta t}f(u,u_t)`$, $`\delta \{0;1\}`$ | $`_a`$, $`J_{ab}`$, $`2_t\delta x_a_a`$ |
| $`2`$ | $`e^uh(u_t)`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`2_ux_a_a`$ |
| $`3`$ | $`|u|^{2\delta }h(u_t)`$, $`\delta 2`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`2t_t+2u_u+\delta x_a_a`$ |
| $`4`$ | $`h(u_t)`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`_u`$, $`D`$ |
| $`5`$ | $`e^{u_t}`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`_u`$, $`D`$, $`x_a_a2t_u`$ |
| $`6`$ | $`|u_t|^\beta `$, $`\beta 0,1,2`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`_u`$, $`D`$, $`(\beta 2)x_a_a2u_u`$ |
| $`7`$ | $`u_t^2`$ | $`2g_{\mu \nu }c^\mu (u)x_\nu x_\varkappa _\varkappa c^\varkappa (u)g_{\mu \nu }x_\mu x_\nu _\varkappa +g_{\mu \nu }b^{\mu \varkappa }(u)x_\nu _\varkappa +`$ $`+d(u)x_\varkappa _\varkappa +a^\varkappa (u)_\varkappa +\eta (u)_u`$ |
| $`8`$ | $`\epsilon _2u_t^2+\epsilon _1`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`_u`$, $`D`$, $`J_{ua}=u_a+\epsilon _1x_a_u`$, $`J_{ta}=t_a+\epsilon _2x_a_t`$, $`J_{ut}=u_t\epsilon _1\epsilon _2t_u`$, $`K_a=2x_aDs^2_a`$, $`K_u=2uD+\epsilon _1s^2_u`$, $`K_t=2tD+\epsilon _2s^2_t`$, where $`s^2=x_ax_a\epsilon _1u^2\epsilon _2t^2`$ |
| $`9`$ | $`\epsilon _2e^uu_t^2+\epsilon _1`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`t_t+2_u`$, $`(t^24\epsilon _1\epsilon _2e^u)_t+4t_u`$ |
| $`10`$ | $`\mathrm{cos}^2uu_t^2+1`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`\mathrm{cos}t\mathrm{tan}u_t\mathrm{sin}t_u`$, $`\mathrm{sin}t\mathrm{tan}u_t+\mathrm{cos}t_u`$ |
| $`11`$ | $`\pm (\mathrm{cos}^2uu_t^21)`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`\mathrm{cosh}t\mathrm{tan}u_t+\mathrm{sinh}t_u`$, $`\mathrm{sinh}t\mathrm{tan}u_t+\mathrm{cosh}t_u`$ |
| $`12`$ | $`\mathrm{cosh}^2uu_t^2+1`$ | $`_a`$, $`J_{ab}`$, $`_t`$, $`\mathrm{cosh}t\mathrm{tanh}u_t\mathrm{sinh}t_u`$, $`\mathrm{sinh}t\mathrm{tanh}u_t\mathrm{cosh}t_u`$ |

6. Proof. The complete proof is very cumbersome. We omit technical calculations, and some interim results are given as lemmas, and without proofs.

Let us denote the maximal invariance algebra in the Lie sense for the equation (1) as $`A^{\mathrm{max}}`$. In investigation of determining equations (2)– (4) two essentially different cases arise: $`F_{u_tu_tu_t}0`$ and $`F_{u_tu_tu_t}=0`$.

If $`𝑭_{𝒖_𝒕𝒖_𝒕𝒖_𝒕}\mathbf{}\mathrm{𝟎}`$, then it follows from (3) and (4), that $`\eta _a=\xi _a^0=\xi _t^a=\xi _u^a=0`$, $`\xi _{1a}^1=0`$, thus $`A^{\mathrm{max}}=A^{\mathrm{ker}}+A^{\mathrm{ext}}`$, and $`A^{\mathrm{ext}}\delta x_a_a+\xi ^0(t,u)_t+\eta (t,u)_u`$, where $`\delta `$ is an arbitrary constant, $`\xi `$ and $`\eta `$ are arbitrary smooth functions of the variables $`t`$ and $`u`$. $`A^{\mathrm{ker}}`$ ia an ideal, and $`A^{\mathrm{ext}}`$ is a subalgebra of the algebra $`A^{\mathrm{max}}`$, and $`QA^{\mathrm{ext}}(Q0)\text{:}(\xi ^0,\eta )(0,0).`$ The dimension of $`A^{\mathrm{ext}}`$ determines dimension of an extension of the algebra $`A^{\mathrm{max}}`$.

If $`dimA^{\mathrm{ext}}>0`$, then we may choose any non-zero operator $`Q`$ from $`A^{\mathrm{ext}}`$. By equivalence transformations (7) it is always possible to reduce it to the form $`Q=2_t\delta x_a_a`$, where $`\delta \{0;1\}`$, after that, having solved equations (4) with $`\xi ^0=2`$, $`\xi ^a=\delta x_a`$, $`\eta =0`$ with respect to $`F`$, we obtain the first case of extension of the algebra $`A^{\mathrm{max}}`$.

Let $`dimA^{\mathrm{ext}}>1`$. Then $`A^{\mathrm{ext}}`$ has to contain an operator $`Q`$ with $`\delta =0`$. By equivalence transformations (7) we can reduce the operator $`Q`$ to the operator $`_t`$. In the following we will assume that $`_tA^{\mathrm{ext}}`$, whence $`F_t=0`$, i.e. $`F=F(u,u_t)`$, and we will say that there is an additional extension of symmetry, if for such $`F`$ the relation $`dimA^{\mathrm{max}}>dimA^{\mathrm{ker}}+1`$ is satisfied.

Let us assume that $`dimA^{\mathrm{ext}}=2,`$ and

$$A^{\mathrm{ext}}=Q^1=_t,Q^2=\delta x_a_a+\xi ^0(t,u)_t+\eta (t,u)_u\text{and}\eta 0.$$

If the algebra $`A^{\mathrm{ext}}`$ is commutative, then $`\xi _t^0=0`$, $`\eta _t=0`$, and by means of the transformation (8) the operator $`Q^2`$ may be reduced to the form $`\stackrel{~}{Q}^2=2_u\stackrel{~}{\delta }x_a_a`$, where $`\stackrel{~}{\delta }\{0;1\}`$ (the operator $`Q^1`$ does not change in this case). Let us substitute the coefficients of the operator $`\stackrel{~}{Q}^2`$ into the equation (4) and integrate it with respect to $`F`$: $`F=e^{\stackrel{~}{\delta }u}h(u_t)`$, where $`h`$ is a smooth function of the variable $`u_t`$. If $`\stackrel{~}{\delta }=1`$, and $`h`$ is not fixed, then the algebra $`A^{\mathrm{ext}}`$ is truly two-dimensional (the second case of an extension).

If the algebra $`A^{\mathrm{ext}}`$ is non-commutative, we can assume, that $`[Q^1,Q^2]=2Q^1`$, whence $`\xi _t^0=2`$, $`\eta _t=0`$, and by the transformation (8) the operator $`Q^2`$ can be reduced to the form $`\stackrel{~}{Q}^2=2t_t+2u_u+\delta x_a_a`$. Let us substitute coefficients of the operator $`\stackrel{~}{Q}^2`$ to the equation (4) and integrate it with respect to $`F`$: $`F=|u|^{2\delta }h(u_t)`$, where $`h`$ is a smooth function of the variable $`u_t`$. If $`\delta 2`$, and $`h`$ is not specified, then the algebra $`A^{\mathrm{ext}}`$ is truly two-dimensional (the third case of an extension).

Let us assume that $`dimA^{\mathrm{ext}}=3,`$ and

$$A^{\mathrm{ext}}=Q^i=\delta _ix_a_a+\xi ^{0i}(t,u)_t+\eta ^i(t,u)_u,i=1,2,Q^3=_t,$$

the functions $`\eta ^1`$ and $`\eta ^2`$ are linearly independent, and the equation $`\eta ^1(\text{4}.2)\eta ^2(\text{4}.1)`$ is an identity with respect to $`F`$ (here $`(\text{4}.i)`$ is the equation, obtained from (4) by substitution of coefficients of the operator $`Q^i`$). By virtue of the latter condition the coefficients of the operators $`Q^1`$ and $`Q^2`$ have to satisfy the following equations:

$$\begin{array}{c}\eta ^1\eta _t^2=\eta ^2\eta _t^1,\eta ^1\left(\eta _u^2\delta _2\right)=\eta ^2\left(\eta _u^1\delta _1\right),\hfill \\ \eta ^1\left(\eta _u^2\xi _t^{02}\right)=\eta ^2\left(\eta _u^1\xi _t^{01}\right),\eta ^1\xi _u^{02}=\eta ^2\xi _u^{01},\hfill \end{array}$$

from which it follows that $`(\delta _1,\delta _2)(0,0)`$, $`\eta _t^i=0`$, $`\xi _t^{0i}=\delta _i`$, $`i=1,2`$, $`\eta ^1=(\delta _1\delta _2\phi )/\phi _u`$, $`\eta ^2=\phi \eta ^1`$ for a certain function $`\phi =\phi (u)const.`$ Let us substitute, if necessary, the operators $`Q^1`$ and $`Q^2`$ by their linear combinations $`\stackrel{~}{Q}^1`$ and $`\stackrel{~}{Q}^2`$ so $`\stackrel{~}{\delta }_1=1`$, $`\stackrel{~}{\delta }_2=0`$. By the transformation (8) the operators $`\stackrel{~}{Q}^1`$ and $`\stackrel{~}{Q}^2`$ are reduced to the form $`\widehat{Q}^1=t_t+u_u+\delta x_a_a`$, $`\widehat{Q}^2=_u`$. After substitution of the coefficients of the operators $`\widehat{Q}^1`$ and $`\widehat{Q}^2`$ into the equation (4) we obtain one equation for the function $`F`$: $`F_u=0`$, whence $`F=h(u_t)`$ where $`h`$ is a smooth function of the variable $`u_t`$. If $`h`$ is not specified then the algebra $`A^{\mathrm{ext}}`$ is truly three- dimensional (the fourth case of extension).

Lemma 1. In all other cases when $`dimA^{\mathrm{ext}}>1`$, the function $`F`$ satisfies the equation

$$\left(Au_t^2+Bu_t+C\right)F_{u_t}=\left(2Au_t+D\right)F,$$
(9)
where $`A`$, $`B`$, $`C`$ and $`D`$ are arbitrary smooth functions of the variable $`u`$. The transformation (8) is an equivalence transformation at the set of equations of the form (9).

When we integrate the equations (9), four non-equivalent cases arise where $`F_{u_tu_tu_t}0`$.

1. $`F=\alpha (u)e^{u_t}`$ where $`\alpha >0`$. In this case an additional extension of the symmetry is possible only for $`\alpha =A|u|^\nu e^{\mu u}`$, where $`A,\nu ,\mu =const`$, $`A>0`$. By means of scale transformations by variables $`x_a`$ we can get $`A=1`$. If $`\nu =0`$, then we can put $`\mu =0`$ (if it is not so, it is sufficient to perform the transformation $`\stackrel{~}{t}=e^{\mu t}/\mu `$, $`\stackrel{~}{u}=e^{\mu t}(u+2t2/\mu )`$; here and in the following variables $`x_a`$ are not transformed if it is not specifically mentioned). Thus we obtain the fifth case of an extension. For $`\mu \nu 0`$ an additional symmetry extension will arise only for $`\nu =2`$, but then it is possible to put $`\mu =0`$ (the respective equivalence transformation is $`\stackrel{~}{t}=e^{\mu t}/\mu `$, $`\stackrel{~}{u}=e^{\mu t}u`$). If $`\nu 0`$ and $`\mu =0`$, then the equation (1) has the same symmetry algebra as in the more general third case of an extension.

2. $`F=\alpha (u)|u_t+\delta |^{\beta (u)},`$ where $`\delta \{0;1\}`$, $`\alpha >0`$. For additional symmetry extension it is necessary $`\beta =const`$ (and $`\beta \{0;1;2\}`$, or otherwise $`F_{u_tu_tu_t}=0`$). Then any equation (1) with $`\delta =0`$ is equivalent to the same equation in which in addition $`\alpha =1`$ (the sixth case of an extension). The equation (1) with $`\delta 0`$ is also reduced to the same case of extension, if $`(\alpha ^{}\alpha )^{}=0`$. Really, then $`\alpha =Ae^{\mu u}`$ and we can put $`A=1`$ (due to scal transformations by the variables $`x_a`$), $`\mu =0`$ and $`\delta =0`$ (the respective equivalence transformation when $`\mu 0`$ has the form $`\stackrel{~}{t}=e^{\mu \delta t/\beta }\beta /(\mu \delta )`$, $`\stackrel{~}{u}=e^{\mu (u+\delta t)/(\beta 2)}(\beta 2)/\mu `$, and when $`\mu =0`$ it has the form $`\stackrel{~}{t}=t`$, $`\stackrel{~}{u}=u+\delta t`$). If $`\delta 0`$ and $`(\alpha ^{}\alpha )^{}0`$, then we will not obtain new cases of symmetry extension.

3. $`F=\alpha (u)e^{\beta (u)\mathrm{tan}^1u_t},`$ where $`\alpha >0`$, $`\beta 0`$. An additional extension of the symmetry exists only on the conditions $`\beta =const`$ and $`\alpha `$ being an exponential or a power function, but it is not wider than extensions for more general cases 2–4 from the Table 1, and for this reason there is not need to list this case separately.

4. $`F=\alpha (u)u_t^2e^{u_t^1},`$ where $`\alpha >0`$. The equation with such second part is equivalent to the equation with $`F=e^{u_t}`$. The equivalence is determined by the transformation (of the type (7)) $`\stackrel{~}{t}=u`$, $`\stackrel{~}{u}=t+\mathrm{ln}\alpha (u)𝑑u`$.

We completed consideration of the case $`F_{u_tu_tu_t}0`$.

Let in the following $`𝑭_{𝒖_𝒕𝒖_𝒕𝒖_𝒕}\mathbf{=}\mathrm{𝟎}`$, whence

$$F=A(t,u)u_t^2+B(t,u)u_t+C(t,u),$$
(10)
where $`A`$, $`B`$, $`C`$ are smooth functions of the variables $`t`$ and $`u`$. Splitting the equations (3) and (4) by the variable $`u_t`$, in addition to (2) we obtain the following system of determining equations:

$$\xi _a^0=A\xi _t^a\frac{1}{2}B\xi _u^a,\eta _a=C\xi _u^a\frac{1}{2}B\xi _t^a,$$
(11)
$$\begin{array}{c}A_t\xi ^0+A_u\eta +B\xi _u^0=2A\left(\xi _t^0\xi _1^1\right),\hfill \\ C_t\xi ^0+C_u\eta +B\eta _t=2C\left(\eta _u\xi _1^1\right),\hfill \\ B_t\xi ^0+B_u\eta +2A\eta _t+2C\xi _u^0=B\left(\eta _u+\xi _t^02\xi _1^1\right).\hfill \end{array}$$
(12)
The equivalence group of the set of equations (1) with second parts quadratic by $`u_t`$ coincides with the general equivalence group of the equations (1). By means of the transformation (7) the coefficients $`A`$, $`B`$ and $`C`$ are changed as follows:

$$\begin{array}{c}\delta ^2\stackrel{~}{A}=A\zeta _t^2B\zeta _t\zeta _u+C\zeta _u^2,\hfill \\ \delta ^2\stackrel{~}{B}=B\left(\zeta _t\phi _u+\zeta _u\phi _t\right)2A\zeta _t\phi _t2C\zeta _u\phi _u,\hfill \\ \delta ^2\stackrel{~}{C}=A\phi _t^2B\phi _t\phi _u+C\phi _u^2,\hfill \end{array}$$
(13)
The condition $`B^24AC=0`$ is invariant with respect to transformations of the type (13), that are generated by equivalence transformations (7), and for this reason it is a classifying condition. If $`B^24AC=0`$, the we can put $`A=1`$, $`B=0`$, whence $`C=0`$ the seventh case of an extension. In the following $`B^24AC0`$.

Lemma 2. On the condition $`B^24AC0`$ it follows from the equations (2) and (11) that coefficients of any symmetry operator of the equation (1) have the following form:

$$\begin{array}{c}\xi ^a=2\gamma _bx_bx_a\gamma _ax_bx_b+\sigma _{ab}x_b+\beta x_a+\alpha ^a,\hfill \\ \xi ^0=\frac{1}{2}\left(A\beta _t\frac{1}{2}B\beta _u\right)x_ax_a+\left(A\alpha _t^a\frac{1}{2}B\alpha _u^a\right)x_a+\alpha ^0,\hfill \\ \eta =\frac{1}{2}\left(C\beta _u\frac{1}{2}B\beta _t\right)x_ax_a+\left(C\alpha _u^a\frac{1}{2}B\alpha _t^a\right)x_a+\alpha ^4,\hfill \end{array}$$
(14)
where $`\gamma _a`$, $`\sigma _{ab}`$ are constants, $`\beta `$, $`\alpha ^0`$, $`\alpha ^a`$, $`\alpha ^4`$ are smooth functions of the variables $`t`$ and $`u`$.

Let us substitute the expressions (14) into the system (12) and split it by the variables $`x_a`$. As a result we obtain $`n+2`$ systems of the same structure

$$\begin{array}{c}H^1A_t+H^2A_u+H_u^1B=2\left(H_t^1\lambda \right)A,\hfill \\ H^1C_t+H^2C_u+H_t^2B=2\left(H_u^2\lambda \right)C,\hfill \\ H^1B_t+H^2B_u+2H_t^2A+2H_u^1C=\left(H_t^1+H_u^22\lambda \right)B,\hfill \end{array}$$
(15)
where $`H^1`$, $`H^2`$ and $`\lambda `$ take the following values:

$$\begin{array}{ccc}H^1=A\beta _t\frac{1}{2}B\beta _u,\hfill & H^2=C\beta _u\frac{1}{2}B\beta _t,\hfill & \lambda =0;\hfill \\ H^1=A\alpha _t^a\frac{1}{2}B\alpha _u^a,\hfill & H^2=C\alpha _u^a\frac{1}{2}B\alpha _t^a,\hfill & \lambda =2\gamma _a;\hfill \\ H^1=\alpha ^0,\hfill & H^2=\alpha ^4,\hfill & \lambda =\beta .\hfill \end{array}$$
(16)
Investigating the systems of determining equations (15), (16), we obtain the following statements.

Lemma 3. If $`dimA^{\mathrm{max}}>dimA^{\mathrm{ker}}`$, then the equation (1), (10) is equivalent to the same equation in which in addition $`A=e^{\delta t}\widehat{A}\left(u\right)`$, $`B=e^{\delta t}\widehat{B}\left(u\right)`$ and $`C=e^{\delta t}\widehat{C}\left(u\right)`$, and for which $`A^{\mathrm{max}}2_t\delta x_a_a`$ (with functions $`\widehat{A}\left(u\right)`$, $`\widehat{B}\left(u\right)`$, $`\widehat{C}\left(u\right)`$ being not specified, $`A^{\mathrm{max}}=_a,J_{ab},\mathrm{\hspace{0.25em}2}_t\delta x_a_a`$, and it is the first case of an extension).

Lemma 4. If $`dimA^{\mathrm{max}}>dimA^{\mathrm{ker}}+1`$, then the equation (1), (10) is equivalent to the same equation in which in addition $`A_t=0`$, $`B_t=0`$ and $`C_t=0`$.

It follows from the Lemmas 3 and 4 that for completion of the classification it is necessary to investigate the equation (1) with the second parts of the form $`F=A(u)u_t^2+B(u)u_t+C(u)`$, where $`A`$, $`B`$ and $`C`$ are smooth functions of the variable $`u`$, $`B^24AC0`$. The equivalence group of the set of these equations coincides with the general equivalence group of the equations (1), whose second parts do not depend on $`t`$. Under the action of the transformation (8) the coefficients $`A`$, $`B`$ and $`C`$ change as follows:

$$\begin{array}{c}\delta ^2\stackrel{~}{A}=A\widehat{\delta }^2B\widehat{\delta }\zeta _u+C\zeta _u^2,\delta ^2\stackrel{~}{B}=B\widehat{\delta }\phi _u2C\zeta _u\phi _u,\delta ^2\stackrel{~}{C}=C\phi _u^2,\hfill \end{array}$$

and we can assume additionally that $`B=0`$, $`C=\epsilon _1=\pm 1`$. Thus in the following $`F=A(u)u_t^2+\epsilon _1`$ where $`A0`$.

The condition $`A_u=0`$ gives the eighth case of a symmetry extension.

The case when $`\mu =const`$: $`(u+\mu )A_u+2A=0`$, is reduced to the case $`A_u=0`$ by means of the equivalence transformation $`\stackrel{~}{t}=(u+\mu )\mathrm{sinh}t`$, $`\stackrel{~}{u}=(u+\mu )\mathrm{cosh}t`$, if $`\epsilon _1A<0`$, or $`\stackrel{~}{t}=(u+\mu )\mathrm{sin}t`$, $`\stackrel{~}{u}=(u+\mu )\mathrm{cos}t`$, if $`\epsilon _1A>0`$.

The case when $`\mu ,\nu =const`$ $`(\nu 0,2)`$: $`(u+\mu )A_u+\nu A=0`$, is reduced to a more general case 3 from the Table 1 by means of the equivalence transformation $`\stackrel{~}{t}=t`$, $`\stackrel{~}{u}=|u+\mu |^{1\nu /2}`$.

The condition $`(A_u/A)_u=0`$, $`A_u0`$ gives the ninth case of a symmetry extension.

Let the function $`A`$ satisfy the equation $`A_u/A=\nu A+\mu `$ for some non-zero constants $`\mu `$ and $`\nu `$. Any solution of this equation depending on values of the constants $`\mu `$ and $`\nu `$ is equivalent (by translations and scale transformations by the variable $`u`$) to one of the following functions:

$$\frac{\epsilon _2}{\mathrm{cosh}^2u},\frac{\epsilon _2}{\mathrm{sinh}^2u},\frac{\epsilon _2}{\mathrm{cos}^2u},$$
(17)
where $`\epsilon _2=\pm 1,`$ $`(\epsilon _1,\epsilon _2)(1,1)`$. Let $`\epsilon _0=1`$ for the first two functions, and $`\epsilon _0=1`$ for the third function.

Lemma 5. The equation (1) with the second parts $`A\left(u\right)u_t^2+\epsilon _1`$ and $`\stackrel{~}{A}\left(u\right)u_t^2+\stackrel{~}{\epsilon }_1`$, where $`A`$ and $`\stackrel{~}{A}`$ are chosen from the set of functions (17), are equivalent if and only if $`\epsilon _0\stackrel{~}{\epsilon }_0=\epsilon _1\stackrel{~}{\epsilon }_1=\epsilon _2\stackrel{~}{\epsilon }_2`$.

By virtue of Lemma 5 the cases 10–12 of the Table 1 exhaust all possible non- equivalent equations from this class.

For any other functions $`A=A(u)`$ there would be no extension of the symmetry for the equation (1).

Non-equivalence of the all cases of extensions adduced in the Table 1, where it was not proved directly by means of application of equivalence transformations, follows explicitly from non-isomorphness of the respective maximal symmetry algebrae, in particular, from the fact they have different dimensions.

7. Conclusion. We presented a complete solution of group classification problem for a class of PDEs with derivatives entering nonlinearly. New invariant equations obtained in this paper are very interesting also from the point of view of finding new conditional symmetries  for higher order equations, if these new equations are used as differential constraints. In particular, this method works for nonlinear wave equations, and that will be the subject of further papers.


text:
# Optical SETI: A Spectroscopic Search for Laser Emission from Nearby Stars¹

## 1 Introduction

Most searches for signals from extraterrestrial intelligent civilizations (SETI) have concentrated on radio wavelengths (e.g. Tarter 2001, Cullers 2000, Werthimer et al. 2000, Leigh and Horowitz 2000). However, a preferred wavelength for SETI work has proved elusive. The motivation for this project stemmed from the suggestion of Schwartz and Townes (1961) that advanced civilizations might use optical or near–infrared lasers as an efficient means of interstellar communication (Townes 1998, Lampton 2000). The advantages of optical wavelengths over radio include the higher bandwidth and the narrower diffraction–limited beam size to promote private communication. In addition, the technology needed to produce lasers with detectable powers is not far from current Earth technology (§2.3).

Optical SETI (OSETI) projects are currently designed to detect highly energetic, pulsed lasers and continuous–wave lasers. Such efforts include dedicated searches, piggyback searches, and data mining. Many OSETI projects are designed to detect pulsed lasers using fast photometry (Howard et al. 2000, Bhathal 2000, Werthimer et al. 2001). Two or three photodetectors are employed to search for laser pulses having a duration of $``$1 ns in broadband visible light. A high–powered nanosecond pulsed laser will outshine its host star during its brief pulse because at most one photon will be received (for a typical 1-m class telescope) from the star during the 1 ns duration.

Alternately, interstellar optical communication may occur via nearly continuous laser emission. In that case, lasers of even moderate power, tens of kilowatts, can outshine the host star when focussed by large apertures during transmission (see §2.3). Here we report a search for ultra narrow emission lines such as predicted by the Schwartz–Townes suggestion.

## 2 Optical SETI with High Resolution Spectroscopy

We examine high resolution spectra of 577 F, G, K, and M stars, most of which reside within 50 pc. We have taken $``$12,500 spectra using the HIRES echelle spectrometer on the 10–m Keck 1 telescope. All target stars are listed in Nidever et al. (2002). The stellar spectra were obtained during the period 1997–2001 in a search for extrasolar planets using relative radial velocities. We are indebted to Drs. R. P. Butler and S. S. Vogt for their contributions in obtaining these Keck spectra. The multiple spectra from 4000 Å – 5000 Å are compared to detect sharp differences such as those caused by nonnatural emission lines, as described in the following subsections.

### 2.1 High Resolution Echelle Spectra

The acquired echelle spectra contain 34 spectral orders, covering wavelengths 3800 Å– 6100 Å with resolution $`\lambda /\mathrm{\Delta }\lambda `$ = 70,000 (see Vogt et al. 1994, Vogt et al. 2000, Butler et al. 2000). Not all orders are used here because the starlight is sent through a glass cell filled with iodine vapor before entering the spectrometer for use in setting the wavelength scale in the planet search. The iodine lines contaminate the stellar spectrum at wavelengths greater than 4980 Å. The remaining usable orders, 72 through 89, contain 1000 Å ranging from 3980 Å – 4980 Å. This region includes Balmer $`\beta `$ and $`\gamma `$, but not the CaII H&K lines. The pixel width is given by $`\lambda /\mathrm{\Delta }\lambda `$ = 145,000 , corresponding to 2.1 km s$`^1`$per pixel or $``$0.032 Å per pixel. The actual spectral resolution is only 70,000 because the instrumental profile has a width of 2.3 pixels.

The recorded two–dimensional raw spectra are reduced in the usual way to one–dimensional spectra. The bias is subtracted from all images including both the stellar spectra and the flat–field exposures. Each spectral order recorded by the CCD has a width proportional to the diameter of the stellar seeing disk along the length of the slit. The typical seeing at Keck/HIRES of 0.7–1.2 arcsec yields orders of width $``$4 CCD rows (each “row” being the on–chip sum of neighboring original CCD rows). Cosmic rays are removed by searching for 7–$`\sigma `$ sharp intensity spikes at one row compared to the intensity in a neighboring row (at the same wavelength). We are indebted to Jason Wright for this cosmic–ray removal algorithm. The stellar spectra are divided by the flat–field exposures. The two–dimensional spectral orders are compressed to one dimension by summing the recorded counts within the $``$10 rows centered on each order. We thank Eric Nielsen for the careful reduction of CCD images. Figure 1 shows the results of this process. The signal–to–noise ratio per pixel in orders 72 through 89 of a typical 1-D spectrum is $``$200.

### 2.2 Overview of the Detection Scheme

We search for laser emission lines in the reduced spectra by comparing each spectrum to a reference spectrum of the same star, searching for statistically significant differences. Over time, a star’s spectrum should remain nearly constant except for minor fluctuations due to photon noise and photospheric changes. A large difference at one wavelength could signify the presence of a laser signal that had been detected during one observation but not during the other.

We would fail to detect an emission line if it occured at the same wavelength in both the observed and reference spectra. A laser fired continuously from the rest frame of the stellar photosphere would go undetected. However, such a continuous laser would be detected in both spectra at different wavelengths if the laser’s velocity vector were to change between observations, e.g. due to orbital motion around its host star. The resulting Doppler shift permit detections in the difference spectrum.

For each star, a series of algorithms is applied to all of its spectra. First a “reference” spectrum is chosen. Then, one by one, the other “test” spectra are aligned in wavelength with the reference for comparison (see §3.2). Candidate laser lines are found in the test spectra by subtracting the reference spectrum and searching for positive differences, or “spikes.” One test spectrum is then used as a secondary reference to check for spikes in the primary reference spectrum. The spikes are subjected to a selection criterion based on their width, to ensure that their shapes are consistent with the instrumental profile of HIRES (see §3.3). Spikes that are narrower than the instrumental profile are likely to be caused by cosmic rays hitting the CCD that eluded detection in the CCD reduction process described in §2.1.

A monochromatic laser line would have an intrinsic width narrower than the resolution of the spectrometer. However, as the laser light passes through the HIRES optical system, it would be smeared out by the HIRES point spread function (PSF), giving the detected emission line an apparent distribution in wavelength characteristic of the HIRES PSF. Thus candidate laser lines can be distinguished from most cosmic rays, most CCD flaws, and stellar emission lines by their predictable spread in wavelength that would be merely the HIRES PSF, as described in detail by Valenti et al. (1995). Notably, emission lines having widths that are narrower than possible from a star (with its thermal, collisional, turbulent, and rotational broadening mechanisms) would stand out as possibly “nonastrophysical” in origin. All of the resulting candidate emission lines of any width, excluding the apparent cosmic rays hits, from the spectra of each star are stored for further analysis as described in §4. In sections §3.1, §3.2, and §3.3, we describe in detail the processes by which we searched for emission lines.

### 2.3 Detection Limits for Laser Emission

To be detected, a laser signal must compete against the light of its host star. A laser beam carries three advantages over starlight. First, the laser light is concentrated in a narrow cone of tiny solid angle relative to starlight, which is omnidirectional. Second, a monochromatic laser line falls into one resolution element, with a width of $``$2.3 pixels, while starlight, composed of a wide range ($``$4000 Å) of wavelengths, is spread out into roughly 50,000 resolution elements. Lastly, the laser needs only to outshine the photon noise of the stellar flux, not the flux itself, within that resolution element since we are looking at differences in the spectra. To estimate these advantages of the laser signal, we consider the solid angles of the starlight and a diffraction–limited laser beam.

$`\mathrm{\Omega }_{\mathrm{star}}`$ $`=`$ $`4\pi `$ (1)
$`\mathrm{\Omega }_{\mathrm{laser}}`$ $``$ $`\left({\displaystyle \frac{\lambda }{a}}\right)^2`$ (2)
Here $`\lambda `$ is the wavelength of the laser light and $`a`$ is the aperture size of the laser transmitter. The fraction of starlight contained within the laser beam is then given by

$$f_{\mathrm{star}}=\frac{1}{4\pi }\left(\frac{\lambda }{a}\right)^2$$
(3)
We note that the intrinsic width of a laser line will be much narrower than the HIRES resolution of $``$0.075 Å, thereby placing most all of the laser light into a single resolution element. The stellar spectral energy distribution may be approximated by a blackbody curve. For a Solar–type star, the SED has a width of $``$4000 Å, yielding the number of resolution elements, $`N_{\mathrm{star}}`$, into which the stellar luminosity falls:

$$N_{\mathrm{star}}\frac{4000\text{ Å}}{0.075\text{ Å}}=5.3\times 10^4$$
(4)
The fraction of starlight contained within one resolution element is therefore given by

$$f_{\mathrm{star},\text{ }\mathrm{\Delta }\lambda }=\frac{1}{4\pi }\left(\frac{\lambda }{a}\right)^2\frac{1}{5.3\times 10^4}$$
(5)
The laser light needs only to outshine the stellar photon noise since we are detecting differences in a star’s spectrum over time. The photon noise relative to the flux itself, $`\sigma `$, for one spectrum is given by 1/(S/N) where S/N is the signal to noise ratio per resolution element. Since the laser must compete with the photon noise of two spectra in one resolution element, we have

$$\sigma \frac{\sqrt{2}}{\text{S/N}}$$
(6)
Finally we find that the laser power necessary for an $`n\sigma `$ detection is

$$P_{\mathrm{laser}}=\frac{1}{4\pi }\left(\frac{\lambda }{a}\right)^2\frac{1}{5.3\times 10^4}\frac{n\sqrt{2}}{\text{S/N}}L_{\mathrm{star}}$$
(7)
where $`L_{\mathrm{star}}`$ is the luminosity of the star. As a representative case, we take $`\lambda `$=4500 Å, S/N=200, n=6, $`L_{\mathrm{star}}=L_{\mathrm{sun}}=3.8\times 10^{26}`$ W, and assume Keck to Keck communication so $`a`$=10 m, which gives

$$P_{\mathrm{laser}}=50\text{ kW}$$
(8)
Thus for this representative case, one can detect a 50 kW laser.

The previous calculation assumes that the laser is on continuously. If, however, a nanosecond pulsed laser signal were transmitted, Equations 7 and 8 would need to be multiplied by the factor $`\frac{\mathrm{\Delta }t_{\mathrm{exp}}}{N_{\mathrm{pulses}}\times 10^9\text{ s}}`$ where $`\mathrm{\Delta }t_{\mathrm{exp}}`$ is the exposure time of the observation and $`N_{\mathrm{pulses}}`$ is the number of pulses received in that time. Using a typical exposure time of 10 minutes and assuming 1 pulse is received, the minimum necessary power of a detectable nanosecond pulsed laser is $`P_{\mathrm{pulsed}}=3\times 10^{16}`$ W. Thus, the minimum necessary energy of a detectable nanosecond pulsed laser is 30 MJ.

For comparison, the most powerful continuous laser on Earth is the Free Electron Laser (FEL) developed at Jefferson Laboratory in Virginia. The FEL has reached an average power of 1.7 kW and the laser is expected to achieve average powers up to $``$10 kW with an upgrade (Jefferson Lab 2001). Lawrence Livermore National Laboratory produced the world’s most powerful pulsed laser, the Petawatt laser. During its 440 femtosecond pulse, it delivered a power greater than 1015 W (LLNL 2001).

We can also determine a minimum detectable flux, $`F_{\mathrm{laser}}`$,

$$F_{\mathrm{laser}}=\frac{P_{\mathrm{laser}}}{A_{\mathrm{laser}}}$$
(9)
where $`A_{\mathrm{laser}}`$ is the area of the laser beam at the Earth and given by

$$A_{\mathrm{laser}}=\frac{\pi }{4}\left(\frac{\lambda }{a}\right)^2d^2,$$
(10)
where $`d`$ is the distance from the laser to Earth. The flux limit can now be expressed as

$$F_{\mathrm{laser}}=\frac{4}{\pi }\frac{1}{5.3\times 10^4}\frac{n\sqrt{2}}{\text{S/N}}\frac{L_{\mathrm{star}}}{4\pi d^2},$$
(11)
where $`L_{\mathrm{star}}/4\pi d^2`$ is the flux of the star. If we take $`d=100`$ light years with the previously adopted values, the minimum flux of a detectable laser is

$$F_{\mathrm{laser}}=3.4\times 10^{17}\text{ W/m}\text{2}$$
(12)
Thus, we can detect monochromatic laser light from nearby stars if the flux is at least $`3.4\times 10^{17}`$ W/m2 at Earth.

## 3 Spectral Analysis

A detailed description of the algorithms developed to analyze the spectra is given in this section.

### 3.1 Choosing Reference Spectra

Two reference spectra are required for each star. The primary reference spectrum is compared to all the other test spectra. The second is compared to the primary reference spectrum. Since only positive differences are considered, any candidate laser lines in the first reference would not be flagged during its comparison with the test spectra. A second reference, one of the previous test spectra, is used to account for this possibility.

The two reference spectra are chosen as follows. Almost every star has at least one “template” spectrum. In the planet–search project, the stellar radial velocities are measured with respect to the template spectra. These observations are of good quality with high signal–to–noise ratios. If there are multiple–template spectra for a star, the two with the highest signal to noise are chosen as the reference spectra. If there is only one template, it is used as the primary reference spectrum and another observation with high signal to noise is chosen as the second reference. Some stars which are new to the extrasolar planet search catalog do not yet have templates. For these stars, the two observations with the highest signal to noise are chosen as the references.

### 3.2 Aligning Spectra

Aligning a pair of spectra precisely in both wavelength and flux is important. If the spectra are exactly the same, except for minor fluctuations due to photon noise, the fractional standard deviation, $`\sigma `$, of the difference spectrum is expected to be

$$\sigma =\sqrt{\frac{1}{\text{(S/N)}_1^2}+\frac{1}{\text{(S/N)}_2^2}}$$
(13)
where (S/N)1 and (S/N)2 are the signal–to–noise ratios of the two spectra. For two spectra with signal–to–noise ratios of 300, $`\sigma `$ can, in principle, be as low as 0.005. It is important to align the two spectra well enough such that the observed error is as close to the theoretical error as possible in order to detect the smallest differences.

We line up a pair of spectra order by order. Various normalization and shifting algorithms are applied to both the reference and test spectra, as follows. The first step is to put the orders on roughly the same flux scale by dividing by the continuum. The continuum is found by first partitioning an order into 10 bins. Each bin is assigned a temporary continuum at the 85th percentile in flux. A third–order polynomial is fit to these individual continua to normalize the entire spectral order.

Each spectrum is Doppler–shifted to correct for the Earth’s orbital motion, which causes a maximum wavelength shift of $``$1 Å corresponding to $``$30 pixels. We determine the Doppler shift of the test spectrum relative to the reference spectrum by minimizing $`\chi ^2`$, and then shift the test spectrum so that the two are aligned. The final shifts are precise to $``$0.01 pixel, and fractional pixel displacements of the spectra are accomplished with a spline interpolation. This shifting is first carried out on segments of length 280 pixels, and then it is done on a 40–pixel scale (“chunks”). If any spectrum contains cosmic rays (or laser emission lines), it is difficult to align. Such features are removed temporarily (replaced with a spline interpolation) to carry out the alignment. Operating on chunks of 40 pixels allows us to account for the wavelength dependence of the Doppler shift and to correct any slopes in the flux with wavelength.

The resulting (temporarily spike–free) reference and test spectra are subtracted to yield a difference spectrum, from which the standard deviation is determined and recorded as the “noise” for each spectral chunk. Finally the reference and test spectra are subtracted with all spikes included, leaving a difference chunk which can be searched for laser–line candidates. Figure 2 shows a chunk of difference spectrum both with (top) and without (bottom) flux spikes of unidentified origin.

### 3.3  Criteria for Laser Line Candidates

A laser beam would travel through the telescope and the HIRES spectrometer optics just as starlight would. Moreover, the laser beam travels through the same column of Earth’s atmosphere, acquiring the same wavefront distortions and seeing profile as does the star. Therefore, an emission line from a laser would span the whole width of a spectral order (perpendicular to dispersion) as would the stellar spectrum on the CCD. While the intrinsic width of the laser line in wavelength would be smaller than the resolution of HIRES, it would be broadened by the point–spread function of HIRES. The B1 slit used in HIRES projects to a width of 2.0 pixels on the CCD. The instrumental profile of HIRES from both its optics and the slit width is known from careful measurements (Valenti et al. 1995) and from simple Gaussian fitting of emission lines from a thorium-argon lamp. The instrumental profile with our HIRES setup has a width of 2.3 pixels (FWHM). Thus, we expect any laser line to have an observed width of 2.3 pixels.

We search for candidate emission lines in the stellar spectra as follows. Each 40–pixel chunk of difference spectrum has a corresponding fractional photon noise, $`\sigma `$, associated with it as described in §3.2. We pass through all the chunks and note spikes having heights above $`6\sigma `$. This threshold is chosen to exclude essentially all noise but detect a signal from a laser with the lowest detectable power.

We also require that candidate laser lines exhibit a width that is at least as broad as the HIRES instrumental profile, with FWHM=2.3 pixels, allowing us to discriminate against (i.e., reject) cosmic–ray hits. Most cosmic rays produce electrons in a single pixel or occasionally in two neighboring pixels, clearly narrower than the HIRES PSF. To implement this “width” criterion, we require that the shorter of the two adjacent pixels on either side of each 6–$`\sigma `$ spike contain at least 30% of the photons in that central spike. Examination of the known HIRES PSF and thorium comparison lines (see Figure 3) shows that this width criterion would be met by almost any emission line that had been broadened by the HIRES PSF. Any spike in the difference spectrum that does not meet both the 6–$`\sigma `$ height criterion and this width criterion is rejected from future consideration as they are probably cosmic rays. Spikes that do meet both criteria are stored in a database that includes the location and characteristics of that spike, for later statistical analysis (§4). We note that emission lines having widths greater than the PSF will be detected and stored as candidate signals. Thus our survey would detect intermittent emission lines of any origin, laser or astrophysical, if they are sufficiently strong.

## 4 Statistical Analysis of Candidates

The spectral analysis described in §3 yields a database of surviving candidate emission lines in each of the $``$12,500 Keck spectra of our 577 target stars. A typical star has $``$13 surviving candidate emission lines in its $``$20 spectra. For each candidate, the database contains the star, date of observation, the wavelength, the peak number of photons, and the width of the candidate line. Detailed examination shows that most (if not all) of the surviving candidate emission lines are simply cosmic ray hits that we had failed to reject.

To discriminate a true emission line from these cosmic rays, we apply a statistical analysis on the collection of emission–line candidates from all spectra (typically 20) for a given star. We search for an overabundance of candidates (“hits”) in the vicinity of any particular wavelength. In contrast, cosmic rays would hit at random locations on the CCD, and would not favor any particular wavelength. A clustering of hits near one wavelength from all the spectra of one star could indicate laser emission.

We note that the Earth’s orbital motion will not cause a Doppler shift of the laser line in our analysis, as we have already forceably shifted the reference and test stellar spectra onto the same wavelength scale (§3). However, acceleration of the emitter relative to the star (such as due to orbital motion) during the typical 4 years of observations would cause the laser line to be Doppler–shifted to neighboring wavelengths.

We make a grid representing the CCD with the spectral orders partitioned into approximately 30 segments of roughly 100 pixels each. The segments overlap one another by half of their size to account for the possibility of an abundance of hits split into two neighboring segments. We record the number of emission–line candidates in each segment.

For each star, we plot the number of 100–pixel segments against the number of hits in the segments. Most segments have no hits at all, and only a few have one hit. These infrequent hits (none or one hit out of $``$20 observations in a wavelength segment) are consistent with random cosmic ray hits. To this histogram, we fit a Poisson distribution (see Figure 4), and determine the probabilities of finding an excess number of hits, $`n`$, in any wavelength segment. A star exhibiting an excess number of hits in a segment, i.e. a low probability of the events occurring randomly, would stimulate further investigation as a potential emission line.

The fit of a Poisson distribution to the histogram of hits is done as follows. We fit the low–$`n`$ (0, 1, 2) portion of the histogram and neglect the high–$`n`$ tail. This approach avoids contaminating the Poisson mean, $`\mu `$, with actual emission lines in the high–$`n`$ tail. The first empty bin in the histogram, $`n_0`$, is determined and an estimate of the expected mean number of hits per segment, $`\mu `$, is calculated neglecting those segments with $`nn_0`$.

The Poisson probability of finding $`n`$ cosmic ray hits in any given wavelength segment is given by (Taylor 1982):

$$P_\mu (n)=e^\mu \frac{\mu ^n}{n!}$$
(14)
Using the estimated $`\mu `$, we calculate $`P_\mu (n)`$ for $`n=0,1,\mathrm{},n_01`$, and check if any of the probabilities are less than 1%. If so, we exclude those bins from the fit as well if $`n_0>2`$. The fit is then applied to bins from $`n=0`$ to some $`n_{\mathrm{max}}`$ hits.

The estimated mean number of hits per segment may not necessarily yield a Poisson distribution that gives the best fit to the observed histogram. The best Poisson fit may be better determined by finding the lowest $`\chi ^2`$ given by

$$\chi ^2=\underset{n=0}{\overset{n_{\mathrm{max}}}{}}\frac{(O_nE_n)^2}{E_n}$$
(15)
where $`O_n`$ is the observed number of segments containing $`n`$ candidate hits and $`E_n`$ is the expected number with $`n`$ hits, based on a Poisson distribution. The only free parameter is $`\mu `$, the mean number of hits per segment. We vary $`\mu `$ to find the value corresponding to the minimum of $`\chi ^2`$. The resulting Poisson probabilities are multiplied by the total number of segments and adopted as the best fit to the observed histogram.

We now determine the probabilities of finding $`n`$ hits in any given segment, extrapolating beyond $`n_{\mathrm{max}}`$. Segments having an excessive number of hits with probabilities less than 0.01% are noted as serious candidate emission lines, as they occur in excess numbers above the background cosmic ray hits. Since there are about 500 wavelength segments on our grid of the CCD for each star, we expect to flag one in 20 stars spuriously, i.e. due to a chance coincidence of cosmic rays.

We find that indeed 45 stars were flagged as having a wavelength segment with an excess number of hits at the 0.01% false-alarm level. Each of the surviving flagged candidate laser lines, of which there were 187, were examined by eye in the raw CCD images.

## 5 An Apparent Emission–Line Detection

Our search for intermittent emission lines revealed an apparent detection for the star HD29528 (K0V). Three separate spectra obtained on different dates contained strong candidates at the same wavelength. The hits were noted in both of two overlapping spectral segments for this star. The probability that these three apparent emission lines were merely spurious cosmic rays occurring within a segment of length $``$200 km s-1 is easily calculated from the Poisson occurrence of cosmic rays (§4).

Figure 4 shows the observed histogram of the number of segments versus number of hits along with the best fitting Poisson distribution. The probability of randomly obtaining 3 hits in a segment for this star is 1.23 $`\times `$ 10-7 or 1 in 8 million.

Figure 5 shows these candidate laser lines in the reduced spectra of the three observations. The aligned reference and test spectra are shown and twice the difference spectra are below them. The candidate laser lines are in spectral order 83 and the central peaks are in pixel 1436 with respect to the reference spectrum.

We examined the raw spectra to see if the candidate laser lines span the whole width of the spectral order, indicating whether the light passed through the telescope and spectrometer. Figure 6 shows the raw CCD image of the spectrum of HD29528 obtained on JD 2451793.1, centered on the first candidate laser line. The reference spectrum is shown below it. Beneath the reference is the raw difference spectrum making it easier to see the apparent laser line. The other two candidate laser lines found in observations obtained on JD 2451882.98 and on JD 2451884.08 are not shown, but those apparent emission lines also span the width of the spectral order. Because all three hits occur at the same pixel (and hence wavelength) with respect to the star’s reference spectrum, and because they span the entire width of the spectral order, the hits are unlikely to be cosmic rays. These are surviving laser line candidates thus far in the analysis.

To further investigate these candidate emission lines in HD29528, we compared the spectrum of HD29528 to the solar spectrum using “The Solar Spectrum 2935 Å to 8770 Å” (Moore et al. 1966). The rest wavelengths of the lines are 4307.31 Å and there are no listed absorption lines due to the Earth’s atmosphere nor anomalous Solar emissions lines at this wavelength.

We further examined the spectrum of a rapidly rotating B–type star taken on the same night as the reference spectrum for HD29528. The B star spectrum should have no prominent spectral features. A strange drop in flux of 20% is seen in the spectrum, exactly at pixel 1436 in order 83, precisely where the hits were found with respect to the reference spectrum of HD29528 (see Figure 7). This clearly spurious feature in the B star strongly suggests that our candidate laser lines are merely artifacts caused by a flaw in the CCD.

Indeed, a close examination of the raw CCD image of the reference spectrum of HD29528 shows the flaw in the CCD (see Figure 8). This flaw is certainly the source of the false detection.

As a side note to this false detection of an emission line, one may calculate the necessary power of a laser that could have produced a “signal” similar to that detected. In Eqn 7 we assume a transmitter aperture, $`a`$=10 meters. The other variables are known. The wavelength of the 13$`\sigma `$ “signal” is 4307 Å, and the signal–to–noise ratio of the spectra is $``$120. The luminosity of a K0 star is 0.42 $`L_{\mathrm{sun}}`$ (Carroll and Ostlie 1996). Using these values in equations 7, we find $`P_{\mathrm{laser}}=68`$ kW. This would be the required power of a laser to produce the spurious signal in Figure 6. This test supports our expectation that lasers having power of 50 kW would have been detected.

## 6 Results

For each of the 45 stars that exhibited an excess of candidate emission lines (above the cosmic ray background noise) we examined the original raw CCD images from HIRES. Approximately 200 raw spectra were examined to follow up on the 45 flagged stars. Stellar emission lines were detected in the spectra of two flare stars, HIP5643 (M4.5) and HIP92403 (M3.5). Both stars exhibited varying emission of the Balmer line H$`\beta `$.

Interestingly, many candidate emission lines were located along the steep slopes of absorption lines in the reduced spectra. However, the corresponding emission spikes were not apparent in the raw CCD image. Such emission spikes are apparently neither cosmic rays nor emission lines. Instead, it is likely that these spikes located in absorption lines are caused simply by tiny changes in the instrumental profile of HIRES from one observation to the next. The changes in the instrumental profile could result from different atmospheric seeing of the star image at the slit (making the slit effectively a different size) or from small focus changes in HIRES itself. The steep sides of absorption lines have local intensities that are very sensitive to the instrumental profile. It is also possible that the spikes located in absorption lines are residual numerical artifacts from the spectral alignment process since interpolation routines normally give the largest errors on steep slopes. Thus these candidate emission spikes located within absorption lines were rejected as spurious.

All other candidate emission lines were visually inspected and found to be inconsistent with actual emission lines. Some were multiple, surviving cosmic rays that occurred, by chance, near the same wavelength. Some emission spikes were caused by bad pixels or flaws in the CCD. No lines consistent with laser emission were found in any of our stars.

## 7 Discussion

We would have detected time–varying emission lines from any of the 577 main sequence stars above a threshold of a few percent of the continuum flux. Emission lines of constant intensity would also have been detected had they Doppler–shifted relative to the star during the $``$4 yr of observations. Lines having intrinsic widths from arbitrarily narrow to tens of km s$`^1`$would have been detected. No narrow lines consistent with laser emission were detected. However, naturally occurring emission lines were found in the spectra of two stars.

Emission lines in the spectra of Solar–type stars are rare between 4000 Å – 5000 Å. Such lines can arise occasionally in the chromospheres or coronae of the most magnetically active stars. In some ($``$20%) F, G, K–type main sequence stars, the cores of the Balmer absorption lines are “filled in” by a few percent due to chromospheric emission (Herbig 1985, Robinson, Cram, & Giampapa 1990). We would have detected temporal variation in that chromospheric filling, if it exceeded a few percent of the continuum level. No such variations were detected in any of our F, G, K–type main sequence stars. This nondetection indicates that such chromospheric variation in the cores of the strong Balmer lines, notably H$`\beta `$ and H$`\gamma `$, did not exceed a few percent of the continuum level in these stars. Flare stars, most of spectral type M0 or later, more commonly exhibit emission of H$`\beta `$. Indeed, we detected time–varying emission of H$`\beta `$ in two M–type flare stars.

Time–varying emission lines from stars may also arise from infalling volatile material and comets, as is observed for $`\beta `$ Pic (see for example Beust et al. 1998). Our search would have detected such astrophysical emission lines if they appeared or varied by more than a few percent of the continuum flux of the star. All such astrophysical origins for emissions lines would yield line widths greater than our resolution of 5 km s$`^1`$due to the usual thermal and turbulent broadening, as well as possible rotational and collisional broadening. Hence such lines would be resolved. Nonetheless, no emission lines due to these effects were found in the 577 stars at our threshold of a few percent of continuum flux.

The entrance slit of the HIRES spectrometer is rectangular with dimensions 0.57 $`\times `$ 3.5 arcsec. At the typical distance of our stars of $``$50 pc, this slit corresponds to 30–175 AU at the star. Thus, sources of emission lines residing within 30 AU of our target stars would have been included in our spectra. Of course, the emission beam must be pointed toward the Keck telescope for detection. Extinction from interstellar dust is negligible toward these target stars due to their proximity ($`d<`$100 pc).

The motivation for this search stemmed from the suggestion of Schwartz and Townes (1961) that advanced civilizations might communicate from one colony to another within a planetary system, or from one planetary system to another, by using optical lasers as the carrier. The advantages of optical communication over radio include its higher bandwidth and its smaller diffraction beam to promote private communication.

A nondetection of extraterrestrial intelligence carries no value unless it excludes a plausible model of life in the universe. One might imagine the following model. The Milky Way Galaxy has an age of $``$10 Gyr, while only 4.6 Gyr were required to spawn a technological species on Earth. Approximately half of the stars in the Galactic disk are older than the Sun, and they formed with comparable amounts of heavy elements. Approximately 50% of these old disk stars do not have a stellar companion within 100 AU (Duquennoy and Mayor 1991), thereby permitting stable planetary orbits. Therefore, $``$25% of the stars in the Galactic disk have requisite characteristics of adequate age, chemical composition, and dynamical quiescence to serve as sites for the development of technological species. The anthropocentric nature of these stellar criteria are apparent, and thus serve merely as a useful guide for a model.

Coincidentally, the target stars used here are systematically middle–aged or older and void of close companions. They were selected for the Doppler planet search at Keck from among stars in the Solar neighborhood. They were chosen to have ages greater than 2 Gyr, as judged from the CaII H&K chromospheric emission line, because of the resulting photospheric stability of older stars. Moreover the target stars have no known stellar companions within 2 arcsec (typically $``$50 AU) to prevent any companion’s light from entering the spectrometer slit (but see Vogt et al. 2002 and Nidever et al. 2002 for a few recently detected, but negligibly bright, close companions). Moreover, spectroscopic analyses of the stellar spectra show that 90% of the target stars contain amounts of heavy elements that are within a factor of 2 of the Solar value (\[Fe/H\] = -0.3 to +0.3; D.A.Fischer, private communication).

Thus, the majority of the 577 target stars for this SETI search are older than the Sun, have roughly Solar chemical composition, and have no perturbing stellar companions within 50 AU. This stellar sample is clearly biased. Thus our observations test a limited model of the Galaxy in which older, single stars harbor technological civilizations.

The speculative model posits that some fraction of these target stars contain civilizations that sent probes to, or established colonies on, planets, moons or other platforms either in their host planetary system or in other systems. This model of behavior, albeit anthropocentric, stems from the tendencies of Homo Sapiens toward exploration. To communicate with such probes or colonies over distances of Astronomical Units to parsecs, some fraction of these civilizations might be modeled as using electromagnetic waves, including optical lasers.

We model the solid angle subtended by laser beams as follows. Some fraction, $`f_{\mathrm{tech}}`$, of the total number of target stars, $`N_{\mathrm{stars}}`$=577, harbor technological civilizations. On average each such civilization emits some number of laser beams, Nbeams, perhaps in arbitrary directions. We consider that each beam subtends an average solid angle $`\mathrm{\Omega }_{\mathrm{beam}}`$.

The total solid angle, $`\mathrm{\Omega }_{\mathrm{Total}}`$, subtended by all arbitrarily oriented beams is

$$\mathrm{\Omega }_{\mathrm{Total}}=N_{\mathrm{stars}}f_{\mathrm{tech}}N_{\mathrm{beams}}\mathrm{\Omega }_{\mathrm{beam}}$$
(16)
We ignore here details of overlap of the beams, their actual distribution of power and beam size, and any purposeful aiming toward or away from us. Any properties of the lasers that are constructed with the Earth as a consideration could drastically change the probability of interception of such beams.

In our survey, we have monitored $`N_{\mathrm{stars}}`$= 577 stars. The nondetection obtained in the survey here suggests that the total solid angle subtended by lasers emitting 50 kW or more is constrained as $`\mathrm{\Omega }_{\mathrm{Total}}<4\pi `$ sr .

Thus for arbitrarily directed beams, we constrain the model as

$$f_{\mathrm{tech}}N_{\mathrm{beams}}\mathrm{\Omega }_{\mathrm{beam}}<4\pi /577$$
(17)
The nondetection reported here suggests lasers emitting detectable power must satisfy the above constraint for the combination of $`f_{\mathrm{tech}}`$, $`N_{\mathrm{lasers}}`$, $`\mathrm{\Omega }_{\mathrm{beam}}`$. Apparently, advanced civilizations communicating by kilowatt optical lasers are not so common as to fill 4$`\pi `$ sr.

As a touchstone, we consider the case of the diffraction–limited laser with an aperture of 10-m (Keck–to–Keck), giving a solid angle for one beam, $`\mathrm{\Omega }_{\mathrm{beam}}10^{15}`$ sr. For transmitters of such aperture, a laser emitting 50 kW would be detectable. For this case, each star would have to emit $`10^{13}`$ arbitrarily oriented beams on average, each with power $`P`$ 50 kW, in order that one would likely be oriented, serendipitously, toward us.

Clearly, the observed nondetection implies that such numerous, narrow beams, however implausible, do not exist in sufficient numbers to permit likely detection. These observations offer a poor constraint on lasers of arbitrary orientation, as their solid angles are simply too small in this model. One might consider a Galactic model consisting of wider laser beams, enhancing the chances of interception. In that case, the present nondetection would impose a constraint of fewer lasers emanating from the target stars, but would require proportionally higher power per laser for detection at 6–$`\sigma `$ above the stellar continuum.

More meaningfully, the nondetection reported here suggests that none of the 577 target stars harbors a civilization that has purposefully directed a laser with $`P`$50 kW toward Earth (for the case of a 10–meter diffraction–limited transmitting aperture). This modest detection threshold of 50 kW highlights the ease with which an advanced civilization could signal us, if desired. But such is apparently not the case.

## 8 Acknowledgements

We thank the SETI Institute for generous funding of this work. We wish to thank R.Paul Butler and Steve Vogt for their endless hours in obtaining the Keck/HIRES spectra. S.Vogt’s design of HIRES made this work possible. We thank Jason Wright for the new CCD reduction algorithm of Keck spectra that efficiently removed cosmic ray hits. We thank Eric Nielsen for reducing the 12,500 CCD images to one–dimensional spectra. We thank Debra Fischer for chemical–abundance measurements of the target stars. We also thank D.Nidever, D.Werthimer, J.Tarter, F.Drake, M.Lampton, B.Walp and referee, William D. Cochran for many useful comments. This work was carried out by A.R. as partial fulfillment of the of the requirements toward a Master’s Degree in physics at San Francisco State University.