I read in a book about optical fibers that the different spectral components of a light pulse transmitted in the fiber propagate with different velocities due to a
wavelength dependent refractive index. Can someone explain that? Why is that silica refractive index depends on the wavelength/frequency of the wave?
[Physics] Wavelength-dependent refractive index
electromagnetic-radiationopticsrefractionwavelength
Related Solutions
For an intuitive answer to your first question, please see both my answer and and that of @theSkinEffect about the Effective Refractive Index. However, a more complete answer concerning the modes of a fiber is found in these lecture notes.
For your second question, what is the dimensions of the waveguide (i.e. diameter)? What is the wavelength of interest? All of these factors will determine how many and what the effective indices are for your problem. For the sake of demonstration, I will assume you are using Corning SMF-28 which has a radius of 8.3$\mu m$ and that your wavelength of interest is 1.55$\mu m$ which are very typical for the present optical communications in the telecom industry.
To answer your question about how to calculate the effective indices, please download the free statistical software package R and drop the following code into it. I normally use Python, but since this program numerically calculates the $\beta_{pm}$, I liked the numerical solver better in R. Before you use this program, please understand better the theory of modes in a cyllindrical waveguide such as fiber because there are some assumptions made in the equations I used for the characteristics equations and you need to understand and change the order of the Bessel functions to get the correct answer.
# Written in R by Ansebbian0
# Revisions have been made to correct the mode angle
# Please be aware that the mode angle is give from the propogation direction,
# NOT from the incident angle. This can be changed by using asin instead of acos
library(pracma)
deg2rad = function(deg) {
return((pi * deg) / 180)
}
rad2deg = function(rad) {
return((180 * rad) / pi)
}
# Characteristic Equation for TE mode beta0
TE_eq <- function(beta0,k0,n1,n2,h){
kapppa <- sqrt( (n2*k0)^2 - beta0^2) # lateral wavevector
kapppamax <- sqrt(k0**2 * (n2**2 - n1**2))
gammma <- sqrt(kapppamax**2 - kapppa**2)
te0 <- (besselJ(kapppa*h, 1) / (kapppa*besselJ(kapppa*h, 0))) + (besselK(gammma*h, 1, expon.scaled=FALSE) / (gammma*besselK(gammma*h, 0, expon.scaled=FALSE)))
return(te0)
}
# Characteristic Equation for TM mode beta0
TM_eq <- function(beta0,k0,n1,n2,h){
kapppa <- sqrt( (n2*k0)^2 - beta0^2)
kapppamax <- sqrt(k0**2 * (n2**2 - n1**2))
gammma <- sqrt(kapppamax**2 - kapppa**2)
tm0 <- ((besselJ(kapppa*h, 1) / (kapppa*besselJ(kapppa*h, 0)))*(k0*n2)**2) + ((besselK(gammma*h, 1, expon.scaled=FALSE) / (gammma*besselK(gammma*h, 0, expon.scaled=FALSE)))*(k0*n1)**2)
return(tm0)
}
# calculation of betas (set seq to small number for high resolution)
beta0_array <- function(lambda,n1,n2,h){
k0 <<- 2*pi/lambda # wavenumber in free space
beta0 <- seq(n1*k0, n2*k0, 1000) # linspace of beta0s from smallest to largest(see pollock)
beta0 <- beta0[-length(beta0)] # removes the last element of beta0
return(beta0)
}
# calculates the crossing points where there are possible zeros
pol_zeros <- function(polarization, beta0){
intervals <- (polarization>=0)-(polarization<0) # gives a matrix of 1 (te0>=0) and -1 (te0<0)
differences = diff(intervals) # returns the difference between each variable
izeros <- which(differences<0) # a neg. return from diff will mean a crossing from + to -
X0 <- cbind(beta0[izeros],beta0[izeros+1]) # makes a Nx2 matrix of the two crossing variables
return(X0)
}
# Numerically solves the characteristic equation [n_eff = beta0_m / k0]
modes <- function(lambda, h, n1, n2){
beta0 <- beta0_array(lambda,n1,n2,h)
te0 <- TE_eq(beta0,k0,n1,n2,h) # calls the characteristic equation for beta0 of te equation
tm0 <- TM_eq(beta0,k0,n1,n2,h) # calls the characteristic equation for beta0 of tm equation
# n_eff for TE and TM numerical calculation
te_Xs <- pol_zeros(te0, beta0)
nzeros_te <- dim(te_Xs)[1]
tm_Xs <- pol_zeros(tm0, beta0)
nzeros_tm <- dim(tm_Xs)[1] # finds how many rows of crossings there were
n_eff<-matrix(data=NA, nrow=max(nzeros_te, nzeros_tm), ncol=10) # initializes a matrix for the n_eff
colnames(n_eff) <- c('te_ne_eff', 'te_mode_beta', 'prop_dir_te_mode_angle', 'tm_ne_eff', 'tm_mode_beta', 'prop_dir_tm_mode_angle','te_kappa', 'tm_kappa', 'te_mode', 'tm_mode')
for(i in 1:nzeros_te){ # feeds crossing pairs into a numerical root solver
n_eff[i,1] <- (fzero(function(x){TE_eq(x,k0,n1,n2,h)}, te_Xs[i,])$x)/k0 # fzero is a root solver in pracma
n_eff[i,2] <- n_eff[i,1]*k0
n_eff[i,3] <- rad2deg(acos(n_eff[i,1]/n2))
n_eff[i,7] <- sqrt( (n2*k0)^2 - n_eff[i,2]^2)
n_eff[i,9] <- (nzeros_te + 1) - i
}
for(i in 1:nzeros_tm){ # feeds crossing pairs into a numerical root solver
n_eff[i,4] <- (fzero(function(x){TM_eq(x,k0,n1,n2,h)}, tm_Xs[i,])$x)/k0 # fzero is a root solver in pracma
n_eff[i,5] <- n_eff[i,4]*k0
n_eff[i,6] <- rad2deg(acos(n_eff[i,4]/n2))
n_eff[i,8] <- sqrt( (n2*k0)^2 - n_eff[i,5]^2)
n_eff[i,10] <- (nzeros_tm + 1) - i
}
return(data.frame(n_eff)) #, nTM, TEparam, TMparam)
}
modes(lambda=1.55e-6, h = 8.3e-6, n1=1.4271, n2=1.4446)
This is the output:
The output should be self explanatory, but if you have questions, let me know and I will tell you about how to use it. This program is not valid for the LP modes because those are composite modes. This is only valid for TE and TM. Look at Buck, Fundamentals of Optical Fibers to find the equation for the LP modes and then just replace that into the TE_eq or TM_eq function and you'll get the results.
PS: One small note. Before the above program works correctly, you need to type this in R:
install.packages("pracma")
This is the numerical solvers library which you need.
Suppose we define the following $\zeta = \ln{\varepsilon}$ and $\xi = \ln{\mu}$, where $\varepsilon$ and $\mu$ are the permittivity and permeability, respectively. In a system with no sources (i.e., $\mathbf{j} = 0$ and $\rho_{c} = 0$), then we know that $\nabla \cdot \mathbf{D} = 0$, where $\mathbf{D} = \varepsilon \ \mathbf{E}$ and $\mathbf{B} = \mu \ \mathbf{H}$. After a little vector calculus we can show that: $$ \nabla \cdot \mathbf{E} = - \nabla \zeta \cdot \mathbf{E} \tag{0} $$ Using this and some manipulation of Faraday's law and Ampêre's law, we can show that the general differential equation in terms of electric fields only is given by: $$ \left( \mu \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \nabla^{2} \right) \mathbf{E} = \left( \mathbf{E} \cdot \nabla \right) \nabla \zeta + \left( \nabla \zeta \cdot \nabla \right) \mathbf{E} + \nabla \left( \zeta + \xi \right) \times \left( \nabla \times \mathbf{E} \right) \tag{1} $$
We can get a tiny amount of reprieve from this by assuming that the permeability is that of free space, i.e., $\nabla \xi = 0$. If we further argue that the only direction in which gradients matter is along $\hat{x}$ and that the incident wave vector, $\mathbf{k}$, is parallel to this, then we can further simplify Equation 1 to: $$ \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) \mathbf{E} = \left( \mathbf{E} \cdot \nabla \right) \zeta' \hat{x} + \left( \zeta' \frac{ \partial }{ \partial x } \right) \mathbf{E} + \zeta' \hat{x} \times \left( \nabla \times \mathbf{E} \right) \tag{2} $$ where $\zeta' = \tfrac{ \partial \zeta }{ \partial x }$.
After some more manipulation, we can break this up into components to show that: $$ \begin{align} \text{x : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{x} & = E_{x} \zeta'' + \zeta' \frac{ \partial E_{x} }{ \partial x } \tag{2a} \\ \text{y : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{y} & = 0 \tag{2b} \\ \text{z : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{z} & = 0 \tag{2c} \end{align} $$
In the limit where the incident wave is entirely transverse, then $E_{x} = 0$ and the x-component (Equation 2a) is entirely zero.
Next you assume that $\mathbf{E} = \mathbf{E}_{o}\left( x \right) e^{i \omega t}$, where $\omega$ is the frequency of the incident wave. Then there will be incident, reflected, and transmitted contributions to the total field at any given point (well the transmitted is always zero in the first medium, of course). Any incident and transmitted contributions with have $\mathbf{k} \cdot \mathbf{x} > 0$ while reflected waves will satisfy $\mathbf{k} \cdot \mathbf{x} < 0$. You define the ratio of the reflected-to-incident fields (well impedances would be more appropriate) to get the coefficient of reflection.
Simpler Approach
A much simpler approach is to know where to look for the answer to these types of questions. As I mentioned in the comments, there has been a ton of work on this very topic (i.e., spatially dependent index of refraction) done for the ionosphere. If we look at, for instance, Roettger [1980] we find a nice, convenient equation for the reflection coefficient, $R$, as a function of the index of refraction, given by:
$$
R = \int \ dx \ \frac{ 1 }{ 2 \ n\left( x \right) } \frac{ \partial n\left( x \right) }{ \partial x } \ e^{-i \ k \ x} \tag{3}
$$
There is no analytical expression for $R$ for your specific index of refraction. However, numerical integration is not difficult if one knows the values for $d$ and $\epsilon$. Note that if we do a Taylor expansion for small $\epsilon$, then the integrand (not including the exponential) is proportional to cosine, to first order in $\epsilon$ (cosine times sine if we go to second order).
References
- Gossard, E.E. "Refractive index variance and its height distribution in different air masses," Radio Sci. 12(1), pp. 89-105, doi:10.1029/RS012i001p00089, 1977.
- Roettger, J. "Reflection and scattering of VHF radar signals from atmospheric refractivity structures," Radio Sci. 15(2), pp. 259-276, doi:10.1029/RS015i002p00259, 1980.
- Roettger, J. and C.H. Liu "Partial reflection and scattering of VHF radar signals from the clear atmosphere," Geophys. Res. Lett. 5(5), pp. 357-360, doi:10.1029/GL005i005p00357, 1978.
Best Answer
The fundamental reason for the wavelength dependance of refractive index ($n$), in fact the fundamental description of refraction itself, is the domain of quantum field theory and is beyond my understanding. Hopefully somebody else can provide an answer on that subject.
However, I can state that it isn't just silica that has a wavelength dependent $n$. In fact, every material has some wavelength dependence, and this property is called dispersion. In optical materials, the dispersion curve is very well approximated by the Sellmeier Equation: $$ n^2(\lambda) = 1 + \sum_k \frac{B_k \lambda^2}{\lambda^2 - C_k} $$
usually taken to $k=3$, where $B_k$ and $C_k$ are measured experimentally. As far as I know this equation is not derived from theory; it is completely empirical.