## Abstract

We derive an explicit formula for the singular part of the fundamental matrix of crystal optics. It consists of a singularity remaining fixed at the origin *x*=0, of delta terms located on the positively curved parts of the wave surface, the well-known Fresnel surface and of a Cauchy principal value distribution on the negatively curved part of the wave surface.

## 1. Introduction and notation

The theory of partial differential equations has evolved in close concurrence with specific examples coming from physics. A prominent case in the history of science was the propagation of light in crystals investigated by such great physicists as Huygens, Fresnel, Biot, Brewster and Laplace. In particular, the prediction of *conical refraction* in biaxial crystals by Hamilton in 1832 and its experimental verification shortly thereafter by Lloyd fascinated physicists and mathematicians alike, triggering a long-lasting search for a mathematical framework encompassing and explaining such phenomena.

This search was advanced by scientists like Navier, Cauchy, Green, Stokes, Lamé, Kovalevskaya, Volterra, Grünwald, Zeilon, Herglotz, etc., cf. the historical account in Gårding (1989), culminating in the modern theory of *hyperbolic differential operators*, initiated by Hadamard, and elaborated by, among others, Petrovsky, Riesz, Gårding and Hörmander. This theory allows to analyse the singularities and lacunas of solutions of hyperbolic operators by means of *localizations* and *microlocal analysis*, i.e. *wavefront sets*, see Petrovsky (1945), Hörmander (1969, 1983*a*,*b*), Atiyah *et al.* (1970, 1973) and Gårding (1984).

Parallel to this qualitative study, several methods were devised for the calculation, or representation by integrals, of *fundamental solutions* of differential operators by Laplace, Poisson, Fourier, Stokes, Zeilon, Fredholm, Herglotz and others. Placed on a sound mathematical basis by Schwartz’s invention of distribution theory around 1950, such integral representations of fundamental solutions were further generalized by Petrovsky, Leray, Gel’fand, Shilov, etc., and became known under the name of *Herglotz–Petrovsky formulae*, cf. the historical sketch in Wagner (2004). Let us mention that explicit formulae for fundamental solutions of irreducible hyperbolic cubic and quartic operators in three and four dimensions were found only recently (see Wagner 1999, 2005, 2006, Ortner & Wagner 2008).

In the present paper, we investigate the system of crystal optics with the purpose of representing explicitly the singular terms in its fundamental matrix. Strictly speaking, one could deduce the main theorem in §10 by combining results from various sources dispersed over the literature. We aim here, however, at giving a unified approach even at the expense of reproducing some earlier calculations.

Let me shortly survey the setup of this article. We first recall, in §2, Maxwell’s equations and the resulting 3×3 systems of second-order partial differential equations, which govern the electric and the magnetic field vectors, respectively, in crystal optics, cf. Courant & Hilbert (1962). As in Gårding (1989), we focus on the system for the magnetic field; the equations for the electric field, see, e.g. Burridge & Qian (2006), could be treated analogously. In §§3 and 4, we first briefly review part of the theory of Atiyah, Bott and Gårding on the propagation of singularities for solutions of hyperbolic equations, and we then specify what this theory implies in the context of crystal optics. This approach, which largely parallels Esser (1987) and Gårding (1989), leads to the investigation of Fresnel’s surface and of the Hamiltonian circles on it.

In §5, we re-derive Grünwald’s formula (5.3) for the fundamental matrix of a strictly hyperbolic system in **R**^{4} (see Grünwald 1904; Ortner & Wagner 2004). In contrast to crystal acoustics, one of the characteristic speeds vanishes in crystal optics, and this implies that the slowness surface decomposes into the plane at infinity *τ*=0 and the quartic Fresnel surface. With respect to Grünwald’s formula, the plane portion of the slowness surface contributes to the fundamental matrix *E* the so-called *static term*, which yields the simplest of the three singular terms in *E*. We calculate it in §6.

The remaining singular parts of *E* originate from the discontinuities arising in the Abelian integrals implicit in Grünwald’s formula when (*t*,*x*) approaches the wave surface. We discuss these Abelian integrals in §7, and we show, in §8, how the Gaussian curvature *κ* on the slowness surface is connected with the jumps of the Abelian integrals under consideration. This parallels work by Buchwald, Duff and Burridge on the sharp waves in anisotropic elasticity (see Buchwald 1959; Duff 1960; Burridge 1967). Finally, in §9, we re-derive Darboux’s formula for the curvature *κ* on Fresnel’s surface, and in §10, we state the final result in which we put together all the information on the singular terms obtained in the paper.

Let us introduce some notation. As usual, Euclidean space is written as **R**^{n} or, if necessary, in order to indicate the variable to which we refer; **S**^{n−1} denotes the unit sphere in **R**^{n}, i.e. **S**^{n−1}:={*x*∈**R**^{n}:|*x*|=1}, with the surface measure d*σ*. In physical space-time **R**^{4}, we use the conjugate pair of variables (*τ*,*ξ*)=(*τ*,*ξ*_{1},*ξ*_{2},*ξ*_{3}) and (*t*,*x*)=(*t*,*x*_{1},*x*_{2},*x*_{3}). In matrix products, which are indicated by a dot, *x* is understood as a column vector, and the raised letter *T* means matrix transposition, such that *x*^{T}⋅*x* is the square of the modulus |*x*| of *x*, while *x*⋅*x*^{T} is a 3×3 matrix. *I*_{l} denotes the *l*×*l* unit matrix, and *A*^{ad} is the adjoint matrix of the *l*×*l* matrix *A*, i.e. The Heaviside function is denoted by *Y* .

We consider only differential operators with constant coefficients and use as differentiation symbols

We employ the standard notation for the distribution spaces the duals of the spaces of ‘test functions’ and of ‘rapidly decreasing functions’, respectively (see Schwartz 1966; Hörmander 1983*a*). By supp *T*, sing supp *T*, sing supp_{A} *T*, we denote the support, the singular support and the analytical singular support of respectively, i.e. the complements of the sets where *T* vanishes, is infinitely differentiable, or is real-analytical, respectively (see Hörmander (1983*a*, §§2.2 and 8.4) and Friedlander & Joshi (1998, definitions 1.4.1 and 8.6.1)). The composition of a distribution in one variable with a submersive function—as, e.g. in —is defined as in Friedlander & Joshi (1998, §7.2, p. 81) or Hörmander (1983*a*, theorem 6.1.2, p. 134).

We use the Fourier transform in the form this being extended to by continuity. (Here and also elsewhere, the Euclidean inner product (*y*,*η*)↦*yη* is simply expressed by juxtaposition.) When the Fourier transform is ‘partial’, say with respect to the space variables, we use notation as .

## 2. The system of crystal optics

Maxwell’s equations read as follows in the Gauß unit system:
and
where and denote the electric and magnetic field, and the dielectric displacement and the magnetic induction, respectively; and *ϱ* denote the densities of current and of charge, respectively, and *c* is the speed of light in vacuum (cf. Marion 1965, (4.18a–d), p. 111; Schwinger *et al.* 1998, (4.60), p. 42).

