\section{Introduction}
\label{sec:intro}

The increasing maturity of space-based gravitational wave detector 
concepts~\cite{NASA,ESA} has in part motivated considerable work in the 
last fifteen years on self-consistent calculations of extreme-mass-ratio 
inspirals (EMRIs).  Such a system consists of a small compact object of mass 
$\mu \simeq 1-10 M_{\odot}$ (e.g., neutron star or black hole) moving on a 
decaying orbit about, and ultimately into, a supermassive black hole of 
mass $M \sim 10^5 - 10^9 M_{\odot} \gg \mu$.  Irrespective of when a detector 
might launch, there is also simply considerable theoretical interest in the 
problem of motion of a point mass in a background geometry in general 
relativity, influenced by its own self-force~\cite{Poisson_2011}.

The extreme mass-ratio lends itself to use of black hole perturbation theory.  
In the limit of $\mu \rightarrow 0$ the small mass orbits on a geodesic of
the massive black hole background with constants of motion.  At next order 
the small mass draws up a small perturbation in the metric, which results 
in gravitational radiation fluxing to infinity and down the horizon of the 
massive black hole.  The metric perturbation (MP) also acts back on the small 
body locally (through a self-force), giving rise to dissipative effects that 
cause the orbit to decay and to small conservative corrections to the motion.  
The perturbation problem is singular in several respects~\cite{Pound_2010}, 
with a divergence in the MP at the particle location making 
the motion correction also formally divergent.  A general understanding of 
how to treat the self-force (i.e., regularize it) in an arbitrary spacetime 
was given by Mino, Sasaki, and Tanaka~\cite{MiSaTa} and Quinn and 
Wald~\cite{QuWa}.  Practical procedures for regularizing the self-force in 
numerical calculations followed (e.g.,~\cite{BO_2000}).

The physical retarded field $p^{\text{ret}}_{\mu\nu}$ can be
conveniently split into regular ($R$) and singular ($S$) parts,
as first introduced by Detweiler and Whiting~\cite{DW_2003}.
The advantage of this split is that while the singular contribution
to the MP $p^{S}_{\mu\nu}$ satisfies the inhomogeneous field equations,
it does not contribute at all to the self-force.
On the other hand, the regular contribution, $p^{R}_{\mu\nu}$, is
a smooth, homogeneous solution to the field equations and,
through a projected gradient, is entirely 
responsible for the self-force.  Indeed, when interpreted this way,
the regular field can be thought of as an external field which sources
the deviation from geodesic motion on the background metric 
$g_{\mu \nu}$.  The motion of the particle is then geodesic on the spacetime
$g_{\mu \nu} + p^{R}_{\mu \nu}$.  The singular part of the MP is 
calculated analytically in 
Lorenz gauge, and an expansion provides the regularization 
parameters~\cite{BO_2000}.  The singular part is then subtracted from the 
full retarded field mode by mode in a spherical harmonic expansion, allowing 
the difference ($p_{\mu \nu}^{R}$) to converge.  While in 
principle~\cite{BO_2001,Barack_2009} the full retarded field could be 
calculated in a variety of gauges, in practice most 
calculations~\cite{BS_2007,BS_2009,BS_2010} have also used Lorenz gauge 
to find $p^{\text{ret}}_{\mu\nu}$.

We are developing a set of techniques and assembling a computer code to 
calculate with high accuracy the first-order MPs from a small compact 
object in a generic orbit about a Schwarzschild black hole.  The need for 
high accuracy is related to a set of arguments that have been made for years 
that EMRIs should be calculated through second-order in perturbation 
theory~\cite{Poisson_2002,Rosenthal_2006,Pound_2009,Galley_2011,Poisson_2011}.
One particular argument centers on calculating the phase evolution of an 
EMRI and using it in interpreting data from a detector.  For a small mass ratio 
$\epsilon = \mu/M$ (and in the absence of transient 
resonances~\cite{FH_2010} which may occur during EMRI evolution on a Kerr 
black hole), we expect that as an EMRI evolves through a detector passband 
the gravitational waveform will accumulate (schematically) a phase of
\begin{equation}
\Phi = \kappa_1 \frac{1}{\epsilon} + \kappa_2 \epsilon^0 
+ \kappa_3 \epsilon^1 + \cdots ,
\end{equation}
where the $\kappa$'s are coefficients of order unity that depend upon, among 
other things, the lower and upper limits in frequency of the detector 
response.  The first term reflects the dissipative effects of the first-order 
self-force in spurring a decay of the orbit.  The second term would result 
from second order in perturbation theory.  For example, with 
$\epsilon = 10^{-6}$, an EMRI might be observed to accumulate a total phase 
of $\Phi \sim 10^6$.  For matched filtering purposes we might need to compute 
the phase to an accuracy $\delta\Phi \lesssim 0.1$, and thus a fractional 
error of $\lesssim 10^{-7}$.  However, the error in phase in using the 
first-order calculation alone is $\sim \mathcal{O}(1)$.  Hence, the need for 
a second approximation to take full advantage of a detector output.  
\emph{But there is a corollary to this argument.}  We cannot possibly hope 
to make use of a second-order calculation if we have not already computed 
the first-order self-force to a relative accuracy \emph{much} better than 
$\mathcal{O}(\epsilon)$.  The requirement might be at least several orders 
of magnitude better than $\mathcal{O}(\epsilon)$ to make a second order 
calculation worthwhile (say $10^{-9}$ to $10^{-8}$ in the example).  
Furthermore, given that computation of the self-force is a numerically 
subtractive procedure, the first-order pre-regularization field contributions 
likely need to be known even more accurately (perhaps $10^{-11}$ to 
$10^{-10}$).

An accurate calculation strongly suggests use of Fourier decomposition and 
frequency domain (FD) methods, to gain the benefit of integrating ordinary 
differential equations.  Ultimately we are interested in the time-dependent 
self-force and must transfer back to the time domain (TD).  For that step, a 
lynchpin of the effort has been use of the recently developed method of 
extended homogeneous solutions (EHS)~\cite{BOS}, which allows partial Fourier 
series sums for the perturbations to avoid the Gibbs phenomenon and to 
converge exponentially even at the location of the point mass and despite 
loss of differentiability there.  This conclusion has guided others as 
well~\cite{akcay_2011,warburton_2012}.

Like other recent calculations, we want to determine the first-order 
MP in Lorenz gauge.  However, our approach is indirect.  
In an earlier paper~\cite{HE_2010}, we calculated 
radiative modes ($\ell\ge 2$) of the MP in Regge-Wheeler (RW) gauge by 
applying EHS to solutions of the master equations of the 
Regge-Wheeler-Zerilli (RWZ) formalism and then determining the metric from 
the master functions.  The new aspect of that work was being able to use FD 
techniques and nevertheless determine the metric amplitudes in the TD with 
accuracy right up to the location of the particle $r=r_p(t)$, an essential 
requirement for computing the self-force accurately.  EHS works by 
recognizing that the solution of a master equation with a moving singular 
source (in the RWZ case the source has both a delta function term and a 
derivative of delta function term) is a weak solution of the form 
$\Psi(t,r) = \Psi^{+}(t,r)\, \th \left[ r-r_p(t)\right] 
+ \Psi^{-}(t,r)\, \th \left[ r_p(t)-r \right]$, 
where $\Psi^{+}(t,r)$ and $\Psi^{-}(t,r)$ are differentiable solutions to the 
source-free master equation.  $\Psi^{+}$ and $\Psi^{-}$ are in turn obtained
as Fourier sums of properly-normalized Fourier-harmonic modes that solve 
the source-free master equation in the FD.  Since the functions in the 
separate Fourier sums are smooth everywhere, the lack of differentiability 
of $\Psi$ is captured entirely by the $\theta$ functions.  We can in turn then 
calculate the MPs from $\Psi$.  The result, however, is in RW gauge.  This 
paper, and a subsequent one, address the calculation of the infinitesimal 
gauge generator that transforms $p^{\text{RW}}_{\mu\nu}$ in RW gauge to its 
counterpart $p^{\text{L}}_{\mu\nu}$ in Lorenz gauge.  

This paper is restricted to finding the odd-parity part of the gauge 
generator, which for each spherical harmonic mode has an amplitude that is 
a solution to a single inhomogeneous wave equation.  While the wave equation 
is simple to express, what is more challenging is to find a way to 
generalize the underlying idea behind EHS to equations with 
non-compact source terms.  A substantial part of this paper is devoted to 
laying out two new analytic/numerical methods (\emph{method of partial 
annihilators} (PA) and \emph{method of extended particular solutions} (EPS)) 
we have developed for solving differential equations of this type.  
We present results from each of these methods, including a comparison of the 
two, showing that they agree to a high accuracy.
These techniques will play a central role in a subsequent paper where we 
present the more involved procedure, based on analysis by Sago, Nakano, 
and Sasaki~\cite{SNS_2003}, for determining the even-parity parts of 
the RW-to-Lorenz gauge generator. The Sago, Nakano, Sasaki approach takes
the gauge generator equations for even-parity, which are most naturally
expressed as a set of three, coupled equations for the vector spherical 
harmonic amplitudes, and transforms them to an altered set of equations in 
terms of different amplitudes.  The resulting equations are a hierarchical 
set of second-order equations, in which solutions to preceding steps in the 
hierarchy form source terms for subsequent steps. The result is a system 
of equations that, while containing more steps (amplitudes), lends itself to 
the immediate application of the techniques developed here.

