text:
# Horizons in 1+1 Dimensional Dilaton Gravity Coupled to Matter

Abstract

We study the local behaviour of static solutions of a general 1+1 dimensional dilaton gravity theory coupled to scalar fields and Abelian gauge fields near horizons. This type of model includes in particular reductions of higher dimensional theories invariant under a sufficiently large isometry group. The solution near the horizon can in general be obtained solving a system of integral equations or in favourable cases in the form of a convergent series in the dilaton field.

1. In the last decade 1+1 dimensional models of dilaton gravity coupled to scalar matter fields were found to describe some important properties of higher-dimensional black holesaaa The best studied example is the string inspired dilaton gravity of CGHS . For a review see .. The connection between high and low dimensions has been exploited in different contexts of gravity and string theory - symmetry reduction, compactification, holographic principle, AdS/CFT correspondence, duality, etc. (see e.g.  and references therein). An important case is given by spherically symmetric solutions of higher dimensional theories. In dimension 1+1, these models are usually not integrable, but a generalization of the CGHS and Jackiw – Teitelboim integrable models was introduced in . For its application to describing black hole evolution see e.g. . The most general class of integrable models of dilaton gravity coupled to scalar fields was recently proposed in . A further reduction from 1+1 to 0+1 dimension can be performed for time-independent solutions. The resulting 0+1 dimensional dilaton gravity models are frequently integrable in closed form (see e.g. , , ), although generically this is not the case. For example, static, spherical black holes carrying Abelian gauge fields are described by integrable models, while non-Abelian gauge fields  generically lead to non-integrable models (for an exception see ).

In general the considered 1+1 dimensional models have no time-independent globally regular solutions, but they have static black hole solutions. The most important properties of black holes (e.g. their mass, charge and their thermodynamics) are encoded in the structure of their horizon. For this reason, local properties of black hole horizons are of significant physical interest and merit some special investigation. Although in principle horizons are global concepts, not completely accessible by local considerations, this is not the case for Killing horizons of static black holes and this is what we shall consider henceforth. We will show how one can find the analytic behaviour of the geometry and of the matter fields near the horizon for the class of models considered. This does not require integrability of the equations, however the results are valid only locally.

2. We consider a fairly general DG described by the Lagrangian density

$$=\sqrt{g}[UR(g)+V+Wg^{ij}\varphi _i\varphi _j+\underset{k=1}{\overset{K}{}}X^{(k)}F_{ij}^{(k)}F^{(k)ij}+Y+\underset{n=1}{\overset{N}{}}Z^{(n)}g^{ij}\psi _i^{(n)}\psi _j^{(n)}],$$
(1)
where $`F_{ij}^{(k)}=_iA_j^{(k)}_jA_i^{(k)}`$, $`\varphi _i=_i\varphi `$, $`\psi _i=_i\psi `$, $`i,j=0,1`$. $`U,V,W`$ are arbitrary functions of $`\varphi `$ and $`X,Y,Z`$ – arbitrary functions of $`\varphi `$ and $`\psi `$. As is well known (and has been very simply and explicitly demonstrated in Refs. ), the pure dilaton gravity coupled to Abelian gauge fields $`F_{ij}`$ (when $`_\psi X0`$, $`YZ0`$) is explicitly integrable for arbitrary functions $`U(\varphi )`$, $`V(\varphi )`$, $`W(\varphi )`$, $`X(\varphi )`$. The general solutions are described by functions of one variable.

It is convenient to write the solutions in the conformal metric

$$ds^2=4f(u,v)dudv.$$
(2)
Then for the static solution we have

$$\varphi (u,v)=\varphi (\tau ),f(u,v)=h(\tau )a^{}(u)b^{}(v),\tau a(u)+b(v),$$
(3)
i.e. the metric is also defined by the function of one variable $`\tau `$, up to conformal transformations $`uA(u)`$, $`bB(v)`$. To further simplify the formulae we assume $`U(\varphi )\varphi `$. achievable by a suitable field redefinition in a domain, where $`U`$ is monotonous. In addition we apply a Weyl transformation to make $`W(\varphi )`$ zero, requiring $`U`$ (resp. $`\varphi `$) to be nonzero. Then  one finds

$$h(\tau )=MN(\varphi ),$$
(4)
where bbbNote that the equation of motion for $`F_{ij}^{(k)}`$ can be solved, because in dimension 1+1 the Maxwell fields do not propagate. Thus, their effect can be included in the effective action as a charge dependent piece in the potentials $`V(\varphi )`$ or $`Y(\varphi ,\psi )`$.

$$N(\varphi )=𝑑\varphi [V(\varphi )+\underset{k}{}2Q_k^2/X^{(k)}(\varphi )]$$
(5)
and the dependence of $`\varphi `$ on $`\tau `$ can be found from

$$𝑑\varphi [MN(\varphi )]^1=\tau \tau _0.$$
(6)
Here $`M`$, $`Q_k`$ and $`\tau _0`$ are integrals of motion and the solutions of the gauge field equations have the form $`F_{ab}^{(k)}=2Q_kh(\tau )/X^{(k)}\varphi )`$ (we consider only electric charges $`Q_k`$, magnetic ones being treated analogously).

Note that the static solution represents the general solution of pure dilaton gravity in terms of one free massless field $`\tau =a(u)+b(v)`$. So, in this case, we might reduce the 1+1 dimensional theory to 0+1 dimension, solve the equations for the variables $`g_{ij}(\tau )`$, $`\varphi (\tau )`$, $`A(\tau )`$, and then take $`\tau =a(u)+b(v)`$ with arbitrary functions $`a(u)`$, $`b(v)`$.

Also note that this general static solution has at least one horizon, because the equation $`MN(\varphi )=0`$ has always at least one root $`\varphi _0`$ for some value of $`M`$.

For some special functions $`U(\varphi )`$, $`V(\varphi )`$, $`W(\varphi )`$ and $`X(\varphi )`$ the model describes the spherical solution of the Einstein – Maxwell theory in any space-time dimension. Let the Lagrangian density in $`d`$ dimensions be

$$=\sqrt{G}[R(G)\frac{1}{4}\mathrm{e}^{\alpha \psi }F_{[2]}^2\frac{1}{2}(\psi )^2].$$
(7)
With $`\alpha =0`$, this is the Einstein-Maxwell theory minimally coupled to a scalar field $`\psi `$. In string theory this represents a consistent truncation of some $`d`$-dimensional supergravities (in this case instead of the antisymmetric Maxwell field two-form $`F_{[2]}`$ one has to introduce $`n`$-forms $`F_{[n]}`$, corresponding to the gauge potentials $`A_{[n1]}`$, see e.g. ). Two-dimensional dilaton-gravity theories may be obtained from (7) by different sorts of dimensional reductions. For illustrative purposes, it is sufficient to mention the spherical reduction, when the $`d`$-dimensional metric is

$$ds^2=g_{ij}dx^idx^j+\mathrm{e}^{4\nu \varphi }d\mathrm{\Omega }_n^2,$$
(8)
where $`\nu =1/n`$, $`n=d2`$, $`x^0=t`$, $`x^1=r`$; $`g_{ij}`$ and $`\varphi `$ depend only on $`t`$ and $`r`$; $`d\mathrm{\Omega }_n^2`$ is the metric on the $`n`$-dimensional sphere. By integrating out the angle variables in the action (7) we may get the 1+1 dimensional action (1) with

$`U`$ $`=`$ $`\mathrm{e}^{2\varphi },V=n(n1)\mathrm{e}^{2\varphi +4\nu \varphi },W=4(1\nu )\mathrm{e}^{2\varphi },`$
$`X`$ $`=`$ $`\beta \mathrm{e}^{2\varphi +\alpha \psi },Y=0,Z=\gamma \mathrm{e}^{2\varphi },`$ (9)
where $`\beta `$ and $`\gamma `$ are some positive constants depending on the dimension $`d`$.

