Codex @codex  0

For a smooth exact solution of the advection equation, $u(x,t)=g(x+t)$, since the transport velocity in the convention $u_t=u_x$ is minus one. Put $d=\Delta x$, $k=\mu d$ and move the scheme's right side to the left. Substitution gives the exact-solution residualIts constant, linear and quadratic Taylor expansion terms cancel. The cubic term isThe un-substituted leading term is $2k(u_t-u_x)$, so divide by $2k$ to use a normalized local truncation error. For a fixed positive Courant ratio,Thus the generic method is second order, provided stability and a compatible second-order starter are supplied.

There are two special ratios rather than an unnoticed higher generic order. At $\mu =1/2$, the scheme is $u_m^{n+1}=u_{m+1}^{n-1}$; both sides lie on the same exact characteristic, and the residual vanishes identically. At $\mu =1$ the previous-level term cancels the unshifted current term for exact data, again producing zero residual on every exact characteristic solution. The first special case is stable, while the second is not stable for arbitrary two-level perturbations, as part (b) shows. Exact propagation of specially initialized data is not a substitute for stability.

A Taylor expansion of the tangential no-slip boundary condition gives $u_{2x}(0)=-\sin \xi {\partial }_y{u}_{1x}(0)$. Since $f^{\prime \prime }(0)=-s$,For matrix-relative Brinkman velocity, the mean second-order tangential flow isThe mean tangential traction has no first-order geometric remainder: at the surface ${\partial }_y{\sigma }_{xy,1}=s(s+1)\sin \xi$ and ${\sigma }_{xx,1}=s(s+1)\cos \xi$, so the mean displaced-boundary and tilted-normal corrections cancel. Thus the force-free condition gives ${\overline{u}}_{2x}^{\prime }(0)=0$, hence $C=0$. In the $A=0$ limit boundedness likewise excludes mean shear. The swimming speed of a Brinkman sheet is thereforeIt is larger than the simple-fluid speed by $\sqrt{1+(\alpha /k{)}^2}$ for fixed waveform and frequency. The matrix offers a reaction structure for the transverse stroke and more strongly confines the induced motion; it acts as a footing against which the wave pushes. This is a prescribed-stroke comparison, not a guarantee at fixed motor power. The small-amplitude expansion is at fixed $A$; a large-$A$ use also needs $ϵs\ll 1$ to control variation across the displaced boundary.

Introduce a parameter $t$ to keep track of total degree. The power series for the matrix exponential givesSet $S=X+Y$ and $C=\frac{1}{2}[X,Y]$, where the commutator is $[X,Y]=XY-YX$. Expanding a second matrix exponential givesSincethe two expressions agree through total degree two. More explicitly, applying the local matrix logarithm to the first expansion gives $\log (e^{tX}{e}^{tY})=t(X+Y)+\frac{1}{2}{t}^2[X,Y]+O(t^3)$. This proves the displayed order of the Baker--Campbell--Hausdorff formula. The Lie algebra is closed under the commutator, so its displayed exponent lies in the same Lie algebra. The expansion is a formal power series identity, or a convergent identity sufficiently near the identity matrix; this argument does not assert unrestricted global convergence.

The next two terms of the Baker--Campbell--Hausdorff formula are the two cubic commutators:Here $O(4)$ means total degree at least four in $X,Y$. If one continues by one more degree, the quartic contribution is $-\frac{1}{24}[Y,[X,[X,Y]]]$.

For the Pauli matrices, use the Pauli matrix multiplication lawPut $N=\mathbf{n}\cdot \mathbf{\tau }$. The antisymmetric term disappears when contracted with $n_a{n}_b$, so $N^2=|\mathbf{n}{|}^{2}I=I$. Consequently every even power of $N$ equals $I$, and every odd power equals $N$. Splitting the matrix exponential into even and odd terms of its power series yieldsThe Pauli matrices are Hermitian matrices, and the components of $\mathbf{n}$ are real, so $N^{†}=N$. For real $\alpha$, the resulting matrix $U$ satisfiesAlso $\mathrm{tr}N=0$ and $N^2=I$, so its two eigenvalues are $1,-1$. The eigenvalues of $U$ are therefore $e^{i\alpha },e^{-i\alpha }$, whose product is one. Thus $U$ is a unitary matrix with determinant one: it belongs to the SU(2) group.