One might ask, why two new methods?  In part, partial annihilators is the 
easier of the two methods to implement but requires that the partial 
annihilator operator be found.  In some applications that may be difficult.  
In contrast, EPS is straightforward if somewhat more involved in terms of the 
number of steps.  Ultimately, the main advantage of finding two methods is
that they provide a powerful check on each other and confirmation that the 
solution has been obtained.  We demonstrate this comparison in 
Sec.~\ref{sec:results}.

The discussion above begs the question, why the need to transform to Lorenz 
gauge?  In part we know that the MP amplitudes in Lorenz gauge are $C^0$ in 
the TD at the particle location.  As we discussed in Ref.~\cite{HE_2010} the 
MP amplitudes in RW gauge are one or two orders of continuity worse behaved 
(i.e., some amplitudes are $C^{-1}$, or discontinuous, and some have a 
delta function term at $r=r_p(t)$).  Fig.~\ref{fig:hthrRW} demonstrates this 
problem graphically for the $\ell = 2, m = 1$ (odd parity) RW amplitudes 
$h_{t}^{2,1}$ and $h_{r}^{2,1}$.  The insets show discontinuities in the MP 
amplitudes at the particle's location.  The plots also show the $\propto r$ 
growth in the wave pulse amplitude as $r \rightarrow \infty$, reflecting 
the fact that RW gauge is not asymptotically flat~\cite{TC_1967}.  
Transformation to Lorenz gauge not only improves the behavior of the modes 
at $r=r_p(t)$, it also removes the non-asymptotically-flat behavior seen in
the radiative modes of RW gauge.

\begin{figure}
\vspace{-.5cm}
\includegraphics[scale=1]{h_t_and_f_h_r_RW}
\vspace{-.5cm}
\caption{
\label{fig:hthrRW}
The $\ell = 2$, $m = 1$ mode of the Regge-Wheeler gauge MP amplitudes 
$h_{t}^{\text{RW}}$ and $h_{r}^{\text{RW}}$.  The orbit is parameterized by 
eccentricity $e = 0.764124$ and semi-latus rectum $p = 8.75455$.  The plots 
show the real and imaginary parts of the MP amplitudes at time $t = 93.58$ 
(where $t = 0$ is at the periapsis).  Dotted vertical lines indicate limits 
of the source libration region.  Insets show the discontinuities at the 
location of the particle.  The lack of asymptotic flatness is evident as 
$r\rightarrow \infty$.  We plot $f h^{\text{RW}}_{r}$ so that the wave 
behavior near the horizon can be seen.}
\vspace{-.5cm}
\end{figure}

Lastly, we note that for non-radiative ($\ell = 0,1$) modes, the RWZ formalism
breaks down.  Various researchers have used different gauges to solve the 
Einstein equations for these modes.  Considering generic motion on a 
Schwarzschild 
background, Zerilli \cite{Zerilli_1970} solved for these modes analytically
in his own gauge which exhibits $C^{-1}$ behavior for certain modes.  For
circular orbits Detweiler and Poisson \cite{DP_2003} showed how to transform 
Zerilli's solutions to Lorenz gauge.  Barack and Sago \cite{BS_2010} solved 
the Lorenz gauge field equations directly for these modes (and higher 
radiative modes).  The non-radiative modes provide a crucial contribution 
to the conservative piece of the self-force.  In this paper, though, we are 
solely concerned with transforming our RW gauge ($\ell \ge 2$) MP amplitudes 
to Lorenz gauge.  However, we note that while RW gauge is not defined 
for $\ell < 2$, there is no such restriction on the gauge transformation.  
Indeed, the work presented here can be directly extended to handle the 
transformations of non-radiative mode solutions from other gauges to Lorenz.

Throughout this paper we use the sign conventions and notation of 
Misner, Thorne, and Wheeler \cite{MTW_1973} and use units in which 
$c = G = 1$.  The background geometry is a non-rotating black hole which is
described in terms of Schwarzschild coordinates 
$x^{\mu} = (t,r,\theta, \varphi )$.

\section{Formalism}
\label{sec:formalism}

Consider the motion of a small compact object of mass $\mu$ in orbit about 
a static black hole of mass $M$, with $\mu/M \ll 1$.  The small body gives 
rise to a perturbation $p_{\mu\nu}$ in the metric relative to the 
Schwarzschild background,
\begin{equation}
ds^2 = g_{\mu\nu} dx^{\mu} dx^{\nu} = 
-f(r) dt^2 + f(r)^{-1} dr^2 
+ r^2 \left( d\theta^2 + \sin^2\theta d\varphi^2 \right) ,
\end{equation}
where $f(r) = 1 - 2M/r$.  We are concerned here only with the first-order 
part $p^{(1)}_{\mu\nu}$ of the MP in an expansion in powers of $\mu/M$ and 
accordingly simply set $p_{\mu\nu} = p^{(1)}_{\mu\nu}$.  The full 
non-stationary metric is then 
$\text{g}_{\mu\nu} = g_{\mu\nu} + p_{\mu\nu}(t,r,\theta,\varphi)$.  It often 
proves convenient to work with the trace-reversed MP
\begin{equation}
\Bar{p}_{\mu\nu} = p_{\mu\nu} 
- \tfrac{1}{2} g_{\mu\nu} p_{\alpha\beta} g^{\alpha\beta} .
\end{equation}
The Einstein field equations can be expanded about the background geometry.
Once the stress-energy tensor $T^{(0)}_{\a\b}$ associated with the motion of 
the small compact object has been specified,
the first-order linear equations
\begin{equation}
\label{eq:linearEFE}
G^{(1)}_{\a\b}(\Bar{p}_{\mu\nu}) = 8 \pi T^{(0)}_{\a\b} ,
\end{equation}
can be solved to determine $\Bar{p}_{\mu\nu}$.  The small body is 
approximated as a point particle and the stress tensor is that of the 
particle moving on a geodesic of the Schwarzschild black hole (i.e., 
zeroth-order approximation).  

In this work, we are interested in bound eccentric 
motion.  The orbit is parameterized by a pair of constants, which can 
alternately be taken to be Darwin's \cite{Darwin} eccentricity $e$ and 
dimensionless semi-latus rectum $p$, or the bounding values of radial 
motion $r_{\text{min}}$ and $r_{\text{max}}$, or the specific energy 
$\cal{E}$ and angular momentum $\cal{L}$.  The first integrals of motion 
are integrated using Darwin's curve parameter $\chi$ and are converted to 
functions of the form $r=r_p(t)$, $\varphi = \varphi_p(t)$, and 
$\tau = \tau_p(t)$ (with $\theta = \pi/2$).  The orbit 
has two fundamental frequencies with
a rate $\Omega_r$ associated with radial libration and a mean rate of 
azimuthal advance $\Omega_{\varphi}$.  See Darwin 
\cite{Darwin} for the orbital integration and Cutler et al.~\cite{CKP_1994} 
for first application to black hole perturbations from bound motion. 

\subsection{Harmonic decomposition of the field equations}

The usual route to solving the field equations (\ref{eq:linearEFE})
on a Schwarzschild background is to
introduce tensor spherical harmonics and decompose the equations into 
individual angular harmonic modes.  This approach was followed by Regge and 
Wheeler \cite{RW_1957} and Zerilli \cite{Zerilli_1970} in solving for 
certain master functions that represent the odd- and even-parity parts, 
respectively, of the gravitational field.  The MP is then
derived from these master functions in Regge-Wheeler gauge.  
In Ref.~\cite{HE_2010},
 we used a variant of this approach, with slightly different versions 
of the two master functions, and combined it with a new analytic/numerical 
method for finding convergent solutions near the location of the point 
particle.  

In that paper we followed Martel and Poisson \cite{MP_2005} in using their 
definitions of the tensor spherical harmonics.  We recap those definitions 
here.  The unit two-sphere $\mathcal{S}^2$ is covered by coordinates 
$(\theta, \varphi)$.  Upper-case Latin indices $A$, $B$, etc.~denote these 
two angular coordinates and associated tensor components.  
The coordinates $(t,r)$ cover the submanifold $\mathcal{M}^2$.  For these 
coordinates and tensor components we use lower-case Latin indices $a$, $b$, 
etc.  The full Schwarzschild spacetime is 
$\mathcal{M} = \mathcal{M}^2 \times \mathcal{S}^2$.  The usual scalar 
spherical harmonic functions are $Y^{\ell m}(\theta,\varphi)$.  From 
these we define even-parity $Y_A^{\ell m} = D_A Y^{\ell m}$ and
odd-parity $X_A^{\ell m} = - {\varepsilon_A}^B \, D_B Y^{\ell m}$ vector 
harmonics, where $D_A$ is the covariant derivative on $\mathcal{S}^2$.  The 
metric on the unit sphere is $\Omega_{AB} = \text{diag}[1, \sin^2\theta]$ and 
the Levi-Civita tensor is $\varepsilon_{AB}$.  Both are compatible with 
$D_A$: $D_C \, \Omega_{AB} = D_C \, \varepsilon_{AB} = 0$.  Martel and 
Poisson define the two even-parity tensor spherical harmonics as 
$\Omega_{AB} Y^{\ell m}$ and
$Y_{AB}^{\ell m}=\left[D_A D_B +\tfrac{1}{2}\ell (\ell +1)\right] Y^{\ell m}$,
where the latter is trace-free and thus differs from that used by Regge and 
Wheeler ($D_{A} D_{B} Y^{\ell m}$)~\cite{RW_1957}.  They take the odd-parity 
tensor harmonic to be
$X_{AB}^{\ell m} = - \tfrac{1}{2} 
\left( {\varepsilon_A}^C D_B + {\varepsilon_B}^C D_A \right) D_C Y^{\ell m}$,
which differs by a minus sign from that of~\cite{RW_1957}.