For anisotropic materials, these equations are supplemented by
where denotes the dielectric tensor (which can be assumed to be diagonal after a suitable rotation of the coordinate axes) and *μ* the magnetic permeability (which is assumed isotropic). We consider only homogeneous media, i.e. media for which *c*,*μ*,*ϵ*_{j} are constants.

By inserting equation (III) into equations (I) and (II) we obtain two uncoupled 3×3 systems of equations for the field vectors and
and
(for (IV) cf. Courant & Hilbert 1962, ch. VI, §3a, equation (8), p. 603; for (V) cf. Gårding 1989, (10), p. 363, 377). System (V) has the form
2.1
where *B*(*ξ*)*x*=−*ξ*×*A*(*ξ*×*x*), *A*=(*c*^{2}/*μ*) *ϵ*^{−1}, *x*∈**R**^{3}. Explicitly, *B*(*ξ*) is given by the symmetric matrix
2.2
where we have set *a*_{j}=*c*^{2}/(*μϵ*_{j}), *j*=1,2,3.

Similarly, system (IV) has the form
with *B*_{0}=*ϵ* and *B*_{2}(*ξ*)=(*c*^{2}/*μ*)(|*ξ*|^{2}*I*_{3}−*ξ*⋅*ξ*^{T}). In the following, we shall consider only the *biaxial case*, which means that the dielectric constants *ϵ*_{j} are pairwise different. By suitably numbering the coordinate axes, we can therefore assume that
2.3

Since *B*_{0},*B*_{2} are symmetric, Grünwald’s formula (see §5 below) applies to system (IV) for as well (cf. Ortner & Wagner 2004, p. 325), but for simplicity, we shall restrict our analysis in the following to the system *P*(*D*) for .

## 3. A brief review of the Atiyah–Bott–Gårding theory

Let *P*(*D*) denote here a general *l*×*l* matrix of linear constant coefficient differential operators on **R**^{n}, *D*=−i∂, and let us fix a direction *N*∈**R**^{n}\{0}. The system *P*(*D*) is called *hyperbolic* in the direction *N* if is a scalar hyperbolic operator with respect to *N*, i.e.
3.1
and
3.2
(cf. Atiyah *et al.* (1970, p. 128), Hörmander (1983*b*, Def. 12.3.3, p. 112) and Gårding (1984, p. 221)).

In particular, if *Q*(*D*) is homogeneous and satisfies condition (3.1), then condition (3.2) is equivalent to the condition
3.3
see Hörmander (1983*b*, theorem 12.4.3, p. 113).

Let us consider first a scalar operator *R*(*D*) which is hyperbolic in the direction *N*. The *hyperbolicity cone* *Γ*(*R*(*D*)) is the connectivity component containing *N* in the set of all *η*∈**R**^{n}\{0} with respect to which *R*(*D*) is hyperbolic. Actually, *Γ*(*R*(*D*)) is an open convex cone and it coincides with the connectivity component containing *N* of {*η*;*R*^{pr}(*η*)≠0} (see Hörmander (1983*b*), corollary 12.4.5, p. 114). Furthermore, *R*(*D*) has one and only one fundamental solution *E*(*R*(*D*)) with support in the half-space *H*_{N}={*y*∈**R**^{n}; *Ny*≥0} (see Hörmander (1983*b*), theorem 12.5.1, p. 120). The support of *E*(*R*(*D*)) is actually contained in the *propagation cone*

For *η*∈**R**^{n}, *R*_{η}(*ζ*) denotes the *localization* of *R* at infinity in the direction *η*∈**R**^{n}, i.e. the lowest (non-vanishing) coefficient with respect to *s* in the MacLaurin series of *s*^{m}*R*(*ζ*+*η*/*s*), . The operators *R*_{η}(*D*) are again hyperbolic. The set is called the *wavefront surface* of *R*(*D*). One of the central results of Atiyah, Bott and Gårding is the following series of inclusions:
3.4
see Atiyah *et al.* (1970, theorem 4.10, p. 144; theorem 7.24, p. 177). In the physically relevant cases, in particular if *n*≤4, all the inclusions in (3.4) are identities, cf. Atiyah *et al.* (1973, theorem 7.7, p. 175).

If *Ξ*={*η*∈**R**^{n}; *R*^{pr}(*η*)=0} is the *slowness surface* of *R*(*D*) and *η*∈*Ξ* fulfils (*R*^{pr})′(*η*)≠0 (and *η* is thus a regular point on *Ξ*), then *R*_{η} is a first-order polynomial and *K*(*R*_{η}(*D*)) is the half-ray **R**⋅(*R*^{pr})′(*η*)∩*H*. Hence, *Ξ**∩*H*⊂sing supp *E*(*R*(*D*)) if *Ξ** denotes the algebraic hypersurface dual to *Ξ*, the so-called *wave surface*, i.e. *Ξ** is the closure of {*t*(*R*^{pr})′(*η*); *t*∈**R**, *η*∈*Ξ*}. (Here, we suppose that *R*^{pr} does not contain multiple factors.) We shall say that *conical refraction* occurs if *Ξ**∩*H* is a proper subset of sing supp *E*(*R*(*D*)). If the inclusions in (3.4) are identities (as is usually the case), then conical refraction occurs if the set
3.5
is not already contained in *Ξ**.

Let us mention that the name ‘conical refraction’ actually refers to the occurrence of the corresponding physical phenomenon in the reflection of light in biaxial crystals. This phenomenon was predicted on mathematical grounds in 1832 by Hamilton (1931), cf. also our determination of the wavefront surface in the next section, and verified experimentally in the same year by H. Lloyd (see Gårding (1989), pp. 360–362).

Let us now consider a hyperbolic *l*×*l* system *P*(*D*) as in the beginning of this section. As in the scalar case, *P*(*D*) has one and only one *fundamental matrix* , i.e. *P*(*D*)*E*=*I*_{l}*δ*, with support in the half-space *H*_{N} (see Atiyah *et al.* (1970), lemma 3.1, p. 127). In fact, we can express *E*(*P*(*D*)) by means of the adjoint system *P*^{ad}(*D*) and the fundamental solution *E*(*Q*(*D*)) of the determinant operator
3.6
cf. Hörmander (1969, Section 3.8, p. 94) and Schwartz (1966, p. 140).

An extension of the method of localization to the matrix case has been sketched in Esser (1987) and recently been elaborated in Ortner & Wagner (2009*b*). We shall not go into the details here. Let us instead just state the following inclusion, which trivially follows from equation (3.6), and which yields an upper estimate for the singular support of *E*(*P*(*D*)):
3.7

## 4. The Fresnel surface

