text:
# Delay Equations and Radiation Damping

## 1 Introduction

In the treatment of the motion of *extended* bodies in classical field theory, the derivation of radiation reaction forces is based upon certain expansions of the retarded field potentials in powers of the retardation . The resulting local equations of motion involve derivatives of the acceleration and generally suffer from the existence of unphysical runaway solutions. Under certain model circumstances, we trace the origin of these problems to the expansion of the functions of the retarded arguments resulting in the replacement of the original nonlocal delay equations of motion by local higher-derivative equations that exhibit runaway solutions. In this general context, the properties of the delay equations that appear in classical field theory were first studied by L. Bel . Although our approach is rather general, for the sake of concreteness we discuss physical situations involving only the gravitational interaction.

Consider, for instance, inspiraling compact binaries that are expected to be promising sources of gravitational radiation. For a binary that is comprised of two compact objects—neutron stars or black holes—with, say, approximately equal masses $`m`$ and $`m^{}`$ in nearly circular orbits about each other, the relative orbital radius decays because orbital energy is emitted in the form of gravitational radiation. The dynamics of a usual binary system can be adequately described using the post-Newtonian approximation scheme that is valid in case the gravitational field is everywhere ‘weak’ and the motion is slow, that is $`vc`$, where $`v`$ is the characteristic orbital speed and $`c`$ is the speed of light. Although Einstein’s equations have a hyperbolic character associated with the retarded gravitational interaction, the standard post-Newtonian approximation scheme of general relativity deals with functions of instantaneous coordinate time $`t`$ rather than the retarded time $`t_r=tr/c`$, where (for the binary system) $`r`$ is the effective distance between the bodies (approximately the relative orbital radius). The gravitational potentials, which are originally functions of the retarded time, are expanded in Taylor series about $`t`$ using the effective small parameter in this expansion that can be written as $`\omega _br/c=v/c`$, where $`\omega _b=2\pi /P_b`$ is the relative orbital frequency and $`P_b`$ is the binary period. Because the gravitational waves emitted by the binary have an effective frequency of $`2\omega _b`$ and wavelength $`\lambda _bcP_b/2`$, the small parameter in the expansion can be reduced to the ratio $`\pi r/\lambda _b`$. Due to the observational fact that in typical astronomical systems $`v/c1`$, the first few terms of such an expansion can be used to derive the post-Newtonian equations of motion that describe, for instance, the orbital evolution of the binary pulsars discovered by Hulse and Taylor . The post-Newtonian equations of motion of binary stars are similar to the Abraham-Lorentz form in electrodynamics  but, because of the tensorial character of the gravitational field, these equations involve not only the third, but the fourth and fifth derivatives of the stars’ positions with respect to time as well. Schematically, the equation of the relative orbital motion reads

$$\ddot{𝐫}=𝐅_0(𝐫)+c^2𝐅_2(𝐫,\dot{𝐫},\ddot{𝐫})+c^4𝐅_4(𝐫,\dot{𝐫},\ddot{𝐫},𝐫^{(3)},𝐫^{(4)})+c^5𝐅_5(𝐫,\dot{𝐫},\ddot{𝐫},𝐫^{(3)},𝐫^{(4)},𝐫^{(5)}),$$
(1)
where $`𝐫`$ is the radius vector connecting the stars, the overdot denotes differentiation with respect to time, $`𝐫^{(n)}:=d^n𝐫/dt^n`$ and per reduced mass $`𝐅_0(𝐫)=(G(m+m^{})/r^3)𝐫`$ is the Newtonian force, $`𝐅_2`$ is the post-Newtonian force, $`𝐅_4`$ is the post-post-Newtonian force and $`𝐅_5`$ is the gravitational radiation reaction force responsible for the decay of the orbital period ($`\dot{P}_b<0`$) associated with the emission of gravitational waves by the binary system. In the quadrupole approximation under consideration here, the gravitational waves carry away energy and angular momentum, but not linear momentum ; therefore, the total momentum of the binary system is conserved and this fact is responsible for the absence of a force term proportional to $`c^3`$ in (1). Moreover, all tidal, spin-orbit, and spin-spin interactions are neglected in (1); the only parameters in equation (1) are the masses and the separation between the centers of mass of the members of the binary system. We mention that relativistic hydrodynamical (Euler) equations similar to system (1) have been derived to describe the motion of the fluid elements of the stars .

The higher time-derivative equations of the form (1) cannot be used directly to predict the dynamical evolution of a physical system because of the existence of so-called runaway modes that have been much discussed in the literature on electrodynamics but not in connection with astrophysical problems involving gravitational radiation reaction and the calculation of templates of the gravitational waves emitted by coalescing binaries . In analogy with electrodynamics, the existence of these runaway modes suggests that the truncated equations of the form (1) may not correctly predict the qualitative behavior of the solutions of the original true dynamical delay-type equations  that involve the retarded time $`t_r=tr/c`$. Moreover, the existence of runaway modes can cause serious difficulties for numerical integration in addition to the problems associated with the inaccuracies inherent in the approximation of higher-order derivatives by finite differences .

In case the post-Newtonian expansion parameter $`r/\lambda _b`$ is sufficiently small, we will provide in the following section a theoretical basis for eliminating the runaway solutions by replacing system (1) with a new model that is a system of second-order ordinary differential equations. Within this theory, high-order vector differential equations like equation (1) are not the desired approximate equations of motion, and they should not be used for numerical integration. Rather, system (1) must be viewed merely as an intermediate step in the derivation of the physically correct, second-order model equation with no runaway solutions that faithfully approximates the dynamics of the underlying delay-type equation. For illustration purposes, we apply this approach in section 3 to a discussion of one-dimensional gravitational dynamics of a two-body system. As expected, the reduced model predicts the correct motion of the binary except possibly when $`r/\lambda _b`$ is not small and residual terms that have been neglected in equation (1) start to play a significant role. Indeed, for an inspiraling binary system, the effective delay $`\omega _br/c`$ increases to some noticeable finite value as the system approaches coalescence. Motivated by this physical scenario, we introduce a simple model involving variable delay in section 4 that can be expressed as a Duffing-type differential-delay equation. This model is then analyzed to show some specific behavior of this delay equation that is not predicted by expansion in powers of the delay. Finally, section 5 contains a discussion of our results.