Using these definitions the MP is split into even- and 
odd-parity parts.  Each parity is decomposed into sums over $\ell$ and $m$ 
of their respective harmonics.  For even-parity there are seven amplitudes 
that Martel and Poisson define, which can be related to those of 
\cite{RW_1957} and \cite{Zerilli_1970}.  They are $h_{tt} = f H_{0}$, 
$h_{tr} = H_{1}$, $h_{rr} = H_{2} / f$, $j_{t} = h_{0}$, $j_{r} = h_{1}$, 
$G_{\rm here} = G_{\text{RW}}$, and $K_{\rm here} = K_{\text{RW}} 
- \ell (\ell + 1) G / 2$.  (Note: here and on many occasions later in the 
text we suppress the spherical harmonic scripts $\ell$ and $m$ for brevity 
when no confusion should arise.)  Then the even-parity MP is
\begin{align}
\label{eq:MPEven}
p_{ab} &= \sum_{\ell, m} h_{ab}^{\ell m} Y^{\ell m}, 
& p_{aB} &= \sum_{\ell, m} j_a^{\ell m} Y_B^{\ell m}, 
& p_{AB} &= r^2 \sum_{\ell, m}\Big( K^{\ell m} \O_{AB} Y^{\ell m} 
+  G^{\ell m} Y^{\ell m}_{AB} \Big).
\end{align}
For odd-parity there are three amplitudes, $ h_{t} = h_0 $, $ h_{r} = h_1 $, 
and $ h_{2}^{\rm here} = -h^{\text{RW}}_2 $, which are equivalent to the 
original definitions up to sign.  The odd-parity MP is then
\begin{align}
\label{eq:MPOdd}
p_{aB} = \sum_{\ell, m} h_a^{\ell m} X_B^{\ell m},  
\quad \quad
p_{AB} = \sum_{\ell, m} h_2^{\ell m} X^{\ell m}_{AB} ,
\end{align}
along with the fact that $p_{ab} = 0$.  For the balance of 
this paper we are only concerned with odd-parity.

Regge-Wheeler gauge places the (algebraic) condition on the metric that
$h_{2}^{\text{RW}} = 0$.  In this gauge the odd-parity field equations 
become
\begin{align}
\begin{split}
\label{eq:oddfieldeqns}
 - \partial_t \partial_r  h_r^{\text{RW}} + \partial_r^2  h_t^{\text{RW}} 
- \frac{2}{r} \partial_t  h_r^{\text{RW}} 
- \frac{\ell (\ell + 1) r - 4M}{r^3 f}  h_t^{\text{RW}} &= P^t, \\
\partial_t^2  h_r^{\text{RW}} - \partial_t \partial_r  h_t^{\text{RW}} 
+ \frac{2}{r} \partial_t  h_t^{\text{RW}} 
+ \frac{ (\ell + 2) (\ell - 1) f }{r^2}  h_r^{\text{RW}} &= P^r, \\
 -\frac{1}{f} \partial_t  h_t^{\text{RW}} + f \partial_r  h_r^{\text{RW}} 
+ \frac{2M}{r^2}  h_r^{\text{RW}} &= P ,
\end{split}
\end{align}
where the source amplitudes $P^t$, $P^r$, and $P$ (for each $\ell$ 
and $m$) are odd-parity projections of the stress tensor,
\begin{align}
P^a (t,r)
\equiv  \frac{16 \pi r^2}{\ell (\ell +1)}
\int T^{aB} X_{B}^{*} \, d\O ,
\quad \quad
P (t,r)  \equiv  
16 \pi r^4 \frac{(\ell - 2)!}{(\ell + 2)!} 
\int T^{AB} X^{*}_{AB}\, d\O .
\label{eq:sourceTerms}
\end{align}
We use an asterisk to denote complex conjugation.
The source amplitudes in turn satisfy the contracted Bianchi identity
\begin{equation}
\frac{\partial P^t}{\partial t} + \frac{\partial P^r}{\partial r} 
+ \frac{2}{r} P^r - \frac{(\ell -1)(\ell +2)}{r^2} P = 0 ,
\end{equation}
and the stress tensor itself is taken to be that of a particle in 
geodesic motion on the background geometry.

While in principle the coupled equations (\ref{eq:oddfieldeqns}) 
might be solved to yield the 
metric in RW gauge, the usual approach involves defining and using one of 
several odd-parity master functions and solving a lone wave equation (master 
equation) for this function.  The odd-parity part of the metric then is 
derived from the master function.  An equivalent master function 
representation is used for even-parity in RW gauge.  In our previous paper 
\cite{HE_2010} we used the odd-parity Cunningham-Price-Moncrief (CPM) 
function~\cite{CPM_1978},
 which we refer to here as $\Psi_{o}$.  This gauge-invariant master 
function is defined in Regge-Wheeler gauge by
\begin{equation}
\label{eq:masterOdd}
\Psi_{o}(t,r) \equiv \frac{2 r}{(\ell -1)(\ell +2)}
\left[ \partial_r h_t^{\text{RW}} - \partial_t h_r^{\text{RW}} 
- \frac{2}{r} h_{t}^{\text{RW}} \right] .
\end{equation}
In terms of the tortoise coordinate $r_{*} = r + 2 M \ln (r/2M -1)$, 
$\Psi_{o}$ satisfies the wave equation
\begin{equation}
\label{eq:waveCPM}
\mathcal{W}_{2} \Psi_{o} (t,r) = S_{o} ,
\end{equation}
where ${\cal W}_{2}$ is the spin-2 Regge-Wheeler operator,
a particular case of the spin-$s$ operator
\begin{equation}
\label{eq:WsDef}
\mathcal{W}_{s} = -\frac{\partial^2}{\partial t^2}  + \frac{\partial^2}{\partial r_*^2}
- f \left[ \frac{\ell \l \ell + 1 \r}{r^2} + \frac{2 M (1 - s^2)}{r^3}\right] .
\end{equation}
Later, we will have need for the Fourier transform 
of this operator ($\partial_{t} \to -i \o$), which
we will denote ${\cal L}_{s}$.
The source term for the master equation involves a combination of moments
of the stress tensor,
\begin{equation}
S_{o} (t,r) \equiv \frac{ 2 r f}{(\ell -1)(\ell +2)}
\left[\frac{1}{f} \partial_t P^{r} + f \partial_r P^{t} + \frac{2M}{r^{2}} P^{t} \right] 
= \tilde G_{o} (t)\, \d \left[ r - r_{p}(t) \right] +
\tilde F_{o} (t)\, \d' \left[ r - r_{p}(t) \right] ,
\label{eq:SOdd}
\end{equation}
and is a distribution (see~\cite{HE_2010} for details).
Once the CPM master function is known the MP amplitudes in 
RW gauge can be reconstructed via the expressions
\begin{align}
\label{eq:ReconstructOdd}
h_t^{\text{RW}} (t,r) = \frac{f}{2} \partial_r \l r \Psi_{o} \r 
- \frac{r^2 f}{(\ell -1)(\ell +2)} P^t,
\quad \quad
h_r^{\text{RW}} (t,r) = \frac{r}{2 f} \partial_t \Psi_{o} 
+ \frac{r^2}{(\ell -1)(\ell +2) f} P^r.
\end{align}
Their numerical determination in the time domain with a convergent 
and accurate behavior everywhere including the vicinity of the moving 
particle was the subject of our previous paper.

Because we have reason to consider it in what follows, it is worthwhile 
noting that the original master function of Regge and Wheeler, 
$\Psi_{\text{RW}}$, is not the CPM master function we use.  They are related 
by
\begin{equation}
\label{eq:CPMRW}
\Psi_{\text{RW}} (t,r) = \frac{f}{r} h_r^{\text{RW}} (t,r) = 
\frac{1}{2} \partial_t \Psi_{o} + \frac{r}{(\ell -1)(\ell +2)} P^r .
\end{equation}
The RW master function satisfies an almost identical wave equation,
\begin{equation}
\label{eq:RWDE}
\mathcal{W}_{2} \Psi_{\text{RW}} (t,r) = S_{\text{RW}}  ,
\end{equation}
with the only difference being the source term
\begin{equation}
S_{\text{RW}} (t,r) \equiv \frac{f}{r}
\left[-P^{r} + f \partial_r P - \frac{2}{r} \l 1 - \frac{3M}{r} \r P \right] .
\label{eq:SRW}
\end{equation}

\subsection{Gauge transformations}

The exact form of the field equations (\ref{eq:linearEFE}) will depend upon 
specifying a gauge.  
As mentioned in the Introduction, two frequent choices are Regge-Wheeler (RW)
gauge and Lorenz (L) gauge.  The small gauge generator $\Xi^{\mu}$ that 
transforms the coordinates 
$x^{\mu}_{\text{L}} = x^{\mu}_{\text{RW}} + \Xi^{\mu}$ between the two 
gauges is on the same order of magnitude as the MP, that is 
$|\Xi_{\mu}| \sim |\Bar{p}_{\mu \nu}| \ll 1$.  The MP then transforms as
\begin{equation}
\label{eq:MPTrans}
\Bar{p}_{\mu \nu}^{\text{L}} = \Bar{p}_{\mu \nu}^{\text{RW}} 
- \Xi_{\mu | \nu} - \Xi_{\nu | \mu} + g_{\mu\nu} {\Xi^{\alpha}}_{| \alpha} ,
\end{equation}
where stroke ${| \mu}$ indicates covariant differentiation with respect
to the background metric.  Lorenz gauge requires the following condition on 
the MP, 
\begin{equation}
{{\Bar p}_{\mu \nu}^{\text{L}}}^{|\nu} = 0 .
\end{equation}
Using this condition in Eq.~(\ref{eq:MPTrans}) then provides a wave 
equation that must be satisfied by the gauge generator,
\begin{equation}
\label{eq:gaugeGen}
{\Xi_{\mu|\nu}}^{\nu} = {\Bar{p}_{\mu \nu}^{\text{RW}}}^{|\nu} .
\end{equation}
A gauge generator that satisfies this
equation is unique only up to some $\Xi^{\prime}_{\mu}$ that satisfies 
the homogeneous version of (\ref{eq:gaugeGen}).  Specifying the initial data
and boundary values (if any) removes the residual gauge freedom and fully 
determines the gauge.

