# Reflection and transmission of waves in a medium of nonuniform density

{???}
Hello all,

Apologies in advance for the text-wall; this is a rather involved question.

I am trying to compute the effective transmission coefficient for a medium of non-uniform refractive index. For simplicity I am assuming the slab has thickness ##d##, that ##n(0)=1##, and that ##n(d)=n## (where ##n## is a constant).

I am familiar with the normal incidence reflection and transmission coefficients, ##R=(n_1-n_2)/(n_1+n_2)##, for light incident from a region of index ##n_1## upon a region of index ##n_2##. I also know how to derive these coefficients from the continuity and smoothness conditions of the electric and magnetic fields at the boundary.

One can equally well consider the problem of a wave on a string of nonuniform density. My question is largely of a mathematical nature.

At first, there does not appear to be a clear problem with applying these same conditions to "infinitesimal discontinuities," i.e. ##n(x)## to the left of the boundary and ##n(x)+dn(x)=n(x+dx)## to the right of the boundary. Now our reflection coefficient ##dR(x)## is infinitesimal, and the transmission coefficient ##T(x)=1+dR(x)##.

I'm aware that the principles leading to the derivation of the wave equation do make assumptions about the uniformity of the medium - this is most clearly illustrated by the wave-on-a-string example. Assuming a non-uniform string breaks the derivation. However, the fact of the matter is an answer exists, and my current approach (of the three I've tried) seems to be the most sensible.

I treat the medium like an infinite sequence of media with "infinitesimal discontinuities" that have infinitesimal reflection coefficients, and transmission coefficients which differ from unity by an infinitesimal.

Much like integrals are the limits of sums of terms whose values tend to zero as the number of terms tends to infinity, what I have for the effective transmission coefficient (assuming no back-and-forth reflections) is a product of factors whose values tend to unity as the number of factors tends to infinity. (This is akin to an exponential of an integral of logarithms.)

This leads to an effective transmission coefficient of ##1/\sqrt{n}## for the non-uniform medium, times ##2n/(n+1)## for the final transmission.

My question is about corrections to this. If the incident wave is reflected back at ##x_1##, then reflected forward at ##x_2##, there will be a doubly-infinitesimal contribution. We can choose any pair of coordinates such that ##0<x_1<d## and ##0<x_2<x_1##, meaning we will have to integrate over a two-dimensional region. This gives a finite and non-infinitesimal correction to the above.

Hence my current predicament. Here are the approaches I've attempted, in decreasing order of preference:

1. Using a "hopscotch"-style boundary conditions approach, I can generalize the relations ##E_t(d+x)=TE_i(d+(n_t/n_i)x)## and ##E_r(d-x)=RE_i(d+x)##. This gives me some nasty recursion relations in the general case, that simplify in the event of infinitesimal layers. However, they simplify to the same expression used to get the na\"ive estimate above, without considering multiple reflections. I would like to use this method, if possible, but I would greatly appreciate some help with figuring out how to get the NLO contribution.
2. Computing the multiple integrals analytically is impossible except for a very specific parametrization of the refractive index, ##n(x)=(1-x/L)^{-1}## for some choice of length scale ##L##. Even for this choice of refractive index, this method requires integration by parts and leads to a dynamic mess at every order, and while I most trust the results it gives me, I am least confident in my ability to generalize to the correction for ##2n## reflections.
3. Using the BCs that are used to derive the reflection and transmission coefficients in the first place, I get equations that don't even agree with my verified results for a uniform slab of thickness ##d## (for which ##T_\mathrm{eff}=T_i^*T_i(1-2R_0^2\cos(nkd)+R_0^4)^{-\frac{1}{2}}##, where ##R_0## is the reflection coefficient at the interface, ##n## the refractive index of the slab, and ##k## the wavenumber of the incident light). Moreover, when I attempt to use this method to compute infinitely many contributions, all I get is a telescoping product, the result of which is ##T=1##. This method would appear to require the least work, but seeing as it does not reproduce either transmission coefficient above, I am reticent to trust it.