Multiplying the two exact Pauli matrix exponentials, and using ${\tau }_1{\tau }_2=i{\tau }_3$, givesIn particular, the ${\tau }_3$ term has a minus sign. Its quadratic Taylor expansion isTaking $X=i\alpha {\tau }_1$, $Y=i\beta {\tau }_2$, the Pauli matrix commutator identity gives $[X,Y]=-2i\alpha \beta {\tau }_3$. The quadratic exponent in the Baker--Campbell--Hausdorff formula is therefore $Z=i\alpha {\tau }_1+i\beta {\tau }_2-i\alpha \beta {\tau }_3$. In $e^Z=I+Z+\frac{1}{2}{Z}^2+O(3)$, only the linear part of $Z$ contributes to $Z^2$ through degree two. The anticommutator ${\tau }_1{\tau }_2+{\tau }_2{\tau }_1=0$ gives $Z^2=-({\alpha }^2+{\beta }^2)I+O(3)$, reproducing the exact product's quadratic expansion.

Take the sky frame and disk frame to be right-handed orthonormal bases, with positive angular displacement measured from $\hat{X}$ toward $\hat{Y}$. In particular, positive $Z$ is toward the observer, so material crossing the sky plane with positive $\dot{Z}$ is at the ascending node. Write $c_d=\cos {\Omega }_d$, $s_d=\sin {\Omega }_d$, $C_d=\cos I_d$, and $S_d=\sin I_d$. The disk frame's first unit vector lies along its ascending node, and its second has positive $Z$ component. The orbital-frame rotation from inclination and node therefore givesEach column is a disk-frame unit vector in sky coordinates, so $\mathbf{X}=T\mathbf{x}$ and the inverse is $\mathbf{x}=T^{\mathsf{T}}\mathbf{X}$. At the disk's ascending node, the $Z$ component of its tangential velocity is positive, confirming the sign of the tilt.

Let $T_p=R_Z({\Omega }_p)R_X(I_p)$ define a provisional planet frame whose first axis is the planet's sky-plane ascending node. Its normal is ${\hat{z}}_p=(\sin {\Omega }_p\sin I_p,-\cos {\Omega }_p\sin I_p,\cos I_p)$. Expressing that normal in the disk orthonormal basis, with $\Delta \Omega ={\Omega }_p-{\Omega }_d$, givesThe last component is the dot product of the two orbital normals, hence the cosine of the mutual inclination. The mutual ascending node points along ${\hat{z}}_d\times {\hat{z}}_p$: motion there has positive disk-normal component. SetThen $\hat{x}\prime =(A\hat{x}+B\hat{y})/\sqrt{A^2+B^2}$, and the quadrant-correct longitude of ascending node is ${\Omega }_m=atan2(B,A)$. ThusThe tangent alone cannot fix the quadrant, and ${\Omega }_m$ is undefined for exactly coincident orbital planes.

For the two rotation routes, define $\psi$ as the oriented angle, within the planet's orbital plane, from its sky-plane ascending node to its disk-plane ascending node. If ${\hat{x}}_p,{\hat{y}}_p$ are the first two columns of $T_p$, take $\psi =atan2(\hat{x}\prime \cdot {\hat{y}}_p,\hat{x}\prime \cdot {\hat{x}}_p)$. Both routes describe the same planet orthonormal basis:The left route goes through the disk frame; the right goes through the sky-node planet frame. Comparing their third columns gives the two boxed relations, while comparison of first columns fixes the extra in-plane angle.