We consider next the spherical harmonic decomposition of the gauge vector. 
Momentarily considering again both even- and odd-parity, $\Xi_{\mu}$ can be
broken down into
\begin{equation}
\Xi_{a} = \sum_{\ell,m} 
\left[ {\delta_a}^t \xi^{\ell m}_t (t,r) 
+ {\delta_a}^r \xi^{\ell m}_r (t,r) \right]
Y_{\ell m},
\quad \quad
\Xi_{A} = \sum_{\ell,m} \Big[ 
\xi^{\ell m}_{e} (t,r) Y_{A}^{\ell m} 
+ \xi^{\ell m}_{o} (t,r) X_{A}^{\ell m} \Big] .
\label{eq:xiExpand}
\end{equation}
There are three even-parity amplitudes and one odd-parity amplitude.  We 
will concern ourselves with determining $\xi_t$, $\xi_r$, and $\xi_e$ in 
a subsequent paper.  In this paper we seek to obtain $\xi_o$.  Substituting 
the decomposition of $\Xi_{\mu}$ into Eq.~(\ref{eq:gaugeGen}), we find after 
a bit of calculation that $\xi_o$ satisfies the differential equation
\begin{equation}
\mathcal{W}_1 \xi_{o} (t,r) = 2 f \Psi_{\text{RW}} + f P .
\label{eq:oddParityGauge}
\end{equation}
Once the gauge generator amplitude is known,
we decompose Eq.~(\ref{eq:MPTrans}) in harmonic amplitudes
and see that the odd-parity MP 
amplitudes are transformed by
\begin{align}
\begin{split}
h_{t}^{\text{L}} (t,r) &= h_{t}^{\text{RW}} - \frac{\partial \xi_{o}}{\partial t}, \\
h_{r}^{\text{L}} (t,r) &= h_{r}^{\text{RW}} - \frac{\partial \xi_{o}}{\partial r} 
+ \frac{2}{r} \xi_{o},  \\
h_{2}^{\text{L}} (t,r) &= - 2 \xi_{o}.
\label{eq:push}
\end{split}
\end{align}

\subsection{Local nature of the metric perturbation
and gauge generator at $r=r_p(t)$}

The RHS of Eq.~(\ref{eq:oddParityGauge}) is singular at the location of the 
particle.  In this sense Eq.~(\ref{eq:oddParityGauge}) 
is very similar to Eq.~(\ref{eq:waveCPM}).
In Ref.~\cite{HE_2010} we examined Eq.~(\ref{eq:waveCPM})
to determine the local behavior of $\Psi_{o}$.
Assuming $\Psi_{o} = \Psi_{o}^{+}\, \th \left[ r - r_{p} (t) \right] +
\Psi_{o}^{-}\, \th \left[ r_{p} (t) - r \right]$, we calculated
jumps in the field, $\llbracket \Psi_{o} \rrbracket_{p}$, and in its
radial derivative, $\llbracket \partial_{r} \Psi_{o} \rrbracket_{p}$.
We use a subscript $p$ to indicate that a function of $r$ is evaluated
at the location of the particle, $r = r_{p}(t)$, becoming a function of time.

Following the same logic here, we postulate a form for 
the gauge amplitude of
$
\xi_{o}
=
\xi_{o}^{+}\, \th \left[ r - r_{p} (t) \right]
+
\xi_{o}^{-}\, \th \left[ r_{p} (t) - r \right].
$
Then, similar analysis to that found in Ref.~\cite{HE_2010}, indicates that
$\xi_{o}$ is $C^{0}$, i.e. $\llbracket \xi_{o} \rrbracket_{p} = 0$.
Further, we find the jump in the first radial derivative is
\begin{equation}
\left\llbracket \partial_{r} \xi_{o} \right\rrbracket_{p} (t)
=
\frac{f_{p}}{f_{p}^{2} - {\dot r_{p}}^{2}}  p .
\end{equation}
Here $p(t)$ comes from the source amplitude $P$.  All three source 
amplitudes are delta distributions with time dependent amplitudes: 
$P = p(t)\, \d [r - r_p(t)]$, $P^t = p^t(t)\, \d [r - r_p(t)]$, and
$P^r = p^r(t)\, \d [r - r_p(t)]$.

Having computed the expected jumps in $\xi_{o}$ and its radial derivative,
we can use Eq.~(\ref{eq:push}) to find the jumps in the Lorenz gauge MP 
amplitudes.  As with the fields $\Psi_{o}$ and $\xi_{o}$, we expect each 
MP amplitude to consist of left and right side differentiable functions 
that are joined at the location of the particle by Heaviside functions.
We calculated the jumps in $h_{t}^{\text{RW}}$ and $h_{r}^{\text{RW}}$ in 
Ref.~\cite{HE_2010}.  The discontinuities in the RW gauge MP amplitudes 
are exactly canceled out by terms arising from the derivatives of $\xi_{o}$ 
and all three Lorenz gauge amplitudes are $C^{0}$ as expected.  The jumps 
in their first derivatives are also analytically computable
(either by examining the jumps in the higher-order derivatives of 
$\xi_{o}$ or more simply by directly analyzing the Lorenz gauge
field equations).  We find
\begin{align}
\label{eq:LGJumps}
\left\llbracket \partial_{r} h_{t} \right\rrbracket_{p} (t)
= \frac{f_{p}^{2}}{f_{p}^{2} - {\dot r_{p}}^{2}}  p^{t},
\quad \quad
\left\llbracket \partial_{r} h_{r} \right\rrbracket_{p} (t)
= - \frac{1}{f_{p}^{2} - {\dot r_{p}}^{2}}  p^{r},
\quad \quad
\left\llbracket \partial_{r} h_{2} \right\rrbracket_{p} (t)
= - \frac{2 f_{p}}{f_{p}^{2} - {\dot r_{p}}^{2}}  p .
\end{align}
We use these expressions later (see Fig.~\ref{fig:convVsN}) 
as a powerful check that we have correctly
solved the gauge transformation equations to high accuracy.

\section{Two EHS-like methods for equations with non-compact sources}
\label{sec:twomethods}

In Ref.~\cite{HE_2010} we solved Eq.~(\ref{eq:waveCPM}) for a variety of 
eccentric orbits.  We used a FD approach to find the Fourier harmonic 
modes of $\Psi_{o}$ and transformed back to the TD using the EHS method.  
The EHS method was first applied to wave equations with delta function 
sources.  It allows TD reconstruction of the spherical harmonic amplitudes 
with exponential convergence, circumventing the Gibbs phenomenon that 
otherwise arises from solving equations with discontinuous or singular 
sources.  Our application of the method also demonstrated it could be 
applied to sources with a derivative of a delta function.

With only a change in spin parameter, Eq.~(\ref{eq:oddParityGauge}) 
has a similar differential operator as Eq.~(\ref{eq:waveCPM}).  Where the 
two equations differ markedly is in their source terms.  While 
the source in Eq.~(\ref{eq:waveCPM}) is point-singular and compact, the
source in Eq.~(\ref{eq:oddParityGauge}) is both distributional and 
non-compact.  We can use the linearity of the equation to split off the 
singular part and split the generator into two parts,
$\xi_{o} = \xi_{o}^{\text{ext}} + \xi_{o}^{\text{sing}}$,
that satisfy separate equations,
\begin{align}
{\cal W}_{1} \xi_{o}^{\text{sing}} (t,r)
&=
 f_{p} p(t)\, \d \left[ r - r_{p} (t) \right], 
\label{eq:xiSing}
 \\
{\cal W}_{1} \xi_{o}^{\text{ext}} (t,r)
&=
2 f \Psi_{\text{RW}} .
\label{eq:xiExtended}
\end{align}
While the former equation can be solved using the EHS method, the latter's 
extended source term is more problematic.  The extended source is both 
non-compact and has a time-dependent discontinuity that moves periodically 
between $r_{\text{min}}$ and $r_{\text{max}}$ as the particle orbits.  In 
this section we present two equivalent methods for solving 
Eq.~(\ref{eq:xiExtended}) using FD methods, both of which provide 
exponential convergence upon returning to the TD.