I don't see another way forward at this point. Could anybody please look over the work I have done? I really just want to know: Whence comes the NLO contribution?

Thanks,
QM

[Post edited by the Mentors at the request of the OP]

Last edited by a moderator:
• Delta2

{???}
There is something called "complex index of refraction" which accounts for decay:
Hi A.T.,

I'm sorry but I don't think this addresses my question. I am familiar with complex refractive indices as they are used for total internal reflection. What I am dealing with is a material with a real, but non-constant refractive index. My assumptions are as follows:

1. We can treat the medium of non-constant refractive indices as a sequence of media of constant refractive indices, the number of which increases and thickness decreases as our approximation gets closer to the physical scenario. This is exactly the same assumption used for Riemann sums and integration, and I don't see any arc-length- or discontinuity-related problems that would cause convergence issues e.g. ##\pi=4##. The sequence of media does indeed approach the nonuniform medium.
2. At the interfaces, there is an infinitesimal difference in refractive index. This leads to an infinitesimal reflection coefficient, and a transmission coefficient that differs from 1 by an infinitesimal.

Are there problems with my assumptions? Because, if not, the three methods I outline should give the same result. (There is another approach, which is simply to sum over the countably infinite number of propagation modes that contribute, and which also agrees with my verified result, but I want to be sure there are no viable alternatives before resorting to this method.)

So, if there are no problems with my assumptions, are there any problems with my methods? I was so certain Method #1 would give a sensible result, but there is nothing involving multiple reflections. I know that modes with multiple reflections MUST be considered because I can calculate them using Method #2 and they are nonzero (albeit suppressed relative to the LO term).

I'm convinced that the effective transmission coefficient ##T## satisfies ##\frac{2\sqrt{n}}{n+1}<T<1##. I just want to know if I have the right idea to find it.

Homework Helper
2022 Award
I am confused as to the physical setup. Is this effectively a stack of very thin flat panes of glass of differing real indices of refraction?
If the index is real then energy of wave is strictly conserved.

Last edited:
{???}
I am confused as to the physical setup. Is this effectively a stack of very thin flat panes of glass of differing real indices of refraction?
If the index is real then energy of wave is strictly conserved.
Hi hutchphd,

Your effective model is correct. I agree that the energy stored in waves is conserved. However, there is a non-trivial division of how much energy is in the transmitted wave versus how much is in the reflected wave. This makes the effective transmission and reflection coefficients nontrivial, although again we must have energy conservation when we consider them together.

All right, you guys have convinced me - this question seems to be confusing.

Instead of optics, let's consider classical mechanics. Let's say there's a wave of wavenumber ##k## on a rope of mass density ##\lambda_0## which encounters a segment of mass density ##\lambda(x)=\lambda_0+x\Delta\lambda=\lambda_0+x(\lambda_1-\lambda_0)/d## beginning at ##x=0## and ending at ##x=d##. For ##x>d## the mass density is again ##\lambda_0##. What is the ratio of the transmitted amplitude to the incident amplitude? Consider as well the infinitesimal contributions which are reflected back one or more times (there are infinitely many of them).

There is https://www.researchgate.net/profile/Fernando_Dall_Agnol/publication/256096820_Wave_propagation_in_a_non-uniform_string/links/00463521beb033eebf000000/Wave-propagation-in-a-non-uniform-string.pdf discussing this exact scenario but I don't understand their approach. I would like to treat this string of nonuniform density as many short segments whose densities differ infinitesimally, and look at reflection and transmission from each "interface" of these "segments." I have ~4 approaches to consider:

