26 views (last 30 days)

Hi,

I want to solve the characteristic equation describing the TE mode in an optical fiber (Light transmission optics, D. Marcuse):

where a is the core radius of the fiber. The key parameter to solve in this equation is the propagation constant β that is related to κ and γ as

Parameters and are constants in the core and cladding of the fiber related to the wavenumber as , being the refractive index of core and cladding, respectively.

I know there are multiple solutions and are related to the zeros of the bessel function in the numerator. I tried to solve this equation using matlab functions solve() and fzero() for the values , , and but I'm stuck. I need some help to find the solution of this equation.

Thanks for your attention.

Carolina Rickenstorff

David Goodmanson
on 18 Feb 2020

Edited: David Goodmanson
on 22 Feb 2020

Hi Carolina,

It's not a real good idea to go with an answer solely on the grounds that 'solve' can find one from an expression in a book. With books, there are definitely issues with the notation, who is using what, so you have to make sure you are solving the right equations. The E and H fields on each side of the boundary r = a are

E: G*J1(k1*a) = alpha*G*H1(k2*a)

H: G*(1/mu1)*k1*J0(k1*a) = alpha*G*(1/mu2)*k2*H0(k2*a)

% not all the constant conversion factors to get from E to H are shown.

% these cancel out when doing equality of H between the two regions.

where E is in the theta direction, around the cylinder, and H is in the z direction, along its length. (There is also a radial component of B that is matched automatically). G is an overall amplitude that doesn't matter for matching, alpha is an unknown constant to get the sizes to match up, and k1 and k2 are

k1 = sqrt(n1^2 k^2 - kz^2) k2 = sqrt(n2^2 k^2 - kz^2)

kz is your beta, which must be the same in both regions since the solution has to have the same z dependence down the length of the cylinder. k1 is your kappa, and k2 is almost your gamma. This is assuming that in your definition of k, k = 2pi / lambda, the wavelength is the free space version, lambda = c / f.

On the right, H1 and H0, which are really denoted H1(+) and H0(+) are certain Hanckel functions. (There is also an H1(-) and H0(-) ).

Dividing one equation by the other yields

mu1*J1(k1*a) / (k1*J0(k1*a)) = mu2*H1(k2*a) / (k2*H0(k2*a))

For a solution, kz has to be such that k1 is real and k2 is pure imaginary. Rewrite this as

k2 = i*g2 where g2 = +sqrt(kz^2 - n2^2 k^2)

and g2 is your gamma.

The fields in the outside region have to decay exponentially for large r, and the bessel functions K have that property. It turns out that

H0(i*g2*a) = -i*(2/pi)*K0(g2*a)

H1(i*g2*a) = -(2/pi)*K1(g2*a)

and after tossing those in we arrive at last at

mu1*J1(k1*a) / (k1*J0(k1*a)) = - mu2*K1(g2*a) / (g2*K0(g2*a))

Here the J's and K's are straight Matlab besselj and besselk.

Note the minus sign on the right. Various definitons of bessel functions include factors of (-)^n where n is the order of the bessel function. That messes with the minus signs.

This reply has gone on for a bit, but I beleive the process is absolutely necessary because of all the ways of defining bessel functions.

The larger the value of k*a, the more modes (roots) there will be. Your value of k*a is 40.54 which is large enough for several modes.

The code below shows where the roots are. I will leave it to you to use fzero or whatever method you choose to find roots.

a = 10e-6;

k = 2*pi/1.55e-6;

mu1 = 1;

n1 = 1.4;

mu2 = 1;

n2 = 1.33;

kzmax = n1*k;

kzmin = n2*k;

kz = linspace(kzmin,kzmax,100000);

kz(1) = [];

kz(end) = [];

k1 = sqrt(n1^2*k^2 - kz.^2);

g2 = sqrt(kz.^2 - n2^2*k^2);

C = mu1*besselj(1,k1*a)./(k1.*besselj(0,k1*a));

% quiet down the plot y values

Cm = mean(abs(C));

C(abs(C)>3*Cm) = nan;

D = -mu2*besselk(1,g2*a)./(g2.*besselk(0,g2*a));

Dm = mean(abs(D));

D(abs(D)>3*Dm) = nan;

figure(1)

plot(kz,C,kz,D)

grid on

kzroot = 5.4617301149e6; % first root from plot

% check, stay with besselh and k2, not besselk and g2

% besselh has an extra input to choose H(+) rather than H(-)

k1root = sqrt(n1^2*k^2 - kzroot^2);

k2root = sqrt(n2^2*k^2 - kzroot^2);