As discussed earlier, an eccentric orbit on Schwarzschild provides two 
fundamental frequencies.  When we Fourier transform 
Eq.~(\ref{eq:xiExtended}), we have a two-fold countably infinite frequency 
spectrum,
\begin{equation}
\o \equiv \o_{mn} = m \O_{\varphi} + n \O_{r}, \quad \quad m,n \in \mathbb{Z}.
\end{equation}
The Fourier transform and standard TD reconstruction of 
$\xi_{o}^{\text{ext}} (t,r)$ is then
\begin{equation}
\tilde \xi_{o}^{\text{ext}} (r) \equiv \frac{1}{T_r} \int_0^{T_r} dt \ \xi_{o}^{\text{ext}} (t,r)
\, e^{i \o t},
\quad \quad
\xi^{\text{ext}}_{o} (t,r) = \sum_{n=-\infty}^{\infty} 
\tilde \xi^{\text{ext}}_{o}(r) \, e^{-i \o t}.
\end{equation}
Note that in addition to the already suppressed indices $\ell$ and $m$,
FD quantities have a third implied index, $n$.
Equivalent expressions are used for the Fourier transforms and 
series representations of the other fields that we consider below.
Note that while we use the standard tilde \~ \ notation with
$\tilde \xi^{\text{ext}}_{o}$ to indicate a FD quantity, 
for other quantities we try to maintain consistency with 
previous literature by changing the base 
symbol.
For the TD function $\Psi_{\text{RW}}$, we
write $R_{\text{RW}}$ in the FD. Similarly, for its  TD 
source term $S_{\text{RW}}$, we write  $Z_{\text{RW}}$ in the FD.

\subsection{First approach: partial annihilator method}
\label{sec:partialannihilator}

Our first method for solving Eq.~(\ref{eq:xiExtended}) is a generalization
of the standard method of annihilators used for solving inhomogeneous 
differential equations.  It hinges on finding an ``annihilator,'' a
differential operator which gives a vanishing result after acting
on the source.  Then, one can act with the annihilator on both sides of the 
differential equation.  What results is a homogeneous differential equation
of higher order.  Our strategy for solving Eq.~(\ref{eq:xiExtended}) is 
essentially the same, except that because the initial source is discontinuous 
the operator that we find does not completely annihilate the RHS but instead 
converts it to a distribution.  Hence, we refer to the operator as 
a \emph{partial annihilator}.

The RHS of Eq.~(\ref{eq:xiExtended}) is well suited to the method
of partial annihilators because the Regge-Wheeler variable satisfies its own 
wave equation with a point-singular source, (\ref{eq:RWDE}).  Therefore, 
upon dividing Eq.~(\ref{eq:xiExtended}) by $f$, we can take ${\cal W}_{2}$ as 
the partial annihilator and act on both sides of the equation
\begin{equation}
f {\cal W}_{2} \left( \frac{1}{f}\, {\cal W}_{1} \xi (t,r) \right)
 = 2 f S_{\text{RW}} (t, r_{p}(t)) 
 = 2 f
  \l \tilde G_{\text{RW}} (t)\, \d \left[ r - r_{p}(t) \right]
 +
 \tilde F_{\text{RW}} (t)\, \d' \left[ r - r_{p}(t) \right] \r .
\label{eq:PATD}
\end{equation}
For simplicity here and in the remainder of this section we drop the 
$^{\text{ext}}_{o}$ tags.  We have multiplied back through by $f$ to 
ensure that the leading-order derivatives have unit coefficients.
This differential equation is now fourth-order, but its source is 
point-singular.  This allow us to solve it using the EHS method, generalized
to fourth-order equations.  The specific form of the source in 
Eq.~(\ref{eq:PATD}) is given by Martel \cite{Martel_2004}, though we 
assume that both $\tilde G_{\text{RW}}$ and $\tilde F_{\text{RW}}$ have
been evaluated at $r = r_{p}(t)$.

We Fourier transform Eq.~(\ref{eq:PATD}) to obtain the FD equation 
\begin{equation}
f {\cal L}_{2} \left( \frac{1}{f}\, {\cal L}_{1} \tilde \xi (r) \right)
 = 2 f Z_{\text{RW}} (r).
\label{eq:PAFD}
\end{equation}
The Fourier transform averages the point source motion in time and produces
$Z_{\text{RW}}(r)$ which has support only within the source libration region 
$r_{\rm min}<r<r_{\rm max}$.

There are four linearly independent homogeneous solutions to
Eq.~(\ref{eq:PAFD}).   Two of these are the solutions to the 
second-order equation ${\cal L}_{1} \tilde \xi = 0$ and we denote them by 
$\tilde{\xi}^{\pm}_{h2}$.  One behaves asymptotically as an outgoing wave 
at infinity while the other is downgoing at the horizon
\begin{equation}
\tilde \xi_{h2}^{-} 
\sim e^{-i \o r_{*}} \quad (r \to 2M), \quad \quad
\tilde \xi_{h2}^{+} 
\sim e^{i \o r_{*}}  \quad (r \to \infty).
\end{equation}
The other two solutions only satisfy the full fourth-order 
equation ${\cal L}_{2} (f^{-1} \, {\cal L}_{1} \tilde \xi ) = 0$.
As such we give them the label ${h4}$ and asymptotic analysis shows that
\begin{equation}
\tilde \xi_{h4}^{-} 
\sim f(r)\, e^{-i \o r_{*}}  \quad (r \to 2M),  \quad \quad
\tilde \xi_{h4}^{+}
 \sim r\, e^{i \o r_{*}}  \quad (r \to \infty) .
\label{eq:homog4th}
\end{equation}
These four solutions form a basis spanning the space of homogeneous
solutions of Eq.~(\ref{eq:PAFD}).  The particular solution will be a 
linear combination of these with variable coefficients
\begin{equation}
\tilde \xi_{p} (r)
 = c_{h2}^{-} (r)\, 
 \tilde \xi_{h2}^{-} (r) 