Unlike the case of pure dilaton gravity, the general solution of the 1+1 dimensional theory with non-constant scalar matter field $`\psi `$ essentially depends on two variables and even simple theories of type (9) may not be integrable. Actually, integrability of the 1+1 dimensional theories is rather an exception than the rule (well-known integrable examples are CGHS, Jackiw-Teitelboim, Liouville and bi-Liouville models (see e.g. ). However, to find the static black hole solutions (and their brane generalizations) one has only to solve the 0+1 dimensional sector of the theory (1) and many of the 0+1 dimensional models emerging from higher dimensions by reductions are explicitly analytically integrable. For example reductions of (7) to 0+1 dimensions lead to integrable systems. However, to describe horizons locally it is not necessary to use any sort of global integrability. Thus in  it was shown that any theory (1) with the potentials $`U`$, $`V`$, $`W`$, $`X`$, $`Y`$, $`Z`$ depending only on $`\varphi `$ does not allow static solution with horizons for non-constant scalar fields, provided that the potentials satisfy some mild restrictions (absence of zeroes and poles).

3. In order to study the local behaviour near horizons let us first rewrite the Lagrangian density for the 0+1 dimensional sector of the theory (1). As already said, we put $`U(\varphi )`$ equal to $`\varphi `$ assuming that $`U(\varphi )`$ has no zeroes and is monotonous and use the Weyl transformation $`g_{ij}\mathrm{\Omega }(\varphi )g_{ij}`$ to make $`W=0`$. This transformation is regular, if $`U(\varphi )`$ has no zeroes. Finally, as explained above, we include the gauge fields $`F_{ij}^{(k)}`$ into the effective potential depending on $`Q_k`$. Then the Lagrangian density (1) may be rewritten as

$$=\sqrt{g}[\varphi R(g)+Y(\varphi ,\psi )+\underset{n=1}{\overset{N}{}}Z^{(n)}(\varphi ,\psi )g^{ij}\psi _i^{(n)}\psi _j^{(n)}].$$
(10)
For reasons of simplicity we consider only the case of a single matter field $`\psi `$, the generalization to any number of scalar fields being straightforward. The static Lagrangian density corresponding to (10) may be written in the form (see )

$$_{st}=\frac{1}{l(\tau )}[\varphi ^{}h^{}h^1+Z\psi _{}^{}{}_{}{}^{2}]+l(\tau )[hY(\varphi ,\psi )],$$
(11)
where $`l(\tau )`$ is a Lagrangian multiplier and the prime denotes the differentiation with respect to $`\tau `$.

To investigate the horizons it is convenient to use the Hamiltonian formulation (although our ‘time’ $`\tau `$ is in reality a space coordinate). Due to the constraint it is sufficient to consider the five first-order equations ($`\overline{Z}Z^1`$)

$`\dot{\varphi }`$ $`=`$ $`\chi (\mathrm{a});\dot{\chi }=hY(\mathrm{b});`$ (12)
$`\dot{\psi }`$ $`=`$ $`\overline{Z}\eta (\mathrm{a});\dot{\eta }={\displaystyle \frac{1}{2}}(hY_\psi +\eta ^2\overline{Z}_\psi )(\mathrm{b});`$ (13)
$`\dot{h}`$ $`=`$ $`{\displaystyle \frac{h}{\chi }}(hY+\overline{Z}\eta ^2).`$ (14)
Here the equations (12a) and (13a) are the definitions of the momentum variables $`\chi `$ and $`\eta `$, (14) is the constraint and the dot denotes now differentiation with respect to the new variable $`T𝑑tl(\tau )`$. The variable $`T`$ coincides with $`\tau `$ in the gauge $`l=1`$.

A horizon is characterized by the vanishing of $`h`$, while $`\varphi `$ and $`\psi `$ remain finite. Obviously this is impossible at any finite value of $`T`$ with $`\chi 0`$, because then $`h`$ would be identically zero. On the other hand, if $`\chi `$ and $`h`$ vanish at some finite value of $`T`$, then $`\chi `$ also vanishes identically. Therefore a horizon can only be obtained for $`T\pm \mathrm{}`$. This means that a horizon corresponds to a fixed point of the ‘dynamical system’ (12) - (14). Hence $`\chi `$ and $`\eta `$ have to tend to zero for $`T\pm \mathrm{}`$. In order to desingularize (14) we put $`h=\chi H`$ and $`\eta =\chi G`$ and use $`\chi `$ as independent variable to get

$`{\displaystyle \frac{d\varphi }{d\chi }}`$ $`=`$ $`{\displaystyle \frac{1}{HY}},{\displaystyle \frac{d\psi }{d\chi }}={\displaystyle \frac{\overline{Z}G}{HY}}`$ (15)
$`{\displaystyle \frac{d(\chi G)}{d\chi }}`$ $`=`$ $`{\displaystyle \frac{Y_\psi }{2Y}}+{\displaystyle \frac{G^2\overline{Z}_\psi \chi }{2HY}},{\displaystyle \frac{dH}{d\chi }}={\displaystyle \frac{\overline{Z}G^2}{Y}}.`$ (16)
In order to obtain a horizon at a finite value $`\varphi _0`$ of $`\varphi `$ we have to assume that $`\chi `$ vanishes there. Obviously the system (15) - (16) is singular at $`\chi =0`$ and a for regular solution the dependent variables have to fulfil further conditions. First we have to require $`Y_0=Y(\varphi _0,\psi _0)0`$, where $`\psi _0=\psi (0)`$. Furthermore for a non-degenrate horizon we need $`H_0=H(0)0`$. We replace the system (15) - (16) by a system of integral equations

$`\varphi `$ $`=`$ $`\varphi _0{\displaystyle _0^\chi }{\displaystyle \frac{d\chi ^{}}{HY}}`$ (17)
$`\psi `$ $`=`$ $`\psi _0{\displaystyle _0^\chi }{\displaystyle \frac{\overline{Z}Gd\chi ^{}}{HY}}`$ (18)
$`H`$ $`=`$ $`H_0+{\displaystyle _0^\chi }{\displaystyle \frac{\overline{Z}G^2}{Y}}𝑑\chi ^{}`$ (19)
$`G`$ $`=`$ $`{\displaystyle \frac{1}{\chi }}{\displaystyle _0^\chi }\left(\chi ^{}{\displaystyle \frac{G^2\overline{Z}_\psi }{2HY}}+{\displaystyle \frac{Y_\psi }{2Y}}\right)𝑑\chi ^{}`$ (20)
incorporating the regular behaviour at $`\chi =0`$. Starting with $`\varphi =\varphi _0,\psi =\psi _0,H=H_0,G=G_0=\frac{Y_\psi (0)}{2Y(0)}`$ we may iterate these equations to obtain a solution defined in a suitably small neighbourhood of $`\chi =0`$ . As is easily seen the solution inherits the smoothness properties of the functions $`Z(\varphi ,\psi )`$ resp. $`Y(\varphi ,\psi )`$. Thus in case the latter are analytic, we obtain a solution that is analytic at the horizon.

A somewhat different approach consists in constructing a power-series expansion in $`\chi `$ of the solution near the horizon. Using the dependent variables of Eqs. (12)-(14) we write the expansion

$`\eta `$ $`=`$ $`\eta _0\chi (1+\eta _1\chi +\eta _2\chi ^2+\mathrm{}),`$
$`h`$ $`=`$ $`h_0\chi (1+h_1\chi +h_2\chi ^2+\mathrm{}),`$
$`\varphi `$ $`=`$ $`\varphi _0+\varphi _1\chi +\varphi _2\chi ^2\mathrm{},`$
$`\psi `$ $`=`$ $`\psi _0+\psi _1\chi +\psi _2\chi ^2+\mathrm{}.`$ (21)
Correspondingly $`Y`$ and $`\overline{Z}`$ have to be expanded in powers of $`\chi `$

$`Y`$ $`=`$ $`Y_0(1+Y_1\chi +Y_2\chi +\mathrm{}),Y_\psi =Y_0^{}(1+Y_1^{}\chi +Y_2^{}\chi ^2+\mathrm{}),`$
$`\overline{Z}`$ $`=`$ $`\overline{Z}_0(1+\overline{Z}_1\chi +\overline{Z}_2\chi +\mathrm{}),\overline{Z}_\psi =\overline{Z}_{0}^{}{}_{}{}^{}(1+\overline{Z}_{1}^{}{}_{}{}^{}\chi +\overline{Z}_{2}^{}{}_{}{}^{}\chi ^2+\mathrm{}),`$ (22)
where $`Y_0`$, $`Y_1^{}`$, $`\overline{Z}_0`$, $`\overline{Z}_0^{}`$ depend only on $`\varphi _0`$, $`\psi _0`$ while $`Y_n`$, $`Y_n^{}`$, $`\overline{Z}_n`$, $`\overline{Z}_n^{}`$ depend on $`\varphi _m`$, $`\psi _m`$ with $`mn`$. In particular $`Y_0=Y(\varphi _0,\psi _0)`$, $`Y_0^{}=Y_\psi (\varphi _0,\psi _0)`$, etc.

Now, substituting the expansions (21), (22) into the equations

$$h^{}=(hY+\eta ^2\overline{Z})(\chi Y)^1,(\mathrm{a})\eta ^{}=(hY_\psi +\eta ^2\overline{Z}_\psi )(2hY)^1(\mathrm{b})$$
(23)
$$\varphi ^{}=\chi (hY)^1,(\mathrm{a})\psi ^{}=\eta \overline{Z}(hY)^1(\mathrm{b}).$$
(24)
we find recurrence relations, which can be solved order by order. The crucial observation is that $`h_0`$, $`\varphi _0`$, $`\psi _0`$ are not constrained by these relations and that $`\eta _0`$ is given by the zeroth order term in Eq. (23b)

$$\eta _0=Y_0^{}/2Y_0=Y_\psi (\varphi _0,\psi _0)/2Y(\varphi _0,\psi _0).$$
(25)
Then from Eqs. (23a), (24) we find $`h_1`$, $`\varphi _1`$ and $`\psi _1`$:

$$h_1=\frac{\eta _0^2\overline{Z}_0}{h_0Y_0}=\frac{Y^2\overline{Z}_0}{4h_0Y_0^3}\varphi _1=\frac{1}{h_0Y_0}\psi _1=\frac{Y_0^{}\overline{Z}_0}{2h_0Y_0^2}.$$
(26)
It is easy to check that this allows us to find $`\eta _1`$ from (23b) and then $`h_2`$, $`\varphi _2`$, $`\psi _2`$. In general, once we know $`\eta _{m1}`$, $`h_m`$, $`\varphi _m`$, $`\psi _m`$ for $`mn`$, we immediately find $`\eta _n`$ and then $`h_{n+1}`$, $`\varphi _{n+1}`$, $`\psi _{n+1}`$.

The convergence of the Taylor series in the analytic case can be inferred from the integral equation approach.

The dependence of the obtained solution on $`T`$ may be found from Eq.(12a) which now gives

$$\varphi _1\mathrm{ln}\chi +2\varphi _2\chi +\frac{3}{2}\varphi _3\chi ^2+\mathrm{}=T\tau _0.$$
(27)
This can be solved iteratively near the horizon. The zeroth approximation is $`\chi _0=\mathrm{exp}(TT_0)/\varphi _1`$ which is exponentially small for $`\varphi _1>0`$ if $`T\mathrm{}`$. Also the corrections decay exponentially near the horizon

$$\chi =\mathrm{e}^{(TT_0)/\varphi _1}(12\frac{\varphi _2}{\varphi _1}\mathrm{e}^{(TT_0)/\varphi _1}+\mathrm{}).$$

4. In the $`0+1`$ dimensional formulation the (Killing) horizon is reached at a singular point corresponding to its (‘unphysical’) bifurcation surface . In order to see that the solutions we have obtained really correspond to a regular horizon we have to approach it along light-like geodesics. This can be achieved using the light-like coordinates $`u`$ resp. $`v`$ together with a suitable space-like coordinate serving as an affine parameter for in- resp. outgoing light rays. For that purpose we rewrite the metric in the form

$$ds^2=4h(\widehat{\tau })du^22dud\widehat{\tau }\mathrm{resp}.ds^2=4h(\widehat{\tau })dv^2+2dvd\widehat{\tau }.$$
(28)
From Eq. (12b) we find $`d\widehat{\tau }=2/Yd\chi `$. Since $`YY_00`$ as $`\chi 0`$ we see that the regularity in $`\chi `$ in fact guarantees the regularity in $`\widehat{\tau }`$. In addition the metrical function $`h`$ has a simple zero in $`\widehat{\tau }`$ as required for a non-degenerate horizon.

5. In summary, we have found the solution near the horizon for a general dilaton gravity coupled to scalar matter. As one can see, the main condition for the existence of the obtained solutions is $`Y(\varphi _0,\psi _0)0`$, which means that the horizon is nondegenerate. The degenerate case can be solved by using a similar technique and will be treated elsewhere. The obtained formulae allow one to derive all geometric and physical characteristics of the horizon. Considering the 0+1 dimensional model as a reduction of certain higher dimensional theory one can also some characteristics of the horizon in the original higher dimensional theory. Note that our approach and results may also be applied to deriving asymptotic behaviour of some nonintegrable cosmological models.

This investigation was partly supported by the Russian Foundation for Basic Research (project no. 00-01-00299).

Appendix

Consider the simplest theory

$$=\frac{1}{l}(\dot{\varphi }\dot{h}h^1+\dot{\psi }^2)+lh2g\psi ^n.$$
$`(\mathrm{A}.1)`$
The Lagrange equations,

$$\ddot{\psi }=ghn\psi ^{n1},(\dot{h}/h)^{}=0,\ddot{\varphi }=2gh\psi ^n,$$
$`(\mathrm{A}.2)`$
give $`h=\mathrm{exp}(\alpha +\beta T)`$ and can easily be solved for $`n=1,2`$.ccc For $`n=2`$ solutions of the equation for $`\psi `$ can be expressed in terms of the Bessel functions of $`T`$.

In the simplest case, $`n=1`$, solving the equations and the constraint we have

$$\psi =\psi _0+c_1T+g\beta ^2\mathrm{e}^{\alpha +\beta T},$$

$$\varphi =\varphi _0+c_1^2\beta ^1T2gc_1\beta ^2\mathrm{e}^{\alpha +\beta T}+2g\beta ^2(2c_1\beta ^1\psi _0)\mathrm{e}^{\alpha +\beta T}g^2(2\beta ^4)^1\mathrm{e}^{2(\alpha +\beta T)}.$$
$`(\mathrm{A}.3)`$
Thus, we have a horizon for $`T\mathrm{}`$ if and only if $`c_1=0`$. Then the expansion of the solution in powers of $`h`$ is finite:

$$\chi =\dot{\varphi }=2gh\beta ^1\psi _0[1+gh(2\psi _0\beta ^2)^1],$$

$$\eta =\dot{\psi }=gh\beta ^1,\psi =\psi _0+gh\beta ^2,$$
$`(\mathrm{A}.4)`$
$$\varphi =\varphi _02gh\psi _0\beta ^2g^2h^2(2\beta ^4)^1.$$

One can rewrite it as an expansion in powers of $`\chi `$ but then it will be an infinite series, whose radius of convergence is $`\chi _c=\psi _{0}^{}{}_{}{}^{2}`$.

Thus, for $`\psi _0^20`$ the radius of convergence is vanishing. If $`\psi _0=0`$, we have a degenerate case in the following sense: $`\chi h^2`$ and $`(\varphi \varphi _0)h^2`$ This means that $`h(\varphi )`$ has a square root singularity at $`\varphi =\varphi _0`$. Such a singularity would require a stronger singularity in the potential $`V(\varphi )`$ ($`V|\varphi \varphi _0|^{\frac{1}{2}}`$) for the pure dilaton gravity and this is unnatural, while the degenerate horizon, $`h(\varphi \varphi _0)^n`$, where $`n`$ is a positive integer, is possible for analytic potentials. Note however that the black hole solution of the $`n=1`$ model (A.1) with $`\psi _0=0`$ is exceptional also for dilaton gravity coupled to matter, because in this case $`Y_0=0`$ and the expansion (21) is impossible (see Eq.(25)).

Of course, we may use a similar expansion in powers of $`h`$ and thus reproduce the exact solution. But finding conditions for the existence of degenerate solutions requires using more detailed information on the potentials $`Y`$, $`Z`$. For example, the model (A.1) with $`n=2`$ does not allow the degeneracy. In this case $`\psi `$ and $`\varphi \varphi _0`$ are convergent series in powers of $`2gh/\beta ^2`$:

$$\psi =\psi _0\underset{n=0}{\overset{\mathrm{}}{}}(2gh\beta ^2)^n(n!)^2,$$

and $`\varphi `$ can be obtained by integrating Eq.(A.2). Thus, they are analytic in $`h`$, and $`(\varphi \varphi _0)h`$, $`(\psi \psi _0)h`$.


format: text

text:
# Numerical methods for solving the time-dependent Maxwell equations1footnote 11footnote 1Invited talk presented at the International Computational Accelerator Physics Conference 2002

## I Introduction

The Maxwell equations describe the evolution of electromagnetic (EM) fields in space and time BornWolf . They apply to a wide range of different physical situations and play an important role in a large number of engineering applications. In many cases, numerical methods are required to solve Maxwell’s equations Taflove ; Kunz . A well-known class of algorithms is based on a method proposed by Yee Yee66 . This finite-difference time-domain (FDTD) approach owes its popularity mainly due to its flexibility and speed while at the same time it is easy to implement Taflove ; Kunz .

A limitation of Yee-based FDTD techniques is that their stability is conditional, depending on the mesh size of the spatial discretization and the time step of the time integration Taflove ; Kunz . Furthermore, in practice, the amount of computational work required to solve the time-dependent Maxwell equations by present FDTD techniques Taflove ; Kunz ; Website ; Zheng1 ; Namiki ; Zheng2 ; Garcia ; Harsh00 ; Kole01 ; Kole02  prohibits applications to a class of important fields such as bioelectromagnetics and VLSI design Taflove ; Gandi ; Houshmand . The basic reason for this is that the time step in the FDTD calculation has to be relatively small in order to maintain stability and a reasonable degree of accuracy in the time integration. Thus, the search for new algorithms that solve the Maxwell equation focuses on removing the conditional stability of FDTD methods and on improving the accuracy/efficiency of the algorithms.

A systematic approach to construct unconditionally stable algorithms is to employ a Suzuki product-formula Suzuki77  to approximate the time evolution operator DeRaedt87 . In the case of EM fields, the latter is the matrix exponential of a skew-symmetric matrix and the approximations take the form of products of orthogonal transformations Kole01 ; Kole02 . The resulting numerical algorithms are unconditionally stable by construction DeRaedt87 ; Smith85 .

The spectral-domain split-operator technique proposed in Ref.Harsh00  is one of the many forms that results from the use of the Lie-Trotter-Suzuki product formulas Suzuki77 . This technique makes use of Fast Fourier Transforms to compute the matrix exponentials of the displacement operators. The choice made in Ref. Harsh00  yields an approximation to the time-evolution operator that is no longer orthogonal and hence unconditional stability is not automatically guaranteed COMMENT1 . In contrast, the methodology proposed in Refs. Kole01 ; Kole02  yields efficient, explicit, unconditionally stable schemes that operate on the EM fields defined on the real space grid only. These algorithms naturally allow for the spatial variations in both the permittivity and the permeability.

The Suzuki product-formula approach also provides a unified framework to construct and analyse other time stepping algorithms DeRaedt87 ; Horvath . To illustrate this point we show that the conventional Yee algorithm and the alternating-direction-implicit (ADI) time-stepping algorithms Zheng1 ; Namiki ; Zheng2 ; Garcia ; Horvath  fit into this framework. Furthermore we propose new variants of the Yee algorithm.

Another route to improve upon the accuracy/efficiency of time-integration schemes is to make use of the Chebyshev polynomial expansions of the matrix exponential TAL-EZER0 ; TAL-EZER ; LEFOR ; Iitaka01 ; SILVER ; LOH . In this paper we also discuss a one-step algorithm, based on Chebyshev polynomials, to solve the time-dependent Maxwell equations for (very) large time steps hdrcheb0 ; hdrcheb1 .

The main purpose of this paper is to review the basic ideas behind the recent developments in numerical algorithms to solve the time-dependent Maxwell equations and to compare the virtues and shortcomings of the different methods. The plan of the paper is as follows. In Section II we briefly discuss the basic physical symmetries of the time-dependent Maxwell equations. The general framework to construct time-integration algorithms is laid out in Section III. We also pay attention to the numerical treatment of the current source term. In Section IV we use the simplest case of the time-dependent Maxwell equations to illustrate how the various algorithms can be implemented. We explicitly show how the conventional Yee algorithm naturally fits into this framework and, by minor modification, construct second-order and fourth-order time-accurate schemes that do not require the use of staggered-in-time fields, nor extra memory to store intermediate results. Then we recall the steps to construct the unconditionally stable algorithms proposed in Ref. Kole01 ; Kole02  and analyse a modification to improve the time-integration accuracy. Finally we discuss the implementation of the ADI and one-step algorithms. A discussion of the results of numerical experiments and our conclusions are given in Section V and VI respectively.

## II Theory

We consider EM fields in linear, isotropic, nondispersive and lossless materials. The time evolution of EM fields in these systems is governed by the time-dependent Maxwell equations BornWolf . Some important physical symmetries of the Maxwell equations can be made explicit by introducing the fields

$$𝐗(t)\sqrt{\mu }𝐇(t)\mathrm{and}𝐘(t)\sqrt{\epsilon }𝐄(t).$$
(1)
Here, $`𝐇(t)=(H_x(𝐫,t),H_y(𝐫,t),H_z(𝐫,t))^T`$ denotes the magnetic and $`𝐄(t)=(E_x(𝐫,t),E_y(𝐫,t),E_z(𝐫,t))^T`$ the electric field vector, while $`\mu =\mu (𝐫)`$ and $`\epsilon =\epsilon (𝐫)`$ denote, respectively, the permeability and the permittivity. In the absence of electric charges, Maxwell’s curl equations Taflove  read

$$\frac{}{t}\left(\begin{array}{c}𝐗(t)\\ 𝐘(t)\end{array}\right)=\left(\begin{array}{c}𝐗(t)\\ 𝐘(t)\end{array}\right)\frac{1}{\sqrt{\epsilon }}\left(\begin{array}{c}0\\ 𝐉(t)\end{array}\right),$$
(2)
where $`𝐉(t)=(J_x(𝐫,t),J_y(𝐫,t),J_z(𝐫,t))^T`$ represents the source of the electric field and $``$ denotes the operator

$$\left(\begin{array}{cc}0& \frac{1}{\sqrt{\mu }}\times \frac{1}{\sqrt{\epsilon }}\\ \frac{1}{\sqrt{\epsilon }}\times \frac{1}{\sqrt{\mu }}& 0\end{array}\right).$$
(3)
Writing $`𝐙(t)=(𝐗(t),𝐘(t))^T`$ it is easy to show that $``$ is skew symmetric, i.e. $`^T=`$, with respect to the inner product $`𝐙(t)|𝐙^{}(t)_V𝐙^T(t)𝐙^{}(t)𝑑𝐫`$, where $`V`$ denotes the system’s volume. In addition to Eq.(2), the EM fields also satisfy $`(\sqrt{\mu }𝐗(t))=0`$ and $`(\sqrt{\epsilon }𝐘(t))=0`$ BornWolf . Throughout this paper we use dimensionless quantities: We measure distances in units of $`\lambda `$ and expresss time and frequency in units of $`\lambda /c`$ and $`c/\lambda `$, respectively.

A numerical algorithm that solves the time-dependent Maxwell equations necessarily involves some discretization procedure of the spatial derivatives in Eq. (2). Ideally, this procedure should not change the basic symmetries of the Maxwell equations. We will not discuss the (important) technicalities of the spatial discretization (we refer the reader to Refs. Taflove ; Kunz ) as this is not essential to the discussion that follows. On a spatial grid Maxwell’s curl equations (2) can be written in the compact form Kole01

$$\frac{}{t}𝚿(t)=H𝚿(t)𝚽(t).$$
(4)
The vector $`𝚿(t)`$ is a representation of $`𝐙(t)`$ on the grid. The matrix $`H`$ is the discrete analogue of the operator (3), and the vector $`𝚽(t)`$ contains all the information on the current source $`𝐉(t)`$. The formal solution of Eq. (4) is given by

$$𝚿(t)=U(t)𝚿(0)_0^tU(tu)𝚽(u)𝑑u,$$
(5)
where

$$U(t)=e^{tH},$$
(6)
denotes the time-evolution matrix. If the discretization procedure preserves the underlying symmetries of the time-dependent Maxwell equations then the matrix $`H`$ is real and skew symmetric Kole01 , implying that $`U(t)`$ is orthogonal WILKINSON . Physically, the orthogonality of $`U(t)`$ implies conservation of energy Kole01 .

## III Time integration algorithms

There are two, closely related, strategies to construct an algorithm for performing the time integration of the time-dependent Maxwell equations defined on the grid Smith85 . The traditional approach is to discretize (with increasing level of sophistication) the derivative with respect to time Smith85 . The other is to approximate the formally exact solution, i.e. the matrix exponential $`U(t)=e^{tH}`$ by some time evolution matrix $`\stackrel{~}{U}(t)`$ Smith85 ; DeRaedt87 . We adopt the latter approach in this paper as it facilitates the construction of algorithms with specific features, such as unconditional stability DeRaedt87 ; Kole01 .

If the approximation $`\stackrel{~}{U}(t)`$ is itself an orthogonal transformation, then $`\stackrel{~}{U}(t)=1`$ where $`X`$ denotes the 2-norm of a vector or matrix $`X`$ WILKINSON . In the absence of source terms (i.e. $`𝚽(t)=\mathrm{𝟎}`$) this implies that $`\stackrel{~}{𝚿}(t)=\stackrel{~}{U}(t)𝚿(0)=𝚿(0)`$, for an arbitrary initial condition $`𝚿(0)`$ and for all times $`t`$ and hence the time integration algorithm defined by $`\stackrel{~}{U}(t)`$ is unconditionally stable by construction Smith85 ; DeRaedt87 .

In the presence of current sources, for general $`\stackrel{~}{U}(t)`$, it follows immediately from Eq.(5) that

$`\stackrel{~}{𝚿}(t)`$ $``$ $`𝚿(t)+\stackrel{~}{ϵ}\left(𝚿(0)+{\displaystyle _0^t}𝚽(u)𝑑u\right),`$ (7)
where $`\stackrel{~}{U}(u)U(u)\stackrel{~}{ϵ}`$ for $`0ut`$ and $`\stackrel{~}{ϵ}`$ is a measure for the accuracy of the approximation $`\stackrel{~}{U}(t)`$.

From Eq.(5) it follows that the EM fields $`𝚿(t)`$ change according to

$$𝚿(t+\tau )=e^{\tau H}𝚿(t)_t^{t+\tau }e^{(t+\tau u)H}𝚽(u)𝑑u.$$
(8)
In the time-stepping approach we approximate the source term in Eq.(8) by the standard 3-point Gauss-Legendre quadrature formula ABRAMOWITZ

$$𝚿(t+\tau )=e^{\tau H}𝚿(t)+\frac{\tau }{2}\underset{i=0}{\overset{2}{}}w_ie^{(1+x_i)\tau H/2}𝚽(t+(1+x_i)\tau /2)+𝒪(\tau ^7),$$
(9)
where $`x_0,x_1,x_2`$ are the zeros of the Legendre polynomial $`P_3(x)=x(5x^23)/2`$ and $`w_i=8/(1x_i^2)(15x_i^23)^2`$ ABRAMOWITZ . In practice we replace $`e^{(1+x_i)\tau H/2}`$ in Eq.(9) by an approximation $`\stackrel{~}{U}((1+x_i)\tau /2)`$.

We now consider three options to construct the approximate time evolution matrix $`\stackrel{~}{U}(t)`$. We exclude from the discussion the exceptional cases for which the matrix exponential $`U(t)=e^{tH}`$ can be calculated directly, as these are usually of little relevance for realistic problems. The first approach yields the conventional Yee algorithm, a higher-order generalization thereof, and the unconditional schemes proposed in Ref.Kole01 . The second option is to use rational approximations to the exponential, yielding the standard ADI methods. Finally, the Chebyshev polynomial approximation to the matrix exponential is used to construct a one-step algorithm.

### III.1 Suzuki product-formula approach

A systematic approach to construct approximations to matrix exponentials is to make use of the Lie-Trotter-Suzuki formula Suzuki77 ; Trotter59

$$e^{tH}=e^{t(H_1+\mathrm{}+H_p)}=\underset{m\mathrm{}}{lim}\left(\underset{i=1}{\overset{p}{}}e^{tH_i/m}\right)^m,$$
(10)
and generalizations thereof Suzuki8591 ; DeRaedt83 . Expression Eq. (10) suggests that

$$U_1(\tau )=e^{\tau H_1}\mathrm{}e^{\tau H_p},$$
(11)
might be a good approximation to $`U(\tau )`$ if $`\tau `$ is sufficiently small. Applied to the case of interest here, if all the $`H_i`$ are real and skew-symmetric $`U_1(\tau )`$ is orthogonal by construction and a numerical scheme based on Eq. (11) will be unconditionally stable. For orthogonal matrices $`U(\tau )`$ and $`U_1(\tau )`$ it can be shown that DeRaedt87

$$U(\tau )U_1(\tau )\frac{\tau ^2}{2}\underset{i<j}{}[H_i,H_j],$$
(12)
where $`[H_i,H_j]=H_iH_jH_jH_i`$ is, in general, non-zero. Relaxing the condition that $`U(\tau )`$ and $`U_1(\tau )`$ are orthogonal matrices changes the $`\tau `$ dependence in Eq. (12) but for small $`\tau `$ the error still vanishes like $`\tau ^2`$ Suzuki8591 . From Eq. (12) it follows that, in general, the Taylor series of $`U(\tau )`$ and $`U_1(\tau )`$ are identical up to first order in $`\tau `$. We will call $`U_1(\tau )`$ a first-order approximation to $`U(\tau )`$.

The product-formula approach provides simple, systematic procedures to improve the accuracy of the approximation to $`U(\tau )`$ without changing its fundamental symmetries. For example the matrix

$$U_2(\tau )=U_1(\tau /2)^TU_1(\tau /2)=e^{\tau H_p/2}\mathrm{}e^{\tau H_1/2}e^{\tau H_1/2}\mathrm{}e^{\tau H_p/2},$$
(13)
is a second-order approximation to $`U(\tau )`$ Suzuki8591 ; DeRaedt83 . If $`U_1(\tau )`$ is orthogonal, so is $`U_2(\tau )`$. For orthogonal $`U_2(\tau )`$ we have

$$U(t)[U_2(t/m)]^mc_2\tau ^2t,$$
(14)
where $`c_2`$ is a positive constant DeRaedt87 .

Suzuki’s fractal decomposition approach Suzuki8591  gives a general method to construct higher-order approximations based on $`U_2(\tau )`$ (or $`U_1(\tau )`$). A particularly useful fourth-order approximation is given by Suzuki8591

$$U_4(\tau )=U_2(a\tau )U_2(a\tau )U_2((14a)\tau )U_2(a\tau )U_2(a\tau ),$$
(15)
where $`a=1/(44^{1/3})`$. The approximations Eqs.(11) and (13), and (15) have proven to be very useful in many applications Suzuki77 ; Chorin78 ; DeRaedt83 ; DeRaedt87 ; Koboyashi94 ; DeRaedt94 ; Rouhi95 ; Shadwick97 ; Krech98 ; Tran98 ; Michielsen98 ; DeRaedt00  and, as we show below, turn out to be equally useful for solving the time-dependent Maxwell equations. As before, for orthogonal $`U_4(\tau )`$ we have DeRaedt87

$$U(t)[U_4(t/m)]^mc_4\tau ^4t,$$
(16)
where $`c_4`$ is a positive constant.

As our numerical results (see below) show, for sufficiently small $`\tau `$, the numerical error of a time integrator vanishes with $`\tau `$ according to the $`\tau `$-dependence of the corresponding rigorous bound, e.g. Eqs.(12),(14), or (16). Our experience shows that if this behavior is not observed, there is a fair chance that the program contains one or more errors.

In practice an efficient implementation of the first-order scheme is all that is needed to construct the higher-order algorithms Eqs.(13) and (15). The crucial step of this approach is to choose the $`H_i`$’s such that the matrix exponentials $`\mathrm{exp}(\tau H_1)`$, …, $`\mathrm{exp}(\tau H_p)`$ can be calculated efficiently. This will turn the formal expressions for $`U_2(\tau )`$ and $`U_4(\tau )`$ into efficient algorithms to solve the time-dependent Maxwell equations.

### III.2 ADI algorithms

Instead of hunting for a decomposition that leads to matrix exponentials $`\mathrm{exp}(\tau H_1)`$, …, $`\mathrm{exp}(\tau H_p)`$ that are easy to compute, one can opt for an algorithm in which each of these exponentials is calculated approximately. In principle this might be beneficial because there is more flexibility in decomposing $`H`$. The standard strategy, preserving the symmetry of $`H_1`$, …, $`H_p`$ is to use rational (Padé) approximations to the exponential Smith85 . For instance, the approximation $`e^x(1+x/2)/(1x/2)`$ with some decompostion $`H=H_1+H_2`$ yields the second-order-accurate ADI algorithm Smith85 ; NumericalRecipes ; Horvath

$$U_2^{ADI}(\tau )=(I\tau H_1/2)^1(I+\tau H_2/2)(I\tau H_2/2)^1(I+\tau H_1/2),$$
(17)
where $`I`$ is the identity matrix. As the subscript indicates, the algorithm (17) is second-order accurate in time. For general skew-symmetrix $`H_1`$ and $`H_2`$, it is easy to show that the algorithm $`U_2^{ADI}(\tau )`$ is unconditionally stable. Following Ref. Horvath  we rearrange factors and obtain

$`\left[U_2^{ADI}(\tau )\right]^m`$ $`=`$ $`(I\tau H_1/2)^1X_2X_1X_2\mathrm{}X_1X_2(I+\tau H_1/2)`$ (18)
$``$ $`(I\tau H_1/2)^1X_2X_1X_2\mathrm{}X_1X_2(I+\tau H_1/2)=(I\tau H_1/2)^1(I+\tau H_1/2).`$
We used the fact that for skew-symmetric $`H_i`$, $`X_i(I\tau H_i/2)^1(I+\tau H_i/2)`$ is orthogonal and that $`X_2X_1X_2\mathrm{}X_1X_2^2=1`$. If $`X`$ is skew-symmetric, it’s eigenvalues are pure imaginary and therefore $`(IX)^1`$ is non-singular. Hence, for any number of time steps $`m`$, $`\left[U_2^{ADI}(\tau )\right]^mC`$, where $`C`$ is some finite positive constant, proving that the $`U_2^{ADI}(\tau )`$ is unconditionally stable in the Lax-Richtmyer sense Smith85 .

The matrix inversions appearing in Eq.(17) suggest that for practical purposes the implicit method $`U_2^{ADI}(\tau )`$ will not be very useful unless $`I\tau H_1/2`$ and $`I\tau H_2/2`$ take special forms that allow efficient matrix inversion Smith85 ; NumericalRecipes .

### III.3 One-step algorithm

The basic idea of this approach is to make use of extremely accurate polynomial approximations to the matrix exponential. First we use the Chebyshev polynomial expansion to approximate $`U(t)`$ and then show how to treat the source term in Eq. (5). We begin by “normalizing” the matrix $`H`$. The eigenvalues of the skew-symmetric matrix $`H`$ are pure imaginary numbers. In practice $`H`$ is sparse so it is easy to compute $`H_1\mathrm{max}_j_i|H_{i,j}|`$. Then, by construction, the eigenvalues of $`BiH/H_1`$ all lie in the interval $`[1,1]`$ WILKINSON . Expanding the initial value $`𝚿(0)`$ in the (unknown) eigenvectors $`𝐛_j`$ of $`B`$, we find from Eq. (5) with $`𝚽(t)0`$:

$$𝚿(t)=e^{izB}𝚿(0)=\underset{j}{}e^{izb_j}𝐛_j𝐛_j|𝚿(0),$$
(19)
where $`z=tH_1`$ and the $`b_j`$ denote the (unknown) eigenvalues of $`B`$. There is no need to know the eigenvalues and eigenvectors of $`B`$ explicitly. We find the Chebyshev polynomial expansion of $`U(t)`$ by computing the expansion coefficients of each of the functions $`e^{izb_j}`$ that appear in Eq. (19). In particular, as $`1b_j1`$, we can use the expansion ABRAMOWITZ  $`e^{izb_j}=J_0(z)+2_{k=1}^{\mathrm{}}i^kJ_k(z)T_k(b_j),`$ where $`J_k(z)`$ is the Bessel function of integer order $`k`$ to write Eq. (19) as

$$𝚿(t)=\left[J_0(z)I+2\underset{k=1}{\overset{\mathrm{}}{}}J_k(z)\widehat{T}_k(B)\right]𝚿(0).$$
(20)
Here $`\widehat{T}_k(B)=i^kT_k(B)`$ is a matrix-valued modified Chebyshev polynomial that is defined by the recursion relations

$$\widehat{T}_0(B)𝚿(0)=𝚿(0),\widehat{T}_1(B)𝚿(0)=iB𝚿(0),$$
(21)
and

$$\widehat{T}_{k+1}(B)𝚿(0)=2iB\widehat{T}_k(B)𝚿(0)+\widehat{T}_{k1}(B)𝚿(0),$$
(22)
for $`k1`$.

As $`\widehat{T}_k(B)1`$ by construction and $`|J_k(z)||z|^k/2^kk!`$ for $`z`$ real ABRAMOWITZ , the resulting error vanishes exponentially fast for sufficiently large $`K`$. Thus, we can obtain an accurate approximation by summing contributions in Eq. (20) with $`kK`$ only. The number $`K`$ is fixed by requiring that $`|J_k(z)|<\kappa `$ for all $`k>K`$. Here, $`\kappa `$ is a control parameter that determines the accuracy of the approximation. For fixed $`\kappa `$, $`K`$ increases linearly with $`z=tH_1`$ (there is no requirement on $`t`$ being small). From numerical analysis it is known that for fixed $`K`$, the Chebyshev polynomial is very nearly the same polynomial as the minimax polynomial NumericalRecipes , i.e. the polynomial of degree $`K`$ that has the smallest maximum deviation from the true function, and is much more accurate than for instance a Taylor expansion of the same degree $`K`$. In Fig.1 we show a plot of $`J_n(z=200)`$ as a function of $`n`$ to illustrate these points. From Fig.1 it is clear that the Chebyshev polynomial expansion will only be useful if $`K`$ lies to the right of the right-most extremum of $`J_n(z=200)`$, i.e. $`K`$ has to be larger than 200 in this example.

We now turn to the treatment of the current source $`𝐉(t)`$. The contribution of the source term to the EM field at time $`t`$ is given by the last term in Eq. (5). For simplicity we only consider the case of a sinusoidal source

$$𝐉(𝐫,t)=\mathrm{\Theta }(t)\mathrm{\Theta }(Tt)𝐬(𝐫)\mathrm{sin}(\mathrm{\Omega }t),$$
(23)
where $`𝐬(𝐫)`$ specifies the spatial distribution and $`\mathrm{\Omega }2\pi f_s`$ the angular frequency of the source. The step functions $`\mathrm{\Theta }(t)`$ and $`\mathrm{\Theta }(Tt)`$ indicate that the source is turned on at $`t=0`$ and is switched off at $`t=T`$. Note that Eq. (23) may be used to compose sources with a more complicated time dependence by a Fourier sine transformation.

The formal expression for the contribution of the sinusoidal source (23) reads

$`{\displaystyle _0^t}e^{(tu)H}𝚽(u)𝑑u`$ $`=`$ $`(\mathrm{\Omega }^2+H^2)^1e^{(tT^{})H}\times (\mathrm{\Omega }e^{T^{}H}\mathrm{\Omega }\mathrm{cos}\mathrm{\Omega }T^{}H\mathrm{sin}\mathrm{\Omega }T^{})𝚵`$ (24)
$``$ $`f(H,t,T^{},\mathrm{\Omega })𝚵,`$
where $`T^{}=\mathrm{min}(t,T)`$ and $`𝚽(u)\mathrm{\Theta }(t)\mathrm{\Theta }(Tt)\mathrm{sin}(\mathrm{\Omega }t)𝚵`$. The vector $`𝚵`$ represents the spatial (time-independent) distribution $`𝐬(𝐫)`$ and has the same dimension as $`𝚿(0)`$. The coefficients of the Chebyshev polynomial expansion of the formal solution (24) are calculated as follows. First we repeat the scaling procedure described above and substitute in Eq. (24) $`H=ixH_1`$, $`t=z/H_1`$, $`T^{}=Z^{}/H_1`$, and $`\mathrm{\Omega }=\omega H_1`$. Then, we compute the (Fast) Fourier Transform with respect to $`x`$ of the function $`f(x,z,Z^{},\omega )`$ (which is non-singular on the interval $`1x1`$). By construction, the Fourier coefficients $`S_k(tH_1)`$ are the coefficients of the Chebyshev polynomial expansion ABRAMOWITZ .

Taking into account all contributions of the source term with $`k`$ smaller than $`K^{}`$ (determined by a procedure similar to the one for $`K`$), the one-step algorithm to compute the EM fields at time $`t`$ reads

$`\widehat{𝚿}(t)`$ $`=`$ $`\left[J_0(tH_1)I+2{\displaystyle \underset{k=1}{\overset{K}{}}}J_k(tH_1)\widehat{T}_k(B)\right]𝚿(0)+\left[S_0(tH_1)I+2{\displaystyle \underset{k=1}{\overset{K^{}}{}}}S_k(tH_1)\widehat{T}_k(B)\right]𝚵.`$ (25)
Note that in this one-step approach the time dependence of the source is taken into account exactly, without actually sampling it.

In a strict sense, the one-step method does not yield an orthogonal approximation. However, for practical purposes it can be viewed as an extremely stable time-integration algorithm because it yields an approximation to the exact time evolution operator $`U(t)=e^{tH}`$ that is exact to nearly machine precision, i.e. in practice the value of $`\stackrel{~}{ϵ}`$ in Eq.(7) is very small. This also implies that within the same precision $`(\mu 𝐇(t))=(\mu 𝐇(t=0))`$ and $`(\epsilon 𝐄(t))=(\epsilon 𝐄(t=0))`$ holds for all times, implying that the numerical scheme will not produce artificial charges during the time integration Taflove ; Kunz .

## IV Implementation

The basic steps in the construction of the product-formula and one-step algorithms are best illustrated by considering the simplest case, i.e. the Maxwell equations of a 1D homogeneous problem. From a conceptual point of view nothing is lost by doing this: the extension to 2D and 3D nonhomogeneous problems is straigthforward, albeit technically non-trivial Kole01 ; Kole02 ; hdrcheb0 ; hdrcheb1 .

We consider a system, infinitely large in the $`y`$ and $`z`$ direction, for which $`\epsilon =1`$ and $`\mu =1`$. Under these conditions, the Maxwell equations reduce to two independent sets of first-order differential equations BornWolf , the transverse electric (TE) mode and the transverse magnetic (TM) mode BornWolf . As the equations of the TE- and TM-mode differ by a sign we can restrict our considerations to the TM-mode only. The magnetic field $`H_y(x,t)`$ and the electric field $`E_z(x,t)`$ of the TM-mode in the 1D cavity of length $`L`$ are solutions of

$`{\displaystyle \frac{}{t}}H_y(x,t)`$ $`=`$ $`{\displaystyle \frac{}{x}}E_z(x,t),`$ (26)
$`{\displaystyle \frac{}{t}}E_z(x,t)`$ $`=`$ $`{\displaystyle \frac{}{x}}H_y(x,t)J_z(x,t),`$ (27)
subject to the boundary condition $`E_z(0,t)=E_z(L,t)=0`$ BornWolf . Note that the divergence of both fields are trivially zero.

Following Yee Yee66 , to discretize Eqs.(26) and (27), it is convenient to assign $`H_y`$ to odd and $`E_z`$ to even numbered lattice sites, as shown in Fig. 2. Using the second-order central-difference approximation to the first derivative with respect to $`x`$, we obtain

$`{\displaystyle \frac{}{t}}H_y(2i+1,t)`$ $`=`$ $`\delta ^1(E_z(2i+2,t)E_z(2i,t)),`$ (28)
$`{\displaystyle \frac{}{t}}E_z(2i,t)`$ $`=`$ $`\delta ^1(H_y(2i+1,t)H_y(2i1,t))J_z(2i,t),`$ (29)
where we have introduced the notation $`A(i,t)=A(x=i\delta /2,t)`$. The integer $`i`$ labels the grid points and $`\delta `$ denotes the distance between two next-nearest neighbors on the lattice (hence the absence of a factor two in the nominator). We define the $`n`$-dimensional vector $`𝚿(t)`$ by

$$𝚿(i,t)=\{\begin{array}{ccc}H_y(i,t),\hfill & i\text{ odd}\hfill & \\ E_z(i,t),\hfill & i\text{ even}\hfill & \end{array}.$$
(30)
The vector $`𝚿(t)`$ contains both the magnetic and the electric field on the lattice points $`i=1,\mathrm{},n`$. The $`i`$-th element of $`𝚿(t)`$ is given by the inner product $`𝚿(i,t)=𝐞_i^T𝚿(t)`$ where $`𝐞_i`$ denotes the $`i`$-th unit vector in the $`n`$-dimensional vector space. Using this notation (which proves most useful for the case of 2D and 3D for which it is rather cumbersome to write down explicit matrix representations), it is easy to show that Eqs.(28) and (29) can be written in the form (4) where the matrix $`H`$ is given by

$$H=\delta ^1\underset{i=1}{\overset{n1}{}}\left(𝐞_i^{}𝐞_{i+1}^T𝐞_{i+1}^{}𝐞_i^T\right)=\left(\begin{array}{ccccc}0& \delta ^1& & & \\ \delta ^1& 0& \delta ^1& & \\ & \mathrm{}& \mathrm{}& \mathrm{}& \\ & & \delta ^1& 0& \delta ^1\\ & & & \delta ^1& 0\end{array}\right).$$
(31)
We immediately see that $`H`$ is sparse and skew-symmetric by construction.

### IV.1 Yee-type algorithms

First we demonstrate that the Yee algorithm fits into the product-formula approach. For the 1D model (31) it is easy to see that one time-step with the Yee algorithm corresponds to the operation

$`U_1^{Yee}(\tau )=(I+\tau A)(I\tau A^T)=e^{\tau A}e^{\tau A^T},`$ (32)
where

$`A`$ $`=`$ $`\delta ^1{\displaystyle \stackrel{}{}_{i=2}^{n1}}\left(𝐞_i^{}𝐞_{i1}^T𝐞_i^{}𝐞_{i+1}^T\right)=\left(\begin{array}{cccccc}0& \delta ^1& 0& 0& 0& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}\\ 0& \delta ^1& 0& \delta ^1& 0& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}\\ 0& 0& 0& \delta ^1& 0& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right),`$ (39)
and we used the arrangements of $`H`$ and $`E`$ fields as defined by Eq.(30). We use the notation $`^{}`$ to indicate that the stride of the summation index is two. Note that since $`A^2=0`$ we have $`e^{\tau A}=1+\tau A`$ exactly. Therefore we recover the time-step operator of the Yee algorithm using the first-order product formula approximation to $`e^{\tau H}`$ and decomposing $`H=AA^T`$. However, the Yee algorithm is second-order, not first order, accurate in time Taflove ; Kunz . This is due to the use of a staggered grid in time Taflove ; Kunz . To perform one time step with the Yee algorithm we need to know the values of $`E_z(t)`$ and $`H_y(t+\tau /2)`$, not $`H_y(t)`$. Another method has to supply the $`H_y`$-field at a time shifted by $`\tau /2`$.

Within the spirit of this approach, we can easily eliminate the staggered-in-time grid at virtually no extra computational cost or progamming effort (if a conventional Yee code is available) by using the second-order product formula

$`U_2^{Yee}(\tau )=e^{\tau A/2}e^{\tau A^T}e^{\tau A/2}=(I+\tau A/2)(I\tau A^T)(I+\tau A/2).`$ (40)
The effect of the last factor is to propagate the $`H_y`$-field by $`\tau /2`$. The middle factor propagates the $`E_z`$-field by $`\tau `$. The first factor again propagates the $`H_y`$ field by $`\tau /2`$. In this scheme all EM fields are to be taken at the same time. The algorithm defined by $`U_2^{Yee}(\tau )`$ is second-order accurate in time by construction DeRaedt87 . Note that $`e^{\tau A/2}`$ is not orthogonal so nothing has been gained in terms of stability. Since $`\left(U_2^{Yee}(\tau )\right)^m=e^{\tau A/2}\left(U_1^{Yee}(\tau )\right)^me^{+\tau A/2}`$, we see that, compared to the original Yee algorithm, the extra computational work is proportional to $`(1+2/m)`$, hence negligible if the number of time steps $`m`$ is large.

According to the general theory outlined in Sec.III, the expression

$`U_4^{Yee}(\tau )=U_2^{Yee}(a\tau )U_2^{Yee}(a\tau )U_2^{Yee}((14a)\tau )U_2^{Yee}(a\tau )U_2^{Yee}(a\tau ),`$ (41)
defines a fourth-order accurate Yee-like scheme, the realization of which requires almost no effort once $`U_2^{Yee}`$ has been implemented. It is easy to see that the above construction of the Yee-like algorithms holds for the much more complicated 2D, and 3D inhomogeneous case as well. Also note that the fourth-order Yee algorithm $`U_4^{Yee}`$ does not require extra storage to hold field values at intermediate times.

### IV.2 Unconditionally stable algorithms

Guided by previous work on Schrödinger and diffusion problems DeRaedt87 ; DeRaedt94 ; Michielsen98 , we split $`H`$ into two parts

$`H_1`$ $`=`$ $`\delta ^1{\displaystyle \stackrel{}{}_{i=1}^{n1}}\left(𝐞_i^{}𝐞_{i+1}^T𝐞_{i+1}^{}𝐞_i^T\right)=\left(\begin{array}{cccccc}0& \delta ^1& 0& 0& 0& \mathrm{}\\ \delta ^1& 0& 0& 0& 0& \mathrm{}\\ 0& 0& 0& \delta ^1& 0& \mathrm{}\\ 0& 0& \delta ^1& 0& 0& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}\\ 0& 0& 0& 0& \delta ^1& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right),`$ (49)
$`H_2`$ $`=`$ $`\delta ^1{\displaystyle \stackrel{}{}_{i=1}^{n2}}\left(𝐞_{i+1}^{}𝐞_{i+2}^T𝐞_{i+2}^{}𝐞_{i+1}^T\right)=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& \mathrm{}\\ 0& 0& \delta ^1& 0& 0& \mathrm{}\\ 0& \delta ^1& 0& 0& 0& \mathrm{}\\ 0& 0& 0& 0& \delta ^1& \mathrm{}\\ 0& 0& 0& \delta ^1& 0& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right).`$ (57)
such that $`H=H_1+H_2`$. In other words we divide the lattice into odd and even numbered cells. Clearly both $`H_1`$ and $`H_2`$ are skew-symmetric block-diagonal matrices, containing one $`1\times 1`$ matrix and $`(n1)/2`$ real, $`2\times 2`$ skew-symmetric matrices. According to the general theory given above, the first-order algorithm is given by

$`\stackrel{~}{U}_1(\tau )`$ $`=`$ $`e^{\tau H_1}e^{\tau H_2}=\left\{{\displaystyle \stackrel{}{}_{i=1}^{n1}}\mathrm{exp}\left[\tau \delta ^1\left(𝐞_i^{}𝐞_{i+1}^T𝐞_{i+1}^{}𝐞_i^T\right)\right]\right\}\left\{{\displaystyle \stackrel{}{}_{i=1}^{n2}}\mathrm{exp}\left[\tau \delta ^1\left(𝐞_{i+1}^{}𝐞_{i+2}^T𝐞_{i+2}^{}𝐞_{i+1}^T\right)\right]\right\}.`$ (58)
To derive Eq.(58) we used the block-diagonal structure of $`H_1`$ and $`H_2`$ (see Eqs.(49) and (57)) and obtained an exact expression for $`\stackrel{~}{U}_1(\tau )`$ in terms of an ordered product of matrix exponentials: the order of the matrix exponentials between each pair of curly brackets is irrelevant as these matrices commute with each other. Each of these matrix exponentials only operates on a pair of elements of $`𝚿(t)`$ and leaves other elements intact. The indices of each of these pairs are given by the subscripts of $`𝐞`$ and $`𝐞^T`$. From Eq.(58) it is clear what a program should do: Make loops over $`i`$ with stride 2. For each $`i`$ pick a pair of elements from $`\mathrm{\Psi }(t)`$ according to the subscripts of $`𝐞`$ and $`𝐞^T`$, compute (or fetch from memory) the elements of the plane rotation (see Eq. (67)), perform the plane rotation, i.e. multiply the $`2\times 2`$ matrices and the vectors of length two, and overwrite the same two elements. As the matrix exponential of a block-diagonal matrix is equal to the block-diagonal matrix of the matrix exponentials of the individual blocks, the numerical calculation of $`e^{\tau H_1}`$ (or $`e^{\tau H_2}`$) reduces to the calculation of $`(n1)/2`$ matrix exponentials of $`2\times 2`$ matrices. The matrix exponential of a typical $`2\times 2`$ matrix appearing in $`e^{\tau H_1}`$ or $`e^{\tau H_2}`$ is given by

$`\mathrm{exp}\left[\alpha \left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\right]\left(\begin{array}{c}\mathrm{\Psi }(i,t)\\ \mathrm{\Psi }(j,t)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{cos}\alpha & \mathrm{sin}\alpha \\ \mathrm{sin}\alpha & \mathrm{cos}\alpha \end{array}\right)\left(\begin{array}{c}\mathrm{\Psi }(i,t)\\ \mathrm{\Psi }(j,t)\end{array}\right).`$ (67)
Using the algorithm to compute Eq.(58), it is easy to construct the unconditionally stable, higher-order algorithms $`\stackrel{~}{U}_2(\tau )`$ and $`\stackrel{~}{U}_4(\tau )`$, see Eq.(13) and Eq.(15).

Obviously, the decomposition into $`H_1`$ Eq.(49) and $`H_2`$ Eq.(57) yields the most simple real-space algorithm. It is not difficult to imagine that a better but slightly more complicated algorithm can be constructed by using blocks of $`3\times 3`$ instead of $`2\times 2`$ matrices. Thus we are lead to consider the decomposition

$`H_3`$ $`=`$ $`\delta ^1{\displaystyle \stackrel{\prime \prime }{}_{i=1}^{n2}}\left(𝐞_i^{}𝐞_{i+1}^T+𝐞_{i+1}^{}𝐞_{i+2}^T𝐞_{i+1}^{}𝐞_i^T𝐞_{i+2}^{}𝐞_{i+1}^T\right)=\left(\begin{array}{cccccc}0& \delta ^1& 0& 0& 0& \mathrm{}\\ \delta ^1& 0& \delta ^1& 0& 0& \mathrm{}\\ 0& \delta ^1& 0& 0& 0& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}\\ 0& 0& 0& 0& \delta ^1& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right),`$ (75)
$`H_4`$ $`=`$ $`\delta ^1{\displaystyle \stackrel{\prime \prime }{}_{i=1}^{n4}}\left(𝐞_{i+2}^{}𝐞_{i+3}^T+𝐞_{i+3}^{}𝐞_{i+4}^T𝐞_{i+3}^{}𝐞_{i+2}^T𝐞_{i+4}^{}𝐞_{i+3}^T\right)=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}\\ 0& 0& 0& \delta ^1& 0& \mathrm{}\\ 0& 0& \delta ^1& 0& \delta ^1& \mathrm{}\\ 0& 0& 0& \delta ^1& 0& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right).`$ (83)
where the double prime indicated that the stride of the index $`i`$ is three. Obviously both $`H_3`$ and $`H_4`$ are skew-symmetric block-diagonal matrices, build from the $`3\times 3`$ skew-symmetric matrix

$`B`$ $`=`$ $`\left(\begin{array}{ccc}0& \delta ^1& 0\\ \delta ^1& 0& \delta ^1\\ 0& \delta ^1& 0\end{array}\right).`$ (87)
As $`B^3=2B`$ we have

$`e^{\tau B}`$ $`=`$ $`1+sB+cB^2=\left(\begin{array}{ccc}1c& s& c\\ s& 12c& s\\ c& s& 1c\end{array}\right),`$ (91)
where $`s=\mathrm{sin}(\sqrt{2}\tau )`$ and $`c=\mathrm{sin}^2(\tau /\sqrt{2})`$. In practice, using the $`3\times 3`$ instead of the $`2\times 2`$ decomposition is marginally more difficult. We will denote the corresponding second and fourth-order algorithm by $`U_2^{3\times 3}`$ and $`U_4^{3\times 3}`$ respectively.

### IV.3 ADI algorithm

For the tri-diagonal matrix (31), the ADI algorithm reduces to the Cranck-Nicholson method NumericalRecipes . The tri-diagonal structure of the matrix $`H`$ permits the calculation of $`(I\tau H/2)^1𝚿`$ in $`𝒪(n)`$ operations by standard linear algebra methods NumericalRecipes .

### IV.4 One-step algorithm

The one-step algorithm is based on the recursion (see Eqs.(21) and (22))

$$𝚿_{k+1}=\frac{2H}{H_1}𝚿_k+𝚿_{k1}.$$
(92)
Thus, the explicit form Eq.(31) is all we need to implement the matrix-vector operation (i.e. $`𝚿^{}H𝚿`$) that enters Eq.(92).

The coefficients $`J_k(z)`$ and $`S_k(t)`$ (see Eq.(25)) should be calculated to high precision. Using the recursion relation of the Bessel functions, all $`K`$ coefficients can be obtained with $`𝒪(K)`$ arithmetic operations NumericalRecipes . The numbers $`S_k(t)`$ can be calculated in $`𝒪(K\mathrm{log}K)`$ by standard Fast Fourier transformation techniques. Clearly both computations are a neglible fraction of the total computational cost for solving the Maxwell equations.

Performing one time step amounts to repeatedly using recursion (22) to obtain $`\widehat{T}_k(B)𝚿(0)`$ for $`k=2,\mathrm{},K`$, multiply the elements of this vector by $`J_k(z)`$ (or $`S_k(z)`$) and add all contributions. This procedure requires storage for two vectors of the same length as $`𝚿(0)`$ and some code to multiply such a vector by the sparse matrix $`H`$. The result of performing one time step yields the solution at time $`t`$, hence the name one-step algorithm. In contrast to what Eqs. (21) and (22) might suggest, the algorithm does not require the use of complex arithmetic.

In the sequel, the caret $`\widehat{}`$ on top of a symbol indicates that the results have been obtained by means of the one-step algorithm.

## V Numerical experiments

Except for the conventional Yee algorithm, all algorithms discussed in this paper operate on the vector of fields defined at the same time $`t`$. We use the one-step algorithm (with a time step $`\tau /2`$) to compute $`E_z(\tau /2)`$ and $`H_y(\tau /2)`$. Then we use $`E_z(0)`$ and $`H_y(\tau /2)`$ as the initial values for the Yee algorithm. In the presence of a current source, there are some ambiguities with this procedure as it is not obvious how to treat the source term in Eq.(9). In order to permit a comparison of the final result of the conventional Yee algorithm with those of the other methods, we use the one-step algorithm once more to shift the time of the $`H_y`$ field by $`\tau /2`$. This procedure to prepare the initial and to analyse the final state of the Yee algorithm does in fact make the results of the Yee algorithm look a little more accurate than they would be if the exact data of the $`\tau /2`$-shifted fields are not available. Thus, the results on the errors of the conventional Yee algorithm presented in this paper give a too optimistic view on the accuracy of this algorithm but we nevertheless adopt the above procedure to make a quantitative comparison between the various algorithms.

We define the error of the solution $`\stackrel{~}{𝚿}(t)`$ for the wave form by $`\stackrel{~}{𝚿}(t)\widehat{𝚿}(t)/\widehat{𝚿}(t)`$ where $`\widehat{𝚿}(t)`$ is the vector of EM fields obtained by the one-step algorithm. Thereby we have already assumed that the one-step algorithm yields the exact (within numerical precision) results but this has to be demonstrated of course. A comparison of the results of an unconditionally stable algorithm, e.g. $`\stackrel{~}{U}_4`$ with those of the one-step algorithm is sufficient to show that within rounding errors the latter yields the exact answer. Using the triangle inequality

$$𝚿(t)\widehat{𝚿}(t)𝚿(t)\stackrel{~}{𝚿}(t)+\stackrel{~}{𝚿}(t)\widehat{𝚿}(t),$$
(93)
and the rigorous bound

$$𝚿(t)\stackrel{~}{𝚿}(t)c_4\tau ^4t\left(𝚿(0)+_0^t𝐉(u)𝑑u\right),$$
(94)
we can be confident that the one-step algorithm yields the numerically exact answer if i) Eq.(94) is not violated and ii) if $`\stackrel{~}{𝚿}(t)\widehat{𝚿}(t)`$ vanishes like $`\tau ^4`$.

