A previous presentation defined a generalized Lambert function of two arguments as the implicit solution of the nonlinear transcendental equation

$W(x,y) =xe−W (x,y) +yeW (x,y)$

Since it is defined and evaluated in the same ways as the ordinary Lambert function, it is an as equally valid analytic function as the latter. A numerical implementation is available in Math.

This double Lambert function can be used to solve a variety of transcendental equations involving trigonometric functions. Consider one relating a function to its hyperbolic cosine:

$w=xcoshw =x2 (ew +e−w ) → w=W( x2 , x2)$

The solution for a hyperbolic sine differs by a sign on an argument:

$w=xsinhw =x2 (ew -e−w ) → w=W( −x2 , x2)$

A linear combination of these two functions on the right-hand side is equally soluble:

$w=xcoshw +ysinhw =x+y2 ew +x-y2 e−w w=W( x-y2 , x+y2 )$

Here is a visualization of the solutions for the two single hyperbolics on the complex plane:

The branch of the solution is set with the slider. Note that the principal branch for a choice of hyperbolic sine is identically zero.

Repeat this for equations containing circular functions. First the cosine:

$w=xcosw =x2 (eiw +e−iw ) → w=−iW( ix2 , ix2)$

The solution for a circular sine again differs by a sign on an argument:

$w=xsinw =x2 (eiw -e−iw ) → w=−iW( −ix2 , ix2)$

A linear combination of these two functions on the right-hand side is again soluble:

$w=xcosw +ysinw =x+y2 eiw +x-y2 e−iw w=−iW( i(x -y)2 , i(x +y)2 )$

Here is a visualization of the solutions for the two single circulars on the complex plane:

These complex branches have the same form as those for the single hyperbolics, but rotated ninety degrees clockwise and the parts interchanged. The principal branch for a choice of circular sine is again identically zero.

Unfortunately the double Lambert function does not appear to be helpful with equations containing the other trigonometric functions, such as $w=xtanhw$ , because of the relative distribution of exponentials. Such equations appear to require another new function.

A few comments are in order concerning the appearance of the complex branches depicted here versus the double Lambert function in general. For the latter as function of either variable one at a time the asymptotic behavior is the logarithm of a square root. The branches will have the familiar corkscrew appearance of the logarithm overall.

The solutions above, however, have balanced arguments that do not just cancel out asymptotically. This leads to behavior on the principal branches of the cosine solutions akin to those of the inverse hyperbolic and circular tangents, with an overall change in elevation of π/2 in moving from one limit to the other.

Higher branches of these solutions clearly do not have the simple behavior of higher branches of the inverse trigonometric functions, which are just displacements in elevation. Their asymptotic values share some of the simple behavior, starting from ±π/2 on the principal branch and being spaced apart on higher branches by π when nonzero.

Kepler’s equation relating mean anomaly to eccentric anomaly

$ωt=τ =u-εsinu$

has been a known quantity in astronomy for some four hundred years. It has no exact solution in terms of standard special functions. While there are a variety of series approximations to the solution, they generally converge very slowly. It is simpler to perform numerical inversion of the equation using a method initiated by Newton himself. For real values of the independent variables this is easily accomplished:

The solution consists of a two parts: constant linear growth and variation around it of the same period as the sine function.

If the solution $u\left(\tau ,\epsilon \right)$ of the equation is treated as a function of complex variables, it is exactly soluble in terms of the double Lambert function. First explicitly exhibit the exponentials in the sine function,

$τ=u +iε2 (eiu -e−iu )$

then shift the variable with $w=\tau -u$ :

$iw =iε2 (eiτ -iw -e−iτ +iw) iw=−ε2 eiτ e−iw +ε2 e−iτ eiw$

The solution of this equation in terms of the double Lambert function is

$w=−iW( −ε2 eiτ , ε2 e−iτ )$

and the exact solution to the Kepler equation is

$u(τ,ε) =τ +iW( −ε2 eiτ , ε2 e−iτ )$

Comparing this solution to those given at the outset, the exponentials in the arguments clearly arise from the presence of the linear term.

As a check of the method, repeat the interactive graphic above along with the real part of the principal branch of this new solution in purple:

Et voilà: the solution works as expected! Note that the imaginary part of this solution is zero on the real axis. The real part of the function with naturally increase in steps of 2π from the linear contribution.

Since the principal branch matches the numerical solution along the real axis, one would expect it to provide a solution on the entire complex plane. One would also expect that specifying the index of the double Lambert function would provide other branches of this complex function:

$un (τ,ε) =τ +iWn( −ε2 eiτ , ε2 e−iτ )$

The function as written has two parts: a linear increase in mean anomaly and a contribution manifestly periodic in the real part of the mean anomaly. For purely real input parameters, this implies that this term will cycle between forms proportional to the solutions above for the single hyperbolic sine and the single circular cosine. This leads to complex branches that are not what one might expect.