## 2 Delay equations with small delays

The delay-type equations of motion with retarded arguments are usually too complicated for mathematical analysis; therefore, we limit our discussion in this section to equations with constant delays. Although this is an unrealistic restriction in general, we note that for an astrophysical binary system consisting of compact point-like neutron stars or black holes moving around each other along circular orbits, the delay is almost constant. In fact, a close approximation to this delay is the ratio $`r/c`$, where $`r`$, the radius of the relative orbit, is changing very slowly due to the emission of gravitational energy in the form of gravitational waves.

Taking into account the last remark, let us consider a family of delay differential equations of the form

$$\dot{x}(t)=F(x(t\tau ),x(t)),$$
(2)
where $`\tau `$ is viewed as a real dimensionless parameter and $`x`$ is a variable in $`^n`$; intuitively, the constant delay $`\tau `$ corresponds in effect to $`\omega _br/c`$. The members of this family are examples of a more general and widely studied class called retarded functional differential equations (see ).

Using the delay equation (2) as an abstraction for the retarded-time model that is supposed to be approximated by a system of the form (1), we will discuss an approach for extracting the ‘correct’ dynamical equations of motion from system (1) that eliminates the runaway solutions.

Our approach assumes the existence of an attractor for the underlying delay-type equation. We will rely on the work of Bel  for (numerical) evidence in favor of the existence of attractors in the retarded equations of motion with space-dependent delays that appear in electrodynamics; but, we know of no mathematical proof for the existence of attractors for these equations or for the similar delay-type equations of astrophysics. Indeed, the proof of the existence of attractors for delay-type equations with space-dependent delays remains a challenging mathematical problem of physical significance. For the delay equation (2), however, if $`|\tau |`$ is sufficiently small, then the corresponding member of the family (2) has a global $`n`$-dimensional attractor such that the restriction of the delay equation to this attractor is equivalent to a first-order system $`𝒮_A`$ of ordinary differential equations. We will eventually outline a proof of this result. But, let us first discuss our approach to eliminating the runaway solutions.

The solutions of the delay equation (2) approach the attractor exponentially fast; therefore, the system $`𝒮_A`$ on this attractor determines (asymptotically) the true dynamical behavior of the system, hence we consider it to be the ‘correct’ physical model. On the other hand, it is easy to see that if equation (2) is expanded in the small parameter $`\tau `$ and truncated at some order $`N`$, then an $`N`$th-order ordinary differential system $`𝒮_N`$ akin to system (1) is obtained such that the coefficient of the $`N`$th-order time derivative of $`x`$ contains the factor $`\tau ^N`$ and is therefore singular in the limit as $`\tau 0`$.

For $`\tau >0`$ and sufficiently small, the high-order differential equation $`𝒮_N`$ has an equivalent first-order system $`𝒮`$ that has an $`n`$-dimensional (invariant) slow manifold. Moreover, this slow manifold has corresponding stable and unstable manifolds; in effect, the first-order system has (physical) solutions that are asymptotically attracted to the slow manifold and (unphysical or runaway) solutions that are asymptotically repelled from the slow manifold. The restriction of the first-order singularly perturbed system of differential equations to its slow manifold is of course an $`n`$-dimensional first-order system of ordinary differential equations $`𝒮_S`$ on this $`n`$-dimensional manifold. Our main result states that *in appropriate local coordinates, the system $`𝒮_S`$ on the slow manifold agrees to order $`N`$ in $`\tau `$ with the first-order system $`𝒮_A`$ on the global attractor of the underlying delay differential equation; therefore, the system $`𝒮_S`$, which can be obtained directly from the high-order differential equation $`𝒮_N`$, is a faithful approximation of the ‘correct’ physical model.* Generalizing to the Abraham-Lorentz type equation (1) (analogous to $`𝒮_N`$), the correct physical model is obtained as the system of ordinary differential equations (analogous to $`𝒮_S`$) that determines the motion on the slow manifold of a corresponding first-order system that is viewed as being singularly perturbed relative to the small parameter $`r/\lambda _b`$.

While there is evidence that the mathematical assertions in the scenario just proposed are valid, some of these assertions have not yet been rigorously justified in full generality, even for the case of fixed delays. In the remainder of this section we will provide some evidence, in the case of fixed delays, for the existence of a global attractor and for the claim that the dynamical system on this attractor is well approximated by the dynamical systems on the slow manifolds of singularly perturbed first-order systems obtained by truncations of the expansion of the delay equation in powers of the delay.

Our approach for the elimination of runaway solutions is equivalent to the procedure of iterative reduction (also called order reduction) that is often used to eliminate runaway solutions by means of the evaluation of the higher time-derivative terms in equations like (1) by the repeated substitution of the equations of motion and the subsequent reduction of the resulting equation to one of the second order (cf. ). Thus, our approach provides a theoretical framework for the rigorous justification of iterative reduction (cf. ), a procedure that has been justified so far by physical intuition.

For *conservative* higher time-derivative systems, the order reduction procedure has been investigated within the frameworks of Lagrangian and Hamiltonian dynamics by a number of authors (see  and the references cited therein). In particular, it can be shown that under suitable conditions—relevant, for instance, to Euler-Lagrange equations analogous to equation (1) truncated at the fourth order—higher-derivative Lagrangians can be iteratively reduced by redefinitions of position variables .