+ c_{h2}^{+} (r)\, 
\tilde  \xi_{h2}^{+} (r)
+ c_{h4}^{-} (r)\, 
\tilde  \xi_{h4}^{-} (r) 
+ c_{h4}^{+} (r)\, 
\tilde  \xi_{h4}^{+} (r).
\end{equation}
The four normalization functions $c_{h2/h4}^{\pm} (r)$ are fixed by the 
method of variation of parameters, which entails solving the equations
\begin{equation}
\label{eq:varOfPar}
\frac{d c^{\pm}_{h2/h4}}{d r_{*}} 
= 2 f Z_{\text{RW}} (r) \frac{W^{\pm}_{h2/h4} (r)}{W(r)}.
\end{equation}
Here $W(r)$ is the Wronskian and $W^{\pm}_{h2/h4} (r)$ is the ``modified 
Wronskian'' (Cramer's rule), which is the Wronskian with the column 
corresponding to the $\xi_{h2/h4}^{\pm} (r)$ homogeneous solution replaced 
by the column vector $(0, 0, 0, 1)$.  Note that because the differential 
operator in Eq.~(\ref{eq:PATD})
is written in terms of $r_{*}$, the derivatives within the Wronskian must
also be taken with respect to $r_{*}$ and the LHS of Eq.~(\ref{eq:varOfPar})
is a derivative taken with respect to $r_{*}$.
For the two ``$+$'' equations, the integral form of Eq.~(\ref{eq:varOfPar}) is 
(we change the variable of integration to $r$ and see the factor of $f$
cancel)
\begin{equation}
c^{+}_{h2/h4} (r) 
= 2 \int_{r_{\rm min}}^{r}
\bigg[ 
\frac{1}{T_{r}}
\int_{0}^{T_{r}}  \Big( 
\tilde G_{\text{RW}} (t)\, \d \left[ r' - r_{p}(t) \right] 
+\tilde F_{\text{RW}} (t)\, \d' \left[ r' - r_{p}(t) \right]
 \Big) e^{i \o t} dt \bigg] 
 \frac{W^{+}_{h2/h4} (r')}{W(r')} dr'.
\end{equation}
Likewise, for the two ``$-$'' equations (note the change on the limits of 
integration),
\begin{equation}
c^{-}_{h2/h4} (r) 
= 2 \int_{r}^{r_{\rm max}}
\bigg[ 
\frac{1}{T_{r}}
\int_{0}^{T_{r}}  \Big( 
\tilde G_{\text{RW}} (t)\, \d \left[ r' - r_{p}(t) \right] 
+\tilde F_{\text{RW}} (t)\, \d' \left[ r' - r_{p}(t) \right]
 \Big) e^{i \o t} dt \bigg] 
 \frac{W^{-}_{h2/h4} (r')}{W(r')} dr'.
\end{equation}
The EHS method requires knowing only the terminal values of the four 
functions $c^{\pm}_{h2/h4}(r)$, i.e.,
$C^{+}_{h2/h4} = c^{+}_{h2/h4}(r_{\text{max}})$ and 
$C^{-}_{h2/h4} = c^{-}_{h2/h4}(r_{\text{min}})$.  Switching the order of 
integration and integrating by parts, we find
\begin{equation}
C^{\pm}_{h2/h4} = 
\frac{2}{T_{r}}
\int_{0}^{T_{r}}  
\Bigg\{ 
\tilde G_{\text{RW}} (t)  
\frac{W^{\pm}_{h2/h4} (r_{p})}{W(r_{p})} 
- \tilde F_{\text{RW}} (t) 
\left[ 
-  \frac{W^{\pm}_{h2/h4} (r_{p} )}{W(r_{p})^{2}} \partial_{r} W (r_{p})
+ \frac{ \partial_{r} W^{\pm}_{h2/h4} (r_{p})}{W (r_{p})} \right]
\Bigg\} e^{i \o t} dt.
\end{equation}
At this point we define the EHS in the FD to be 
\begin{equation}
\tilde  \xi_{h}^{-} (r) \equiv C^{-}_{h2}\tilde  \xi_{h2}^{-} (r) 
+ C^{-}_{h4} \tilde \xi_{h4}^{-} (r), 
\quad \quad
\tilde \xi_{h}^{+} (r) \equiv C^{+}_{h2} \tilde \xi_{h2}^{+} (r) 
+ C^{+}_{h4} \tilde \xi_{h4}^{+} (r), 
\end{equation}
and the EHS in the TD are defined by the Fourier
sums (recall the suppressed $\ell, m, n$ indices)
\begin{equation}
\xi^{\pm} (t,r) \equiv \sum_{n} \tilde  \xi_{h}^{\pm} (r) e^{-i \o t}.
\end{equation}
The extension of these solutions to $r=r_{p}(t)$ then gives the desired
solution to Eq.~(\ref{eq:xiExtended}),
\begin{equation}
\xi_{o}^{\text{ext}} (t,r) =
\xi^{+} (t,r)\, \th \left[ r - r_{p} (t) \right]
+
\xi^{-} (t,r)\, \th \left[ r_{p} (t) - r \right].
\label{eq:xiPartAnn}
\end{equation}

 
\subsection{Second approach: method of extended particular solutions}
\label{sec:epsmethod}

Now we look for a solution to Eq.~(\ref{eq:xiExtended}) 
that does not require a partial annihilator.
In the FD the equation transforms to
\begin{equation}
{\cal L}_{1} \tilde \xi_{o}^{\text{ext}} (r)
= 2 f R_{\text{RW}}.
\label{eq:xiFD}
\end{equation}
Again, for notational simplicity we drop the $^{\rm ext}_{o}$ tags
for the remainder of this section.  In the end we want solutions to
Eq.~(\ref{eq:xiFD}) that allow us to form an exponentially converging 
solution to Eq.~(\ref{eq:xiExtended}) when we transfer to the TD.  This 
will require a new technique which we call
\emph{extended particular solutions}, and is closely analogous to the
EHS method.  First, though, we 
consider how to get the correct causal solution to Eq.~(\ref{eq:xiFD})
from a ``standard'' approach.

In the subsequent sections we make a distinction between quantities with
$^{\infty}$ and $^{H}$ tags which designate functions computed
from a ``standard'' RHS source (Eq.~(\ref{eq:RStdRW}) below) 
and those quantities with
$^{+}$ and $^{-}$ tags which designate functions computed 
from an ``extended'' RHS source (Eq.~(\ref{eq:R_EHS}) below). 
Because the homogeneous solutions do not depend on the source, we always
tag them with $^{+}$ or $^{-}$.  We distinguish between
particular and homogeneous solutions by using the respective subscripts
$_{p}$  and $_{h}$.

\vspace{-.2cm}
\subsubsection{Finding standard FD solutions with causal boundary conditions}

By examining the source, $2f R_{\text{RW}}$, and the differential operator,
${\cal L}_{1}$, we can obtain asymptotic and Taylor expansions of the 
particular solution $\tilde \xi_{p}$ near infinity and the horizon, 
respectively.  The expansions are useful numerically but for our purposes 
here we need only consider the leading asymptotic dependence.
(See App.~\ref{sec:asympExp} for discussion of the asymptotic expansion 
($r \rightarrow \infty$) of $\tilde \xi_{p}$ and how it couples to the 
expansion of $R_{\text{RW}}$.)

Consider first the spatial infinity side.  Let the RW function have an 
asymptotic amplitude $C^{+}_{\text{RW}}$, so 
$R_{\text{RW}} = C^{+}_{\text{RW}} e^{ i \o r_{*}}$ as $ r_{*} \to  \infty$.
We then make the ansatz that 
$\tilde \xi_{p} = C^{\infty}_{p} r e^{i \o r_{*}}$ as $r_{*} \to \infty$.  
Using an asymptotic approximation to Eq.~(\ref{eq:xiFD}) we find
\begin{align}
\l \frac{d^{2}}{d r_{*}^{2}} + \o^{2} \r 
\l C^{\infty}_{p} r e^{ i \o r_{*}} \r
= 2 C^{+}_{\text{RW}}  e^{ i \o r_{*}}  \quad \quad
\Rightarrow
\quad\q
C^{\infty}_{p} 
= \frac{1}{i \o} C^{+}_{\text{RW}}.
\label{eq:xiCPlus}
\end{align}
Therefore, the asymptotic form of $\tilde \xi_{p}^{\infty}$ is 
\begin{equation}
\tilde \xi_{p}^{\infty} 
= - \frac{i}{\o} C^{+}_{\text{RW}} r  e^{i \o r_{*}}, \quad \quad 
r \to  \infty.
\label{eq:xipInfty}
\end{equation}
Next we consider the horizon side.  The RW function is asymptotically
$R_{\text{RW}} = C^{-}_{\text{RW}} e^{-i\o r_{*}}$ as $r_{*} \to -\infty$.
In this case we expect the particular solution to behave as
$\tilde \xi_{p} = C^{H}_{p}\, f e^{-i \o r_{*}}$ as $r_{*} \to -\infty$. 
Again, acting with the near-horizon leading parts of the differential 
operator we find
\begin{align}
\l \frac{d^{2}}{d r_{*}^{2}} + \o^{2} \r 
\l C^{H}_{p} f e^{ -i \o r_{*}} \r
= 2 f C^{-}_{\text{RW}} e^{ -i \o r_{*}}  \quad 
\Rightarrow
\quad
C^{H}_{p} 
=  2 \l \frac{1}{4M^{2}}- \frac{i \o}{M} \r^{-1} C^{-}_{\text{RW}}.
\label{eq:xiCMinus}
\end{align}
Therefore, the near-horizon form of $\tilde \xi_{p}^{H}$ is
\begin{equation}
\tilde \xi_{p}^{H} 
= 
2 \l \frac{1}{4M^{2}}- \frac{i \o}{M} \r^{-1} C^{-}_{\text{RW}}
 f  e^{-i \o r_{*}}, \quad \quad 
r \to 2M.
\label{eq:xipH}
\end{equation}

We can use Eqs.~(\ref{eq:xipInfty}) and (\ref{eq:xipH}) to set boundary 
conditions (B.C.'s) for two separate integrations of the inhomogeneous 
differential equation, (\ref{eq:xiFD}) (yielding two different particular 
solutions that differ by some homogeneous solution).  This integration 
requires the source $R_{\text{RW}}$, which is itself the solution to the 
differential equation
\begin{equation}
{\cal L}_{2} R_{\text{RW}} (r)  = Z_{\text{RW}} (r).
\label{eq:RWeq}
\end{equation}
We find it by variation of parameters, which yields
\begin{equation}
R^{\text{std}}_{\text{RW}} (r) 
= c^{+}_{\text{RW}} (r) \hat R^{+} (r) + c^{-}_{\text{RW}} (r) \hat R^{-} (r),
\label{eq:RStdRW}
\end{equation}
where $\hat R^{\pm} (r)$ are unit-normalized homogeneous solutions to 
Eq.~(\ref{eq:RWeq}).  Note that 
$c_{\text{RW}}^{+} (r \ge r_{\text{max}}) = C^{+}_{\text{RW}}$
and 
$c_{\text{RW}}^{-} (r \le r_{\text{min}}) = C^{-}_{\text{RW}}$.  Furthermore, 
the solution $R^{\text{std}}_{\text{RW}} (r)$ is \emph{not} the same as the 
EHS to Eq.~(\ref{eq:RWeq}).

Having solved Eq.~(\ref{eq:RWeq}) for the source term in 
Eq.~(\ref{eq:RStdRW}), and determined the 
B.C.'s, we are ready to solve Eq.~(\ref{eq:xiFD}).  The idea is to find the 
two different particular solutions, neither of which has the proper causal 
behavior, and then correct for the acausality by adding appropriate
homogeneous solutions.  The result of this process is a 
solution to Eq.~(\ref{eq:xiFD}) with causal behavior on both sides.  The 
details follow in a series of steps. 

\begin{enumerate}
\item 
{\sc Solve for the particular solution from the spatial-infinity side
(see Fig.~\ref{fig:xiPInfty}).}

We set a B.C. to Eq.~(\ref{eq:xiFD}) using 
Eq.~(\ref{eq:xipInfty}) and integrate  to large negative
$r_{*}$ using Eq.~(\ref{eq:RStdRW}) as the source.
Although the starting B.C.
specified no homogeneous contribution, homogeneous
solutions on the horizon side will be excited.
See the left side of Fig.~\ref{fig:xiPInfty}. 
The particular solution integrated from 
the spatial infinity side has the asymptotic behavior
\begin{equation}
\tilde \xi_{p}^{\infty} 
=
\left\{
\begin{array}{ll}
C^{\infty}_{p} r e^{i \o r_{*}}, 
& r_{*} \to +\infty, \\
C^{H}_{p} f e^{-i \o r_{*}} 
+ \k^{+} e^{i \o r_{*}} + \k^{-} e^{-i \o r_{*}},
&
 r_{*} \to -\infty .
\end{array}
\right.
\label{eq:xiPInftyOnH}
\end{equation}  
The term with coefficient $C^{H}_{p}$ is the part directly dependent on the 
source and it is sub-dominant in comparison to the homogeneous solutions.
The coefficients $\k^{\pm}$ are to-be-determined.
Importantly, $\k^{+} e^{i \o r_{*}}$ is an acausal term (upgoing from the 
past horizon).

\begin{figure}[h!]
\includegraphics[scale=.95]{xi_p_Plus}
\caption{
\label{fig:xiPInfty}
Integration from large $r_{*}$ of the particular solution, 
$\tilde \xi_{p}^{\infty}$.
Dotted lines indicate the source libration region.
}
\vspace{-.3cm}
\end{figure}

\item {\sc Solve for the particular solution from the horizon side
(see Fig.~\ref{fig:xiPH}).}

We set a B.C. to Eq.~(\ref{eq:xiFD}) using Eq.~(\ref{eq:xipH})
and integrate to large positive
$r_{*}$ using Eq.~(\ref{eq:RStdRW}) as the source.
 Although the starting B.C.
specified no homogeneous contribution, homogeneous
solutions on the spatial infinity side will be excited.
The effect can be seen on the right side of Fig.~\ref{fig:xiPH}.  
The particular solution dominates, but the constant offset between real
and imaginary parts indicates the presence of 
asymptotically-constant-amplitude homogeneous terms.  Analysis shows that 
the particular solution integrated from the horizon will behave as
\begin{equation}
\tilde \xi_{p}^{H} 
=
\left\{
\begin{array}{ll}
C^{H}_{p} f e^{ -i \o r_{*}},
& r_{*} \to -\infty, \\
C^{\infty}_{p} r e^{i \o r_{*}} 
+ \lambda^{-} e^{-i \o r_{*}} + \lambda^{+} e^{i \o r_{*}},
&
 r_{*} \to +\infty.
\end{array}
\right.
\label{eq:xiPHOnInfty}
\end{equation}  
The coefficients $\lambda^{\pm}$ are to-be-determined and again we find an 
acausal term ($\lambda^{-} e^{-i \o r_{*}}$), which in this case is ingoing 
from past null infinity.

\begin{figure}[h!]
\includegraphics[scale=.95]{xi_p_Minus}
\caption{
\label{fig:xiPH}
Integration from large negative
$r_{*}$ of the particular solution, $\tilde \xi_{p}^{H}$.
Dotted lines indicate the source libration region.
}
\vspace{-.3cm}
\end{figure}

\item {\sc Solve for the homogeneous solution from the spatial infinity
side.}

Next, we set an outgoing B.C. to the homogeneous version of 
Eq.~(\ref{eq:xiFD}) at large positive $r_{*}$.
We integrate to solve the scattering problem for
reflection and transmission amplitudes $R^{+}$ and $T^{+}$
\cite{Chandra_1983}.  In the terminology of Gal'tsov \cite{Galtsov}
this is an ``up'' mode,
\begin{equation}
\tilde \xi_{h}^{+} 
=
\left\{
\begin{array}{ll}
T^{+} e^{i \o r_{*}},
& r_{*} \to +\infty, \\
R^{+} e^{-i \o r_{*}} +  e^{i \o r_{*}},
&
 r_{*} \to -\infty.
\end{array}
\right.
\label{eq:xiHPlus}
\end{equation}  
Scaled appropriately, this solution can be added to
Eq.~(\ref{eq:xiPInftyOnH}) to remove its acausality.

\item {\sc Solve for the homogeneous solution from the horizon side.}

We set a downgoing B.C. to the homogeneous version of 
Eq.~(\ref{eq:xiFD}) at large negative $r_{*}$ and
integrate to solve the scattering problem for reflection and transmission 
amplitudes $R^{-}$ and $T^{-}$.  This is an ``in'' mode,
\begin{equation}
\tilde \xi_{h}^{-} 
=
\left\{
\begin{array}{ll}
 T^{-} e^{-i \o r_{*}},
& r_{*} \to -\infty, \\
R^{-} e^{i \o r_{*}} +  e^{-i \o r_{*}},
&
 r_{*} \to +\infty.
\end{array}
\right.
\label{eq:xiHMinus}
\end{equation}  
Scaled appropriately, this solution can be added to
Eq.~(\ref{eq:xiPHOnInfty}) to remove its acausality.

\item {\sc Resolve the acausality in the particular solutions 
(see Fig.~\ref{fig:xi_p_Std}).}
\label{resolveStep}

The acausal piece in Eq.~(\ref{eq:xiPInftyOnH}) is
$\k^{+} e^{i \o r_{*}}$.
Using Eq.~(\ref{eq:xiHPlus}), we can remove this by subtracting 
$\k^{+} \tilde \xi^{+}_{h}$,
\begin{equation}
\tilde \xi^{\infty}_{p} - \k^{+} \tilde \xi^{+}_{h}
=
\left\{
\begin{array}{ll}
C^{\infty}_{p} r e^{i \o r_{*}} 
- \k^{+} T^{+}  e^{i \o r_{*}},
& r_{*} \to +\infty, \\
C^{H}_{p} f e^{-i \o r_{*}} 
+ \k^{-} e^{-i \o r_{*}}
- \k^{+} R^{+} e^{-i \o r_{*}},
&
 r_{*} \to -\infty.
\end{array}
\label{eq:xiPA}
\right.
\end{equation}
The acausal piece in Eq.~(\ref{eq:xiPHOnInfty}) is
$\lambda^{-} e^{-i \o r_{*}}$.
Using Eq.~(\ref{eq:xiHMinus}), we can remove this by subtracting 
$\lambda^{-} \tilde \xi^{-}_{h}$,
\begin{equation}
\tilde \xi^{H}_{p} - \lambda^{-} \tilde \xi^{-}_{h}
=
\left\{
\begin{array}{ll}
C^{H}_{p} f e^{ -i \o r_{*}}
- \lambda^{-} T^{-} e^{ -i \o r_{*}},
& r_{*} \to -\infty, \\
C^{\infty}_{p} r e^{i \o r_{*}} 
 + \lambda^{+} e^{i \o r_{*}}
- \lambda^{-} R^{-} e^{ i \o r_{*}} ,
&
 r_{*} \to +\infty.
\end{array}
\right.
\label{eq:xiPB}
\end{equation}  
Eqs.~(\ref{eq:xiPA}) and (\ref{eq:xiPB}) are both solutions
to Eq.~(\ref{eq:xiFD}) and both satisfy the 
causal nature of the problem.  Therefore they must be equal.
In order to form them, we must know $\k^{+}$ and $\lambda^{-}$.
We find them by setting Eqs.~(\ref{eq:xiPA}) and (\ref{eq:xiPB}) 
and their first derivatives equal at any point,  
\begin{align}
\tilde \xi^{H}_{p} - \lambda^{-} \tilde \xi^{-}_{h} &=
\tilde \xi^{\infty}_{p} - \k^{+} \tilde \xi^{+}_{h},
\label{eq:XiEq}
 \\
\partial_{r_{*}} \tilde \xi^{H}_{p} - \lambda^{-} \partial_{r_{*}} \tilde \xi^{-}_{h} &=
\partial_{r_{*}} \tilde \xi^{\infty}_{p} - \k^{+} \partial_{r_{*}} \tilde \xi^{+}_{h} .
\label{eq:diffXiEq}
\end{align}
We solve these equations for $\k^{+}$ and $\lambda^{-}$ and 
form  $\tilde \xi^{H}_{p} - \lambda^{-} \tilde \xi^{-}_{h}$ and 
$\tilde \xi^{\infty}_{p} - \k^{+} \tilde \xi^{+}_{h}$, which are equivalent.
In principle one could pick any point and expect the same result.  In 
practice, slight numerical differences occur.  In fact, we use several points 
to determine these constants and use the discrepancies that are found as a 
measure of the order of magnitude of the error.  Cumulative numerical error 
in the solutions is addressed in Sec.~\ref{sec:results}.

\begin{figure}[h!]
\includegraphics[scale=1]{xi_ph_Minus}
\caption{
\label{fig:xi_p_Std}
Causally correct solution to Eq.~(\ref{eq:xiFD}),
$\tilde \xi^{\text{std}}$.
Dotted lines indicate the source libration region.}
\end{figure}
\end{enumerate}

The function $\tilde \xi^{\text{std}}_{p} 
= \tilde \xi^{H}_{p} - \lambda^{-} \tilde \xi^{-}_{h}
= \tilde \xi^{\infty}_{p} - \k^{+} \tilde \xi^{+}_{h}$ represents the 
standard solution to Eq.~(\ref{eq:xiFD}).  If the TD source were 
differentiable, we would be able to find the corresponding TD solution via 
an exponentially converging Fourier synthesis.  However, the source in this 
case is non-differentiable and we need an EHS-like trick to complete the
method.

\subsubsection{Restoring exponential convergence with extended particular 
solutions}

The EHS of the Regge-Wheeler equation Eq.~(\ref{eq:RWeq}) are found by 
taking the constants $C^{\pm}_{\text{RW}}$ and scaling the unit-normalized 
homogeneous solutions 
\begin{equation}
R^{\pm}_{\text{RW}} (r) \equiv C^{\pm}_{\text{RW}} \hat R^{\pm} (r) .
\label{eq:R_EHS}
\end{equation}
These solutions are defined for all $r>2M$.  In like fashion we seek to find 
FD EPS of Eq.~(\ref{eq:xiFD}) and denote these by $\tilde \xi^{\pm}$.  We 
first find $\tilde \xi^{\pm}_{p}$ by separately integrating 
Eq.~(\ref{eq:xiFD}) with the modified source terms $R^{\pm}_{\text{RW}}$.  
The solutions are each made to match the exterior behavior of 
$\tilde \xi^{\text{std}}_{p}$ by adding the correctly scaled homogeneous 
solutions found in Step \ref{resolveStep} above.  We then define
\begin{equation}
\tilde \xi^{+} \equiv \tilde \xi_{p}^{+} - \k^{+} \tilde \xi_{h}^{+},
\quad \quad
\tilde \xi^{-} \equiv \tilde \xi_{p}^{-} - \lambda^{-} \tilde \xi_{h}^{-}.
\end{equation}
See Fig.~\ref{fig:xi_EPS}, which contrasts Fig.~\ref{fig:xi_p_Std}
in the source region.

\begin{figure}[h!]
\includegraphics[scale=1]{xi_EPS.pdf}
\caption{
\label{fig:xi_EPS}
Causally correct extended particular solution,
$\tilde \xi^{\pm}$.  Note 
the difference within the libration region (shown in dotted lines), 
between $\tilde \xi^{\pm}$ and 
$\tilde \xi^{\text{std}}$ (Fig.~\ref{fig:xi_p_Std}).
The difference in TD convergence between these two solutions is
shown later in Fig.~\ref{fig:convVsX}.}
\end{figure}

These FD EPS can be transferred to the TD via Fourier series
\begin{equation}
\xi^{\pm} (t,r) \equiv \sum_{n} \tilde \xi^{\pm} (r) e^{- i \o t}.
\end{equation}
The solution to Eq.~(\ref{eq:xiExtended}) is then the weak solution,
\begin{equation}
\xi_{o}^{\text{ext}} (t,r) =
 \xi^{+} (t,r)\, \th \left[ r - r_{p} (t) \right]
+
\xi^{-} (t,r)\, \th \left[ r_{p} (t) - r \right].
\label{eq:EPS}
\end{equation}
The support for this claim has three legs.  Firstly, the same arguments about 
EHS, based on analytic continuation, made by Barack, Ori, and Sago in 
Ref.~\cite{BOS} appear to apply in extension to Eq.~(\ref{eq:xiExtended}) 
as well.  Secondly, we demonstrate existence \emph{numerically} by integrating 
the equation, with causal boundary conditions, and checking that the jump 
conditions (internal boundary conditions) at the particle are satisfied.  
One then appeals to the linearity of the equation to establish uniqueness.  
Finally, we have an independent numerical solution found through the method 
of partial annihilators and given in Eq.~(\ref{eq:xiPartAnn}).  We have 
confirmed that the two methods give entirely consistent solutions.  These 
results are covered in detail in the next section.

\begin{figure}[h!]
\includegraphics[scale=1]{xi_odd}
\caption{
\label{fig:xiOdd}
The $\ell = 2, m = 1$ mode of the odd-parity RW-to-Lorenz gauge generator 
amplitude $\xi_{o}$.  The Lorenz gauge MP amplitude $h_{2}$ differs from 
$\xi_{o}$ only by a factor of $-2$.  Dotted lines indicate the region of 
libration.  Orbital parameters are given in the text.  While the field 
$h_{2}$ grows asymptotically, upon transforming to an orthonormal frame it 
would contribute a term that falls off as $1/r$.}
\end{figure}

\section{Results}
\label{sec:results}

The methods of the previous section allow us to transform odd-parity 
solutions of the first-order Einstein equations from RW to Lorenz gauge.  
As an example, we consider an orbit with eccentricity $e = 0.764124$ and 
semi-latus rectum $p = 8.75455$.  This orbit was used in 
Fig.~\ref{fig:hthrRW} where we showed the RW amplitudes $h_{t}$ and $h_{r}$ 
for $\ell = 2$ and $m = 1$ ($h_{2} = 0$).  The MPs can be evaluated at any 
time but we chose to display results at $t = 93.58$ (where $t = 0$ is at 
the periapsis).  The RW modes are discontinuous at $r=r_p(t)$ and lack 
asymptotic flatness.  The gauge generator to go from RW to Lorenz gauge 
is computed for this same orbit and at the same time in the TD.  It is 
used to obtain the MPs in Lorenz gauge using Eq.~(\ref{eq:push}).
Fig.~\ref{fig:xiOdd} shows the $\ell=2$, $m=1$ amplitude of the gauge 
generator itself, which differs from $h_{2}$ in Lorenz gauge only by a factor 
of $-2$.  Fig.~\ref{fig:hthrL} shows the Lorenz gauge metric amplitudes 
$h^{\text{L}}_{t}$ and $h^{\text{L}}_{r}$ for the same mode.  The MPs are 
now $C^0$ at $r=r_p(t)$ and are asymptotically flat.

\begin{figure}[h!]
\includegraphics[scale=1]{h_t_and_f_h_r_L}
\caption{
\label{fig:hthrL}
The $\ell = 2, m = 1$ mode of
the Lorenz gauge MP amplitudes $h^{\text{L}}_{t}$
and $h^{\text{L}}_{r}$.
Dotted lines indicate the region of libration.
Orbital parameters are given in the text.
Note (comparing to Fig.~\ref{fig:hthrRW}) 
the discontinuity at the location of the particle has vanished 
and the wave no longer grows asymptotically.
We plot $f h^{\text{L}}_{r}$ so we can see the wave behavior
near the horizon.}
\end{figure}

Of key importance to our method is the exponential convergence of the TD 
solutions.  We can first consider self-convergence of the modes for all 
$r$.  As an example we choose an orbit with $e = 0.188917$ and 
$p = 7.50478$ at time $t = 96.44$.  In Fig.~\ref{fig:convVsX} we show the 
self-convergence of $\xi_{o}(r)$ for a set of partial Fourier sums over 
$n$ from $-N \le n \le +N$ for various $N$.  The right panel of this figure 
shows exponential self-convergence of the EPS method as a function of $r$, 
including at the particle.  This result is in contrast with the left panel 
which shows that the standard method is only algebraically convergent in 
the source libration region.  Note that the convergence is initially 
exponential before becoming algebraic around an error level of $10^{-4}$.  
This transition is due to the equation for $\xi_{o}$ having singular and
extended parts (see Eqs.~(\ref{eq:xiSing}) and (\ref{eq:xiExtended})).  
We find the singular part using EHS, which converges exponentially.  This 
part of the solution dominates the self-convergence in the left panel at 
first.  Eventually, the lack of differentiability of the extended source 
and the use (for comparison) of the standard Fourier series for that part 
of the field manifests itself.  The appearance of Gibbs behavior stalls 
the convergence in the libration region.

\begin{figure}
\includegraphics[scale=1]{convVsX.pdf}
\caption{
\label{fig:convVsX}
Self convergence of the $\ell = 2$, $m = 1$ mode of $\xi_{o}$.  We show 
results for the standard method in the left panel and our (equivalent)
partial annihilators and EPS on the right.  The orbital parameters are 
$e = 0.188917$ and $p = 7.50478$, and the fields are calculated at 
time $t = 96.44$.}
\end{figure}

Beyond self-convergence, we can check absolute convergence to the 
analytically-known jump conditions.  This test can be applied to both the 
MP amplitudes and their first radial derivatives.  As shown in the left 
panel of Fig.~\ref{fig:convVsN}, we find exponential convergence to the 
analytically computed values in Eq.~(\ref{eq:LGJumps}), in this case using 
the partial annihilator method.  Here the orbit is the more eccentric one 
with $e = 0.764124$ and $p = 8.75455$.  Each partial Fourier sum ranges 
over all harmonics from $-N \le n \le +N$ for different values of $N$ as 
seen on the horizontal axis.  The jump conditions are time dependent and 
thus we compare our results at several moments in time (in this case at 
20 points) throughout the orbit.  The left panel plots the maximum error 
encountered in each quantity throughout an orbit.  Since the Lorenz gauge 
amplitudes are all $C^{0}$, we plot absolute convergence for the jumps in 
the amplitudes themselves (which are expected to converge to 0) and relative 
convergence for the jumps in the $r$ derivatives of the amplitudes.  The 
convergence appears to bottom out around $10^{-12}$ to $10^{-11}$.

In the right panel of Fig.~\ref{fig:convVsN} we compare the accuracy of 
the EPS and PA methods.  For the jumps in Lorenz gauge MP amplitudes and 
their radial derivatives we show the relative error between the two methods 
as a function of time throughout one orbit.  The same high eccentricity 
orbit is used, though to compare the two methods the partial Fourier 
sums were fixed and taken to range over $-85 \le n \le 106$.  The two 
methods agree with each other to the level of $10^{-12}$ to $10^{-10}$.

\begin{figure}
\includegraphics[scale=1]{compAndConv.pdf}
\caption{
\label{fig:convVsN}
Convergence of the jump conditions as a measure of solution error.  In the 
left panel the partial annihilator method was used to compute the gauge 
generator and MPs for a high eccentricity orbit with $e = 0.764124$ and 
$p = 8.75455$.  Partial Fourier sums over $n$ are computed with 
$-N \le n \le +N$ and for various $N$.  Exponential convergence is exhibited 
in the various sums.  See text for further discussion.  In the right panel 
a comparison is made of discrepancies between the EPS method and the PA 
method as a function of time about the orbit.  See further discussion in the 
text.}
\end{figure}

Although we have only displayed results in this paper for the 
$\ell = 2, m = 1$ mode, we have run the code on many different modes and 
for different orbits.  We have no difficulty in computing the gauge 
transformation from RW to Lorenz for odd-parity modes with high accuracy.  
It now remains for us to apply these methods to the even-parity part of 
the gauge transformation, a somewhat more involved procedure.  We will 
turn to that issue in a subsequent paper.  

\section{Conclusion}

This paper is the first of two on the transformation of metric perturbations 
from Regge-Wheeler gauge to Lorenz gauge.  This first paper was confined 
to treating the odd-parity part of the MPs and devoted much of the 
discussion to the development of two new analytic/numerical methods for 
using frequency domain methods to find accurate solutions in the time 
domain.  The follow-on paper will be primarily devoted to discussing the 
analytic problem of finding the even-parity part of the gauge transformation, 
and will draw upon the numerical methods which we have detailed here.

\acknowledgments

The authors thank Chad Galley and Ian Hinder for helpful discussions.  We 
also appreciate suggestions made by the referee.  CRE acknowledges support 
from the Bahnson Fund at the University of North Carolina--Chapel Hill.

