# Calculating reflection of light when the refractive index changes continuously

Suppose you have 2 materials, one with refractive index $$n_1$$ and the other with refractive index $$n_2$$, and a plane-wave coming from the first material hits the interface with an incident angle of $$0^\circ$$.

Fresnel tells us that the reflected power will be $$r=\frac{n_1-n_2}{n_1+n_2} \Rightarrow R=\left(\frac{n_1-n_2}{n_1+n_2}\right)^2$$

Now, if you have a set of materials each with thickness $$d_i$$ and refractive index $$n_i$$, you can use the simple formula above and interference to calculate the net reflection and net transmission (for example, one can multiply the matrices associated with the transmission and reflection of each interface and the matrices propagating in each material).

I'm having trouble with a similar but different problem. I have this optical fiber with refractive index $$n_0$$, and at some point along the fiber the refractive index changes periodically and continuously: $$n(x)=n_0+\epsilon\sin\left(\frac{2\pi}{d}x\right), \; \epsilon\ll1$$

where the wavelength $$\lambda$$ is not negligible compared to the period $$d$$. After $$N$$ periods, the refractive index returns to $$n_0$$.

The question is: How do I calculate the net reflection and transmission of such grating? The refractive index varies continuously, not in discrete steps for which I can use Fresnel's equations.

Thanks a lot!

• You might consider using variational calculus. This would certainly allow you to solve the ray tracing problem for continuous index of refraction; whether or not it is useful for the problem you have in mind is not clear, but it seems like an obvious thing to look up. Jan 21, 2017 at 15:48
• You may find section 3 of this paper a good place to start. They are trying to solve the inverse problem - but start by solving a problem similar to the one you are asking about. Jan 21, 2017 at 15:56
• @ZeroTheHero that's not my problem :/ . My problem is completely one dimensional. Jan 21, 2017 at 18:58
• @Floris thank you but unfortunately it didn't have the answer, do you have any tips for searching relative info? Jan 21, 2017 at 19:02
• @honeste_vivere I am aware of this equation :) Sep 28, 2021 at 16:02

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  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.

Electromagnetic wave equation for $$\vec E$$ in an isotropic medium with $$\mu=\mathrm{const}$$ and no free currents nor charges reads

$$-\nabla^2\vec E+\mu\varepsilon\partial_t^2\vec E=-\nabla\left(\vec E\cdot\frac{\nabla\varepsilon}{\varepsilon}\right).\tag1$$

(It can be derived from Maxwell's equations similarly to how Wikipedia does it for microscopic wave equations.)

If $$\varepsilon$$ only varies along $$\hat x$$, and a plane wave incident from empty space also propagates along $$\hat x$$, then further evolution should also be symmetric along $$\hat y$$ and $$\hat z$$. Thus, the dot product on the RHS of $$(1)$$ vanishes, and we get the usual wave equation

$$-\nabla^2\vec E+\mu\varepsilon\partial_t^2\vec E=0,\tag2$$

from which we can obtain the equation for eigenmodes by setting $$\vec E=\vec E_0e^{i\omega t}$$:

$$-\frac{c^2}{n^2}\nabla^2\vec E_0=\omega^2\vec E_0.\tag3$$

For a wave polarized along $$\hat z$$, we have $$E_x=E_y\equiv0$$, and we are left with a single ODE for $$E_z$$:

$$-\frac{c^2}{n^2}\frac{\mathrm d^2 E_{0z}}{\mathrm d x^2}=\omega^2 E_{0z}.\tag4$$

This one is generally not solvable analytically. With a periodic medium, though, the domain for numerical computation can be reduced to a single period by Bloch theorem (for some frequencies your Bloch wavenumber may become imaginary, that's OK for band-gap solutions). Then you can match your incident & reflected plane waves with this set of solutions by requiring appropriate interface conditions: continuity of $$E_{0z}$$ and its first derivative (even if $$n(x)$$ is discontinuous at the interface); these conditions can be derived directly from $$(4)$$.

After you do the matching described above, your reflected wave will have amplitude $$r$$, from which your reflection coefficient is $$R=|r|^2$$.

• Block theorem will work if I had infinite periods. When I use infinite periods I will have for some frequencies R=1 and for others R=0. Since I'm interested in the case of finite $N$, the reflection coefficient will change continuously - that's exactly what I'm looking for :) . Maybe I didn't understand your answer, can you please expand how this is compatible with finite periods? Sep 27, 2021 at 9:48
• @OfekGillon remember how we solve the problem with finite-depth potential barrier in quantum mechanics? We use piecewise wavefunction, where each piece is a superposition of solutions for a free particle. Here is exactly the same situation: you find general solutions for infinite space in each subdomain of your problem (vacuum, fiber, vacuum), and then find the coefficients that let you satisfy the interface conditions. Here you just need to find the two linearly independent solutions for the middle part numerically. Sep 27, 2021 at 12:40

I will outline two approaches. Both involve numerical solution.

Method 1. Use the transmission matrix concept, but for a very large number of small steps through the material. For each small step $$\delta z$$ find the matrix based on the local change of index, and multiply it onto the matrix you have so far.

Method 2. Solution of wave equation. You have $$\frac{\partial^2 E}{\partial t^2} - \frac{c^2}{n^2} \frac{\partial^2 E}{\partial z^2} = 0$$ (I think that's ok even when $$n$$ is a function of $$z$$, but to be honest I am not completely sure and you would need to look into this before proceeding. Born and Wolf is a good resource.) Assuming this differential equation is the right one (or replacing it with another one which is right), then all you need to do is use a standard numerical partial differential equation solver. Matlab has them, so does Python, and many other languages.

I wrote this answer just off the top of my head. I would not trust it until I had done some further working out, and I would probably want to set up both methods and thus use them to check each other.