Returning to the delay equation (2), we note that it has an infinite-dimensional state space of initial conditions. For example, if the delay $`\tau `$ is a fixed positive number, then the natural state space of initial conditions is the infinite-dimensional vector space of continuous $`^n`$-valued functions on the interval $`[\tau ,0]`$. This space endowed with the supremum norm is a Banach space that we will denote by $`𝒞`$. Note that for an arbitrary continuous $`^n`$-valued function $`\gamma `$ defined on the interval $`[t\tau ,t]`$, the function $`\gamma _t`$ given by $`\gamma _t(\vartheta )=\gamma (t+\vartheta )`$ is in $`𝒞`$. Under the assumption that $`F`$ is a smooth function and $`\varphi 𝒞`$, there is a unique continuous solution $`y`$ of the corresponding delay equation in the family (2) such that $`y`$ is uniquely defined for $`t\tau `$ and $`y_0=\varphi `$ (see, for example, ). The state of the system at time $`t>0`$ is defined to be the function $`y_t`$ in $`𝒞`$.

To see that there is an attractor for the family (2) in case $`\tau `$ is sufficiently small, it is convenient to introduce the fast time $`s:=t/\tau `$, valid for $`\tau 0`$, so that with $`y(s):=x(t)`$ the family (2) takes the form

$$y^{}(s)=\tau F(y(s1),y(s))$$
(3)
and each member of this family, parametrized by $`\tau `$, has the same state space—the continuous functions on the interval $`[1,0]`$. For each $`\tau 0`$ the delay equation (3) is equivalent to the corresponding member of the family (2). For $`\tau =0`$ the corresponding differential equations are not equivalent, but this is of no consequence because we are only interested in the solutions of the family (2) for $`\tau 0`$. By viewing the unperturbed system (3), namely $`y^{}(s)=0`$, as a delay equation with unit delay, it is clear that the solution with initial state $`\varphi `$ is given by $`y=\varphi `$ on the interval $`[1,0]`$ and by the constant $`y=\varphi (0)`$ for $`t0`$. The initial state in $`𝒞`$ thus evolves at time $`t=1`$ to its final constant state, the function defined on the interval $`[1,0]`$ with the constant value $`\varphi (0)`$. Thus, we conclude that the $`n`$-dimensional space of constant functions on $`[1,0]`$ is a global attractor for the delay equation $`y^{}(s)=0`$. Moreover, the convergence to this attractor is faster than any exponential (the solution reaches the attractor in finite time), and the dynamical system on this attractor is given by the ordinary differential equation $`y^{}(s)=0`$. If $`F`$ is appropriately bounded and $`\tau `$ is sufficiently small, then, because the contraction rate to the attractor is exponentially fast, the attractor persists in the family (3) in analogy with the persistence of attractors in finite-dimensional dynamical systems; in fact, each corresponding member of the family (3) has an $`n`$-dimensional attractor in the state space $`𝒞`$ and the restriction of the dynamical system to this attractor is an ordinary differential equation. In particular, the family (3) has a corresponding family of invariant manifolds (that is, manifolds consisting of a union of solutions) that depend smoothly on the parameter $`\tau `$.

To identify the dynamical system on an attractor of a delay equation, let us suppose that the delay equation (2) has a family of $`n`$-dimensional invariant manifolds parametrized by $`\tau `$. Moreover, let $`\xi `$ denote the local coordinate on these invariant manifolds, and let $`x(t,\xi ,\tau )`$ denote the solution with the initial condition $`x(0,\xi ,\tau )=\xi `$ on the invariant manifold corresponding to the parameter value $`\tau `$. Because these solutions satisfy the delay equation (2), we have that

$$\dot{x}(t,\xi ,\tau )=F(x(t\tau ,\xi ,\tau ),x(t,\xi ,\tau ));$$
(4)
therefore, the generator of the dynamical system on the attractor is the vector field

$$X(\xi ,\tau ):=\frac{}{t}x(t,\xi ,\tau )|_{t=0}=F(x(\tau ,\xi ,\tau ),\xi ).$$
(5)
Under our assumption that this vector field is analytic in $`\tau `$, it can be expanded as a Taylor series about $`\tau =0`$ by differentiating the function $`F(x(\tau ,\xi ,\tau ),\xi )`$ with respect to $`\tau `$. To this end, we note that the partial derivatives of $`x(\tau ,\xi ,\tau )`$ with respect to its first argument can be evaluated using equation (4), and partial derivatives with respect to its third argument vanish at $`\tau =0`$ since $`x(0,\xi ,\tau )=\xi `$; moreover, its mixed partial derivatives can be evaluated by differentiation of equation (4) with respect to $`\tau `$.

As a concrete and instructive example of the construction of the dynamical system on an attractor, let us consider a simple case of equation (2) by replacing it with $`\dot{x}(t)=\widehat{f}(x(t\tau ))`$ where the function $`\widehat{f}:`$ is the scalar linear function given by $`\widehat{f}(x)=ax`$ so that the associated family of delay equations is

$$\dot{x}(t)=ax(t\tau ).$$
(6)
In this case, there is a corresponding family of solutions given by

$$x(t,\xi ,\tau )=e^{\lambda (\tau )t}\xi ,$$
(7)
where $`\lambda (\tau )`$ is the unique *real* root of the equation $`\lambda =a\mathrm{exp}(\lambda \tau )`$, a fact that is easily checked by direct substitution of equation (7) into the delay equation (6).

The dynamical system on the invariant manifold is generated by the family of vector fields $`X(\xi ,\tau )=\widehat{f}(x(\tau ,\xi ,\tau ))=ax(\tau ,\xi ,\tau )=a\mathrm{exp}(\lambda (\tau )\tau )\xi =\lambda (\tau )\xi `$. By an application of the Lagrange inversion formula , the Taylor series expansion of $`X(\xi ,\tau )`$ about $`\tau =0`$ is

$$X(\xi ,\tau )=\xi \underset{n=1}{\overset{\mathrm{}}{}}\frac{(1)^{n1}n^{n1}a^n}{n!}\tau ^{n1},$$
(8)
and its radius of convergence is $`\tau ^{}=(|a|e)^1`$. The qualitative behavior of solutions of the system (6) for small $`\tau `$ is clear: all solutions are attracted to a one-dimensional invariant manifold on which the dynamical system is the linear ordinary differential equation $`\dot{x}=X(x,\tau )`$. For example, if $`a<0`$ and $`\tau `$ is sufficiently small, then all solutions are attracted to the trivial solution $`x(t)0`$.