In Fig.3 we show a typical result of a one-step calculation on a grid of $`n=5001`$ sites with $`\delta =0.1`$ (corresponding to a physical length of $`250.05`$), and a current source placed at $`i=2500`$ to eliminate possible artifacts of the boundaries. The frequency of the source is set to one ($`f_s=1`$) and the number of periods the source radiates is set to six (i.e. $`T=6`$). In Table 1 (Fig.4) we present results for the errors, as obtained by repeating the simulation shown in Fig.3 using eight different FDTD methods. In Tables 2 and 3 (Figs.5 and 6 respectively) we shown similar results but instead of using a current source, a random wave form (Table 2) and Gaussian wave packet (Table 3) was taken as the initial condition.

From the data in Tables 1, 2 and 3 we conclude that the error of algorithm $`\stackrel{~}{U}_4`$ vanishes like $`\tau ^4`$, demonstrating that the one-step algorithm yields the numerically exact result (see Eqs.93 and 94). The results presented in Tables 2 and 3 have been obtained by using a vector of initial values that is normalized to one, i.e. $`𝚿(0)=1`$. As $`\widehat{𝚿}(t)=1`$ for $`0t100`$ to at least 10 digits, $`\stackrel{~}{𝚿}(t)\widehat{𝚿}(t)/\widehat{𝚿}(t)=\stackrel{~}{𝚿}(t)\widehat{𝚿}(t)`$ for all entries in Tables 2 and 3. The high precision of the one-step algorithm also allows us to use it for genuine time stepping with arbitrarily large time steps, this in spite of the fact that strictly speaking, the one-step algorithm is not unconditionally stable.