Begin with the behavior of the solution as a function of eccentricity. When the arguments of the double Lambert function are large in absolute value, one can use the asymptotic approximation

$Wn (x,y) ≈lnx −y +nπi$

which becomes exact in the same limit. For the solution to the Kepler equation this naively implies

$limε→∞ un (τ,ε) ≈τ+i[ ln−ε2 eiτ −ε2 e−iτ +nπi] ≈−nπ$

so that the branches of the solution would appear to be uniquely labeled by the asymptotic value of its real part. Here is how the branches actually appear as a function of eccentricity on the complex plane:

The mean anomaly is entered in this graphic as a string of the sort one normally uses to express complex numbers. The same will apply to input boxes below whose default value is such a string.

If the value of a function on the real axis is taken to be part of its principal branch, then the shifted region of each sheet of the real part is not an error. This can be verified by plotting the defining equation and varying the eccentricity up and down its real axis:

For all settings of mean anomaly other than multiples of π, there is a change in value of π in traversing the real ε-axis, which is precisely the behavior seen in the previous graphic. Note that there will always be a real solution to the Kepler equation on the principal branch for real parameters. Other real intersections that appear belong to adjacent branches.

The location in terms of the argument of the eccentricity where the shift to a different asymptotic sheet occurs depends on the real part of the mean anomaly. This can be seen by setting that value to multiples of π/2, leading to either flat sheets with no jump or those with a jump precisely in the middle. This behavior is due to the fact that the real part of the mean anomaly determines the relative contributions of the solutions above for the single hyperbolic sine and the single circular cosine.

Now consider the complex branches as a function of mean anomaly on the complex plane. Since the behavior in the real direction is periodic, its range is limited to ±4π:

The jumps in the imaginary part are rather unpleasant, but again not an error. They arise because of how the complex sheets of the double Lambert function are organized. Since this latter function has the well-behaved logarithmic asymptotic expression given above, the set of its solutions is divided into groups closest in their imaginary parts to this asymptotic behavior.

For the complex Kepler function, the extra factor of the imaginary unit means that its set of solutions is divided by proximity in their real parts to its own asymptotic behavior. This can be understood by writing the complex Kepler equation as

$τ-u+εsinu =0$

and visualizing the real (blue) and imaginary (red) parts of the resulting contours along with a couple helpful indicators:

The gray dot here is the value corresponding to replacing the full double Lambert function in the solution to the complex Kepler equation with its logarithmic asymptotic expression. Due to the fact that a numerical exponential operates modulus 2πi on its input, the value can be offset from the naive limit above by a factor of π but is clearly visually in the proper location on the real axis. The set of intersections representing solutions to the complex Kepler equation are organized by their real parts relative to this value.

The purple line is located at the real part of the solution to the complex Kepler equation using the full double Lambert function. It indicates which of the two nearby intersections are actually taken to be part of the complex branch labeled by n. As the imaginary part of the mean anomaly is altered, this line passes through either an upper or a lower intersection, and when the choice changes there will be discontinuities in the visualized complex sheet. Such jumps occur for the double Lambert function itself but are not as distracting to the eye. Perhaps this is a result of the magnification resulting from an exponential function.

One could remove the discontinuities occurring as a function of the mean anomaly by arbitrarily choosing to use either the upper or lower intersection for the entire sheet. This would lead to a nice visualization here but introduces odd behavior as a function of eccentricity. Since there is a solid line of reasoning as to how the complex branches of the double Lambert function are organized, it is perhaps better just to accept these odd looking branches for what they are.

As the imaginary part of the mean anomaly approaches infinity, the first argument of the double Lambert function goes to zero. This means it can be replaced with an ordinary Lambert function of negative argument and negative sign:

$limτI →∞ u(τ,ε) =τ -iW( −ε2 e−iτ )$

Replacing the Lambert function with its known approximation

$W(z) ≈lnz-lnlnz$

$limτI →∞ u(τ,ε) ≈τ -iln( −ε2 e−iτ ) +iln[ln( −ε2 e−iτ )] limτI →∞ u(τ,ε) =−iln( −ε2 ) +iln[ τI -iτR +ln(−ε2 )]$

so that the overall behavior of the imaginary part is a logarithm. As the imaginary part of the mean anomaly approaches negative infinity, the second argument of the double Lambert function goes to zero. This means it can be replaced with an ordinary Lambert function:

$limτI →−∞ u(τ,ε) =τ +iW( −ε2 eiτ )$

Using the approximation gives

$limτI →−∞ u(τ,ε) ≈τ +iln( −ε2 eiτ ) -iln[ln( −ε2 eiτ )] limτI →−∞ u(τ,ε) =iln( −ε2 ) -iln[ −τI +iτR +ln(−ε2 )]$

and the overall behavior of the imaginary part is again a logarithm. Oddly enough as the integer labeling the branch is altered in their visualization the imaginary part changes sign depending on whether n is even or odd. This behavior is not captured by the two limits just given, but is readily apparent from the contours of the complex Kepler equation.