Let us now turn to the standard approach in physics where an underlying delay equation is expanded in powers of the delay to obtain an ordinary differential equation of motion. To illustrate this, let us consider a special scalar case of the delay equation (2) given by

$$\dot{x}=f(x(t\tau ))+g(x(t)),$$
(9)
and let us suppose, in analogy with the true dynamical delay-type equations that might arise in theories of electromagnetism and gravitation, that the true equation of motion for some process is the delay equation (9). The result of expanding equation (9) to order $`\tau ^2`$ is the second-order differential equation (an analogue of equation (1))

$$\dot{x}=f(x)+g(x)\tau f^{}(x)\dot{x}+\frac{\tau ^2}{2}[f^{\prime \prime }(x)\dot{x}^2+f^{}(x)\ddot{x}],$$
(10)
where a prime denotes differentiation with respect to $`x`$. Although we only write the second-order expansion, we note that the coefficient of the $`N`$th-order time derivative of $`x`$ in the $`N`$th-order expansion is $`(\tau ^N/N!)f^{}(x)`$. Hence, the corresponding $`N`$th-order ordinary differential equation is singular in the limit as $`\tau 0`$. Also, if we assume that system (9) has a (smooth) family of attractors parametrized by $`\tau `$, then the corresponding family of vector fields generating the dynamical systems on these attractors is given to second order in $`\tau `$ by

$`X(x,\tau )`$ $`=`$ $`f(x)+g(x)\tau f^{}(x)(f(x)+g(x))`$ (11)
$`+{\displaystyle \frac{1}{2}}\tau ^2\{f^{\prime \prime }(x)[f(x)+g(x)]^2`$
$`+f^{}(x)[3f^{}(x)+g^{}(x)][f(x)+g(x)]\}+O(\tau ^3).`$
The ‘correct’ model (that is, the dynamical system on the attractor in the original delay equation) can be obtained by treating the expanded and truncated system akin to system (10) as a singular perturbation problem, which can be analyzed using Fenichel’s geometric theory of singular perturbations . A basic result of this theory states that if an $`N`$th-order singular perturbation problem with small parameter $`\tau `$ is recast as a first-order (‘fast’) system and the corresponding unperturbed system has an invariant manifold that satisfies certain conditions (normal hyperbolicity), then for sufficiently small $`\tau `$ each member of the family of perturbed first-order systems has an invariant slow manifold. The dynamical system on this slow manifold for the perturbed first-order family obtained from the $`N`$th-order truncation of the delay equation is the desired faithful approximation to the correct model. For example, let us recast the second-order ordinary differential equation (10) as the first-order singular perturbation problem

$`\dot{x}`$ $`=`$ $`u,`$
$`\tau ^2\dot{u}`$ $`=`$ $`(f^{}(x))^1\left[2(1+\tau f^{}(x))u2f(x)2g(x)\tau ^2f^{\prime \prime }(x)u^2\right].`$ (12)
Using Fenichel’s theory, it is easy to show that each member of this family, corresponding to a sufficiently small value of $`|\tau |`$, has a slow manifold. Also, it is possible to prove that the family of vector fields on these manifolds agrees to order $`\tau ^2`$ with the family (11) of vector fields on the attractor in the state space of the underlying family of delay equations (9). We note that these results are valid for the vector case of delay equation (9) as well.

For the delay equation (9), and also for more general families of delay equations where the delay is viewed as a small parameter, we conjecture that the slow vector field, for an appropriately defined first-order system that is equivalent to the $`N`$th-order truncation of the expansion of the family in powers of the small delay, agrees to order $`N`$ with the vector field on the attractor in the state space of the original delay equation. We have just mentioned that this conjecture is true for the delay equation (9) in case $`N=2`$. It can be shown that the conjecture is true in general for the linear delay equation $`\dot{x}(t)=Ax(t\tau )`$, where $`x`$ is a variable in $`^n`$ and $`A`$ is a nonsingular $`n\times n`$ matrix .

As we have already discussed, singular equations of motion like system (2) generally have unphysical runaway solutions. To eliminate these solutions and leave only the physical solutions, the singular system must be replaced by the dynamical system on the corresponding slow manifold. In effect, the truncated equations obtained from the underlying delay equation after expansion in the small delay must be replaced by the system obtained using iterative reduction; this system is equivalent to the dynamical system on the slow manifold. Without this replacement, the appearance of spurious runaway modes in inevitable, and their existence will cause overflows in numerical simulations.

## 3 Gravitational radiation damping

To illustrate the singular perturbation procedure described in section 2 as a method for the elimination of runaway solutions, we examine an application of this approach to a one-dimensional Abraham-Lorentz equation of the form (1).

Let us consider an ideal linear quadrupole oscillator (that is, two masses $`m`$ and $`m^{}`$ connected by a spring of negligible mass), where the only source of damping is the gravitational radiation reaction force associated with the emission of gravitational radiation due to the variable quadrupole moment of the system. A model for the (dimensionless) relative position $`z`$ of these particles, with gravitational radiation damping included, is the fifth-order ordinary differential equation

$$\mu z\frac{d^5z^2}{dt^5}+\frac{d^2z}{dt^2}+z=1,$$
(13)
where the small parameter is given by $`\mu =4G\mu _0\mathrm{}_0^2\omega _0^3/(15c^5)`$, $`\mu _0`$ is the reduced mass ($`\mu _0^1=m^1+m^1`$), $`\mathrm{}_0`$ is the spatial scale parameter and $`\omega _0`$ is the frequency of the ideal linear oscillator such that $`\omega _0^1`$ is the temporal scale parameter. In equation (13), we have neglected the Newtonian gravitational interaction between $`m`$ and $`m^{}`$ as well as all relativistic effects except for radiation damping. Let us note that this oscillator has an equilibrium solution $`z(t)1`$. According to our general approach described in section 2, this system can be treated as a singular perturbation problem, where the physically correct dynamical system would be the system defined on the slow manifold of a corresponding, appropriately chosen, first-order system.