Let us specialize the above to the system *P*(*D*) in equation (2.1). Hence, *η*=(*τ*,*ξ*)∈**R**^{4} and *P*(*τ*,*ξ*)=*B*(*ξ*)−*I*_{3}*τ*^{2}, where *B*(*ξ*)=−*ξ*×*A*(*ξ*×*x*). Owing to *B*(*ξ*)*ξ*=0, we have and from the relations
(where *a*_{j+3}:=*a*_{j}, *j*=1,2,3), we conclude
with
4.1
cf. Courant & Hilbert (1962, (13), p. 606), Gårding (1989, (119), p. 363) and Liess (1993, (6.7.6), p. 272).

Due to *a*_{j}>0, *j*=1,2,3, Jacobi’s criterion implies that *B*(*ξ*) in equation (2.2) is positive semi-definite. Hence, its eigenvalues *λ*_{1}(*ξ*)=0,*λ*_{2}(*ξ*),*λ*_{3}(*ξ*) are non-negative and the polynomials have the real roots 0 (twice), Thus the conditions (3.1) and (3.3) are satisfied with *N*=(1,0,0,0), which means that *P*(*D*),*Q*(*D*) and *R*(*D*) are hyperbolic with respect to the time variable.

Let us investigate the slowness surface *Ξ*={(*τ*,*ξ*)∈**R**^{4}; *R*(*τ*,*ξ*)=0} of the hyperbolic operator *R*(*D*). The projective surface corresponding to the cone *Ξ* as well as its affine representation *X*={*ξ*∈**R**^{3}; *R*(1,*ξ*)=0} are called *Fresnel’s surface*. Note that *X* is the union of the two sheets where either one of the two non-zero eigenvalues *λ*_{2,3}(*ξ*) of the matrix *B*(*ξ*) equals 1. Therefore, *X* is non-singular apart from the points where *λ*_{2} and *λ*_{3} coincide, i.e. where *B*(*ξ*)−*I*_{3} has rank one. For such *ξ*, all vectors orthogonal to *ξ* must be eigenvectors of *B*(*ξ*), and inserting (*ξ*_{2},−*ξ*_{1},0) yields *ξ*_{1}*ξ*_{2}*ξ*_{3}(*a*_{1}−*a*_{2})=0. Owing to the inequalities in (2.3), we conclude that the singular points of *X* must lie on one of the three coordinate planes. The rank one condition then readily yields that *X* has exactly four singular points, which—due to our ordering of the constants *a*_{j} in (2.3)—are given by
4.2
cf. Liess (1993, p. 273).

Hence *X* is homeomorphic to two disjoint spheres glued together at the four singular points of *X*, which two by two are pairwise opposite and span the ‘optical axes’, see figure 1, which is taken from Fladt & Baur (1975, §6, p. 346) by courtesy of Springer. In this figure, the outer sheet of *X* and two of the four singular points are plainly visible. (Note that the *x*-axis in figure 1 corresponds to the *ξ*_{2}-axis in our notation.)

Let us next calculate the wavefront surface *W*(*R*(*D*)), which—according to §3—coincides with the singular support of the fundamental solution *E*(*R*(*D*)).

We first determine the wave surface *Ξ** in a similar way as in Poincaré (1889), Esser (1987, (6.a), p. 203) and Gårding (1989, pp. 359, 360). If (*τ*,*ξ*)∈**R**^{4} such that *τ*^{2}≠*a*_{j}|*ξ*|^{2}, *j*=1,2,3, then an easy calculation shows that *R*(*τ*,*ξ*)=0 is equivalent to
4.3
Therefore, *Ξ** is the closure of the cone {(*t*,*x*)=(∂*f*/∂*τ*,∇*f*)(*τ*,*ξ*); *f*(*τ*,*ξ*)=0}.

We set *η*_{j}=*τ*^{2}−*a*_{j}|*ξ*|^{2} and *y*_{j}=*a*_{j}*t*^{2}−|*x*|^{2}. Then
Because of
*f*(*τ*,*ξ*)=0 implies

Using *f*(*τ*,*ξ*)=0 again we conclude from this
and
Hence
4.4
and
which is equivalent to
Therefore, *Ξ** is the closure of the cone

The same conversion which led from the polynomial *R* in equation (4.1) to equation (4.3) then yields that *Ξ** is the zero set of the polynomial
4.5
Hence, the wave surface *Ξ** is again a Fresnel surface and has the same appearance as *Ξ* depicted in figure 1.

As explained in §3, the *wavefront surface* *W*(*R*(*D*))=sing supp *E*(*R*(*D*)) is obtained from the *wave surface* *Ξ** by adjoining to *Ξ**∩*H* the union *C* of the propagation cones *K*(*R*_{η}(*D*)), where *η* runs through the singular points on *Ξ*, cf. equation (3.5).

If we set *η*=(1,*η*_{1},0,*η*_{3}) with
and
in accordance with equation (4.2), then
and therefore the localization of *R* at infinity in the direction *η* is given by Upon differentiating *R* in equation (4.1) and setting *ζ*=(*τ*,*ξ*), *ξ*∈**R**^{3}, a straight-forward calculation yields
cf. Liess (1991, p. 75).

The symbol of the operator *R*_{η}(*D*) becomes more transparent when expressed in orthonormal coordinates *y*_{1},*y*_{2},*y*_{3}, such that *y*_{1} runs in the direction of (*a*_{3}−*a*_{2})*η*_{1}*x*_{1}+(*a*_{1}−*a*_{2})*η*_{3}*x*_{3}, *y*_{2}=*x*_{2} and If, furthermore,
then
and hence
4.6
is a shifted *circular* wave operator in two space dimensions.

Starting from the fundamental solution of the planar wave operator in **R**^{4}, i.e.
and using a linear transformation, cf. Wagner (1984, Satz 4, p. 10), we obtain from equation (4.6)
Therefore, the propagation cone of *R*_{η}(*D*) is given by
4.7

Hence, the affine representation of the wavefront surface of *R*(*D*) consists of the Fresnel surface *X** augmented by four circular discs which are bounded by the so-called *Hamiltonian circles*, cf. Ludwig (1961, p. 117) and Liess (1991, p. 68). Owing to equation (4.7), these circles lie in the planes where i.e. *η*_{1}*x*_{1}+*η*_{3}*x*_{3}+1=0, they have the centres *y*_{1}=−*c*/2, *y*_{2}=0, , i.e. *x*=−(1/2)(*η*_{1}(*a*_{2}+*a*_{3}),0,*η*_{3}(*a*_{1}+*a*_{2})), and the diameter *c*, cf. Esser (1987, p. 206).

Since the circular cone *R*_{η}(*ζ*)=0 is the tangent variety of *R* at the singular point *η*, the dual cone, which coincides with the boundary of the propagation cone *K*(*R*_{η}(*D*)), is contained in *Ξ**. Hence the Hamiltonian circles, two of which are highlighted in figure 1, belong to *X**. Furthermore, because *X*^{**}=*X*, the planes tangential to *X** at points of the Hamiltonian circles correspond to the four singular points of *X* given in equation (4.2), i.e. they coincide with the planes *η*_{1}*x*_{1}+*η*_{3}*x*_{3}+1=0 containing the circular discs {*x*∈**R**^{3}; (1,*x*)∈*K*(*R*_{η}(*D*))}. Therefore, these four circular discs lie on the convex completion of the outer sheet of *X** (cf. also Esser (1987), pp. 206–209).