alpha = besselj(1,k1root*a)/besselh(1,1,k2root*a);

% E and H should agree on each side

E1 = besselj(1, k1root*a)

E2 = alpha*besselh(1,1,k2root*a)

% there are constant conversion factors to get from E to H, not shown here.

% these do not affect an equality check of H between regions 1and 2

H1 = (1/mu1)*k1root*besselj(0, k1root*a)

H2 = alpha*(1/mu2)*k2root*besselh(0,1,k2root*a)

David Goodmanson
on 27 Feb 2020

Thank you, Walter, for your comment.

Carolina,

Looks like you got this thing solved. Since you asked,

yes the largest possible order here is five, the number of roots shown on the plot. (Sometimes people start counting at 0 instead of 1, but five altogether). The limits

kzmax = n1*k; kzmin = n2*k;

are there for plotting purposes that's true, but also because outside of that range either k1 goes imaginary or (in the convention you are using), (i*k2) goes real. Either of those produces a bad nonphysical solution.

Hope you don't mind but I have a couple more comments. First of all, the code uses w*mu0 to find Ephi which is correct numerically but I think it better illustrates the physics to use

mu0*w = mu0*c*k = 120pi*k = 377*k

where mu0*c = 377 ohms is the impedence of free space. Then

Ephi1 =-1i.*b.*(120*pi)*(k/k1root).*bessellj(...

i.e., E [V/m] is a resitance times a dimensionless quantity (k/k1root) times H [A/m].

Using ind_roots=find(abs(C-D)<=0.5e-10); will work most of the time but is not the most reliable of methods. Multiple roots make it a bit inconvenient but a zero finder such as fzero is preferable.

In the definition of kz, 500000 points is way more than too much. It's only there because I temporarily upped the number in order to inspect the transition point r = a and then forgot to put it back to something better. Maybe 1000 or whatever you think, makes sense.

Sign in to comment.

Carolina Rickenstorff
on 28 Feb 2020

Edited: Carolina Rickenstorff
on 28 Feb 2020

Sign in to comment.

David Goodmanson
on 2 Mar 2020

Hi Carolina,

do you have a reference for where you read that dielectric fibers have no cutoff frequency? I would be interested to see who is maintaining that. If you run your code and look at the plot of C vs D, you can see that D is always negative. D is a function of k2*a, and allowed values of k2*a depend on the details of this problem, but you can use any positive value you like for k2*a and D is always negative. That means for a solution,

C = besselj(1,k1*a)./(k1.*besselj(0,k1*a))

has to be negative. But C is only negative if J0 and J1 have opposite sign. If you plot these out, both start out positive and fiirst have opposite sign when J0 goes negative at k1*a > 2.4048, the first root of J0. The actual solution of the problem lies higher than that, but at the very least k1*a > 2.4048. Then from

k1 = sqrt(n1^2*k^2 - kz.^2);

if follows that n1*k*a > k1*a > 2.4048

and since w = c*k there is a cutoff frequency. This exact cutoff frequency turns out to be larger, but this demo does does show the existence of a cutoff frequency.

For the transmitted power, you have to integrate the Poynting vector ExH over the area of the fiber. (Let's say you don't calculate the transmitted power in the cladding since it is not really available for use). This means doing an integral of

(amplitudes)*J1(k1*r)^2 * 2*pi*r dr which is over the area of the fiber. Fortunately that integral is doable analytically and Matlab even has it. Try

syms r k1 a

int(2*pi*r*besselj(1,k1*r)^2,0,a)

David Goodmanson
on 3 Mar 2020

Hi Carolina,

thanks very much for sending for the chapter from the book. Thanks also for the compliment, but I think I have reached the end of my knowledge of optical fibers. The solution I gave was related to m = 1 on p. 44, although his angular dependence is different. My impression is that he is playing fast and loose with some of the equations and I don't know that his boundary conditions are even correct. And he does not show any expressions for what the E and B fields actually are. Balanced against this impression is the fact that it's a Springer book in its second edition, so how incorrect could it be?

Everyone agrees that there is a transverse 01 mode with no lower frequency cutoff, which unfortunately my answer does not apply to (although I think the answer has a correct description of certain other modes). I have not yet found an expression for the fields for 01, so it looks like a trip to the university library is on the agenda.

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 8 Comments

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798542

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798542

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798565

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798565

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798566

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798566

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798572

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798572

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798573

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798573

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798578

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798578

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798649

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798649

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798654

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/506014-te-mode-trascendental-equation-problem#comment_798654

Sign in to comment.