We emphasize that the fifth-order differential equation (13) is not the correct physical model; for example, to specify a solution, the initial position and its first four time derivatives must be given. Even with the obvious choice for these initial conditions—that is, the initial conditions for a sinusoidal oscillation—a numerical integration shows that such solutions do not oscillate; rather, they are divergent.

To obtain a system with the expected dynamical behavior of an under-damped oscillator, the differential equation (13) must be replaced by its restriction to an appropriate slow manifold. We will not carry out the complete reduction procedure here . We note, however, that the system matrix for the linearization of system (13) at the steady state solution $`z=1`$ has five distinct eigenvalues that are given to lowest order in the small parameter by

$$(2\mu )^{1/3},(2\mu )^{1/3}(\frac{1}{2}\pm i\frac{\sqrt{3}}{2}),\mu \pm i.$$

For small $`\mu `$, the first three eigenvalues are ‘fast’ and the last two are ‘slow’. This suggests that the nonlinear system has a two-dimensional slow manifold. In fact, in accordance with our general scheme, the restriction of the dynamical system to this invariant manifold is a second-order system that gives the correct post-Newtonian dynamics. In this case, the dynamical system on the slow manifold to first order in the small parameter is given by the second-order differential equation

$$\ddot{z}+32\mu z(z\frac{15}{16})\dot{z}+z=1.$$
(14)
A unique solution of this equation is obtained by specifying only the initial relative position and velocity of the oscillating masses. For $`z`$ near the equilibrium state $`z=1`$, the expected dynamics for the radiating system is revealed: the relative motion is an under-damped oscillator. Numerical integration of this equation using standard algorithms is stable and produces the expected result.

The iterative reduction procedure can also be used to obtain equation (14) from equation (13). Even our simple example illustrates the necessity of reducing the higher-order equations of motion involving radiation reaction before numerical integration. For the more realistic hydrodynamic equations that include conservative post-Newtonian terms as well as radiation reaction, the corresponding Euler equation must involve these forces in the *reduced* form, that is, they should contain at most the position and the velocity of the fluid element .

## 4 Delay equations with sufficiently large delays

It is important to point out that the reduction procedure described in sections 2 and 3 cannot in general be expected to produce a good approximation to the true dynamics for ‘large’ delays.

As a simple but revealing example, let us reconsider the scalar linear delay equation $`\dot{x}(t)=ax(t\tau )`$ with $`a>0`$. For small $`|\tau |`$, we have already shown that all orbits in the state space are attracted to a one-dimensional attractor on which the dynamical system is given by the vector field (8) with $`aa`$. For $`|\tau |`$ less than the radius of convergence of this series $`\tau ^{}=(|a|e)^1`$, the correct dynamical behavior of the delay equation is predicted by this vector field. Because, in this case, the zeroth-order approximation $`\dot{x}=ax`$ already has a hyperbolic structure (that is, all solutions are attracted to the rest point at the origin exponentially fast), even the zeroth order approximation determines the qualitative dynamics for these values of $`\tau `$. By inspection of this delay equation, it might seem natural to conclude that the fixed delay $`\tau `$ does not influence the behavior for sufficiently large $`t`$ and the approximation $`\dot{x}=ax`$ remains valid for all fixed delays. This is not true. For instance, if $`\tau =\pi /(2a)`$, then the delay equation has the two-parameter family of exact solutions

$$tc_1\mathrm{cos}at+c_2\mathrm{sin}at.$$
(15)
Therefore, the qualitative behavior of the delay equation $`\dot{x}(t)=ax(t\pi /(2a))`$ is certainly not predicted by the ordinary differential equation $`\dot{x}=ax`$, or by the corrections to this equation within the radius of convergence of the slow vector field. The transition of the dynamical behavior of this delay equation from a stable rest point to a periodic regime as $`\tau `$ increases is easily seen to be the result of a degenerate Hopf bifurcation . Indeed, we recall that $`x(t)=\mathrm{exp}(\lambda (\tau )t)`$ is a solution of the delay equation under consideration if $`\lambda `$ is a solution of the characteristic equation $`\lambda =a\mathrm{exp}(\lambda \tau )`$. For $`\tau <\tau ^{}`$, the solutions of this equation have negative real parts and all such solutions are therefore attracted to the zero solution. If $`\tau =\pi /(2a)`$, then the characteristic equation has a pair of pure imaginary roots that give rise to the two-parameter family of periodic solutions (15). For $`\tau >\pi /(2a)`$, the characteristic equation has roots with positive real parts; therefore, there are solutions that grow without bound. Nevertheless, for these values of $`\tau `$, the delay equation has an attractor. In fact, for $`\pi /(2a)\tau <\pi /(2a)+2\pi /a`$, there is a two-dimensional attractor and the dynamical system on the attractor has the form $`\ddot{x}2\theta \dot{x}+(\theta ^2+\phi ^2)x=0`$ corresponding to the roots $`\lambda _\pm =\theta (\tau )\pm i\phi (\tau )`$ of the characteristic equation $`\lambda =a\mathrm{exp}(\lambda \tau )`$ with positive real parts. As $`\tau `$ increases further, the dimension of the attractor increases discontinuously by two at each $`\tau =\pi /(2a)+2N\pi /a`$, where $`N=2,3,4\mathrm{}`$.

The Hopf bifurcation for delay equations with constant delays has been studied in detail. For instance, a more sophisticated analysis (see, for example, \[25, p. 341\]) shows that $`\tau =\pi /2`$ is a supercritical Hopf bifurcation value for the nonlinear scalar delay equation

$$\dot{x}(t)=[1+x(t)]x(t\tau ).$$

Moreover, this system has a nontrivial periodic orbit for each $`\tau >\pi /2`$.