The data in Tables 1, 2 and 3 suggests that there does not seem to be a significant difference between the conventional Yee algorithm and its variant $`U_2^{Yee}`$ but in fact there is. The time evolution matrix corresponding to the Yee and the $`U_2^{Yee}`$ algorithm is not orthogonal. Therefore the energy of the electromagnetic field ($`𝚿(t)^2`$) is not conserved. Furthermore, in the conventional Yee algorithm, the $`E`$ and $`H`$-fields are time-shifted by $`\tau /2`$ with respect to each other. These artifacts of are less prominent if we use $`U_2^{Yee}`$ instead of the Yee algorithm. In Fig.7 we show results of the time evolution of the total energy of the EM field, for a system of n=97 sites, a mesh size $`\delta =0.1`$ and a time step of $`\tau =0.01`$. A normalized ($`𝚿(0)^2=1`$) random initial condition was used. Furthermore, for this type of application, it makes no sense to invoke the procedure described at the beginning of this section to time-shift one of the EM-fields by $`\tau /2`$: as this operation has to be performed at each time step and is computationally expensive (because it is numerically exact), we could as well use the same numerically exact procedure for time stepping itself.

For the Yee algorithm, the fluctuations of the energy are a factor of ten larger than in the case of the $`U_2^{Yee}`$ algorithm. As expected on theoretical grounds, the $`\stackrel{~}{U}_4`$ algorithm (dotted, horizontal line) exactly conserves the energy. The fact that $`U_2^{Yee}`$ conserves EM-field energy much better than the Yee algorithm also has a considerable impact on the quality of the eigenmode distribution. The latter is obtained by Fourier transformation of $`𝚿^T(t)𝚿(0)`$ (see ref.Kole01  for more details). In Fig.8 we show the low-frequency part of the eigenmode distribution of the same system as the one of Fig.7. It is obvious that there is a significant improvement in the quality of the spectrum if we use $`U_2^{Yee}`$ instead of the Yee algorithm but for this application $`\stackrel{~}{U}_4`$ performs much better than the Yee-type algorithms.