1. Go directly to the continuum limit. I take ##T=1\pm dR## and compute the infinite product of transmission coefficients over ##0<x<d##. Then we have the LO contribution. The NLO contribution has two positions defining the reflections, which we integrate over 2D space to get a nonzero, finite answer. Repeat an infinite number of times.
2. Play "hopscotch" and treat the left-/right-propagating wave at a given position as the sum of one reflected and one transmitted wave. This approach avoids computing infinite sums. I have worked out the general case for however many nodes separated by an arbitrary distance, but taking the limit as the separation distance becomes infinitesimal, I don't see any way to get the NLO contribution discussed in #1. Honestly, I would definitely prefer this method if I can just figure out what's going wrong with my differentials. Since I want an infinite product, I'm throwing away all terms which have the product of two or more differentials in them. If someone could take a look at my earlier attachment, that would be great.
3. Just use continuity/smoothness across each boundary. This is how the coefficients are traditionally derived. However, trying this does not seem to agree with my earlier results.
4. Compute the infinite sums over propagation modes explicitly. I've done this for 2-4 interfaces and checked them for consistency. However, with an infinite number of interfaces, requiring infinite sums of infinite sums of ..., this looks pretty daunting. I'm sure it can be done, but even by my standards this seems pretty nasty.

Would it help elicit more feedback if I actually did all three of these methods to completion (not just to the point where I see a disagreement/problem)? I just want to know what I'm doing wrong, such that essential features are not present.

Thanks,
QM

Gold Member
I'm not sure I understand your approach, and the document you attached has tons of equations and not enough words to understand what you are doing. In any case, I see two straightforward ways to solve this using basic tools that you probably already know. The first is analytical and more complicated due to special functions, and the second is numerical and pretty simple to setup. For me, the numerical approach would be much faster...

