MATLAB Answers

TE mode trascendental equation problem

31 views (last 30 days)
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


Show 5 older comments
Star Strider
Star Strider on 17 Feb 2020
I usually use fsolve if fzero has problems:
X = fsolve(F,x0)
however even if I only use one of the ‘x0’ values, this throws:
Error using trustnleqn (line 28)
Objective function is returning undefined values at initial point. FSOLVE cannot
So your ‘x0’ is apparently not appropriate. I am not familiar with what you are doing, so I cannot help you choose more appropriate values that are consistent with what you want to solve. You need to experiment with those.
Carolina  Rickenstorff
Carolina Rickenstorff on 17 Feb 2020
It is very weird that there are books as Fiber Optics (F. Mitschke) where the same equation is written slightly different and I have no problem solving it.
%%%%Te mode
syms u w
%solution=solve((1/u).*besselj(0,u)./besselj(1,u)==(1/w).*besselk(0,w)./besselk(1,w), w);
solution=solve(u.*besselj(1,u)./besselj(0,u)==w.*besselk(1,w)./besselk(0,w), w);
Star Strider
Star Strider on 17 Feb 2020
In that situation, use the version you can solve!

Sign in to comment.

Accepted Answer

David Goodmanson
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;
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)


Show 7 older comments
Carolina  Rickenstorff
Carolina Rickenstorff on 25 Feb 2020
Dear David,
Here is my final program for the TE mode fundamental mode (last kzroot). I returned to Hankel functions despite they produce little artifacts. The field plots have smooth transitions now. Thanks to point out the fact that the biggest kzroot is the one that produces the fundamental mode. In my particular case the biggest radial order of the TE mode is 5 isn´t?.
clear all
close all
%%%%Optical parameters
c = 3e8;
k = 2*pi/1.55e-6;
w = k.*c;
mu0 = 1.26e-6;
%%%%Fiber parameters inside the nucleus (subscript 1) and outside the nucleus (subscript 2)
mu1 = 1; n1 = 1.4;
mu2 = 1; n2 = 1.33;
%%%%Fiber radius
a = 10e-6;
%%%%The variable to solve is the propagation constant kz
%%%%kz range, n2*k < kz < n1*k
kzmax = n1*k;
kzmin = n2*k;
kz = linspace(kzmin,kzmax,500000);
kz(1) = [];
kz(end) = [];
%%%%TE equation mode in an optical fiber (D. Marcuse)
%%%%LHS evaluated at r=a
%%%%k1 variable (D. Marcuse)
k1 = sqrt(n1^2*k^2 - kz.^2);
%%%%Let's define the ratio
C = 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;
%%%%RHS evaluated at r=a
%%%%k2 variable (D. Marcuse)
k2 = sqrt(kz.^2-n2^2*k^2);
%%%%Let´s define the ratio
D = -1i.*besselh(1,1,1i.*k2*a)./(k2.*besselh(0,1,1i.*k2*a));
%%%%quiet down the plot y values
Dm = mean(abs(D));
D(abs(D)>3*Dm) = nan;
%%%%Trying to find the intersection points of both equations (kzroots) more or less automatically
%%%%Plot the functions inside and outside the nucleus
legend('Inside nucleus','Outside nucleus')
xlabel('kz (1/m)')
hold on
hold off
grid on
%%%%Verification of kzroots values in both sides of equation
kzroot = kzroots(end); % 5.4617301149e6;
k1root = sqrt(n1^2*k^2 - kzroot^2);
k2root = sqrt(kzroot^2-n2^2*k^2);
%%%%Both LHS/RHS are real, DD has a very small imaginary artifact e-23
%%%%For TE mode only coefficients b, d and fields Hz, Hr, Ephi are different to 0
%%%%I choose b=1, then
b=1; %%%%real
d=b.*besselj(0,k1root*a)./besselh(0,1,1i.*k2root*a); %%%%imaginary
%%%%Finally Hz, Hr y Ephi are
r=linspace(0, N-1, N);
%%%%Inside nucleus
Hz1 = b .*besselj(0,k1root.*r(1:N/2+1)); %%%%real
Hr1 = 1i.*b.*(kzroot./k1root).*besselj(1,k1root.*r(1:N/2+1)); %%%%imaginary
Ephi1 =-1i.*b.*(w.*mu0./k1root).*besselj(1,k1root.*r(1:N/2+1)); %%%%imaginary
%%%%Outside nucleus
Hz2 = d .*besselh(0,1,1i.*k2root.*r(N/2+2:end)); %%%%real
Hr2 = d.*(kzroot./k2root).*besselh(1,1,1i.*k2root.*r(N/2+2:end)); %%%%imaginary
Ephi2 = -d.*(w.*mu0./k2root).*besselh(1,1,1i.*k2root.*r(N/2+2:end)); %%%%imaginary
Hz=[Hz1 Hz2];
Hr=[Hr1 Hr2];
Ephi=[Ephi1 Ephi2];
legend('Hz'); ylabel('A/m'); grid on
legend('Hr'); ylabel('A/m'); grid on
xlabel('r/a, a=nucleus radius'); ylabel('V/m')
legend('Ephi'); grid on
With these data I plan to get the total intensity I=abs(Er)^2+abs(Ephi)^2+abs(Ez)^2 and the Poynting vector S=real(ExH*) for to calculate the gradient force Fg, scatering force Fs and absorbing force Fa excerted in a metallic nanoparticle along the radial direction.
Sr = r.*real(Ephi.*conj(Hz)-Ez.*conj(Hphi));
Sphi= r.*real(Ez.*conj(Hr)+Er.*conj(Hz));
Sz = r.*real(Er.*conj(Hphi)-Ephi.*conj(Hr));
Best regards.
Carolina R
Walter Roberson
Walter Roberson on 25 Feb 2020
Thank you David for the nice breakdown!
David Goodmanson
David Goodmanson on 27 Feb 2020
Thank you, Walter, for your comment.
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.

More Answers (2)

Carolina  Rickenstorff
Carolina Rickenstorff on 28 Feb 2020
Edited: Carolina Rickenstorff on 28 Feb 2020
Dear David,
I read that the fundamental modes TE0 and TM0 ideally have no cutoff and they are always present in an optical fibre. One excludes the other or both modes can exist simultaneusly in the fibre? The one that can transport information is TE0 mode or both can be used to codify information?
Another question please. If I want to simulate a 100mW beam in my example then I have to play with factors B, D until the magnitude of Ephi using the last kzroot is equal to sqrt(100e-3) V/m?
I appreciate a lot your help. Thanks for heeding me so far.


Sign in to comment.

David Goodmanson
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


Carolina  Rickenstorff
Carolina Rickenstorff on 2 Mar 2020
Hello David,
I read the issue of cutoff frequency in the book Fiber Optics of Fedor Mitsche Ed. Springer (2016). I attach you the charper 3 where at end of section 3.10 and end of section 3.12 it says so.
Thank you very much for your guidance. I'm glad to have such a good teacher.
Best regards.
Carolina R
David Goodmanson
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.