The delay-type equations of astrophysics generally do not have constant delays. But as we have mentioned, for two coalescing neutron stars with nearly equal masses and on nearly circular orbits, the delays involved are almost constant; in fact, this ‘fast’ periodic motion evolves as a result of radiation damping on a timescale that is much longer than $`P_b`$. During this ‘slow’ evolution, the delay increases as the radius of the binary decreases due to the emission of gravitational waves. Motivated by this astrophysical scenario, we have studied an oscillator model with a time-dependent delay. This example is not intended to be a realistic model, rather it is meant to illustrate some of the bifurcation phenomena that occur in delay equations with time-dependent delays. Our example is the second-order differential-delay equation

$$\ddot{x}(t)+\mathrm{\Omega }^2x(t)+\alpha x(tw)\beta x^3(tw)=0,$$

where $`\alpha `$, $`\beta `$ and $`\mathrm{\Omega }`$ are constant system parameters and $`x`$ is viewed as the state variable of a (Duffing) oscillator with variable delay $`w(t)`$ such that $`\dot{w}(t)+\rho w(t)=\rho \nu `$. Here $`\rho `$ and $`\nu `$ are constants; hence, $`w(t)\nu `$ is an exponentially decreasing or increasing function of time depending on whether $`\rho `$ is positive or negative, respectively. In any case, we have a dynamic delay that is asymptotic to the constant value $`\nu `$. Note that if $`\rho =0`$, the delay is constant; in this case, the corresponding second-order differential equation on the slow manifold (to first order in $`w`$) is given by

$$\ddot{x}+w(3\beta x^2\alpha )\dot{x}+(\alpha +\mathrm{\Omega }^2)x\beta x^3=0,$$
(16)
a form of van der Pol’s equation. In case $`\mathrm{\Omega }0`$, this differential equation typically has a stable limit cycle for $`w>0`$. But for $`\mathrm{\Omega }=0`$ (that is when all forces are retarded), it is easy to prove that no periodic orbits exist and most solutions are unbounded.

A typical plot of $`x`$ versus $`t`$ for system (4) for $`\mathrm{\Omega }^2\alpha >0`$ is given in figure 1, where the delay increases from its initial value $`w(0)=0`$ to $`\nu =10.5`$. The initial response of the system (where the delay is small) is characterized by an oscillation as expected from equation (16), which follows from the expansion of equation (4) to first order in $`w1`$. But as $`w`$ increases, the qualitative behavior of the system is affected by three additional bifurcations not accounted for by equation (16). At the third bifurcation, the stable oscillation disappears. Additional bifurcations of the same type occur if $`\nu `$ is set to a larger value. Numerical experiments suggest that these bifurcations are not Hopf bifurcations; instead, they are ‘center bifurcations’, where at some parameter value there is a rest point of center type and one of the periodic orbits surrounding this rest point continues to exist as the parameter is changed. The family (16) with the parameter values as in figure 1 has a bifurcation of this type as $`w`$ increases through $`w=0`$.

The behavior depicted in figure 1 is suggested by an analysis of the roots of the characteristic equation,

$$\zeta ^2+\mathrm{\Omega }^2+\alpha e^{\zeta w}=0,$$

for the linearization of the delay equation (4). The bifurcation points (corresponding to the existence of centers) are given by

$$w_{\mathrm{}}=\mathrm{}\pi (\mathrm{\Omega }^2+\alpha \mathrm{cos}\mathrm{}\pi )^{1/2},$$
(17)
where $`\mathrm{}`$ is a non-negative integer. These are the values of $`w`$ such that the characteristic equation has pure imaginary roots. A computation shows that if $`\mathrm{}`$ is even, then as $`w`$ increases a pair of pure imaginary roots crosses the imaginary axis into the right half-plane, and if $`\mathrm{}`$ is odd, then the roots cross into the left half-plane. Under the assumption that the bifurcations are supercritical, a stable limit cycle appears after the bifurcation in the first case; in the second case, a stable limit cycle disappears. For the parameter values used to obtain figure 1, the bifurcation values computed from equation (17) are (approximately) $`0`$, $`3.3`$, $`6.0`$, $`9.9`$ for $`\mathrm{}=0,1,2,3`$ such that $`w_{\mathrm{}}\nu `$. At $`w_0=0`$ a limit cycle appears, at $`w_13.3`$ the limit cycle disappears, and so on. Thus, these bifurcations account for the appearance and disappearance of oscillations in figure 1. We note that a similar sequence of bifurcations occurs whenever $`\mathrm{\Omega }^2\alpha >0`$. On the other hand, if $`\mathrm{\Omega }^2\alpha <0`$ (for example if $`\mathrm{\Omega }=0`$), then all bifurcation points correspond to roots of the characteristic equation crossing into the right half-plane. In this case, the bifurcations can be subcritical. Indeed, for $`\mathrm{\Omega }=0`$, numerical simulations indicate that no limit cycle appears. As a result, solutions starting near the unstable rest point become unbounded.

The slow dynamical system, obtained by reduction from a truncation of an expansion of a delay equation in powers of the delay, approximates the dynamics on the global attractor of the delay equation as long as the delay is sufficiently small; but, as our examples show, the ordinary differential equations obtained by expansion, truncation, and reduction *cannot* be used in general to predict the correct dynamical behavior for sufficiently large delays. We have mentioned, for example, that the dimension of the attractor of a family of delay equations, parametrized by the delay, can increase in dimension so that that the corresponding slow vector field is no longer defined on the attractor. But this is not the only possible scenario for the appearance of new attractor dynamics; for example, the attractor could cease to exist or be a manifold for some values of the delay.

The Abraham-Lorentz type equation (1) can be used, after reduction to a slow manifold, to predict the relative orbital motion of a relativistic binary system in the regime where the delay is sufficiently small. The size of the maximum allowed delay would have to be computed on a case-by-case basis using the explicit form of the delay equation that models the dynamics of a coalescing pair of neutron stars. The results of this section show that for sufficiently large delay the attractor does not in general correspond to the slow manifold. The question remains whether such a divergence of behaviors could ever occur in the case of retarded equations of classical field theory. This is an interesting open problem.

## 5 Discussion