Table 1 suggests that $`\stackrel{~}{U}_2`$ is the least efficient of the five FDTD methods: It uses more arithmetic operations than the Yee algorithm and yields errors that are larger than those of the Yee algorithm. However, this conclusion is biased by the choice of the model problem and does not generalize. If the initial EM field distribution is random then, for sufficiently small $`\tau `$, algorithm $`\stackrel{~}{U}_2`$ is more accurate than the two second-order accurate Yee algorithms, as is clear from the data in Table 2 explains . Also in this case, for the largest $`\tau `$ in Table 2, the Yee and $`U_2^{Yee}`$ algorithm are operating at the point of instability, signaled by the fact that the norm of $`𝚿(t)`$ grows rapidly, resulting in errors that are very large. From Tables 1 and 2 one might conclude that the decomposition that generates Yee-type algorithms yields the least accurate approximations to the time evolution operator, although the difference is not really significant, but, as Table 3 shows, this conclusion would be wrong. If the initial state is a Gaussian wave packet that is fairly broad, the Yee-type algorithms are much more accurate than the unconditionally stable algorithms employed in this paper. From the data in Tables 1, 2 and 3 we conclude that there is no good reason to use the ADI algorithm (even disregarding the fact that it is slower than the other second-order methods). In general the $`\stackrel{~}{U}_2^{3\times 3}`$ ($`\stackrel{~}{U}_4^{3\times 3}`$) algorithm performs a little better than $`\stackrel{~}{U}_2`$ ($`\stackrel{~}{U}_4`$) but the gain is marginal. In contrast to the numerical data presented in Ref.Horvath , for all algorithms the data of Tables 1,2 and 3 clearly agree with the theoretically expected behavior of the error as a function of $`\tau `$ if $`\tau `$ is small enough exception .