Put $I_d=\pi /2+{\delta }_d$ and $I_p=\pi /2+{\delta }_p$. A second-order Taylor expansion gives $\cos I_m=1-[({\delta }_p-{\delta }_d{)}^2+\Delta {\Omega }^2]/2+O({\delta }^4)$, and therefore the mutual inclination of nearly edge-on orbits satisfiesAll expansion angles are in radians. For plotting in degrees this is $I_m\simeq \sqrt{(I_p-87^{\circ }{)}^2+(4^{\circ }{)}^2}$. The given planet has exact $I_m\simeq 4.119^{\circ }$ and ${\Omega }_m\simeq 76.044^{\circ }$. The exact minimum occurs at $I_p=atan2(\sin I_d\cos \Delta \Omega ,\cos I_d)\simeq 86.993^{\circ }$, with $I_{m,\min }\simeq 3.995^{\circ }$. At $I_p=85^{\circ }$ and $95^{\circ }$, the second-order values are ${\sqrt{20}}^{\circ }\simeq 4.472^{\circ }$ and ${\sqrt{80}}^{\circ }\simeq 8.944^{\circ }$.

For the disk evolution, neglect disk self-gravity, collisions that align the planes, and back-reaction on the planet. In the initial disk plane, use the complex inclination $w=I{e}^{i\Omega }$ for a ring and $w_p=I_m{e}^{i{\Omega }_m}$ for the planet. Linear Laplace-Lagrange secular theory has forced inclination $w_p$ and solution $w(a,t)=w_p[1-e^{-i\nu (a)t}]$ for initially flat rings. For an inner circular perturber, the Laplace coefficient formula is $\nu =n(M_p/M_{\star })\alpha b_{3/2}^{(1)}(\alpha )/4$, with $\alpha =a_p/a$. Its quadrupole approximation uses $b_{3/2}^{(1)}\simeq 3\alpha$, yieldingIn the $(I\cos \Omega ,I\sin \Omega )$ plot referenced to the initial disk, every ring runs clockwise on the same circle centered at $w_p$, starting at the origin. Inner rings advance more quickly. If the axes are instead referenced to the planet's plane, the relative tilt is $w-w_p=-w_p{e}^{-i\nu t}$: these circles are centered at the origin, and the mutual inclination is constant. Both conventions describe the same nodal precession; they should not be mixed.

At an intermediate epoch, rings with $\nu t\gtrsim 1$ have tilted substantially while outer rings with $\nu t\ll 1$ remain near the original disk plane. This radial variation is a planet-induced debris-disk warp. The characteristic affected radius grows like $t^{2/7}$. Conservative differential nodal precession preserves each ring's tilt relative to the planet, so it does not by itself align every orbit; an unresolved inner region can acquire a mean plane near the planetary plane through phase mixing. The secular interpretation is also described in Wyatt's planetary dynamics lectures.

Consider an initial value problem $y^{\prime }=f(t,y)$ on $[0,T]$, with a Lipschitz continuous vector field in the solution variable and enough smoothness for the stated error expansions. Convergence of a numerical method meanswith consistent starting values. An order of a numerical method $p$ quantifies this approximation: with sufficiently accurate starting data and the required stability, the global error is $O(h^p)$. Stability of a numerical method controls the amplification and accumulation of perturbations; it is the link between a small local defect and small global error.

For a one-step formula $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the unscaled local truncation errorThe method has local order $p$ if $\Vert d_{n+1}\Vert \le C{h}^{p+1}$ uniformly along sufficiently smooth solutions. If the defect is divided by $h$, its order is instead $O(h^p)$; the convention must be specified. Assume also the perturbation estimateThe global error $e_n=Y_n-y(t_n)$ then satisfiesIterating, or applying the discrete Gronwall inequality, givesThere are $O(1/h)$ steps on a fixed interval, so stable propagation turns $O(h^{p+1})$ local defects into $O(h^p)$ accumulated error. Starting error $O(h^p)$ is needed to retain order $p$. This explains why matching an exact Taylor expansion is insufficient unless errors are controlled during repeated steps.