Let us return now to the system *P*(*D*) of crystal optics given in equations (2.1) and (2.2). If, as above, *Ξ* is the slowness surface of the operator *R*(*D*) in equation (4.1), then the determinant has the slowness surface with the dual This implies *W*(*Q*(*D*))=*W*(*R*(*D*))∪{(*t*,0); *t*≥0} for the wavefront surfaces, and, due to equation (3.7),
4.8
Actually, the inclusion in (4.8) is an equality, as one can see by applying the localization theory for systems in Ortner & Wagner (2009*b*). Summing it up, we obtain that the intersection of the singular support of *E*(*P*(*D*)) with the plane *t*=1 consists of the Fresnel surface *X** given by *R**(1,*x*)=0, *R** as in equation (4.5), of the four circular discs bounded by the Hamiltonian circles, and of the origin *x*=0.

## 5. Grünwald’s formula

For the system *P*(*D*) of crystal optics, a representation of *E*(*P*(*D*)) by integrals over the Fresnel surface was first obtained in Grünwald (1904), see the historic account in Gårding (1989, pp. 377, 378). For general homogeneous systems, such representations later came to be known as the *Herglotz–Petrovsky formulas* (see Herglotz (1926,1928), Petrovsky (1945), Leray (1952), Gel’fand & Shilov (1964, ch. I, 6.1), Atiyah *et al.* (1970, pp. 176, 177), Hörmander (1983*b*, (12.6.10), p. 129), Gårding (1989, theorem 2, p. 375), Ortner & Wagner (2004, §2) and Ortner & Wagner (2008, proposition 1, p. 419)). For completeness, we shall shortly derive Grünwald’s formula.

Let us suppose here, more generally, that is an *l*×*l* system of homogeneous second-order differential operators in **R**^{4} such that *B* is a symmetric *l*×*l* matrix. We assume, furthermore, that *P*(*D*) is hyperbolic, which—according to (3.3)—means that *B*(*ξ*) is positive semi-definite for each *ξ*∈**R**^{3}.

The fundamental matrix *E*=*E*(*P*(*D*)) with support in *t*≥0 is temperate, cf. Ortner & Wagner (1992, proposition 1, p. 530), and partial Fourier transform with respect to the space variables yields the equation This system has the uniquely determined solution and hence we arrive at the representation
5.1
Note that is a convergent power series into which the matrix *B*(*ξ*) can be inserted; the result is a temperate distribution as *B*(*ξ*) is positive semi-definite, and it depends continuously on *t*; therefore, the multiplication with the Heaviside function in equation (5.1) is legitimate.

More explicitly, let 0≤*λ*_{1}(*ξ*)≤⋯≤*λ*_{l}(*ξ*) be the eigenvalues of *B*(*ξ*) (counted with multiplicity), and *e*_{1}(*ξ*),…,*e*_{l}(*ξ*) an orthonormal basis of corresponding eigenvectors. Then
and
5.2

We observe that is even and homogeneous of degree one owing to the homogeneity of *B*(*ξ*). Let us assume first that the eigenvalues *λ*_{j}(*ξ*) are positive and pairwise different for *ξ*≠0. (The second assumption just means that the determinant operator is *strictly hyperbolic*.) Then, *ρ*=*ρ*_{j}(*ξ*) is on **R**^{3}\{0} and hence is a submanifold of **R**^{3}. We can express *ξ*∈**R**^{3}\{0} by *polar coordinates* with respect to *ρ*_{j}, i.e.
are diffeomorphisms for *j*=1,…,*l*.

The maps *F*_{j} become measure-theoretic isomorphisms when *Γ*_{j} is equipped with the positive measure d*μ*, which corresponds to the Leray form *ξ*_{1} d*ξ*_{2}∧d*ξ*_{3}+*ξ*_{2} d*ξ*_{3}∧d*ξ*_{1}+*ξ*_{3} d*ξ*_{1}∧d*ξ*_{2}, and has the measure *ρ*^{2} d*ρ*, i.e.
cf. Wagner (1989, p. 409). Hence, the functions occurring in equation (5.2) can be written as integrals over *Γ*_{j}:
where the distributions *T*_{ω}∈𝒮′(**R**^{3}) are given by
(Here we have used that *ρ*_{j} and *e*_{j}(*ω*)⋅*e*_{j}(*ω*)^{T} are even; also note that the distribution-valued function *ω*↦*e*_{j}(*ω*)⋅*e*_{j}(*ω*)^{T}⋅*T*_{ω} can be integrated over *Γ*_{j} as it is continuous.)

*T*_{ω} is essentially a distribution in one variable, and we can easily determine its inverse Fourier transform. If *ω*=(*ω*_{1},0,0), *ω*_{1}>0, say, then
and thus
Therefore, generally,

The slowness surface *X*={*ω*∈**R**^{3}; *Q*(1,*ω*)=0} of the determinant operator is the disjoint union of the manifolds *Γ*_{j}, *j*=1,…,*l*, since
For *ω*∈*X*, let *e*(*ω*) be an eigenvector of *B*(*ω*) for the eigenvalue 1 and of length 1. Hence, *e*(*ω*)=±*e*_{j}(*ω*) if *ω*∈*Γ*_{j}. If we finally take into account that *X* as well as *e*(*ω*)⋅*e*(*ω*)^{T}, *ω*∈*X* are even, we obtain Grünwald’s formula:
5.3
cf. Grünwald (1904, p. 520), Gårding (1989, (47), p. 378) and Ortner & Wagner (2004, (G), p. 326).

The precise meaning of formula (5.3) is the following: The function
is integrated over *X* and thus yields a function This resulting function is then differentiated and multiplied by *Y* (*t*). Incidentally, we note that *E* is still continuous at 0 as a function of *t*, i.e. This is clear from equation (5.2), but it also follows from equation (5.3), as the substitution *ω*↦−*ω* shows that vanishes.

## 6. The static term and geometric considerations

Formula (5.3) is not directly applicable to the system of crystal optics for two reasons: First, *λ*_{1}(*ξ*), the smallest of the three eigenvalues of the matrix *B*(*ξ*) in equation (2.2) vanishes, and, second, *λ*_{2}(*ξ*) and *λ*_{3}(*ξ*) are not distinct in the four singular points of *X*.

Let us start by addressing the first problem. Note that formula (5.2) is still valid with the proviso that the value of the holomorphic function at the origin is *t*. Owing to *e*_{1}(*ξ*)=*ξ*/|*ξ*| for *ξ*≠0, the summand in equation (5.2) pertaining to *j*=1, the so-called *static term*, is given by
6.1
Indeed, the three-dimensional Laplace operator Δ_{3} possesses one and only one homogeneous fundamental solution, namely −1/(4*π*|*x*|), which therefore has to coincide with .