Usually if a current source is present we have $`𝚿(0)=\mathrm{𝟎}`$. Then the one-step algorithm requires $`K^{}`$ (sparse) matrix-vector operations ($`𝚿^{}M𝚿`$) to compute $`𝚿(t)`$. For a 1D system the standard Yee, $`U_2^{Yee}`$ and $`U_4^{Yee}`$ ,$`\stackrel{~}{U}_2`$ , and $`\stackrel{~}{U}_4`$ algorithms perform (in worst case, without additional optimization), respectively, 1, 3/2, 6, 3/2, and 6 $`M𝚿`$-operations per time step. The one-step algorithm carries out $`K^{}=2103`$ matrix-vector operations $`𝚿^{}M𝚿`$ to complete this simulation. This implies that for all $`\tau <t/K^{}`$, the FDTD algorithms will perform more $`𝚿^{}M𝚿`$ operations than the one-step algorithm. For the data presented in this paper, this is the case if $`\tau <0.05`$ for the Yee algorithm and is always the case for $`\stackrel{~}{U}_4`$ because the latter uses a factor of 6 more $`𝚿^{}M𝚿`$ operations than the Yee algorithm.

## VI Conclusion

The answer to the question which of the algorithms is the most efficient crucially one depends on the accuracy that one finds acceptable. Taking the data of Table 1 as an example we see that if one is satisfied with an error of more than 2%, one could use the Yee algorithm. With $`\tau =0.05`$ it needs 2000 time steps to find the solution at $`t=100`$, close to the $`K^{}=2103`$. Nevertheless we recommend to use the one-step algorithm because then the time-integration error is neglegible. The Yee algorithm is no competition for $`\stackrel{~}{U}_4`$ if one requires an error of less than 1% but then $`\stackrel{~}{U}_4`$ is not nearly as efficient (by a factor of about 6) as the one-step algorithm. Increasing the dimensionality of the problem favors the one-step algorithm hdrcheb0 ; hdrcheb1 . These conclusions seem to be quite general and are in concert with numerical experiments on 1D, 2D and 3D systems hdrcheb1 . A simple theoretical analysis of the $`\tau `$ dependence of the error shows that the one-step algorithm is more efficient than any other FDTD method if we are interested in the EM fields at a particular (large) time only hdrcheb0 ; hdrcheb1 . This may open possibilities to solve problems in computational electrodynamics that are currently intractable. The Yee-like algorithms do not conserve the energy of the EM fields and therefore they are less suited for the calculation of the eigenvalue distributions (density of states), a problem for which the $`\stackrel{~}{U}_4`$ algorithm may be the most efficient of all the algorithms covered in the paper.