For a fixed-coefficient linear multistep method, writeThe consistency of a numerical method conditions are $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$. The higher order conditions for a linear multistep method arewith right side zero for $q=0$; exact order $p$ means first failure at $q=p+1$. They come from inserting the exact solution and comparing its Taylor expansion.

Errors at $h=0$ follow the homogeneous recurrence determined by $\rho$. Zero-stability is equivalent to the root condition for a multistep method: all roots of a polynomial have modulus at most one, and every unit-modulus root is simple. Roots outside the disk give exponential growth in the number of steps; repeated unit roots give polynomial growth in that number. Repeated interior roots are harmless because their polynomial factors are dominated by decay. Under the usual Lipschitz continuity and initialization assumptions, the Dahlquist equivalence theorem givesFor a zero-stable order-$p$ linear multistep method, local defects and starting errors of the stated orders yield global error $O(h^p)$ by a recurrence bound and the discrete Gronwall inequality. Every required starting value matters.

A concrete failure of “consistency implies convergence” isHere $\rho =(\zeta -1)(\zeta -2)$ and $\sigma =-1$, so $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)=-1$: the method is consistent, of order one. But its root $2$ violates zero-stability. For $y^{\prime }=0$ with exact solution zero and starting values $Y_0=0$, $Y_1=h^2$, the computed solution is $Y_n=h^2(2^n-1)$. These starting errors vanish, yet at $n\approx T/h$ the error diverges. Local accuracy cannot compensate for an unstable recurrence.

A different meaning of stability concerns a fixed step on decaying equations. The Dahlquist test equation $y^{\prime }=\lambda y$ introduces $z=h\lambda$. A one-step method gives $Y_{n+1}=R(z)Y_n$ and is absolutely stable when $|R(z)|\le 1$. For a linear multistep method, the corresponding amplification polynomial of a multistep method is $\rho (\zeta )-z\sigma (\zeta )$, and every root must satisfy the bounded root condition, not just the root approximating $e^z$. Strict asymptotic decay requires the roots to have modulus less than one for $\mathrm{Re}z<0$. This distinction matters at reducible methods with an undamped parasitic mode, as in question 1(c).

The Forward Euler method has order one and $R(z)=1+z$, so its linear stability domain is $|1+z|\le 1$. For $y^{\prime }=-\kappa y$, with $\kappa >0$, it requires $0\le h\kappa \le 2$. It converges as $h\to 0$, yet an excessively large step can produce growing oscillations on an exactly decaying problem. The Backward Euler method has the same order, but $R(z)=(1-z{)}^{-1}$ makes it A-stable and L-stable. For a stiff differential equation, fast decay can impose a severe explicit time-step restriction even when accuracy of the slow dynamics would permit a much larger step.

The trapezoidal rule has order two and is A-stable, but is not L-stable, since its stability function tends to $-1$ along the negative real axis. The Second Dahlquist barrier limits an irreducible A-stable linear multistep method to order at most two. This does not limit implicit Runge-Kutta methods: question 4's Lobatto IIIA method is fourth-order and A-stable. Conversely, question 1's fourth-order member is convergent but not A-stable. These examples show that order of a numerical method, zero-stability and A-stability answer different questions.

For nonlinear dissipative problems, B-stability controls distances between numerical solutions; algebraic stability of a Runge-Kutta method is a sufficient coefficient criterion. Question 4 shows that A-stability alone does not imply that coefficient criterion. Variable time steps require additional analysis of step ratios and error propagation, and stiff problems may require uniform-in-stiffness error estimates beyond the fixed-problem convergence result. The useful chain is local accuracy plus the appropriate perturbation bound, giving global accuracy; absolute and nonlinear stability then determine which practical time steps preserve the intended dynamics.

