Programming function containing cut in negative imaginary axis

Question

FINAL UPDATE

I emailed the author of the paper and as it turns out there was a mistake in the equation for sigma. I gave the best answer to pv because they did help answer the problem as stated.

FIRST ATTEMPT I am trying to program a numerical representation of the function below:

,

and the '+'/'-' superscripts indicate the limits as z approaches the branch cut, which lies along the negative imaginary half axis. The H and J are Hankel and Bessel functions. The rest of the variables (n_r, m, R) depend on the geometry of the problem. I wish to plot this function along the negative imaginary half-axis with respect to k. My current code (with pv's addition) is as follows)

import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
from numpy import pi
from scipy.special import jv, iv, kv, jvp, ivp, kvp

m = 11  # Mode number
nr = 2  # Index of refraction
R = 1   # Radius of cylinder
eps = 10e-8


def yv_imcut(n, z):
    return -2/(pi*(1j)**n)*kv(n, -1j*z) + 2*(1j)**n/pi * (0.5j*pi) * iv(n,-1j*z)

def yvp_imcut(n, z):
    return (n/z)*yv_imcut(n,z) - yv_imcut(n+1,z)

def hankel1_imcut(n, z):
    return jv(n, z) + 1j*yv_imcut(n, z)

def h1vp_imcut(n, z):
    return jvp(n, z, 1) + 1j*yvp_imcut(n, z)

# Define the characteristic equation
def Dm(n, z):
    return nr*jvp(n, nr*z, 1) * hankel1_imcut(n, z) - jv(n, nr*z)*h1vp_imcut(n,z)

# Define the cut pole density function
def sigma(k,n):
    return  4*(nr**2 - 1)*jv(n,nr*k*R)/(pi**2 * k * ((Dm(n, k*R-eps).real)**2 + (Dm(n, k*R+eps).imag)**2))

k = np.linspace(-eps*1j, -15j,1000)
y = sigma(k,m)
x = np.linspace(0,15,1000)

plt.plot(x, y.imag)
plt.show()

Here's my plot of sigma.imag along the negative imaginary axis:

Here's what the plot should look like (look at the m = 11 curve on the right):

User pv helped me move the cut of the Hankel function to the negative imaginary half axis, but my plot of sigma is still off. I noted in the paper that it states that sigma is "purely imaginary" (top of page five, first column)

These equations and the graph come from page 4 in this article: http://arxiv.org/pdf/1302.0245v1.pdf

SECOND ATTEMPT

Appendix B of the article states the difference of the Hankel function across the cut as

From this relation one can also find the difference of the first derivative of the Hankel function across the cut:

I wrote a script with these formulae:

def hankel1_minus(n,z):
    return hankel1(n,z) - 4*jv(n,z)

def h1vp_minus(n,z):
    return (n/z)*hankel1_minus(n,z) - hankel1_minus(n+1,z)

def Dm_plus(n, z):
    return nr *jvp(n, nr*z, 1) * hankel1(n, z) - jv(n, nr*z)*h1vp(n,z)

def Dm_minus(n, z):
    return nr *jvp(n, nr*z, 1) * hankel1_minus(n,z) - jv(n, nr*z)*h1vp_minus(n,z)

def sigma(k,n):
    return  4*(nr**2 - 1)*jv(n,nr*k*R)/(pi**2 * k * (Dm_plus(n, k*R) *     Dm_minus(n,k*R)).real)

Plotting this sigma gives the same result as the first method.

Answer 1

The branch cut of H1 in scipy is (-inf, 0) and not in (-1j*inf, 0) as expected in the paper you cite, which explains why you get incorrect results.

The problem can be worked about with some creative use of an argument transformation for Y_nu, as I point out above in the comments.

Let us assume integer order n. We have

hankel1(n, z) = jv(n, z) + 1j*yv(n, z)

jv has no branch cuts (integer order), but yv has. The transformation formula reads

yv(n, 1j*z) = -2/(pi*(1j)**n)*kv(n, z) + 2*(1j)**n/pi * (log(1j*z) - log(z))*iv(n,z)

or, in other words,

yv(n, z) = -2/(pi*(1j)**n)*kv(n, -1j*z) + 2*(1j)**n/pi * (log(z) - log(-1j*z))*iv(n,-1j*z)

kv(n,z) is defined in Scipy to have branch cut at (-inf, 0), and iv(n,z) has no branch cuts (integer order). Apart from the logarithms, the branch cuts on the RHS are therefore in (-1j*inf,0 ), exactly at where we want them to be. The only thing left to do is to choose the branch cut of the logarithm term appropriately.

The correct analytic continuation with branch cut in (-1j*inf, 0) is then

def yv_imcut(n, z):
    return -2/(pi*(1j)**n)*kv(n, -1j*z) + 2*(1j)**n/pi * (0.5j*pi) * iv(n,-1j*z)

This coincides exactly with yv(n, z) in 3 of the 4 quadrants. It has a branch cut in (-1j*inf,0). Moreover, it is obviously an analytic function. Therefore, it is the same as yv, but with a different branch cut choice.

We then have

def hankel1_imcut(n, z):
    return jv(n, z) + 1j*yv_imcut(n, z)

It is obviously then the Hankel function with branch cut in (-1j*inf, 0).

Based on this, you can also work out the derivatives.