The main limitation of the one-step algorithm lies in its mathematical justification. The Chebyshev approach requires that $`H`$ is diagonalizable and that its eigenvalues are real or pure imaginary. The effect of relaxing these conditions on the applicability of the Chebyshev approach is left for future research.

In this paper we have focused entirely on the accuracy of the time integration algorithms, using the most simple discretization of the spatial derivatives. In practice it is straightforward, though technically non-trivial, to treat more sophisticated discretization schemes Taflove ; Kole02  by the methodology reviewed is this paper.

###### Acknowledgements.

H.D.R. and K.M. are grateful to T. Iitaka for drawing our attention to the potential of the Chebyshev method and for illuminating discussions. This work is partially supported by the Dutch ‘Stichting Nationale Computer Faciliteiten’ (NCF), and the EC IST project CANVAD.


text:
# On the definition and observability of the neutrino charge radius.

## 1 Introduction

The neutrino electromagnetic form-factor and the neutrino charge radius (NCR) have constituted an important theoretical puzzle for over than three decades. Since the dawn of the Standard Model (SM) it had been pointed out that one-loop radiative corrections to the scattering amplitude of a charged particle (such as the electron) with a neutrino would generate an effective photon-neutrino vertex $`\mathrm{\Gamma }_{A\nu \overline{\nu }}^\mu `$, which, in turn, would give rise to a non-zero NCR . Traditionally (and, of course, rather heuristically) the NCR has been interpreted as a measure of the “size” of the neutrino $`\nu _i`$ when probed electromagnetically, owing to its classical definition (in the static limit) as the second moment of the spatial neutrino charge density $`\rho _\nu (𝐫)`$, i.e. $`<r_\nu ^2>=e^1𝑑𝐫r^2\rho _\nu (𝐫)`$. Even though in the SM the one-loop computation of the entire $`S`$-matrix element describing the electron-neutrino scattering is conceptually straighforward, the identification of a sub-amplitude, which would serve as the effective $`\mathrm{\Gamma }_{A\nu \overline{\nu }}^\mu `$ has been faced with serious complications, associated with the simultaneous reconciliation of crucial requirements such as gauge-invariance, finiteness, and target-independence (for various references see ). The crux of the problem is that since in non-Abelian gauge theories individual off-shell Green’s functions are in general unphysical, the definition of quantities familiar from scalar theories or QED, such as effective charges and form-factors, is in general problematic. In the case of the NCR, various attempts to define its SM value through the one-loop $`\gamma \nu \nu `$ vertex calculated in the renormalizable ($`R_\xi `$) gauges reveal that the corresponding electromagnetic form-factor depends explicitly on the gauge-fixing parameter $`\xi `$ in a prohibiting way. In particular, even though in the static limit of zero momentum transfer, $`q^20`$, the form-factor becomes independent of $`\xi `$, its first derivative with respect to $`q^2`$, which corresponds to the definition of the NCR, continues to depend on it. Similar (and some times worse) problems occur in the context of other gauges (e.g. unitary gauge). These complications have obscured the entire concept of an NCR, and have casted serious doubts on whether it can be regarded as a genuine physical observable.

In the last two years significant progress has been accomplished in our understanding of the NCR. To begin with, the field-theoretcial difficulties mentioned above have been conclusively settled in , by resorting to the well-defined electroweak gauge-invariant separation of physical amplitudes into effective self-energy, vertex and box sub-amplitudes, implemented by the pinch technique formalism . In addition, in a very recent work , the observable nature of the NCR was established. In particular it was shown that the probe-independent charge radius of the neutrino may be extracted from experiment, at least in principle. This was accomplished by expressing a set of experimental $`\nu _\mu e`$ cross-sections in terms of the NCR and two additional gauge- and renormalization-group-invariant quantities, corresponding to the electroweak effective charge and mixing angle. In what follows we will present a brief summary of these new developments.

## 2 On the PT definition of the NCR

The PT is a diagrammatic method which exploits the underlying symmetries encoded in a physical amplitude such as an $`S`$-matrix element, in order to construct effective Green’s functions with special properties. The aforementioned symmetries, even though they are always present, they are usually concealed by the gauge-fixing procedure. The PT makes them manifest by means of a fixed algorithm, which does not depend on the gauge-fixing scheme one uses in order to quantize the theory, i.e. regardless of the set of Feynman rules used when writing down the $`S`$-matrix element. In particular, the PT uses the elementary Ward identities triggered by the longitudinal momenta appearing inside Feynman diagrams in order to enforce massive cancellations. The realization of these cancellations mixes non-trivially contributions stemming from diagrams of different kinematic nature (propagators, vertices, boxes). Thus, a given physical amplitude is reorganized into sub-amplitudes, which have the same kinematic properties as conventional $`n`$-point functions and, in addition, are endowed with desirable physical properties. Most importantly, they are independent of the gauge-fixing parameter, and satisfy naive (ghost-free, QED-like) Ward identities instead of the usual Slavnov-Taylor identities.

For the specific problem of the NCR, by resorting to the PT one can construct a genuine, target-independent form-factor endowed with the crucial properties of gauge-independence, gauge-invariance and finiteness, for arbitrary values of the momentum transfer. The essential ingredient for accomplishing this was the realization that, when the target fermions involved are considered to be massless, all gauge-dependent contributions stemming from box and vertex-like Feynman diagrams are effectively propagator-like, i.e. have no dependence on the kinematic properties or quantum numbers of the initial and final states, as boxes and vertices normally do. Thus, all gauge dependence stemming from vertices and boxes cancels precisely against the gauge-dependence stemming from the conventional one-loop self-energies. As a result, the remaining gauge-independent structures retain their initial kinematic identity; in particular one can speak in terms of gauge-independent effective boxes, vertices, and self-energies. Thus, once the gauge-dependent pieces have been extracted from the box, the remaining gauge-independent “pure” box should not be considered as a part of the resulting form-factor, which should be entirely determined from the “pure” gauge-independent set of vertex graphs.

The aforementioned propagator-like pieces are extracted by tracking down the action of the longitudinal momenta appearing in the $`S`$-matrix element of $`F\overline{F}\nu \overline{\nu }`$ . Longitudinal momenta originate from the tree-level gauge-boson propagators and tri-linear gauge-boson vertices appearing inside loops. In particular, in the $`R_\xi `$-scheme the gauge-boson propagators have the general form

$$\mathrm{\Delta }_V^{\mu \nu }(k)=\left[g^{\mu \nu }\frac{(1\xi _V)k^\mu k^\nu }{k^2\xi _VM_V^2}\right](k^2M_V^2)^1$$
(1)
where $`V=W,N`$ with $`N=Z,\gamma `$ and $`M_\gamma ^2=0`$; $`k`$ denotes the virtual four-momentum circulating in the loop. Clearly, in the case of $`\mathrm{\Delta }_V^{\mu \nu }(k)`$ the longitudinal momenta are those proportional to $`(1\xi _V)`$. The longitudinal terms arising from the tri-linear vertex may be identified by splitting $`\mathrm{\Gamma }_{\alpha \mu \nu }(q,k,kq)`$ appearing inside the one-loop diagrams into two parts ($`q`$ denotes the physical four-momentum entering into the vertex):

$$\mathrm{\Gamma }_{\alpha \mu \nu }=\mathrm{\Gamma }_{\alpha \mu \nu }^F+\mathrm{\Gamma }_{\alpha \mu \nu }^P,$$
(2)
with

$`\mathrm{\Gamma }_{\alpha \mu \nu }^F`$ $`=`$ $`(2k+q)_\alpha g_{\mu \nu }+2q_\nu g_{\alpha \mu }2q_\mu g_{\alpha \nu },`$
$`\mathrm{\Gamma }_{\alpha \mu \nu }^P`$ $`=`$ $`(k+q)_\nu g_{\alpha \mu }k_\mu g_{\alpha \nu }.`$ (3)
The above decomposition assigns a special role to the $`q`$-leg, and allows $`\mathrm{\Gamma }_{\alpha \mu \nu }^F`$ to satisfy the Ward identity

$$q^\alpha \mathrm{\Gamma }_{\alpha \mu \nu }^F=(k+q)^2g_{\mu \nu }k^2g_{\mu \nu },$$
(4)
All aforementioned longitudinal momenta originating from $`\mathrm{\Delta }_V^{\mu \nu }(k)`$ and $`\mathrm{\Gamma }_{\alpha \mu \nu }^P`$ trigger the following WI when contracted with the appropriate $`\gamma `$ matrix appearing in the various elementary vertices. In the absence of fermion masses

$`\overline{)}kP_L`$ $`=`$ $`(\overline{)}k+\overline{)}p)P_LP_R\overline{)}p,`$ (5)
$`=`$ $`S_F^{}^1(\overline{)}k+\overline{)}p)P_LP_RS_F^1(\overline{)}p)`$
where $`P_{R(L)}=[1+()\gamma _5]/2`$ is the chirality projection operator and $`S_F`$ is the tree-level propagator of the fermion $`F`$; $`F^{}`$ is the isodoublet-partner of the external fermion $`F`$. The result of this contraction is that the term in Eq.(5) proportional to $`S_F^{}^1`$, i.e the inverse of the internal fermion propagator gives rise to a self-energy-like term, whose coupling to the external fermions is proportional to the vertex

$$\mathrm{\Gamma }_{WF\overline{F}}^\mu =i\left(\frac{g_w}{2}\right)\gamma _\mu P_L.$$
(6)
The effective vertex $`\mathrm{\Gamma }_{WF\overline{F}}^\mu `$ can be written as a linear combination of the two tree-level vertices $`\mathrm{\Gamma }_{\gamma F\overline{F}}^\mu `$ and $`\mathrm{\Gamma }_{ZF\overline{F}}^\mu `$ as follows:

$$\mathrm{\Gamma }_{WF\overline{F}}^\mu =\left(\frac{s_w}{2T_z^F}\right)\mathrm{\Gamma }_{\gamma F\overline{F}}^\mu \left(\frac{c_w}{2T_z^F}\right)\mathrm{\Gamma }_{ZF\overline{F}}^\mu $$
(7)
In the above formulas $`Q_F`$ is the electric charge of the fermion $`F`$, and $`T_z^F`$ its $`z`$-component of the weak isospin, with $`c_w=\sqrt{1s_w^2}=M_W/M_Z`$, $`e=g_ws_w`$. The identity of Eq.(7) allows to combine the propagator-like parts emerging from box-diagrams and vertex-diagrams after the application of the WI in Eq.(5) with the conventional self-energy graphs $`\mathrm{\Pi }_{\mu \nu }^{\gamma Z}`$ and $`\mathrm{\Pi }_{\mu \nu }^{ZZ}`$ , by judiciously multiplying the former by $`D_N(q)D_N^1(q)`$.

Finally, from the gauge-independent one-loop proper vertex $`\widehat{\mathrm{\Gamma }}_{A\nu _i\overline{\nu }_i}^\mu `$ constructed using this method one extracts the dimension-full electromagnetic form-factor $`\widehat{F}_{\nu _i}(q^2)`$ as $`\widehat{\mathrm{\Gamma }}_{A\nu _i\overline{\nu }_i}^\mu =ieq^2\widehat{F}_{\nu _i}(q^2)\gamma _\mu (1\gamma _5)`$ . The NCR, to be denoted by $`<r_{\nu _i}^2>`$, is then defined as $`<r_{\nu _i}^2>=6\widehat{F}_{\nu _i}(0)`$, and thus one obtains

$$<r_{\nu _i}^2>=\frac{G_F}{4\sqrt{2}\pi ^2}\left[32\mathrm{log}\left(\frac{m_i^2}{M_W^2}\right)\right],$$
(8)
where $`i=e,\mu ,\tau `$, $`m_i`$ denotes the mass of the charged iso-doublet partner of the neutrino under consideration, and $`G_F`$ is the Fermi constant.

## 3 Measuring the probe-independent NCR

After arriving at a physically menaingful definition for the NCR, the next crucial question is whether the NCR so defined constitutes a genuine physical observable. In the rest of this section we will briefly discuss the method proposed in  for the extraction of the NCR from experiment.

It is important to emphasize that measuring the entire process $`f^\pm \nu f^\pm \nu `$ does not constitute a measurement of the NCR, because by changing the target fermions $`f^\pm `$ one will generally change the answer, thus introducing a target-dependence into a quantity which (supossedly) constitutes an intrinsic property of the neutrino. Instead, what we want to measure is the target-independent Standard Model NCR only, stripped of any target dependent contributions. Specifically, as mentioned above, the PT rearrangement of the $`S`$-matrix makes possible the definition of distinct, physically meaningful sub-amplitudes, one of which, $`\widehat{\mathrm{\Gamma }}_{A\nu _i\overline{\nu _i}}^\mu `$, is finite and directly related to the NCR. However, the remaining sub-amplitudes, such as self-energy, vertex- and box-corrections, even though they do no enter into the definition of the NCR, still contribute numerically to the entire $`S`$-matrix; in fact, some of them combine to form additional physical observables of interest, most notably the effective (running) electroweak charge and mixing angle. Thus, in order to isolate the NCR, one must conceive of a combination of experiments and kinematical conditions, such that all contributions not related to the NCR will be eliminated.