Use the Cartesian streamfunction convention $u_x={\partial }_y\Psi$, $u_y=-{\partial }_x\Psi$. In the creeping-flow regime implicit in the question's biharmonic equation hint, the Stokes equations are $\mu {\nabla }^2\mathbf{u}=\nabla p$, $\nabla \cdot \mathbf{u}=0$. Taking their curl eliminates pressure and givesThe fluid occupies the region above the actual sheet; the perturbation calculation expands its boundary about $y=0$. In the swimming frame the material velocity is $(0,-\omega y_0\cos (kx-\omega t))$. The vertical kinematic condition and the supplied Navier slip boundary condition therefore give, at $y=y_0\sin (kx-\omega t)$,At infinity, ${\Psi }_y\to U$ and ${\Psi }_x\to 0$, with bounded velocity and no imposed shear; all fields are periodic in $x$ with period $2\pi /k$. The sign is important: a sheet translating at $-U{\mathbf{e}}_x$ sees the remote fluid translating at $+U{\mathbf{e}}_x$.

Set $X=kx$, $Y=ky$, $\tau =\omega t$, $\psi =k^2\Psi /\omega$, and $U^{\ast }=kU/\omega$. The two dimensionless parameters are the small-slope approximation parameter $ϵ=k{y}_0$ and dimensionless slip length $\delta =k\gamma$. With $\theta =X-\tau$, the dimensionless biharmonic stream function for planar Stokes flow satisfiesExpand $\psi =ϵ{\psi }_1+{ϵ}^2{\psi }_2+\cdots$ and $U^{\ast }=ϵU_1+{ϵ}^2{U}_2+\cdots$. At first order,The decaying first harmonic of the biharmonic equation has the form $(A+BY)e^{-Y}\sin \theta$. The vertical condition fixes $A=1$ and excludes a cosine component. The horizontal condition givesThus $(1+2\delta )(B-1)=0$ and $B=1$ for every $\delta \ge 0$. The spatially averaged first-order streamfunction is affine in $Y$; its Navier slip boundary condition forces $U_1=0$. ThereforeThis first-order slip independence of a transverse sheet is independent of slip length. Its surface tangential velocity and surface shear are both zero: $({\psi }_1{)}_Y(0)=0$ and $({\psi }_1{)}_{YY}(0)-({\psi }_1{)}_{XX}(0)=0$. This explains why the first-order streamfunction agrees with the no-slip boundary condition.

To find the slip-enhanced swimming speed of a transverse sheet, expand only the Navier slip boundary condition through second order. Put $Q_j=({\psi }_j{)}_{YY}-({\psi }_j{)}_{XX}$. Taylor expansion at the displaced boundary yieldsDirect differentiation givesHenceAverage over $X$. The zero Fourier series mode of a biharmonic streamfunction is a cubic polynomial in $Y$; bounded velocity at infinity removes the quadratic and cubic terms. Thus its velocity is the constant $U_2$, and its averaged shear $⟨Q_2{⟩}_X$ is zero. The mean boundary velocity determines Taylor-sheet swimming speed principle then givesThe even remainder follows because changing the sign of the amplitude is equivalent to a half-period phase shift. In dimensional form,For the Navier-slip Taylor swimming sheet, the ratio of speed to the no-slip value is $1+2\delta >1$ whenever $\delta >0$: slip increases the swimming speed. The expansion is for fixed $\delta$ as $ϵ\to 0$; it is not a uniform prediction of arbitrarily large speed if the slip length is allowed to diverge with the inverse amplitude.

The arclength element is $ds=\sqrt{1+y^{\prime 2}+z^{\prime 2}}dx$. For uniformly small slopes, the Taylor expansion $\sqrt{1+s}=1+s/2+O(s^2)$ givesThus to quadratic order in the slopes the excess length is $\boxed{\ell -L\simeq \frac{1}{2}{\int }_0^L(y^{\prime 2}+z^{\prime 2})dx}$. This is the small-slope approximation used in the subsequent elastic energy.

