Caustics

Quantum Mechanics
                       K. Konishi, G. Paffuti

Study the caustic curves for a uniform gravitational field and for the two body gravitational problem.

ClearAll["Global`*"];
Off[General::spell, General::spell1];

Uniform gravitational field

Let us consider a particle of mass m moving in a uniform gravitational field g. We limit ourselves to a fixed plane, x is the horizontal direction and g is directed toward negative y. y=0 is the ground.

We consider a pencil of trajectories at fixed energy. The starting point is y = (Subscript[y, 0]). The angle φ of starting momentum with respect to x axis is the parameter distinguishing the trajectories.

As it is well known from elementary mechanics (we take the origin at the starting point of the trajectories):

x = v Cos[φ] t ;           y =  v Sin[φ]t - (1/2)g (t^2) .

fx[t_,φ_]:= v Cos[φ] t; fy[t_,φ_]:=v Sin[φ]t - (1/2)g (t^2)

The envelope of the family is defined by

| ({{(∂_t)x, (∂_t)y}, {(∂_φ)x, (∂_φ)y}})| = (∂_t)x (∂_φ)y - (∂_t)y(∂_φ)x = 0 .

envelope = (∂_t)fx[t,φ] (∂_φ)fy[t,φ] - (∂_φ)fx[t,φ] (∂_t)fy[t,φ];

The caustic is obtained by the previous equations, eliminating t and φ (we limit ourselves to the region x > 0 y > 0 ):

Off[Eliminate::"ifun"];
implicitCaustic=Simplify[Eliminate[{x== fx[t,φ],y== fy[t,φ],envelope== 0},{t,φ}],Assumptions-> {x>0,y>0}]

Solve[implicitCaustic,y]

This is a graphical view of the trajectories and of the associated caustic

Block[{v=1.,g=1.,y0=0.},
Show[{ParametricPlot[{fx[t,#1],fy[t,#1]},{t,0,((2(v^2))/g)}]&/@Range[(π/4),(π/2),0.05],
Plot[(((v^4)-(g^2) (x^2))/(2 g (v^2))),{x,0,((v^2)/g)},PlotStyle-> {Red,Thickness[0.0025]} ]}
,PlotRange-> {{0,((v^2)/g)},{0,((v^2)/(2g))}}]]

Let us note that in this case the equation for the caustic could be obtained by the following procedure. The curve is a (singular) solution of Hamilton Jacobi equation, then it represents a possible trajectory. In a constant gravitational field this means that the curve is a parabola. By symmetry consideration its vertex is at x=0 and can be computed by the highest point of a vertical trajectory:  y = (v^2)/2g.   The parabola must be symmetric for x → - x so a second point fixes the curve. For t→∞ the trajectories behave as x = v t, y ∼ -1/2 g (t^2)= -1/2 g x^2/(v^2) then the caustic has the form

y = ((v^2)/(2g))- (1/2)g ((x^2)/(v^2)).

This is exactly what we have found by the explicit computation.

Central field

Let us consider now a central field problem. We can imagine a family of trajectories as the orbits of a rocket fired from the "Earth" in the gravitational field of the "Sun". We will limit again to a planar family.

It is known that finite orbits are ellipses, the Sun at one of the foci. An explicit solution of this problem is first due to Jacobi, we proceed as in the previous case: a singular solution of Hamilton Jacobi equation is an orbit, then in this case an ellipse. By symmetry the major axis lie on the segment earth-sun, let us call it x-axis. The earth has coordinates (-(Subscript[r, 0]),0). The trajectories must be tangent to the caustic, then it is easy to see that one extremal orbit is a rocket moving toward negative x axis, up to a maximum distance where its orbit touch the caustic. A second orbit is given by considering the earth at the aphelion or at the perihelion, the opposite point touch the caustic.

Let v the velocity of the rocket. The (fixed) total energy will be (m = 1)

ϵ = (1/2)(v^2) - (α/(Subscript[r, 0]))

We take ϵ < 0. For the first orbit the maximum distance (Subscript[r, 1]) is (see fig. below)

|E|=  (α/(Subscript[r, 1]))=(α/(Subscript[r, 0]))- (1/2)(v^2);      (Subscript[r, 1])=((2α (Subscript[r, 0]))/(2α-(Subscript[r, 0])(v^2)));   AS=(Subscript[r, 1])=(α/(|E|));

Consider now the second orbit. It is known that the major axis of the ellipse is related to energy by

TB = (α/(|E|))= AS

The caustic ellipse then has a major axis

AB = (AS - TS) + TB = 2 TB - TS = 2 (Subscript[r, 1])- (Subscript[r, 0]).

To fix the eccentricity parameter ϵ we have to determine if B is the perihelion or the aphelion. The distance SB is  a (1 ∓ ϵ) in the two cases, a = AB/2 is the semi-axis.

Once found ϵ we can compute minor axis b and focal distance c by c = ϵ a, b = a (Sqrt[1-(ϵ^2)]), and the ellipse can be written in polar coordinates as

(((b^2)/a)/r) = 1+ϵ Cos[φ];

or in parametric form

x = a(Cos[ξ]-ϵ);   y = b Sin[ξ].

Let us verify our conclusions by solving the Newton equations for the rocket and plotting its orbit, for different inclination of initial velocity.

r0=2;Tmax=20;

plotOrbit[v_,φ_]:= Module[{sol},
sol = NDSolve[{y''[t] == -((y[t])/(((((x[t])^2)+((y[t])^2)))^(3/2))),x''[t] == -((x[t])/(((((x[t])^2)+((y[t])^2)))^(3/2))),
x[0]== -r0,y[0]==0,x'[0] == v Cos[φ],y'[0]==v Sin[φ]},{x,y},{t,0,Tmax}];
ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,Tmax}]]

plotCaustic[v_]:= Module[{AS,AB,SB,a,b,c,ϵ,absene},
absene= 1/r0 - 0.5(v^2);
AS = 1/absene; AB = 2 AS-r0;
a = AB/2;  SB = AS-r0;
ϵ = Abs[1- SB/a];  b = a (Sqrt[1-(ϵ^2)]);
ParametricPlot[{a(Cos[ξ]-ϵ),b Sin[ξ]},{ξ,0,2π},PlotStyle-> {RGBColor[1,0,0],Thickness[.005]}]];

v=0.45;

orbits= plotOrbit[v,#1]&/@Range[0.2,1.99Pi/2,0.25];

caustic = plotCaustic[v];

Show[{orbits,caustic},PlotRange-> All,Ticks-> None]