How to calculate the refracted light path when refraction index continuously increasing?

Question

Suppose an incident light from vacuum ($n_1=1.0$) into some media ($n_2=n_1+\mu\; x^2$) as in the figure below.

How to calculate the refracted light path curve in closed form?

Update:

Try to set up ordinary differential equation for the refracted light path per Snell's law.

Suppose the curve is $y=y(x)$;

Since $n_i \sin\theta_i=\text{constant}=n_1\sin\alpha=\sin\alpha$.

For any point $P:(x_0,y(x_0))$ on the path $y(x)$, we have: $$\tan(\theta_P)=\dfrac{\sin\theta_P}{\cos\theta_P}=y'(x)=\dfrac{\rm{d}y}{\rm{d}x},\quad \text{where }\theta_P \text{ is incident / refracted angle}$$

Since $\theta_P$ is always an acute angle, we have:

$$\dfrac{\sin^2\theta_P}{{1-\sin^2\theta_P}}=y'(x)^2\Rightarrow \sin\theta_P=\dfrac{\pm y'(x)}{\sqrt{1+y'(x)^2}}$$

Clearly $n_P\sin\theta_P=\sin\alpha$, where $n_P=1+\mu x^2$, then we have:

$$\left(1+\mu x^2\right)\dfrac{\pm y'(x)}{\sqrt{1+y'(x)^2}}=\sin\alpha\quad\text{with: y(0)=5} ||y'(0)=\tan\alpha$$

Then it becomes how to solve the ODE with a boundary condition. Can the ODE be solved in closed form?

Answer 1

This may (or may not) lead to the same answer as CuriousOne's suggestion above, but the most appropriate (and the longest) way of attempting a solution would to be to employ the Fermat's principle. The method's nicely described in the link, but in a nutshell, you would be led to a condition of the type $$\delta \int n ds = 0$$ where this $ds$ can be cast in terms of your 2D co-ordinates. Now, substitute for the spatial dependence of $n$ and arrive at $$\delta \int n(x,y) \sqrt{(1+(dy/dx)^2)} dx = 0$$

This is a sort of an ab-initio approach. I won't be surprised if there's a shorter method (maybe CuriousOne's suggestion.)

Answer 2

The derivation of the "equations of motion" for the light ray from the Fermat principle is given in the book "Reflections on Relativity", chapter 8.4 "Refractions on Relativity".

We know that the index of refraction $n$ at a point $(x,y)$ equals $c/v$, where $v$ is the velocity of light at that point. Thus, if we parameterize the path by the equations $x = x(u)$ and $y = y(u)$, the "optical path length" from point $A$ to point $B$ (i.e., the time taken by a light beam to traverse the path) is given by the integral

$$L=\int\limits_A^B n\,\mathrm{d}s=\int\limits_A^Bn\sqrt{\dot{x}^2+\dot{y}^2}\,\mathrm{d}u$$

where dots signify derivatives with respect to the parameter $u$. To make this integral an extremum, let $f$ denote the integrand function

$$f(x,y,\dot{x},\dot{y})=n(x,y)\sqrt{\dot{x}^2+\dot{y}^2}$$

Then the Euler equations (introduced in Section 5.4) are

\begin{align} \frac{\partial n}{\partial x}=\frac{\mathrm d}{\mathrm{d}u}\left(\frac{\partial f}{\partial \dot x}\right) && {} && \frac{\partial n}{\partial y}=\frac{\mathrm d}{\mathrm{d}u}\left(\frac{\partial f}{\partial \dot y}\right) \end{align}

which gives

\begin{align} \frac{\partial n}{\partial x}\sqrt{\dot{x}^2+\dot{y}^2}=\frac{\mathrm d}{\mathrm{d}u}\left[\frac{n\dot x}{\sqrt{\dot{x}^2+\dot{y}^2}}\right] && {} && \frac{\partial n}{\partial y}\sqrt{\dot{x}^2+\dot{y}^2}=\frac{\mathrm d}{\mathrm{d}u}\left[\frac{n\dot y}{\sqrt{\dot{x}^2+\dot{y}^2}}\right] \end{align}

Now, if we define our parameter $u$ as the spatial path length $s$, then we have $\dot{x}^2+\dot{y}^2=1,$ and so the above equations reduce to

$$\frac{\partial n}{\partial x}=\frac{\mathrm d}{\mathrm{d}s}\left(n\frac{\mathrm{d} x}{\mathrm{d} s}\right)\tag{1a}$$ $$\frac{\partial n}{\partial y}=\frac{\mathrm d}{\mathrm{d}s}\left(n\frac{\mathrm{d} y}{\mathrm{d} s}\right)\tag{1b} $$

These are the "equations of motion" for a photon in a heterogeneous medium, as they are usually formulated, in terms of the spatial path parameter $s$.

Now, solving these equations for $x(s)$ and $y(s)$, you'll get your ray curve in your medium $n(x,y)$.