Consider the elastic processes $`f(k_1)\nu (p_1)f(k_2)\nu (p_2)`$ and $`f(k_1)\overline{\nu }(p_1)f(k_2)\overline{\nu }(p_2)`$, where $`f`$ denotes an electrically charged fermion belonging to a different iso-doublet than the neutrino $`\nu `$, in order to eliminate the diagrams mediated by a charged $`W`$-boson. The Mandelstam variables are defined as $`s=(k_1+p_1)^2=(k_2+p_2)^2`$, $`t=q^2=(p_1p_2)^2=(k_1k_2)^2`$, $`u=(k_1p_2)^2=(k_2p_1)^2`$, and $`s+t+u=0`$. In what follows we will restrict ourselves to the limit $`t=q^20`$ of the above amplitudes, assuming that all external (on-shell) fermions are massless. As a result of this special kinematic situation we have the following relations: $`p_1^2=p_2^2=k_1^2=k_2^2=p_1p_2=k_1k_2=0`$ and $`p_1k_1=p_1k_2=p_2k_1=p_2k_2=s/2`$. In the center-of-mass system we have that $`t=2E_\nu E_\nu ^{}(1x)0`$, where $`E_\nu `$ and $`E_\nu ^{}`$ are the energies of the neutrino before and after the scattering, respectively, and $`x\mathrm{cos}\theta _{cm}`$, where $`\theta _{cm}`$ is the scattering angle. Clearly, the condition $`t=0`$ corresponds to the exactly forward amplitude, with $`\theta _{cm}=0`$, $`x=1`$.

At tree-level the amplitude $`f\nu f\nu `$ is mediated by an off-shell $`Z`$-boson, coupled to the fermions by means of the bare vertex $`\mathrm{\Gamma }_{Zf\overline{f}}^\mu =i(g_w/c_w)\gamma ^\mu [v_f+a_f\gamma _5]`$ with $`v_f=s_w^2Q_f\frac{1}{2}T_z^f`$ and $`a_f=\frac{1}{2}T_z^f`$.

At one-loop, the relevant contributions are determined through the PT rearrangement of the amplitude, giving rise to gauge-independent sub-amplitudes. In particular, the one-loop $`AZ`$ self-energy $`\widehat{\mathrm{\Sigma }}_{AZ}^{\mu \nu }(q^2)`$ obtained is transverse, for both the fermionic and the bosonic contributions, i.e. $`\widehat{\mathrm{\Sigma }}_{AZ}^{\mu \nu }(q^2)=(q^2g^{\mu \nu }q^\mu q^\nu )\widehat{\mathrm{\Pi }}_{AZ}(q^2)`$. Since the external currents are conserved, from the $`ZZ`$ self-energy $`\widehat{\mathrm{\Sigma }}_{ZZ}^{\mu \nu }(q^2)`$ we keep only the part proportional to $`g^{\mu \nu }`$, whose dimension-full cofactor will be denoted by $`\widehat{\mathrm{\Sigma }}_{ZZ}(q^2)`$. Furthermore, as is well-known, the one-loop vertex $`\widehat{\mathrm{\Gamma }}_{ZF\overline{F}}^\mu (q,p_1,p_2)`$, with $`F=f`$ or $`F=\nu `$, satisfies a QED-like Ward identity, relating it to the one-loop inverse fermion propagators $`\widehat{\mathrm{\Sigma }}_F`$, i.e $`q_\mu \widehat{\mathrm{\Gamma }}_{ZF\overline{F}}^\mu (q,p_1,p_2)=\widehat{\mathrm{\Sigma }}_F(p_1)\widehat{\mathrm{\Sigma }}_F(p_2)`$. It is then easy to show that, in the limit of $`q^20`$, $`\widehat{\mathrm{\Gamma }}_{ZF\overline{F}}^\mu q^2\gamma ^\mu (c_1+c_2\gamma _5)`$; since it is multiplied by a massive $`Z`$ boson propagator $`(q^2M_Z)^1`$, its contribution to the amplitude vanishes when $`q^20`$. This is to be contrasted with the $`\widehat{\mathrm{\Gamma }}_{A\nu _i\overline{\nu }_i}^\mu `$, which is accompanied by a $`(1/q^2)`$ photon-propagator, thus giving rise to a contact interaction between the target-fermion and the neutrino, described by the NCR.

We next proceed to eliminate the target-dependent box-contributions; to accomplish this we resort to the “neutrino–anti-neutrino” method. The basic observation is that the tree-level amplitudes $`_{\nu f}^{(0)}`$ as well as the part of the one-loop amplitude $`_{\nu f}^{(B)}`$ consisting of the propagator and vertex corrections (the “Born-improved” amplitude), are proportional to $`[\overline{u}_f(k_2)\gamma _\mu (v_f+a_f\gamma _5)u_f(k_1)][\overline{v}(p_1)\gamma _\mu P_Lv(p_2)]`$, and therefore transform differently than the boxes under the replacement $`\nu \overline{\nu }`$. In particular, the coupling of the $`Z`$ boson to a pair of on-shell anti-neutrinos may be written in terms of on-shell neutrinos provided that one changes the chirality projector from $`P_L`$ to $`P_R`$, and supplies a relative minus sign , i.e.

$$\overline{v}(p_1)\gamma _\mu P_Lv(p_2)=\overline{u}(p_2)\gamma _\mu P_Ru(p_1)$$
(9)
Thus, under the above transformation, $`_{\nu f}^{(0)}+_{\nu f}^{(B)}`$ reverse sign once, whereas the box contributions reverse sign twice. These distinct transformation properties allow for the isolation of the box contributions when judicious combinations of the forward differential cross-sections $`(d\sigma _{\nu f}/dx)_{x=1}`$ and $`(d\sigma _{\overline{\nu }f}/dx)_{x=1}`$ are formed. In particular, $`\sigma _{\nu f}^{(+)}(d\sigma _{\nu f}/dx)_{x=1}+(d\sigma _{\overline{\nu }f}/dx)_{x=1}`$ does not contain boxes, whereas the conjugate combination $`\sigma _{\nu f}^{()}(d\sigma _{\nu f}/dx)_{x=1}(d\sigma _{\overline{\nu }f}/dx)_{x=1}`$ isolates the contribution of the boxes.

Finally, a detailed analysis shows that in the kinematic limit we consider, the Bremsstrahlung contribution vanishes, due to a a completely destructive interference between the two relevant diagrams corresponding to the processes $`fA\nu (\overline{\nu })f\nu (\overline{\nu })`$ and $`f\nu (\overline{\nu })fA\nu (\overline{\nu })`$. The absence of such corrections is consistent with the fact that there are no infrared divergent contributions from the (vanishing) vertex $`\widehat{\mathrm{\Gamma }}_{ZF\overline{F}}^\mu `$, to be cancelled against.

$`\sigma _{\nu f}^{(+)}`$ receives contributions from the tree-level exchange of a $`Z`$-boson, the one-loop contributions from the ultraviolet divergent quantities $`\widehat{\mathrm{\Sigma }}_{ZZ}(0)`$ and $`\widehat{\mathrm{\Pi }}^{AZ}(0)`$, and the (finite) NCR, coming from the proper vertex $`\widehat{\mathrm{\Gamma }}_{A\nu _i\overline{\nu }_i}^\mu `$. The first three contributions are universal, i.e. common to all neutrino species, whereas that of the proper vertex $`\widehat{\mathrm{\Gamma }}_{A\nu _i\overline{\nu }_i}^\mu `$ is flavor-dependent.

To proceed, the renormalization of $`\widehat{\mathrm{\Sigma }}_{ZZ}(0)`$ and $`\widehat{\mathrm{\Pi }}_{AZ}(0)`$ must be carried out. It turns out that, by virtue of the Abelian-like Ward-identities enforced after the pinch technique rearrangement , the resulting expressions combine in such a way as to form manifestly renormalization-group invariant combinations . In particular, after carrying out the standard re-diagonalization, two such quantities may be constructed:

$`\overline{R}_Z(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _w}{c_w^2}}\left[q^2M_Z^2+\mathrm{}e\{\widehat{\mathrm{\Sigma }}_{ZZ}(q^2)\}\right]^1`$
$`\overline{s}_w^2(q^2)`$ $`=`$ $`s_w^2\left(1{\displaystyle \frac{c_w}{s_w}}\mathrm{}e\{\widehat{\mathrm{\Pi }}_{AZ}(q^2)\}\right).`$ (10)
where $`\alpha _w=g_w^2/4\pi `$, and $`\mathrm{}e\{\mathrm{}\}`$ denotes the real part. In addition to being renormalization-group invariant, both quantities defined in Eq.(10) are universal (process-independent); $`\overline{R}_Z(q^2)`$ corresponds to the $`Z`$-boson effective charge, while $`\overline{s}_w^2(q^2)`$ corresponds to an effective mixing angle. We emphasize that the renormalized $`\widehat{\mathrm{\Pi }}_{AZ}(0)`$ cannot form part of the NCR, because it fails to form a renormalization-group invariant quantity on its own. Instead, $`\widehat{\mathrm{\Pi }}_{AZ}(0)`$ must be combined with the appropriate tree-level contribution (which evidently does not enter into the definition of the NCR, since it is $`Z`$-mediated) in order to form the effective $`\overline{s}_w^2(q^2)`$ acting on the electron vertex, whereas the finite NCR will be determined from the proper neutrino vertex only.

After recasting $`\sigma _{\nu f}^{(+)}`$ in terms of manifestly renormalization-group invariant building blocks, one may fix $`\nu =\nu _\mu `$, and then consider three different choices for $`f`$: (i) right-handed electrons, $`e_R`$; (ii) left-handed electrons, $`e_L`$, and (iii) neutrinos, $`\nu _i`$ other than $`\nu _\mu `$, i.e. $`i=e,\tau `$. Thus, we arrive at the system

$`\sigma _{\nu _\mu \nu _i}^{(+)}`$ $`=`$ $`s\pi \overline{R}^2(0)`$
$`\sigma _{\nu _\mu e_R}^{(+)}`$ $`=`$ $`s\pi \overline{R}^2(0)\overline{s}_w^4(0)2\lambda s_w^2<r_{\nu _\mu }^2>`$
$`\sigma _{\nu _\mu e_L}^{(+)}`$ $`=`$ $`s\pi \overline{R}^2(0)\left({\displaystyle \frac{1}{2}}\overline{s}_w^2(0)\right)^2`$ (11)
$`+\lambda (12s_w^2)<r_{\nu _\mu }^2>`$
where $`\lambda (2\sqrt{2}/3)s\alpha G_F`$, $`\alpha =e^2/4\pi `$. $`\overline{R}^2(0)`$, $`\overline{s}_w^2(0)`$, and $`<r_{\nu _\mu }^2>`$ are treated as three unknown quantities, to be determined from the above equations. The corresponding solutions are given by $`\overline{s}_w^2(0)=s_w^2\pm \mathrm{\Omega }^{1/2}`$ and

$`<r_{\nu _\mu }^2>`$ $`=`$ $`\lambda ^1[(s_w^2{\displaystyle \frac{1}{4}}\pm \mathrm{\Omega }^{1/2})\sigma _{\nu _\mu \nu _i}^{(+)}`$ (12)
$`+\sigma _{\nu _\mu e_L}^{(+)}\sigma _{\nu _\mu e_R}^{(+)}]`$
where the discriminant $`\mathrm{\Omega }`$ is given by

$$\mathrm{\Omega }=(12s_w^2)\left(\frac{\sigma _{\nu _\mu e_R}^{(+)}}{\sigma _{\nu _\mu \nu _i}^{(+)}}\frac{1}{2}s_w^2\right)+2s_w^2\frac{\sigma _{\nu _\mu e_L}^{(+)}}{\sigma _{\nu _\mu \nu _i}^{(+)}}$$
(13)
and must satisfy $`\mathrm{\Omega }>0`$. The actual sign in front of $`\mathrm{\Omega }`$ may be chosen by requiring that it correctly accounts for the sign of the shift of $`\overline{s}_w^2(0)`$ with respect to $`s_w^2`$ predicted by the theory . To extract the experimental values of the quantities $`\overline{R}^2(0)`$, $`\overline{s}_w^2(0)`$, and $`<r_{\nu _\mu }^2>`$, one must substitute in Eq.(12) and Eq.(13) the experimentally measured values for the differential cross-sections $`\sigma _{\nu _\mu e_R}^{(+)}`$, $`\sigma _{\nu _\mu e_L}^{(+)}`$, and $`\sigma _{\nu _\mu \nu _i}^{(+)}`$. This means that to solve the system one would have to carry out three different pairs of experiments.

In summary, we have seen that the neutrino charge radius can be defined unambiguously by means of the pinch technique, and that its observability has been established.

Acknowledgments: This work has been supported by the Grant AEN-99/0692 of the Spanish CICYT. J.P. thanks the organizers of IMFP2002 for their hospitality at Jaca and the stimulating environment they have created.