Away from 0, ∂_{j}(|*x*|^{−1})=−*x*_{j}|*x*|^{−3}, and this is a locally integrable function. As both sides of the last equation are homogeneous distributions of degree −2, and as *δ* and its derivatives have smaller homogeneities in **R**^{3}, we conclude that the equation ∂_{j}(|*x*|^{−1})=−*x*_{j}|*x*|^{−3} holds in In contrast, the function ∂_{k}∂_{j}(|*x*|^{−1})=(3*x*_{j}*x*_{k}−|*x*|^{2}*δ*_{jk})|*x*|^{−5} is no longer locally integrable. We can associate a distribution to the function (3*x*_{j}*x*_{k}−|*x*|^{2}*δ*_{jk})|*x*|^{−5} by defining the principal value
Again, it follows that
holds for *x*≠0, and hence the difference of these two distributions in which is homogeneous of degree −3, is a multiple of *δ*.

Generally, if *h*(*x*) is in **R**^{n}\{0} and homogeneous of degree 1−*n*, then
as one readily sees by application to a radially symmetric test function, cf. Petersen (1983, theorem 15.8, p. 42). In our case, *h*(*x*)=−*x*_{j}|*x*|^{−3} and
Hence, we conclude from equation (6.1) that
6.2

We mention that formula (6.2) for the static term *E*_{1} was calculated by a different method in Ortner & Wagner (2004, pp. 334, 335), cf. also Wagner (2004, p. 415). The static term for the electric field (governed by system (IV) in §2) was expressed in a very complicated manner in Burridge & Qian (2006, (5.11), p. 78, and (B19), p. 93); a simpler formula for it comparable with equation (6.2) was derived in the study of Ortner & Wagner (2009*a*).

Let us now consider the remaining two, positive eigenvalues *λ*_{2}(*ξ*)≤*λ*_{3}(*ξ*) of *B*(*ξ*) and tackle the second of the necessary adaptations to formula (5.3) articulated at the beginning of this section. The sets *Γ*_{j}={*ω*∈**R**^{3}; *λ*_{j}(*ω*)=1}, *j*=2,3 are submanifolds only outside their intersection *Γ*_{2}∩*Γ*_{3}, which consists of the four singular points given in equation (4.2). The union *Γ*_{2}∪*Γ*_{3} is the Fresnel surface *X*={*ξ*∈**R**^{3}; *R*(1,*ξ*)=0} investigated in §4 and depicted in figure 1. The sets *Γ*_{2} and *Γ*_{3} constitute the outer and the inner sheet of *X*, respectively.

The deduction of formula (5.3) for the two parts in the sum in equation (5.2) pertaining to the indices *j*=2,3 remains practically unchanged with the sole difference that *e*(*ω*)⋅*e*(*ω*)^{T} is discontinuous in the four singular points of *X*, and hence the integral of *e*(*ω*)⋅*e*(*ω*)^{T}⋅*δ*(*t*+*ωx*) over *X* is now one over a distribution-valued function, which is merely bounded in these four points. Thus, we conclude that the fundamental matrix, *E*=*E*(*P*(*D*)) of crystal optics is given by
6.3
with *X*={*ω*∈**R**^{3}; *R*(1,*ω*)=0}, *R* as in equation (4.1).

Because the fundamental matrix *E*(*t*,*x*) is homogeneous of degree −2 and vanishes for negative *t*, it will be enough to consider its values for *t*=1. The complement *U*={*x*∈**R**^{3}; (1,*x*)∉*W*(*R*(*D*))} of the wavefront surface of *R*(*D*) consists of seven connectivity components (cf. figure 2), and we shall successively discuss the form of the integral representation (6.3) in these seven regions. Concerning figure 2, let us remark that the intersection of *X** with the plane *x*_{2}=0 containing the optical axes is given by the equation
and hence this intersection is made up of a circle and an ellipse meeting in the four singular points of *X**.

The component

*J*of*U*containing 0 is bounded by the innermost sheet of*X**. Since*X** is the dual surface of*X*, we obtain that*δ*(1+*ωx*) vanishes for*ω*∈*X*and*x*∈*J*, and hence*E*agrees with the static term*E*_{1}inside*J*, cf. Zeilon (1921, 18., p. 108). Incidentally, we notice that*E*is independent of the constants*a*_{1},*a*_{2},*a*_{3}in the region*J*.According to §3,

*E*vanishes outside the propagation cone*K*(*Q*(*D*))=*K*(*R*(*D*)), which is the convex completion of the wavefront surface*W*(*R*(*D*)). Hence,*E*_{1}coincides with for (*t*,*x*) outside*K*(*R*(*D*)).For

*x*in the connected region between the two sheets of*X**, which region is split into the four areas lying between the circle and the ellipse in the section of figure 2, the distribution*δ*(1+*ωx*) vanishes when*ω*belongs to the inner sheet*Γ*_{3}of*X*and hence 6.4Finally, if

*x*belongs to one of the four regions bounded by the lids of conical refraction and the parts of the outer sheets of*X**, where the Gaussian curvature is negative, then formula (6.4) holds again.

The cases (c) and (d) differ in the following geometric circumstance: For *x* as in case (c), the plane *ϵ*:1+*ωx*=0 intersects *Γ*_{2} in a single curve, whereas *ϵ* intersects *Γ*_{2} in two separate curves if *x* lies in one of the four domains of case (d).

## 7. Abelian integrals

In formula (6.3), the fundamental matrix *E* of system (2.1) is represented by means of a distribution-valued integral over *ω*∈*X*. Let us consider this integral, i.e.
now in more detail. By Fubini’s theorem, the evaluation of *F* on a test function furnishes the following:
where the curves *C*_{t,x}:={*ω*∈*X*; *t*+*ωx*=0} are provided with the positive measure d*ν* induced by d*μ* (see below). This shows that *F*(*t*,*x*) is a locally integrable function given by
7.1

Let us first comment on formula (7.1) from the viewpoint of algebraic geometry. The real curve *C*_{t,x} corresponds to a cycle in the complex projective curve
which, generically, is a compact Riemannian surface of genus 3 (see Griffiths & Harris 1978, p. 220). The measure d*ν* is induced by a holomorphic one-form on *S* (see Griffiths & Harris 1978, p. 221), and equation (7.1) thus represents *F* as an Abelian integral, cf. Zeilon (1921, p. 124).