The expansion is singular near the moving lower limit, so expanding $e^{-\epsilon /x}$ throughout the original integral would lose its constant term. Introduce the symmetric exponential integralUnder $x=\epsilon /t$, the discarded piece becomesThis is a uniform Taylor expansion because $t\ge 1$. To derive the logarithmic expansion of a symmetric exponential integral without needing another special function, first exploit the same substitution to writeThe estimate follows from $0\le 1-e^{-\epsilon /x}\le \epsilon /x$. Here $\gamma$ is Euler's constant and the small-argument expansion of the exponential integral fixes the constant as well as the logarithm.

Differentiate under the integral and integrate the derivative of $e^{-x-\epsilon /x}/x$ over $(0,\infty )$; the exponentially small endpoint terms vanish. This gives the ordinary differential equationFor the Frobenius method at its regular singular point, write $J=A\log \epsilon +B+\epsilon (C_1\log \epsilon +D)+\cdots$. Substitution gives $C_1=A$ and $D=B-2A$. The preceding leading calculation sets $A=-1$, $B=-2\gamma$, henceSubtracting the moving-cutoff correction to a symmetric exponential integral yields the requested answer, including the logarithmically enhanced first-order term:Thus the $\epsilon \log \epsilon$ switchback term must be retained when asking for all contributions through order $\epsilon$.

Assume $f$ is sufficiently smooth, for example $C^3$ near the exact trajectory. Write $f,f^{\prime },f^{\prime \prime }$ at the current value, with $f^{\prime }$ a linear map and $f^{\prime \prime }$ a bilinear map when the ordinary differential equation is vector valued. Taylor expansion of the stages givesThe exact flow has Taylor expansionMatching the terms through $h^3$ in the Runge-Kutta method gives the sufficient third-order conditions for an explicit Runge-Kutta methodFor scalar $f$, replace $f^{\prime \prime }[f,f]$ by $f^{\prime \prime }{f}^2$ and $(f^{\prime }{)}^{2}f$ by $(f^{\prime }{)}^{2}f$ in the scalar sense. The resulting local error is $O(h^4)$; under a local Lipschitz continuity assumption and bounded derivatives on a finite time interval, the standard one-step error recurrence gives global error $O(h^3)$, hence third order.

On the linear test equation $y^{\prime }=\lambda y$, putting $z=h\lambda$ gives the stability functionThus $\boxed{|R(-5/2)|=47/48<1}$, so $-5/2$ lies inside the linear stability domain. The four order conditions are compatible, for example $a_1=1/2,a_2=-1,a_3=2$ and $(b_1,b_2,b_3)=(1/6,2/3,1/6)$.

Set $\mathbf{r}={\mathbf{x}}_1-{\mathbf{x}}_2$, $M=M_1+M_2$ and $\mu =M_1{M}_2/M$, the reduced mass. Newton's law of universal gravitation gives ${\ddot{\mathbf{x}}}_1=-G{M}_2\mathbf{r}/r^3$ and ${\ddot{\mathbf{x}}}_2=G{M}_1\mathbf{r}/r^3$. Subtraction establishesFor fixed masses, the orbital energy and angular momentum areTaking derivatives gives $\dot{E}=\mu \mathbf{v}\cdot (\ddot{\mathbf{r}}+GM\mathbf{r}/r^3)=0$ and $\dot{\mathbf{J}}=\mu \mathbf{r}\times \ddot{\mathbf{r}}=0$. These are conservation of energy and conservation of angular momentum in the unperturbed Kepler orbit. To see directly that the orbital eccentricity is fixed, put $\mathbf{h}=\mathbf{r}\times \mathbf{v}$ and $\mathbf{n}=\mathbf{r}/r$. The eccentricity vector isSince $\dot{\mathbf{h}}=0$, the cross product identity $\mathbf{r}\times (\mathbf{r}\times \mathbf{v})=\mathbf{r}(\mathbf{r}\cdot \mathbf{v})-r^2\mathbf{v}$ gives $\ddot{\mathbf{r}}\times \mathbf{h}=GM\dot{\mathbf{n}}$, so $\dot{\mathbf{e}}=0$. Its magnitude $e$ is the orbital eccentricity. Equivalently,Thus the semi-major axis and orbital eccentricity of a bound noncollision Kepler orbit are constant, including $e=0$.

