I have been making strides on my project with the help of others here but have run into another snag: I am having difficulty in expressing the Lienard-Wiechert potentials with retarded time. For context, my project is to show that using Mathematica to turn the Liennard-Wiechert potentials into fields is simpler than coding the direct field expressions (Jefimenko Equations). Due to such, I need to properly express my potential equations with respect to retarded time.

A fellow user provided this equation for the retarded time:

R1 = {x, y, z};
R2 = {xb[t], yb[t], zb[t]};
    R = R1 - R2;
tr[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] := 
  r /. FindRoot[
    r == t - Sqrt[(x - xb[r])^2 + (y - yb[r])^2 + (z - zb[r])^2], {r, 
     t - Sqrt[(x - xb[t])^2 + (y - yb[t])^2 + (z - zb[t])^2]}];

where R1 is the field point and R2 is the location of the source charge. I appropriated this code into my own equations for the potentials as follows:

Rs = {xs[t], ys[t], zs[t]};(*location of the source*)
Rp = {x, y, z};(*location of the field point*)
    R = Rp - Rs;

tr[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
  r /. FindRoot[
    r == t - Sqrt[(x - xs[r])^2 + (y - ys[r])^2 + (z - zs[r])^2]/c, {r, t - Sqrt[(x - xs[t])^2 +
         (y - ys[t])^2 + (z - zs[t])^2]}];

xs[t_] := (u*Cos[\[Theta]]*t) ;(*x-pos of source*)
ys[t_] := (u*Sin[\[Theta]]*t) ;(*y-pos of source*)
zs[t_] := 0;
v[t_] = D[Rs[t], t];(*velocity of source*)

LWS[Rp_, t_] := (1/4*Pi*\[CurlyEpsilon])*(q*c/
      (Norm[R]*c - Dot[R, v[t]]));(*equation for the Lienard-Wiechert scalar potential*)

LWV[Rp_, t_] := (\[Mu]/4*Pi) (q*c*v[t]/
      (Norm[R]*c -Dot[R, v[t]]));(*equation for the Lienard-Wiechert vector potential*)

I then tried to change my electric field equation to accommodate the retarded time:

LWVdt[t_] := D[LWV[Rp, t] /. Abs@x_^2 :> x^2 /. Abs@y_^2 :> y^2, t];
LWSgrad[t_] := Grad[LWS[Rp, t] /. Abs@x_^2 :> x^2 /. Abs@y_^2 :> y^2, Rp];
LWEField[Rp_, t_] := -LWSgrad@t - LWVdt@t /. t -> tr[t, x, y, z];(*LW electric field equation*)

But when I try to plot it using

VectorPlot3D[Evaluate@LWEField[Rp, t] /. t -> -.25, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]

Mathematica just keeps running and running, never actually producing a vector plot. Am I doing something wrong with the evaluation here? Maybe some strange combination of set and set delayed?

  • $\begingroup$ If I just try something simple, like LWEField[Rp, t] /. t -> -.25 /. {x -> -1, y -> 1, z -> 0}, I get a screen full of messages from FindRoot about non-numeric quantities. Did you provide constant values for things like q ,c, epsilon, etc.? I don't see those anywhere. You also seem to be using theta in places, but your plot is in Cartesian x,y,z, so that's going to be another problem. I recommend just trying a few simple test cases before attempting a plot. Maybe evaluate using Table for a few points first, and examine the output to make sure you get numbers back. $\endgroup$
    – user87932
    Oct 17, 2023 at 3:55
  • $\begingroup$ My apologies, I forgot to provide the constants in the code. I have them set as q = 1, c = 1, \[CurlyEpsilon] = 1, \[Mu] = 1, u = 0.1*c, \[Theta] = 0. I only used theta for xs and ys in order to express spherical coordinates in their Cartesian components. Does it not work the way I have it? When I expressed the potentials without retarded time, the fields plot just fine. $\endgroup$
    – JDRobin
    Oct 17, 2023 at 19:13
  • $\begingroup$ After setting the constants, I still seem to see a problem. After running Evaluate[LWEField[Rp, t] /. t -> -.25], then setting x,y,z to some random values (% /.{x -> -.2, y -> .3, z -> .7}), I don't get a set of numbers. I see things like attempts to take the derivative of a list of numbers vs a number, etc. However, the answer posted below seems to resolve these problems, so perhaps it's best to just study it and see what changes were made to get this to work. $\endgroup$
    – user87932
    Oct 17, 2023 at 19:59
  • $\begingroup$ Understood, thanks for letting me know! I made the necessary changes from below and evaluating LWEField does now give me a vector {x, y, z} where each component is some string of numbers, but when I go to plot that evaluated field, Mathematica just keeps running and running—it is going on three minutes as of submitting this comment. Do some 3D vector plots take a while to show up? Is this unusual? I still have plenty of memory left on my PC as it evaluates. $\endgroup$
    – JDRobin
    Oct 17, 2023 at 21:29
  • $\begingroup$ I don't know. It might be best to ask Alex Trouney, since he made such a plot. Since many plotting routines do adaptive sampling, they might spend time in specific regions to fulfill some goal. One alternative is to use Table to generate a set of values on some uniform grid, then use ListVectorPlot3D and compare the timing vs VectorPlot3D. $\endgroup$
    – user87932
    Oct 17, 2023 at 22:11

1 Answer 1


We take the retarded time expression from my answer here

R1 = {x, y, z};
R2 = {xb[t], yb[t], zb[t]};
    R = R1 - R2;
tr[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] := 
  r /. FindRoot[
    r == t - Sqrt[(x - xb[r])^2 + (y - yb[r])^2 + (z - zb[r])^2], {r, 
     t - Sqrt[(x - xb[t])^2 + (y - yb[t])^2 + (z - zb[t])^2]}];

To use this expression with the Lienard-Wiechert potentials we define electric and magnetic field as follows (here we use CGS system with $c=1,q=1,\varepsilon =\mu=1$):

 Vt = D[R2, t, t];
    V0 = D[R2, t];  
    Ev = ((1 - V0 . V0)*(R - V0*Norm[R]) + 
              Cross[R, Cross[(R - V0*Norm[R]), Vt]])/(Norm[R] - R . V0)^3 /. 
           t -> tr[t, x, y, z];R0 = R /. t -> tr[t, x, y, z]; B = Cross[R0, Ev]/Norm[R0]; 

Let define particle trajectory in a form

xb[t_] := u*Cos[t] ;(*x-pos of source*)
yb[t_] := u*Sin[t] ;(*y-pos of source*)
zb[t_] := 0;

With these definitions we can plot electric and magnetic field at a given point $x=2,y=2,z=2$ vs time

p={Plot[Evaluate[Ev /. {x -> 2, y -> 2, z -> 2}], {t, 0, 2 Pi}, 
  Frame -> True, FrameLabel -> {"t", "E"}, 
  PlotLegends -> {"Ex", "Ey", "Ez"}, PlotRange -> All], 
 Plot[Evaluate[B /. {x -> 2, y -> 2, z -> 2}], {t, 0, 2 Pi}, 
  Frame -> True, FrameLabel -> {"t", "B"}, 
  PlotLegends -> {"Bx", "By", "Bz"}, PlotRange -> All]} 

Figure 1

We also can define retarded time using Nest[] as follows

tr1 = Nest[
   t - Sqrt[(x - xb[#])^2 + (y - yb[#])^2 + (z - zb[#])^2] &, 
   t - Sqrt[(x - xb[t])^2 + (y - yb[t])^2 + (z - zb[t])^2], 6];

With this definition we can directly use Lienard-Wiechert potentials as follows

 LWs = 1/(Sqrt[R . R] - Dot[R, V0]) /. t -> tr1;
LWv = V0/(Sqrt[R . R] - Dot[R, V0]) /. t -> tr1;  (*equations for the Lienard-Wiechert scalar and vector potential*)

And then we can define electric and magnetic fields

E1 = -D[LWv, t] - Grad[LWs, {x, y, z}]; b1 = Curl[LWv, {x, y, z}];


p1 = {Plot[Evaluate[E1 /. {x -> 2, y -> 2, z -> 2}], {t, 0, 2 Pi}, 
   Frame -> True, FrameLabel -> {"t", "E"}, 
   PlotLegends -> {"Ex", "Ey", "Ez"}, PlotRange -> All, 
   PlotStyle -> {{Red, Dashed}, {Blue, Dashed}, {Black, Dashed}}], 
  Plot[Evaluate[b1 /. {x -> 2, y -> 2, z -> 2}], {t, 0, 2 Pi}, 
   Frame -> True, FrameLabel -> {"t", "B"}, 
   PlotLegends -> {"Bx", "By", "Bz"}, PlotRange -> All, 
   PlotStyle -> {{Red, Dashed}, {Blue, Dashed}, {Black, Dashed}}]}

Figure 2 Now we can compare plots p,p1 in one plot as

Table[Show[p[[i]], p1[[i]]], {i, 2}]

Figure 3

As we can see, the match is relatively good, so both approaches are equivalent. To plot vector field in 3D we use, for example

 Evaluate[E1 /. {t -> 2.2}], {x, -2, 2}, {y, -2, 2}, {z, 1, 2}, 
 VectorStyle -> Arrowheads[0]]

Figure 4

Finally we can compare retarded time in a form of tr, tr1, and in an explicit form tr0=t - Sqrt[(x - xb[t])^2 + (y - yb[t])^2 + (z - zb[t])^2] in one plot at $x=2,y=2,z=2$

Plot[{tr[t, 2, 2, 2], tr1 /. {x -> 2, y -> 2, z -> 2}, 
  t - Sqrt[(x - xb[t])^2 + (y - yb[t])^2 + (z - zb[t])^2] /. {x -> 2, 
    y -> 2, z -> 2}}, {t, 0, 4 Pi}, 
 PlotStyle -> {{Blue}, {Red, Dashed}, {Green}}, 
 PlotLegends -> {"FindRoot", "Nest", "Explicit"}, Frame -> True, 
 FrameLabel -> Automatic] 

Figure 5 As we can see from this figure, the retarded times calculated using FindRoot and Nest are close to each other. Whereas the explicit expression is very far from both of them.

  • $\begingroup$ In your visualization for the fields produced by the lienard-wiechert potentials, you evaluated Ev and B instead of Ev1 and B1. Was that intentional? Also, what order should I be introducing the code in? As of right now I start with the constants, then the trajectory (xs, ys, zs), then the positions (Rs and Rp), then the retarded time, then Vt and V0, and finally the potentials and fields. Does my starting order affect whether a plot shows up? Lastly, why do we not set delay R1, R2, Vt, V0, and the field equations? $\endgroup$
    – JDRobin
    Oct 17, 2023 at 20:56
  • $\begingroup$ @JDRobin See updated version for second code with new definition for retarded time. $\endgroup$ Oct 18, 2023 at 5:15
  • $\begingroup$ Thank you for the update. I could not get the potentials to work when using set delay and t -> tr, but it works perfectly fine when I do away with the set delay and function definition and instead use t -> tr1. Would you mind explaining how is it that nesting the expression in tr1 six times in its function is the same thing as finding the root r in tr? I looked at the full equation and do not understand how it lends to a retarded time. $\endgroup$
    – JDRobin
    Oct 19, 2023 at 18:42
  • $\begingroup$ @JDRobin This is a mathematical trick. It relies on the fixed point theorem. $\endgroup$ Oct 20, 2023 at 3:22
  • $\begingroup$ Thanks for letting me know. I looked up the theorem and would like to know what your tolerance was, e.g., 10^-4, 10^-5, etc.? How did you choose the tolerance value? $\endgroup$
    – JDRobin
    Oct 20, 2023 at 14:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.