It is expected that interferometric gravitational wave detectors that are presently under construction will be able to detect signals from massive coalescing binary systems. For the analysis of such forthcoming data, it is important to have theoretically predicted wave forms (‘templates’) for the relevant astrophysical processes. To this end, extensive computations are necessary that need to take gravitational radiation reaction into account . The standard approach leads to higher time-derivative equations that involve runaway modes and inevitably produce incorrect results.

We have determined the source of the difficulty by investigating delay equations, which are essentially nonlocal, and the higher time-derivative equations that are obtained by truncations of the (post-Newtonian) expansions in powers of the delay. For sufficiently small delays, a proper justification is provided for the usual method of replacing terms with higher-derivatives by terms with at most first derivatives using repeated substitution of the equations of motion (‘iterative reduction’). We have shown that in the investigation of the solutions of higher-derivative equations that represent phenomena involving radiation reaction, it is essential to reduce such equations to the corresponding slow manifolds before numerical analysis. Our work suggests, however, that unexpected nonlocal phenomena could occur for sufficiently large delays that cannot be predicted using the local equations of motion even after iterative reduction.

## References


format: text

text:
# Untitled Document

For the Quantum Heisenberg Ferromagnet, Some Conjectured Approximations

Paul Federbush

Department of Mathematics

University of Michigan

Ann Arbor, MI 48109-1109

(pfed@math.lsa.umich.edu)

Abstract

We present some conjectured approximations for spin expectations in a Quantum Heisenberg system. The conjectures are based on numerical experimentation, some theoretical insights and underpinning, and aesthetic value. We hope theoretical developments will follow from these ideas, even leading to a proof of the phase transition (in three dimensions).

We organize this paper into three sections. The first presents preliminary definitions and the conjectures. The second contains a rigorous theoretical development of a useful framework for the system. The final section introduces an “average-field” structure that may lead to an understanding, and hopefully a proof, of the conjectures.

I) The Conjectures.

We consider a lattice, $`\mathrm{\Lambda }`$, and the associated Quantum Heisenberg Hamiltonian

$$H=\underset{ij}{}(I_{ij}1)$$
(1)
where $`I_{ij}`$ interchanges the spins of the two neighboring sites $`i`$ and $`j`$ in the lattice $`\mathrm{\Lambda }`$. We let $`p_i`$ be the projection onto spin up at site $`i`$.

$$p_i=\frac{1}{2}\left(\sigma _{zi}+1\right).$$
(2)
We consider a state $`\psi _0`$ with spin up at sites in $`𝒮_0`$, and spin down at the complementary sites.