Let us next express the measure d*ν* on *C*_{t,x} by means of the arc-length d*s*. From *R*(1,*ξ*)=0 on *X*, we infer that the surface measure d*σ* on *X* is given by d*σ*=|∇*R*|/|∂_{3}*R*| d*ξ*_{1} d*ξ*_{2}; on the other hand, owing to and *ξ*^{T}⋅∇*R*(1,*ξ*)=−∂_{τ}*R*(1,*ξ*) on *X* by Euler, the Leray measure d*μ* is given by d*μ*=|∂_{τ}*R*/∂_{3}*R*| d*ξ*_{1} d*ξ*_{2}. Hence, d*μ*=|∂_{τ}*R*|/|∇*R*| d*σ*. By construction, d*ν*/d*s* equals d*μ*/d*σ* divided by the length of the part of *x* tangential to *X*, i.e. by |*x*×∇*R*/|∇*R*||. Therefore, we obtain d*ν*=|∂_{τ}*R*|/|*x*×∇*R*| d*s* and
7.2
(Let us point out that the measure d*ν* does not merely depend on the curve *C*_{t,x}, but also on |*x*|. In fact, the function *F*(*t*,*x*) is homogeneous of degree −1.) We mention that the integrals analogous to equation (7.2), but corresponding to the electric field have been numerically evaluated in Burridge & Qian (2006, pp. 76, 77).

Formula (7.2) also entails that *F* is finite and as long as the plane *t*+*ξx*=0 does intersect *X* neither tangentially nor in the four singular points of *X* given in equation (4.2). In fact, in this case, the vectors ∇*R*(1,*ω*) and *x* are not parallel at any of the points *ω* of *C*_{t,x} and the denominator in formula (7.2) does not vanish. This condition precisely means that (*t*,*x*) does not belong to the wavefront set *W*(*R*(*D*)), and hence we obtain a second proof of the inclusion relation sing supp *E*⊂*W*(*R*(*D*)); a fact we investigated in the framework of the Atiyah–Bott–Gårding theory in §§3 and 4. We also observe that for fixed *t*. In fact, equation (5.2) implies
and thus as the functions *λ*_{j} being homogeneous of degree −2.

Let us yet evaluate the rank-one matrix *e*(*ξ*)⋅*e*(*ξ*)^{T}, *ξ*∈*X*, in terms of the dual projective variable (*t*,*x*)=∂*R*(1,*ξ*), which yields the tangent plane {*η*∈**R**^{3}; *t*+*ηx*=0} of *X* in *ξ*. First recall that *e*(*ξ*) is an eigenvector of *B*(*ξ*) to the eigenvalue 1. An easy calculation shows that the vector *u* satisfying
belongs to the kernel of the matrix *B*(*ξ*)−*I*_{3} if *R*(1,*ξ*)=0. If we set *η*_{j}=1−*a*_{j}|*ξ*|^{2} and *y*_{j}=*a*_{j}*t*^{2}−|*x*|^{2} as in §4 and employ equation (4.4), we conclude that
and, finally,
7.3

## 8. The sharp waves

As we have seen in §7, the singularities of *E*−*E*_{1} arise from the discontinuities of the integral in equation (7.2) at points (*t*,*x*) on the wavefront surface *W*(*R*(*D*)). The form of these discontinuities depends on how the tangent plane {*η*∈**R**^{3}; *t*+*ηx*=0} intersects *X* in its point of tangency *ξ*∈*X*, i.e. whether the Gaussian curvature of *X* in *ξ* is positive, negative or zero, respectively. For elastic waves in anisotropic solids, a comparable analysis of the discontinuities of *F* in equation (7.2) was carried out in Buchwald (1959), Duff (1960) and Burridge (1967). Owing to the differentiation with respect to *t* in equation (6.3), these discontinuities then yield the singular part of *E*−*E*_{1}, the so-called ‘sharp waves’, cf. Duff (1960).

In a summary fashion, points on *X* with positive and negative Gaussian curvature, respectively, produce delta waves and Cauchy principal values along the wave surface *Ξ**, respectively, whereas the singular points on *X* do not contribute to the singular part of the fundamental matrix *E*, but yield only discontinuities in *E* along the lids of conical refraction.

Let us start by investigating the singularities produced by a point *ξ*∈*X* at which *X* has positive Gaussian curvature. If (*t*,*x*)=±∂*R*(1,*ξ*), *t*>0, and *ξ* belongs to the inner sheet *Γ*_{3} of the wave surface, then *C*_{s,x} consists, for *s*<*t* of two curves one of which, say shrinks to the point *ξ* when *s* approaches *t* from below. Similarly, if *ξ*∈*Γ*_{2}, then is just one curve shrinking to {*ξ*} for *s*↗*t*. The matrix *e*(*ω*)⋅*e*(*ω*)^{T} converges on to *e*(*ξ*)⋅*e*(*ξ*)^{T}, which has been expressed as a function of (*t*,*x*) in equation (7.3). Hence, it remains to calculate the limit of for *s*↗*t*. According to the formula for d*ν* in §7, we obtain

In linear coordinates *y* centered at *ξ*, the polynomial *R*(1,*ω*) can be expressed in the form where
and the positive numbers *κ*_{1},*κ*_{2} are the principal curvatures of *X* in *ξ*. Furthermore, *x*/|*x*| corresponds to (0,0,1)^{T} and coincides in the first approximation with the ellipse for some positive *ϵ*. The parametrization
yields
and hence
8.1
*κ*=*κ*_{1}*κ*_{2} being the Gaussian curvature of *X* at *ξ*. As *F* is homogeneous of degree −1, formula (8.1) generally holds for (*t*,*x*)∈*Ξ**∩*H* if {*η*∈**R**^{3}; *t*+*ηx*=0} is the tangent plane of *X* in *ξ* and *κ*>0.

Similarly, if the Gaussian curvature *κ* of *X* in *ξ* is negative, we have *κ*=*κ*_{1}*κ*_{2} with *κ*_{1}>0, *κ*_{2}<0, and the curves
are hyperbolas. For *ϵ*>0, say, they are parametrized by
and this implies
and hence .

We have to integrate over the intersection of *C*_{ϵ} with a fixed neighbourhood of the origin, e.g. ⃛he one that is determined by the inequality
The range of *φ* is therefore asymptotically of the length for *ϵ*→0, and
8.2
is convergent.

Finally, let (*t*,*x*) belong to *W*(*R*(*D*))\*Ξ**, which means that *t*>0 and *x*/*t* belongs to one of the four circular discs of conical refraction, or, equivalently, that (*t*,*x*) lies in the interior of one of the cones *K*(*R*_{η}(*D*)) given in equation (4.7). Then *C*_{t,x} consists of one of the singular points *ξ* in equation (4.2) and of a curve. For *s* near *t*, *C*_{s,x} consists of two curves, a small loop shrinking to *ξ* for *s*→*t*, and a large one that converges to *C*_{t,x}\{*ξ*} for *s*→*t*. As *x*/*t* is in the *interior* of one of the lids of conical refraction, the angle between ∇*R*(1,*ω*) and *x* cannot converge to 0, when *ω*∈*C*_{s,x} and *s*→*t*. Therefore, the contribution to *F*(*s*,*x*) in equation (7.2) from the small loops in *C*_{s,x} converges to 0 for *s*→*t*, i.e. *s*↦*F*(*s*,*x*) is continuous at *t*.