If you want an analytical solution, then as part of the solution you will need to solve a differential equation. If we assume ##e^{i\omega t}## time dependence (and suppress that time dependence) then the transverse Electric field for a wave in an inhomogeneous dielectric satisfies
$$\nabla^2 \mathbf{E} + k_0^2 \epsilon_r(\mathbf{r})\mathbf{E} = 0$$
where ##k_0 = \omega / c## is the free-space wavenumber and ##\epsilon_r(\mathbf{r})## is the relative permittivity, which as usual is related to the index of refraction ##n## by ## \epsilon_r(\mathbf{r}) = n^2(\mathbf{r})##. You should derive this if you haven't seen it before - it is just a few lines of algebra from Maxwell's equations. You have a simple case, where your variation is only in one direction (let's say the positive ##z## direction) and you can pick an arbitrary polarization for the field. For example, if we only have an ##x## component to the electric field, then the equation is
$$\frac{d^2 E}{dz^2} +k_0^2 \epsilon_r(z) E =0$$
Of course this only holds in the slab (##0\leq z \leq d##). Below the slab (##z<0##) you have your incident wave and the reflected wave, and above the slab (##z>d##) you just have a transmitted wave. You still have boundary conditions at ##z=0## and ##z=d##, but for the field inside the slab you will need to solve the differential equation (which should yield two independent solutions). Of course, if ##E## is in the x direction, then ##H## is in the ##y## direction. If ##\epsilon_r(z)## is linear then the solution will be in terms of Airy functions (which are another way to write Bessel functions of order ##1/3##). See https://dlmf.nist.gov for info on special functions. If ##\epsilon_r## is a different profile then it will be other functions. So this gets messy in a hurry, just because of the special functions. Your final reflection coefficient will be in terms of these special functions.

If you are content with numerical solutions then you numerically solve the case for an N-layer slab, where each slab has a constant index of refraction. Basically you write the form of the fields for each slab, then apply boundary conditions at each boundary. You will then have 2 boundary conditions for each of ##N+1## boundaries, so will solve a ##2N+2 \, \times 2N+2## system of linear equations. Numerically this is easy; analytically it is not. You can start with N=4, or some other small number, compute the answer, then double N and re-compute, etc. The result should converge to the correct value pretty rapidly. One nice thing about this approach is that it is trivial to try different profiles for ##\epsilon_r(z)## once you have the code written. For the analytical approach above you have to solve a different differential equation for each profile.

Jason

• Delta2 and hutchphd
Homework Helper
2022 Award
I agree wholeheartedly with @jasonRF. Sorry, I am also unable to attack your calculation

A few thoughts:
1. This is a pretty standard scattering problem which (as mentioned) has no happy analytic solution for even well-behaved "scattering potentials" .
2. There are a host of approximate methods (WKB, eikonal, perturbation expansions). I feel certain that you are reproducing certain aspects of them
3. It is not clear to me that the slab method will be well-behaved, particularly because of possible "resonances" for layers of wavelength commensurate widths
4. One test measure of the calculation is unitarity (conservation of energy).
5. A numerical method that might work (this may be your "hopscotch "...not clear to me) would be a Monte Carlo approach: tracing many pinball-like paths and adding them (with phases and weights) when they get to each edge. Might work well. Maybe a bad idea...
6. https://en.wikipedia.org/wiki/Eikonal_equation

The nice thing about this problem is that it is well-framed, useful in many realms of physics, and nontrivial. But I think you know that.

• vanhees71, Andy Resnick and jasonRF
Gold Member
I like hutchphd's list. He is absolutely correct in item 3, and I should have known better than to arbitrarily pick N=4 or some small number. Each slab should be a small fraction of a wavelength (in that slab) in order to avoid those effects.

{???}
Hi Jason and Hutch,

Thank you both for your in-depth responses.
I'm not sure I understand your approach, and the document you attached has tons of equations and not enough words to understand what you are doing. In any case, I see two straightforward ways to solve this using basic tools that you probably already know. The first is analytical and more complicated due to special functions, and the second is numerical and pretty simple to setup. For me, the numerical approach would be much faster...

If you want an analytical solution, then as part of the solution you will need to solve a differential equation. If we assume ##e^{i\omega t}## time dependence (and suppress that time dependence) then the transverse Electric field for a wave in an inhomogeneous dielectric satisfies
$$\nabla^2 \mathbf{E} + k_0^2 \epsilon_r(\mathbf{r})\mathbf{E} = 0$$
where ##k_0 = \omega / c## is the free-space wavenumber and ##\epsilon_r(\mathbf{r})## is the relative permittivity, which as usual is related to the index of refraction ##n## by ## \epsilon_r(\mathbf{r}) = n^2(\mathbf{r})##. You should derive this if you haven't seen it before - it is just a few lines of algebra from Maxwell's equations. You have a simple case, where your variation is only in one direction (let's say the positive ##z## direction) and you can pick an arbitrary polarization for the field. For example, if we only have an ##x## component to the electric field, then the equation is
$$\frac{d^2 E}{dz^2} +k_0^2 \epsilon_r(z) E =0$$
Of course this only holds in the slab (##0\leq z \leq d##). Below the slab (##z<0##) you have your incident wave and the reflected wave, and above the slab (##z>d##) you just have a transmitted wave. You still have boundary conditions at ##z=0## and ##z=d##, but for the field inside the slab you will need to solve the differential equation (which should yield two independent solutions). Of course, if ##E## is in the x direction, then ##H## is in the ##y## direction. If ##\epsilon_r(z)## is linear then the solution will be in terms of Airy functions (which are another way to write Bessel functions of order ##1/3##). See https://dlmf.nist.gov for info on special functions. If ##\epsilon_r## is a different profile then it will be other functions. So this gets messy in a hurry, just because of the special functions. Your final reflection coefficient will be in terms of these special functions.

If you are content with numerical solutions then you numerically solve the case for an N-layer slab, where each slab has a constant index of refraction. Basically you write the form of the fields for each slab, then apply boundary conditions at each boundary. You will then have 2 boundary conditions for each of ##N+1## boundaries, so will solve a ##2N+2 \, \times 2N+2## system of linear equations. Numerically this is easy; analytically it is not. You can start with N=4, or some other small number, compute the answer, then double N and re-compute, etc. The result should converge to the correct value pretty rapidly. One nice thing about this approach is that it is trivial to try different profiles for ##\epsilon_r(z)## once you have the code written. For the analytical approach above you have to solve a different differential equation for each profile.

Jason
My main problem with this approach was simply that I was expecting the derivation to break when we did not assume a uniform medium. This was true in my freshman mechanics derivation for transverse waves on a string. However, I found another derivation here which makes no such assumption:
http://www.math.ubc.ca/~feldman/m256/wave.pdf
I would think that there are some non-constant scattering potentials that would give a closed-form solution in elementary functions. Do you know of any that would work?
I agree wholeheartedly with @jasonRF. Sorry, I am also unable to attack your calculation

A few thoughts:
1. This is a pretty standard scattering problem which (as mentioned) has no happy analytic solution for even well-behaved "scattering potentials" .
2. There are a host of approximate methods (WKB, eikonal, perturbation expansions). I feel certain that you are reproducing certain aspects of them
3. It is not clear to me that the slab method will be well-behaved, particularly because of possible "resonances" for layers of wavelength commensurate widths
4. One test measure of the calculation is unitarity (conservation of energy).
5. A numerical method that might work (this may be your "hopscotch "...not clear to me) would be a Monte Carlo approach: tracing many pinball-like paths and adding them (with phases and weights) when they get to each edge. Might work well. Maybe a bad idea...
6. https://en.wikipedia.org/wiki/Eikonal_equation

The nice thing about this problem is that it is well-framed, useful in many realms of physics, and nontrivial. But I think you know that.
You are correct about both of those sentences!

My hopscotch method is essentially using the formula
$$E_{12}(d+x)=T_1E_{01}(d+n_{12}/n_{01}x)\pm R_1E_{21}(d-x)$$
recursively with the arguments shifted.

The fact of the matter is that I am taking a generic width of slab and allowing it to vary arbitrarily. That shouldn't be an issue.

Should the worst come to the worst, I can just numerically solve the ODE and see how well the solution agrees with my successive approximations. However, my maths instinct says this:

The PDE is not solvable in terms of elementary functions if ##n## is an arbitrary profile. Yet the LO and hopscotch methods both yield ##2\sqrt{n}/(n+1)## as a profile-independent result. Thus there must be at least one more contribution beyond the LO, and the hopscotch method fails to converge to the right answer.

The most vexing thing is that with a certain potential, I can get a closed-form elementary expression for the ##\mathrm{N}^n\mathrm{LO}## contribution for every ##n##. It would be interesting to check this potential using the analytical differential equation.

EDIT: I think I may have some luck with Euler ODEs. I'll give that a try.

Thanks again for all your help! This has been very informative and interesting.

Cheers,
QM

Homework Helper
2022 Award
• {???}
{???}
There is a potential called after Poschl-Teller ( a "reflectionless" potential) that is well studied and at least partially soluble.

https://en.wikipedia.org/wiki/Pöschl–Teller_potential

There is a lot of stuff about it if you explore a bit. Enjoy.
A reflectionless potential?!

See, it's posts like yours that make me love PhysicsForums. I've studied a lot of interesting things that branch off of basic physics ideas - the Capstan equation, Maxwell's fish-eye lens, and this problem of late - and it's little gems like what you've just shared that make physics fun!

But just a quick update for you (all): The potential for which I can calculate in closed form the order-by-order contributions is also miraculously a solvable ODE. A surprisingly easy one, in fact - I've solved much harder ODEs - and at least upon first glance the solution seems to satisfy basic sanity checks like "the wave kinda looks wave-like if ##n## is almost constant."

However, I am about to perform one last sanity check that I'm not so confident about. See, for uniform regions/slabs the effective reflection/transmission coefficients don't depend on the wavelength. It seems like this should also be true for my test potential, but I definitely need to check it by hand.

Thanks again,
QM