$$p_i\psi _0=\{\begin{array}{cc}\psi _0,\hfill & i𝒮_0\hfill \\ & \\ 0,\hfill & i𝒮_0\hfill \end{array}$$
(3)
We assume there are $`N`$ spin ups,

$$\mathrm{\#}\{𝒮_0\}=N.$$
(4)
We let $`\varphi _\mu (i)`$ be a solution of the lattice heat equation

$$\frac{}{\mu }\varphi _\mu (i)=(\mathrm{\Delta }\varphi _\mu )(i)$$
(5)
with initial conditions

$$\varphi _0(i)=\{\begin{array}{cc}1,\hfill & i𝒮_0\hfill \\ & \\ 0,\hfill & i𝒮_0\hfill \end{array}$$
(6)
We define

$$\rho _\mu (i)=\frac{\varphi _\mu ^2(i)}{\varphi _\mu ^2(i)+(1\varphi _\mu (i))^2}$$
(7)
and

$$<p_i>_\mu =\frac{e^{\mu H}\psi _0,p_ie^{\mu H}\psi _0}{e^{\mu H}\psi _0,e^{\mu H}\psi _0}.$$
(8)
Conjecture 1.

$$|<p_i>_\mu \rho _\mu (i)|<c_d<1.$$
(9)
Conjecture 2.

$$|<p_i>_\mu \varphi _\mu (i)|<c_d<1.$$
(10)
In one-dimension the corresponding $`c_1`$ may be picked to be $`\frac{1}{2}`$, according to our numerical studies.

Conjecture 3.

$$\underset{\mu 0}{lim}\frac{1}{\mu }|<p_i>_\mu \rho _\mu (i)|=0$$
(11)
the limit taken in $`\mathrm{}^{\mathrm{}}(\mathrm{\Lambda })`$, and convergence independent of $`𝒮_0`$ and $`N`$. In one-dimension the $`\mu `$ in (11) may be replaced by $`\mu ^{2\epsilon }`$.

Contrary to our earlier expectations (as presented in a previous version of this note) the behavior of $`<p_i>_\mu `$ as $`\mu `$ becomes large is not simple. A better approximation than $`\rho _\mu (i)`$ or $`\varphi _\mu (i)`$ when $`\mu `$ is large is realized in $`\varphi _{\frac{\mu }{2}}(i)`$. We present in our next conjecture the result of our numeric study:

Conjecture 4.

For $`\mu 4`$ one has

$$|<p_i>_\mu \varphi _{\frac{\mu }{2}}(i)|<.1.$$
(12)
We have been specific with numbers in (12) to give the flavor of the estimate’s quality. We have some theoretical understanding of the reason $`\varphi _{\frac{\mu }{2}}(i)`$ is a good approximation to $`<p_i>_\mu `$, but we do not discuss it in this note, restricting our attention to $`\rho _\mu (i)`$ and $`\varphi _\mu (i)`$ in later sections.

Implicit in all these estimates is a locality property of $`<p_i>_\mu `$. We state a very weak form of this in the following conjecture.

Conjecture 5.

For any $`\epsilon >0`$, there is an $`L_{\epsilon ,\mu }`$, such that $`<p_i>_\mu `$ is determined within $`\epsilon `$ by knowledge of the spin configuration (as specified in (3)) in a region within distance $`L_{\epsilon ,\mu }`$ of site $`i`$, i.e. changing spins outside this distance cannot effect $`<p_i>_\mu `$ by more than $`\epsilon `$

Conjecture 5 is the most basic of our assertions, and should fit into some very general theoretical framework.

Cave Adfirmationes: Most of the numerical investigation was on a one-dimensional lattice. (But also in periodic two-dimensional sets and three-dimensional sets.)

II) Some Simple Theory.

The Hilbert space of the system is naturally viewed as a direct sum

$$=_0_1_2\mathrm{}_{|\mathrm{\Lambda }|}$$

where in $`_N`$ there are $`N`$ spin ups. We write $`H_N`$ for $`H`$ restricted to $`_N`$. The space $`_N`$ is an invariant subspace of $`H`$, the set of $`N`$ spin waves.

We let $`Q`$ be an operator interchanging spin up and spin down, $`Q`$ a unitary operator commuting with $`H`$. $`Q`$ interchanges $`H_N`$ and $`H_{|\mathrm{\Lambda }|N}`$ as follows

$$Q\left(\underset{i𝒮}{}\left(\begin{array}{c}1\\ 0\end{array}\right)_i\underset{j𝒮}{}\left(\begin{array}{c}0\\ 1\end{array}\right)_j\right)=\underset{i𝒮}{}\left(\begin{array}{c}0\\ 1\end{array}\right)_i\underset{j𝒮}{}\left(\begin{array}{c}1\\ 0\end{array}\right)_j.$$
(13)
Vectors in $`_N`$ are described by symmetric functions on $`N`$ distinct lattice sites. $`f=f(.,..,.)`$ is associated to vectors in $`_N`$ as follows

$$f\underset{i_1,\mathrm{},i_N}{}f(i_1,\mathrm{},i_N)\underset{i\{i_i,\mathrm{},i_N\}}{}\left(\begin{array}{c}1\\ 0\end{array}\right)_i\underset{j\{i_i,\mathrm{},i_N\}}{}\left(\begin{array}{c}0\\ 1\end{array}\right)_j.$$
(14)
The sum in (14) is over distinct indices.

For $`N>M`$ there is a linear map from $`_N`$ to $`_M`$ called $`P_{N,M}`$. Let $`f`$ be in $`_N`$, then we define $`P_{N,M}f`$ in $`_M`$ by

$$(P_{N,M}f)(i_1,\mathrm{},i_M)=\underset{i_{M+1},\mathrm{},i_N}{}f(i_1,\mathrm{},i_M,i_{M+1},\mathrm{},i_N)$$
(15)
$`P_{N,M}`$ commutes with $`H`$ and interlaces $`H_N`$ and $`H_M`$

$$P_{N,M}H_N=H_MP_{N,M}.$$
(16)
If $`2N|\mathrm{\Lambda }|`$ then it is easy to show $`P_{N,M}`$ is onto. The preceding structure is related to the invariance of the system under global rotations.

For $`f`$ in any $`_N`$ we define

$$f_\mu e^{\mu H}f.$$
(17)
Of course

$$e^{\mu H_M}P_{N,M}=P_{N,M}e^{\mu H_N}.$$
(18)
In $`_1`$ , $`f_\mu `$ satisfies the heat equation

$$\frac{}{\mu }f_\mu =H_1f_\mu =\mathrm{\Delta }f_\mu .$$
(19)
So for $`f`$ in $`_N`$ we note the amusing fact that $`P_{N,1}f_\mu `$ satisfies the heat equation.

III) An Average-Field Approximation.

Let $`\psi _0`$ be a state in $`_N`$, sharp in the spins, with spin up at $`i`$ if $`i𝒮_0`$, spin down if $`i𝒮_0`$, $`𝒮_0`$ a set of $`N`$ sites.

$$\psi _0\underset{i𝒮_0}{}\left(\begin{array}{c}1\\ 0\end{array}\right)_i\underset{j𝒮_0}{}\left(\begin{array}{c}0\\ 1\end{array}\right)_j$$
(20)
and

$$\psi _\mu =e^{\mu H_N}\psi _0.$$
(21)
We define

$$\varphi _\mu =NP_{N,1}\psi _\mu $$
(22)
$`\varphi _\mu `$ in $`H_1`$, satisfies the heat equation and

$$\varphi _0(i)=\{\begin{array}{cc}1,\hfill & i𝒮_0\hfill \\ & \\ 0,\hfill & i𝒮_0\hfill \end{array}$$
(23)
We introduce an “average-field”-like “approximation” to $`\psi _\mu `$.

$$\psi _\mu ^{AP}\underset{i}{}\left(\begin{array}{c}\varphi _\mu (i)\\ \\ 1\varphi _\mu (i)\end{array}\right)_i.$$
(24)
This approximation has two nice features.

The approximation is “invariant and $`Q`$”. That is, it is $`Q`$ of the approximation obtained starting with $`Q\psi _0`$ instead of $`\psi _0`$.

The approximation is not sharp in spin wave number. (It does not lie in a single $`_N`$.) But in a reasonable sense it projects using $`\{P_{N,1}\}`$ onto $`\varphi _\mu (i)`$ in $`_1`$, that does satisfy the heat equation.

We note that $`\rho _\mu (i)`$, from equation (7), is given by

$$\rho _\mu (i)=\frac{\psi _\mu ^{AP},p_i\psi _\mu ^{AP}}{\psi _\mu ^{AP},\psi _\mu ^{AP}}.$$
(25)
Thus our approximate wave function, the “average-field” function (24) yields the spin up probabilities of approximation 1, equation (9). The “average-field” wavefunction satisfies the equation

$$\frac{d}{d\mu }\psi _\mu ^{AP}=H\psi _\mu ^{AP}$$
(26)
in the limit of nearly constant $`\varphi _\mu (i)`$. We expect there is some truth to the average field wavefunction; perhaps enough, so that its study (objective genitive) leads towards a proof of the phase transition for magnetization. We note that $`\varphi _\mu (i)`$ and $`\rho _\mu (i)`$ differ by less than .16. Both $`\varphi _\mu (i)`$ and $`\rho _\mu (i)`$ are “invariant under $`Q`$” as approximations.