Let us now investigate the impact of the discontinuities of the integral *F* defined in equation (7.1) on the fundamental matrix *E*. From formula (6.3), we infer that the jumps owing to the parts of *X* with positive Gaussian curvature produce delta terms in *E*. More precisely, if *ξ*∈*X* is dual to (*t*_{0},*x*_{0})∈*Ξ**∩*H*, i.e. the tangent plane to *X* in *ξ* is {*η*∈**R**^{3}; *t*_{0}+*x*_{0}*η*=0}, and if *X* has positive curvature *κ* in *ξ*, then *s*↦*F*(*s*,*x*_{0}) is discontinuous at *t*_{0}, and fulfils, by equation (8.1),

If we resolve the equation *R**(*t*,*x*)=0 defining *Ξ** in a neighbourhood of (*t*_{0},*x*_{0}) in the form *t*=*t*(*x*), then the *jump formula*, see, e.g. Schwartz (1966, ch. II, §3, example 1) or Wagner (2010, (3.1), p. 1188), implies that
is locally integrable near (*t*_{0},*x*_{0}). Here, *e*(*ξ*)⋅*e*(*ξ*)^{T} is considered as a function of (*t*,*x*) by means of equation (7.3), and the Gaussian curvature *κ*=*κ*_{X}(*ξ*) will be expressed by (*t*,*x*) in the next section. Since *δ*(*t*−*t*(*x*))=*δ*(*R**(*t*,*x*))⋅|∂_{t}*R**|, we thus obtain that
8.3
holds near the parts of *Ξ**, where *κ*_{X}(*ξ*)>0, cf. Burridge (1967, (6.8), p. 47).

We observe that *δ*(*R**(*t*,*x*)) is a Radon measure in **R**^{4}, which is well-defined also near the conical singularities of *Ξ**, as one sees by comparison with We shall verify in the next section that is continuous and vanishes at the singular points of *Ξ**. Hence, the product in equation (8.3) is meaningful.

If *ξ*∈*X* and *y*∈*X**={*y*∈**R**^{3}; (1,*y*)∈*Ξ**} are dual, i.e. {*η*∈**R**^{3}; 1+*ηy*=0} is the tangent plane to *X* in *ξ*, then *κ*_{X}(*ξ*)*κ*_{X*}(*y*)=|*ξ*|^{−4}|*y*|^{−4}, cf. §9 below, and hence *κ*_{X}(*ξ*)>0 iff *X** has positive curvature in *x*/*t*. We also point out that, as the prefactor of the delta function in equation (8.3) vanishes if *x*/*t* lies on one of the Hamiltonian circles.

Similarly, near a point (*t*_{0},*x*_{0})∈*Ξ** where the Gaussian curvature *κ*_{X}(*ξ*) is negative, we obtain a Cauchy principal value from the expression (8.2) by differentiation. Hence
8.4
holds near the curved-inward section of *Ξ**, cf. Burridge (1967, (6.9), p. 47). (Here *ξ*∈*X* is a function of (*t*,*x*) such that *ξ* is dual to (*t*,*x*) if (*t*,*x*)∈*Ξ**∩*H*.)

Let us yet give a rigorous justification for the fact that there do not appear any additional singularities in *E* at points (*t*,*x*) for which *t*>0 and *x*/*t* belongs to the union *M* consisting of the Hamiltonian circles *M*_{1} and the four singular points *M*_{2} of *X**. In fact, as for fixed *t* (cf. §7), we can write where 1−*χ*_{ϵ} are the characteristic functions of the sets {*x*∈**R**^{3}; *d*(*x*,*M*)<*ϵ*}. We then observe that *F*⋅∂_{t}*χ*_{ϵ}(*x*/*t*) converges to zero in for *ϵ*↘0. For example, for *M*_{2}, one obtains in suitable coordinates *y* a limit of the form where and Hence, the singular part of *E* is, apart from *E*_{1}, indeed represented by equations (8.3) and (8.4).

In order to obtain completely explicit results for the singular part of *E*, we will re-derive Darboux’s formula for the curvature *κ*_{X}(*ξ*).

## 9. Darboux’s curvature formula

Let us first explain more in detail the notion of *duality* for hypersurfaces, which notion was implicit already in the definition of *X* and *X** in §4. For complex algebraic curves, this subject is beautifully expounded in Griffiths & Harris (1978, p. 263 ff.).

If *V* is a finite-dimensional vector space over **R**, then the corresponding projective space is the set of all one-dimensional subspaces in *V*, i.e.
The projective space **P**(*V* *) is canonically identified with the set of all subspaces of *V* of codimension one and is called the *dual* projective space. If *Ξ*⊂**P**(*V*) is an algebraic hypersurface given as the zero-set of a homogeneous polynomial *R*, i.e. *Ξ*={[*v*]∈**P**(*V*); *R*(*v*)=0}, then the set of all hyperplanes tangential to *Ξ* is an algebraic variety in **P**(*V* *), called the *dual variety* *Ξ**. In particular, for *V* =**R**^{n+1}, we identify **P**(*V*) and **P**(*V* *), and also *Ξ* with its affine representation *X*⊂**R**^{n}, i.e. *Ξ* is the completion of {[1,*ξ*]; *ξ*∈*X*} in **P**(**R**^{n+1}). Then *Ξ** corresponds to *X**⊂**R**^{n} given by the set of *x*∈**R**^{n}, such that {*η*∈**R**^{n}; 1+*xη*=0} is tangential to *X*. We also obtain the *duality map*
where {*η*∈**R**^{n}; 1+*D*(*ξ*)*η*=0} is the hyperplane tangential to *X* at *ξ*. (Strictly speaking, *D* is defined only for non-singular points where the tangent hyperplane does not contain the origin. We shall also assume here that *X** is again a hypersurface.)

If *X* is near *ξ*_{0} and locally parametrized by *ξ*=*ξ*(*s*_{1},…,*s*_{n−1}), then *x*(*s*_{1},…,*s*_{n−1})=*D*(*ξ*(*s*_{1},…,*s*_{n−1})) yields the corresponding part of *X**, and *ξ*^{T}⋅*x*=−1 and ∂_{j}*ξ*^{T}⋅*x*=0 imply *ξ*^{T}⋅∂_{j}*x*=0, or, in other words, *X*^{**}=*X*, cf. Griffiths & Harris (1978, p. 278) for the geometric intuition behind this biduality relation. For *v*_{1},…,*v*_{n−1}∈**R**^{n}, we define the cross-product *w*=*v*_{1}×⋯×*v*_{n−1}∈**R**^{n}, such that for all *y*∈**R**^{n}. Then ∂_{1}*ξ*×⋯×∂_{n−1}*ξ* is parallel to *x*, and, dually, ∂_{1}*x*×⋯×∂_{n−1}*x* is parallel to *ξ*. Hence, we conclude from *ξ*^{T}⋅*x*=−1 that
9.1