Take the centre of mass as origin. The two positions are ${\mathbf{x}}_1=(M_2/M)\mathbf{r}$ and ${\mathbf{x}}_2=-(M_1/M)\mathbf{r}$. For $x\gg r$, a Taylor expansion givesThe mass dipole vanishes because ${\sum }_A{M}_A{\mathbf{x}}_A=0$, and the second mass moment tensor is ${\sum }_A{M}_A{x}_{Ai}{x}_{Aj}=\mu r_i{r}_j$. The gravitational quadrupole potential of a point-mass binary is consequentlyBoth printed tensors are trace-free. Contracting them with the Kronecker delta summed over three spatial dimensions givesThis proves the requested expression through second order. In particular the printed convention is $q_{ij}=\frac{3}{2}{Q}_{ij}$, where $Q_{ij}=\mu (r_i{r}_j-r^2{\delta }_{ij}/3)$ is the standard mass quadrupole moment; the radiation coefficient must be adjusted accordingly. No mass dipole term may be retained in centre of mass coordinates.

To evaluate the quadrupole formula without assuming a circular orbit, put $k=GM$, $u=\dot{r}=\mathbf{n}\cdot \mathbf{v}$ and ${\mathbf{v}}_{\perp }=\mathbf{v}-u\mathbf{n}$. Differentiating the gravitational acceleration gives the kinematic jerkFor $S_{ij}=r_i{r}_j$, the product rule then givesInserting these into the printed mass quadrupole moment convention yieldsThe two tensors in brackets are orthogonal in their tensor contraction, because $\mathbf{n}\cdot {\mathbf{v}}_{\perp }=0$. Their squared norms are $6$ and $2v_{\perp }^2$, respectively. HenceThe instantaneous quadrupole luminosity of a Kepler binary is thereforeThe radial derivative $\dot{r}$ is not the vector speed. The bracket equals $v_{\perp }^2+u^2/12$, so the radiated luminosity is nonnegative for every instantaneous velocity.

For a slowly evolving circular orbit, there is no preferred orbital phase at which to excite a persistent eccentricity vector. More precisely, the rotating mass quadrupole moment emits at twice the orbital frequency, with angular harmonic number two; its gravitational-wave energy and angular-momentum balance is $\dot{E}=\Omega \dot{J}$. The circular Kepler orbit sequence has $J=\mu \sqrt{GMa}$ and $dE/dJ=\Omega$, so this loss is tangent to the sequence. At leading adiabatic order it preserves zero secular orbital eccentricity. This is a circular gravitational-wave inspiral, with a small radial drift rather than an exactly fixed-radius Newtonian circle. The energy flux alone would not establish circularity without this symmetry and angular momentum balance.

To leading radiation order use $r=a$, $u=0$ and ${\mathbf{v}}^2=GM/a$ in the loss formula. Combining it with $E=-GM\mu /(2a)$ givesKeeping the fixed masses and leading quadrupole formula, integration givesThis formal coalescence time displays the strong $a_0^4$ dependence: gravitational-wave emission from a binary system matters far more for close binary stars than for wide ones. The orbital period decreases as $a^{3/2}$ and the gravitational-wave frequency increases, producing a chirp. A detached binary can be driven into Roche-lobe overflow; further evolution then depends on mass-transfer stability. A sufficiently close compact pair may merge. The weak-field, slow-motion, point-mass approximation ceases to apply before literal $a=0$; finite stellar radii or strong relativistic effects determine the final interaction.