If *X* is at *ξ* and *T*_{ξ}*X* denotes the tangent space to *X* at *ξ* in the sense of differential geometry, then the *Weingarten map* originates by differentiating the unit normal *ν*(*ξ*)=*D*(*ξ*)/|*D*(*ξ*)|, and the Gaussian curvature *κ* at *ξ* is, by definition, if the dimension *n*−1 of *X* is even, and else (since, for odd *n*−1, this determinant changes its sign in accordance with a sign change of *ν*). In other words,
From *ν*=*x*/|*x*|, we obtain that ∂_{j}*ν*−(∂_{j}*x*)/|*x*| is parallel to *x* and hence
Therefore, by equation (9.1),
9.2

Finally, if denotes the tangent map to *D* at *ξ* and its determinant with respect to orthonormal bases *e*_{1},…,*e*_{n−1}∈*T*_{ξ}*X* and *f*_{1},…,*f*_{n−1}∈*T*_{x}*X**, such that *x*,*e*_{1},…,*e*_{n−1} and *x*,*f*_{1},…,*f*_{n−1} induce the same orientation in **R**^{n}, then and therefore equation (9.2) implies the formula
9.3
Note that equation (9.3) entails the equation *κ*_{X}(*ξ*)⋅*κ*_{X*}(*x*)=|*ξ*|^{−n−1}|*x*|^{−n−1} we have mentioned in §8.

Next, let us assume that the duality map *x*=*Dξ* is given locally near *ξ*∈*X* as the resolution of the implicit equation *F*(*ξ*,*x*)=0 with being and submersive near (*ξ*,*x*). Then *D* is the restriction to *X* of a function defined in a neighbourhood *U* of *ξ* in **R**^{n}, and
(Note that *ξ*⊥*T*_{x}*X** and *x*,*f*_{1},…,*f*_{n−1} induces the same orientation as −*ξ*,*f*_{1},…,*f*_{n−1}.) If we denote by *F*_{ξ},*F*_{x} the *n*×*n* matrices consisting of the partial derivatives of *F* with respect to the *ξ* and the *x* variables, respectively, then and thus . Therefore, equation (9.3) implies, for odd *n*,
9.4

Let us now apply formula (9.4) to the Fresnel surface studied in §4, i.e. to *X*={*ξ*∈**R**^{3}; *R*(1,*ξ*)=0} with *R* as in equation (4.1). If we set *t*=*τ*=1 in equation (4.4), we obtain the following implicit equation *F*(*ξ*,*x*)=0 for the duality map
If *A* denotes—as in §2—the diagonal matrix with entries *a*_{1},*a*_{2},*a*_{3}, we can write *F* in short-hand as *F*(*ξ*,*x*)=*x*−|*ξ*|^{2}*Ax*+|*x*|^{2}*ξ*−*Aξ*.

For the matrix of partial derivatives of *F* with respect to *ξ*, we obtain
and hence *F*_{ξ}*x*=(*A*+|*x*|^{2}*I*_{3})*x* owing to *ξ*^{T}⋅*x*+1=0. Furthermore, setting *u*=−2(|*x*|^{2}*I*_{3}−*A*)^{−1}*Ax*, we have
Then equations (4.3) and (4.4) imply
and hence

Analogously, we obtain for the derivatives of *F* with respect to *x*
where *v*=2(*I*_{3}−|*ξ*|^{2}*A*)^{−1}*ξ* and *v*^{T}⋅*x*=−2, and hence

Finally, we have
and, due to *ξ*^{T}⋅*v*=0,
If we insert these equations into equation (9.4) and once more use equation (4.4), we arrive at *Darboux’s curvature formula*
9.5
cf. Darboux (1972, note VIII, (82), p. 484) and Liess (1991, (A.1), p. 91). (Of course, equation (9.5) holds only in those points where neither of the coordinates of *ξ* or of *η*, respectively, vanishes.)

For points (*t*,*x*)∈*Ξ**∩*H*, such that *ξ* is dual to (*t*,*x*), we then obtain by homogeneity
If we take into account that *ξ* is given by *ξ*_{j}=∂_{j}*R**/∂_{t}*R** owing to *X*^{**}=*X*, we infer the formula
9.6
In particular, when (*t*,*x*)∈*Ξ**∩*H* converges to a singular point of *Ξ**, then ∂_{2}*R**/*x*_{2} and *t*∂_{1}*R**/*x*_{1}∂_{t}*R**=*tξ*_{1}/*x*_{1} remain bounded, whereas ∂_{3}*R**/*x*_{3} converges to 0. This shows that the Radon measure *δ*(*R**) can be multiplied by its prefactor in equation (8.3).

## 10. The singular terms of the fundamental matrix

Let us now collect all the partial results on *E* found above. In the connectivity component *J*⊂**R**^{4} inside the innermost sheet of the wavefront surface *W*(*R*(*D*)) (figure 2), the fundamental matrix *E* coincides with *E*_{1} in equation (6.2). We observe that *J* is determined by the inequalities *t*>0, *R**>0, ∂_{t}*R**>0. On the other hand, *E* vanishes outside the propagation cone *K*=*K*(*R*(*D*)), which is the convex hull of the Fresnel surface *W*(*R*(*D*)). In *K*\*J*, the formulae (7.3), (8.3), (8.4) and (9.6) yield the singular parts of *E*.

### Theorem 10.1

*Let E be the retarded fundamental matrix (i.e. the one with support in t≥0) of the system of crystal optics* *as in equation (*2.2*). Let the polynomials R, R** *be as in equations (*4.1*) and (*4.5*), respectively, and denote by χ*_{J}*,χ*_{K\J}*=χ*_{K}*−χ*_{J} *the characteristic functions pertaining to the sets J={(t,x)∈***R**^{4}*; t>0,R***>0, ∂*_{t}*R***>0} and K=K(R(D))= convex hull of {(t,x)∈***R**^{4}*; t≥0, R***(t,x)=0}. Let the 3×3 matrix e(ξ)⋅e(ξ)*^{T} *be given as in equation (*7.3*) with y*_{j}*=a*_{j}*t*^{2}*−|x|*^{2}*.*

*Then E differs by a regular distribution* *from* *defined by
*
*(Here* vp *h denotes the principal value defined as the convergent limit* *in the sense of distribution theory.)*

Let us yet observe that the regular part of *E*—which is denoted by *f* in the theorem—is homogeneous of degree −2, as this is true for *E* and It is expressed by derivatives of Abelian integrals through formulae (5.3) and (7.2). A numerical procedure for evaluating such integrals is explained in Burridge & Qian (2006).

## Acknowledgements

For our continuous, fruitful cooperation, I am deeply indebted to Norbert Ortner, to whom this work too owes a lot.

- Received January 21, 2011.
- Accepted March 23, 2011.

- This journal is © 2011 The Royal Society