Errata and bug: about signs for imaginary parts of complex refractive index and refracted angles

Hi, I found that some detailed description in the document https://mitsuba.readthedocs.io/en/latest/src/key_topics/polarization.html#theory seems to contain errata, and the current implementation has a bug in terms of sign flipping in imaginary parts. I'm not sure where this issue should be posted, but I want to discuss it for completeness of Mitsuba 3 documentation and code.

Errata in documentation

What I want to discuss is about the following paragraphs in https://mitsuba.readthedocs.io/en/latest/src/key_topics/polarization.html.

Here, $\eta = n-ki$ comes from $\exp\left(+i\omega t\right)$ rather than $\exp\left(-i \alpha \eta t\right)$ (while difference between my notation $\omega$ and $\alpha\eta$ in Mitsuba documentation is not importance, their signs are very important). Similarly, $\eta = n+ki$ comes from $\exp\left(-i\omega t\right)$ contrast to what rewritten in the last paragraph. The current document has flipped signs, and description of the reasoning of these signs is wrong. In a wave formulation $\exp\left(\pm i \left(\omega t - \boldsymbol\kappa \cdot \mathbf r\right)\right)$, the coefficient $\omega$ of time $t$ should be a positive real number, and the complex coefficient is $\boldsymbol\kappa$.

Issue in code: $\cos\theta_t$ on conducting media

fresnel_polarized in fresnel.h computes complex refracted angles for total reflection between dielectric media and incidence on a conducting medium. While the former is implemented correctly (i.e., consistent to the convention of $\eta = n-ki$), I found that the later one, complex refracted angles on conducting media seem not to be correctly implemented.

import mitsuba as mi
mi.set_variant('scalar_rgb')
print(mi.fresnel_polarized(0.5, 0.5))
print(mi.fresnel_polarized(0.5, 0.5+1e-8j)) # Note that input `0.5+1e-8j` will be converted to `0.5-1e-8j`
                                            # at the inside of the function following Mitsuba 3 convention

([-0.3333333432674408 + 0.9428090453147888i], [-0.9393939971923828 + 0.34283965826034546i], 0.0, [0.5 + 0.0i], [2.0 - 0.0i])
([-0.3333333432674408 - 0.9428090453147888i], [-0.9393939971923828 - 0.34283965826034546i], 0.0, [0.5 - 9.99999993922529e-09i], [2.0 + 3.999999975690116e-08i])

As one effect of inconsistent implementation, here we can see a discontinuous behavior for the first and second entries (a_s and a_p) for the output tuples for very little change of imaginary parts of the refractive index.

The actual derivation for complex angles is written in the next section.

Detailed description

Since this sign convention problem is very very, so much confusing, I want to start from the most fundamental facts and derive results step-by-step to argue the errata.

Plane waves

From Maxwell's equations in a medium with permittivity $\epsilon$, permeability $\mu$, and conductivity $g$, assuming no free charge and $\mathbf j_{\mathrm{free}} = g \mathbf {E}$ yields the following equations:

$$ \nabla^2 \mathbf E\left(\mathbf r,t\right) - \epsilon \mu \frac{\partial^2 \mathbf E}{\partial t^2}\left(\mathbf r,t\right) - g \mu \frac{\partial \mathbf E}{\partial t}\left(\mathbf r,t\right) = 0, \ \nabla^2 \mathbf H\left(\mathbf r,t\right) - \epsilon \mu \frac{\partial^2 \mathbf H}{\partial t^2}\left(\mathbf r,t\right) - g \mu \frac{\partial \mathbf H}{\partial t}\left(\mathbf r,t\right) = 0. $$

It holds independent of choice of conventions yet.

References

[Reitz] Chap. 17-4, Eq. (16-29) in p. 423

[Hecht] Chap. 4.8, Eq. (4.74) in p. 139

Assume a plane wave solution $\mathbf E\left(\mathbf r,t \right) = \mathbf E e^{\pm i\left(\omega t- \boldsymbol{\kappa} \cdot \mathbf r\right)}$. Note that taking the sign $\pm$ yields each of different conventions. To make clear which part depends the conventions and which are not, we keep to write $\pm$ symbols rather than choosing one convention. (Then $\mp$ indicates the opposite sign to the convention symbol $\pm$)

Substituting the plane wave form into the first equation, we obtain that it becomes a solution whenever the following relation between $\omega$ and $\boldsymbol\kappa$ holds.

$$ \epsilon \mu \omega^2 \mp ig\mu \omega - \boldsymbol\kappa\cdot \boldsymbol\kappa = 0. $$

Note

Here, for a complex vector $\mathbf v=\left(v_x, v_y, v_z\right)$, the dot product is defined as $\mathbf v\cdot \mathbf v = v_x^2 + v_y^2+v_z^2$, which may be negative or even complex.

Defining a complex permittivity as $\tilde \epsilon \coloneqq \epsilon\mp \frac{g}{\omega}i $, the above equation can be rewritten as,

$$ \tilde \epsilon \mu \omega^2 -\boldsymbol\kappa\cdot \boldsymbol\kappa = 0. \quad\quad\quad(1) $$

Derivation

In the original differential equation, $\frac{\partial}{\partial t}$ can be substituted by $\pm i\omega$, and $\nabla$ does by $\mp i \boldsymbol\kappa$.

$- \left(\boldsymbol\kappa\cdot\boldsymbol\kappa\right) \mathbf E + \epsilon \mu\omega^2 \mathbf E \mp i g \mu \omega \mathbf E = 0$.

Reference

[Reitz] Chap. 17-4, Eq. (17-45) in p. 424

Two semi-infinite media

Now considering a scenario that a plane wave at Medium 1 with $\epsilon_1\in\mathbb R$ Incidents to Medium 2 with $\tilde \epsilon_2\in \mathbb C$, and a common permeability $\mu$.

$$ \begin{align} \text{incident: }& \mathbf E_0\left(\mathbf r,t\right) &=\hat{\mathbf E}_0 e^{\pm i\left(\omega t-\boldsymbol\kappa_0\cdot \mathbf r\right)} &\ \text{reflected: }& \mathbf E_1\left(\mathbf r,t\right) &=\hat{\mathbf E}_1 e^{\pm i\left(\omega t-\boldsymbol\kappa_1\cdot \mathbf r\right)} &\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad(2) \ \text{transmitted: }& \mathbf E_2\left(\mathbf r,t\right) &=\hat{\mathbf E}_2 e^{\pm i\left(\omega t-\tilde {\boldsymbol\kappa}_2\cdot \mathbf r\right)}&= {\hat{\mathbf E}}_2 e^{\pm \Im {\boldsymbol\kappa_2} \cdot \mathbf{r}} e^{\pm i\left(\omega t- \Re\boldsymbol\kappa\cdot \mathbf r\right)} \end{align} $$

Here, I am writing the three waves with the same temporal frequency $\omega$, which is explained the next next paragraph. While $\boldsymbol\kappa_0$ for the incident wave is a real vector, we do not know ${\boldsymbol\kappa}_1$ and $\tilde{\boldsymbol\kappa}_2$ are real or complex yet, but I will prove that ${\boldsymbol\kappa}_1$ is always real but $\tilde{\boldsymbol\kappa}_2$ may be complex.

Boundary conditions of Maxwell's equations for these three wave functions consist of four equations of linear combinations of them such as $\epsilon1\mathbf E{0}\left(\mathbf r,t\right)\cdot \mathbf n + \epsilon1\mathbf E{1}\left(\mathbf r,t\right)\cdot \mathbf n = \epsilon2\mathbf E{2}\left(\mathbf r,t\right)\cdot \mathbf n $. This can be splitted into two conditions (where both should be satisfied): (1) constraints for $e^{\pm i\left(\omega t- \boldsymbol\kappa\cdot \mathbf r\right)}$, and (2) constraints for amplitudes (including phases) $\hat{\mathbf E}_0$, $\hat{\mathbf E}_1$, and $\hat{\mathbf E}_2$, which yield well known Fresnel equations.

(1) Constraints for $e^{\pm i\left(\omega t-\boldsymbol\kappa\cdot \mathbf r\right)}$

Claim: incident, reflected, and transmitted waves have the same real temporal frequency $\omega$ and may have complex values in waves numbers $\boldsymbol\kappa$.

Note that the boundary conditions should hold for any time and any position on the interface, $z=0$. It first yields the three waves have the identical coefficient of the time dependency at their exponents, so that the temporal frequency $\omega$ for the refracted wave should be real, contrast to what is written in the current Mitsuba 3 documentation. Note that the equations above have already been written using the same $\omega$ based on this fact.

Reference

[Reitz] Chap. 17-4 p. 425 $\mathbf E\left(\mathbf r, t\right) = \left(\hat {\mathbf E}e^{-\mathbf \kappa_i\cdot \mathbf r} \right)e^{-i\left(\omega t - \mathbf \kappa_r \cdot \mathbf r\right)}$

$\mathbf B\left(\mathbf r, t\right) = \left(\hat {\mathbf B}e^{-\mathbf \kappa_i\cdot \mathbf r} \right)e^{-i\left(\omega t - \mathbf \kappa_r \cdot \mathbf r\right)}$

Next, $e^{\mp i \boldsymbol\kappa\cdot \mathbf r}$ for the three waves should be identical at any position at $z=0$. It can be rewritten as follows:

$$ {\boldsymbol \kappa}_0\times \mathbf n = {\boldsymbol \kappa}_1\times \mathbf n = \tilde {\boldsymbol \kappa}_2\times \mathbf n \quad\quad\quad(3) $$

From Equation (3), we observe that all wave vectors $\boldsymbol\kappa_0$, $\boldsymbol\kappa_1$, and $\tilde{\boldsymbol\kappa}_2$ lie on a single plane. Assigning a coordinate system as follows, $\boldsymbol\kappa_0$, $\boldsymbol\kappa_1$, and $\tilde{\boldsymbol\kappa}_2$ only have $x$ and $y$ coordinates. Then they are rewritten using angles $\theta_i$, $\theta_r$, and $\tilde \theta_t$ as follows:

fresnel coordinates

$$ \begin{align} {\boldsymbol\kappa}_0 &= \left(\kappa_0 \sin\theta_i, -\kappa_0\cos\theta_i, 0\right) \ {\boldsymbol\kappa}_1 &= \left(\kappa_1 \sin \theta_r, \kappa_1\cos\theta_r, 0\right) \quad\quad\quad\quad\quad(4)\ \tilde{\boldsymbol\kappa}_2 &= \left(\tilde \kappa_2 \sin\tilde\theta_t, -\tilde\kappa_2\cos\tilde\theta_t, 0\right)\ \end{align}, $$

where $\kappa_0 = \sqrt{\boldsymbol\kappa_0 \cdot\boldsymbol\kappa_0} = \sqrt{\epsilon_1\mu}\omega\ge 0$.

Q. Does any complex 2-vector can be written in this way ($r$ and $\theta$)?

For a complex vector $\tilde{\mathbf v} = \left(\tilde x, \tilde y\right)\in {\mathbb C}^2$, find $\tilde r$ and $\tilde \theta$ such that $\tilde{\mathbf v} = \left(\tilde r\cos \tilde \theta, \tilde r \sin \tilde \theta\right)$ $\tilde x^2 + \tilde y^2 = \tilde r^2 \quad \Rightarrow \quad \tilde r = \left(\tilde x^2 + \tilde y^2\right)^{\frac12}$ Letting $\tilde \theta = \frac1i\ln\frac{\tilde x+ i \tilde y}{\tilde r} + 2k\pi,$ for any integer $k$, $\tilde x=\tilde r\cos\tilde \theta$ and $\tilde y = \tilde r \sin\tilde \theta$ hold.

Carefully considering that powers to the $1/2$ and $\ln$ have multiple values for complex variables,

excluding exceptional cases of $\tilde x^2+\tilde y^2 = 0 $,

there only exists two pairs of solution for $\tilde r\in \mathbb C$ and $\tilde \theta \in \left\lbrace z\in \mathbb C \mid 0\le \Re z< 2\pi \right\rbrace$ where the two pairs can be converted to each other by substituting $\tilde r \leftarrow -\tilde r$ and $\tilde \theta\leftarrow\tilde \theta + \pi$.

From Equations (1) and (3), we easily get $\kappa_1 =\kappa_0$ and $\theta_r = \theta_i$. For $\tilde{\boldsymbol\kappa}_2$, Equation (1) is rewritten as follows.

$$ \tilde \kappa_2^2 = \tilde \epsilon_2 \mu \omega^2 = \frac{\tilde \epsilon_2}{\epsilon_1} \kappa_0^2 $$

Define the relative refractive index $\eta$ as one of square root of $\tilde \epsilon_2 / \epsilon_1$, $\tilde \kappa_2 = \eta \kappa_0$

Claim: $\exp\left(\pm i\omega t\right)$ convention for plane waves yields complex refractive indices of signs $\eta = n \mp k i$, respectively.

I believe the real part of $\tilde \kappa_2$ should be taken as positive since it represents the actual wavelength of the wave. For example, if we consider two media with the same permittivity and zero conductivity for one and infinitesimally small conductivity for the other, as the example code above, it will be natural that two media yield very similar values of wave numbers $\kappa$. Even though I could not found a concrete reference for this reasoning, I provide several examples of references consistent of this claim.

Reference

[Reitz] : $-i\omega t$ convention and $n+ik$

Sec. 17-4 p. 425

$\mathbf E\left(\mathbf r, t\right) = \left(\hat {\mathbf E}e^{-\mathbf \kappa_i\cdot \mathbf r} \right)e^{-i\left(\omega t - \mathbf \kappa_r \cdot \mathbf r\right)}$

$\mathbf B\left(\mathbf r, t\right) = \left(\hat {\mathbf B}e^{-\mathbf \kappa_i\cdot \mathbf r} \right)e^{-i\left(\omega t - \mathbf \kappa_r \cdot \mathbf r\right)}$

In p. 426 Eq. (17-54) $\hat n = n + ik$

[Hecht] : $+i\omega t$ convention and $n-ik$

Sec. 4.8 p. 139 Eq. (4.76) $\tilde n = n_R - in_I$

$\vec{\mathbf E} = \vec{\mathbf E}_0 e^{\left(-\omega n_I y / c\right)}e^{i\omega \left(t-n_R y / c\right)}$

[Mueller] : explicitly stating $\pm i\omega t$ convention yields $n\mp ik$

p. 17 TABLE 4

Alternatives: $e^{-i\omega t}$ $e^{i\omega t}$

Effect on complex refractive index: $n\left(1+i \kappa\right)$ $n\left(1-i \kappa\right)$

Alternatives:	$e^{-i\omega t}$	$e^{i\omega t}$
Effect on complex refractive index:	$n\left(1+i \kappa\right)$	$n\left(1-i \kappa\right)$

Note that taking the sign of $\eta = \left(\tilde \epsilon \mu\right)^{\frac12}$ in this way, the argument (the angle between given complex number and the positive real axis in the complex plane) of $\eta$ belongs to the following range.

$$ \arg\eta \in \begin{cases} \left[-\frac\pi2, 0\right] & e^{+i\omega t} \text{ convention} \ \left[0, \frac\pi2\right] & e^{-i\omega t} \text{ convention} \end{cases} \quad\quad\quad\quad\quad(5) $$

Note that $\tilde \kappa_2 = \eta \kappa_0$ has the argument in the same range as $\eta$.

Now we can get $\tilde \theta_t$. From Equation (3) and (4), we get $\kappa_0\sin\theta_i = \tilde\kappa_2 \sin\tilde \theta_t$. Substituting $\tilde \kappa_2=\eta\kappa_0$ to here, $\sin\tilde\theta_t = \frac1{\eta}\sin\theta_i$. Then $\cos\tilde\theta_t$ becomes one of square roots of $\cos^2\tilde\theta_t=1-\sin^2\tilde\theta_t = 1-\frac1{\eta^2}\sin^2\theta_i$. From Equation (5), we get:

$$ \arg\left(\cos^2\tilde\theta_t\right) \in \begin{cases} \left[-\pi, 0\right] & e^{+i\omega t}$ \text{ convention} \ \left[0, \pi\right] & e^{-i\omega t}\text{ convention} \end{cases} $$

To determine the sign of $\cos\tilde\theta_t$, we should select the sign which makes the refracted wave decays rather than diverge. Substituting Equation (4) into Equation (2) for the refracted wave,

$$ \mathbf E_2\left(\mathbf r,t\right) =\hat{\mathbf E}_2 e^{\pm i\left(\omega t-\tilde {\boldsymbol\kappa}_2\cdot \mathbf r\right)}= \hat{\mathbf E}_2 e^{\mp \Im\left(\tilde\kappa_2\cos\tilde\theta_t\right)y} e^{\pm i\left(\omega t- \tilde\kappa_2\sin\tilde\theta_t x + \Re\left(\tilde\kappa_2\cos\tilde\theta_t\right)y\right)} $$

In our coordinate system, the wave decays to zero at $y=-\infty$. It can be written as $\mp \Im\left(\tilde\kappa_2\cos\tilde\theta_t\right)\cdot \left(-\infty\right)\le 0$, furthermore,

$$ \mathrm{sign}\left(\tilde\kappa_2\cos\tilde\theta_t\right) = \mp 1, $$

for $e^{\pm i\omega t}$ convention, respectively. According ranges of arguments of $\tilde \kappa_2$ and $\cos^2\tilde\theta_t$, which are already obtained, we should take $\cos\tilde\theta_t$ with the argument in the following range.

$$ \arg\left(\cos\tilde\theta_t\right) \in \begin{cases} \left[-\frac\pi2, 0\right] & e^{+i\omega t} \text{ convention}\ \left[0, \frac\pi2\right] & e^{-i\omega t} \text{ convention} \end{cases} $$

Q. Does the sign of $\cos\tilde\theta_t$ depends on the choice of coordinate system?

When we flip the $y$ coordinate, Equation (4) should be rewritten as ${\tilde{\boldsymbol\kappa}}_2=\left(\tilde \kappa_2 \sin\tilde\theta_t, \tilde\kappa_2\cos\tilde\theta_t, 0\right)$. In addition, the refracted wave should converges to zero at $y=\infty$ rather than $y=-\infty$. It yields the same condition for $\cos\tilde\theta_t$ as before.

In summary,

Claim: For both total internal reflection between dielectric media and reflection on a conductor, $\Im \cos\tilde\theta_t$ should have $\mp$ sign for $e^{\pm i\omega t}$ convention, respectively.

Reference

A relevant reference is already cited in the comment in the Mitsuba 3 code.

mitsuba3:include/mitsuba/render/fresnel.h:std::tuple<dr::Complex<Float>, dr::Complex<Float>, Float, Float, Float> fresnel_polarized(Float cos_theta_i, Float eta)

"Stellar Polarimetry" by David Clarke

Note that this reference uses $e^{+i \omega t}$ convention (equivalent to $\cos\left(+\omega t\right)$) same as Mitsuba 3. (Equations in p.371-372)

This reference uses $\cos\theta_t$ as $-i \sqrt{\sin^2 \chi_i - n^2}$ in Equations (A38-39) in p.377 for total internal refelection cases.

However, I am sorry that I could not find any reference for signs of $\Im \cos\theta_t$ on conducting media.

I found that in Mitsuba 3 ($e^{+i\omega t}$ convention) while $\cos\tilde\theta_t$ for total internal reflection is implemented correctly,

https://github.com/mitsuba-renderer/mitsuba3/blob/3013adb4f12a491f8dd37c32bcedf55c7998f9e8/include/mitsuba/render/fresnel.h#L150

however, $\cos\tilde\theta_t$ for conductors is implemented wrong, as best of my knowledge.

https://github.com/mitsuba-renderer/mitsuba3/blob/3013adb4f12a491f8dd37c32bcedf55c7998f9e8/include/mitsuba/render/fresnel.h#L227

The discontinuous complex reflectance change in the code example above is yielded by this implementation issue for conductors.

(2) Constraints for $\hat {\mathbf E}_0$, $\hat {\mathbf E}_1$, and $\hat {\mathbf E}_2$: Fresnel equations

Constraints (boundary conditions) for amplitudes yields the well known Fresnel equation. I have not found any error in the Fresnel equation in Mitsuba 3 itself, I omit detailed description for this part.

References

[Reitz] Reitz, Milford, and Christy, Foundations of Electromagnetic Theory, Forth Edition, Benjamin Cummings
[Hecht] Eugene Hecht, Optics, Fifth Edition, Pearson
[Muller] Rolf H. Muller. Definitions and conventions in ellipsometry. Surface Science, 16:14–33, August 1969. doi:10.1016/0039-6028(69)90003-X.

@tizian this might interest you.

Thanks for letting me know. I'm away for the next 2.5 weeks or so but I'll try to take a closer look once I'm back.

The sign flip given in the example with a small complex component is indeed very suspicious. I think a first step would be to reproduce the plots from this section in the documentation. I do recall running similar experiments where I basically reproduced the dielectric behaviour with the complex case (and kappa ~= 0). But I don't see them in the document so probably didn't include them at the time.

It also has been a while since I've been thinking intensely about the various signs of the polarization code. It might take me a bit to refamiliarize myself with it 🙂

Thank you for the response. Following the suggestion, I also reproduced Figures 12 and 14 in the documentation.

import numpy as np
import matplotlib.pyplot as plt

def reflectance(cos_i, cos_t, eta):
    rs = (cos_i - eta*cos_t) / (cos_i + eta*cos_t)
    rp = (eta*cos_i - cos_t) / (eta*cos_i + cos_t)
    return rs, rp

theta_i = np.arange(91, dtype=np.float64)
cos_i = np.cos(np.deg2rad(theta_i))
sinsq_i = 1 - cos_i**2

def show_delta(eta):
    cos_t = np.sqrt(1+0j - sinsq_i / (eta**2))
    cos_t_ipos = np.where(cos_t.imag < 0, -cos_t, cos_t)
    cos_t_ineg = np.where(cos_t.imag > 0, -cos_t, cos_t)

    fig, axes = plt.subplots(3, 2, figsize=(8, 10))
    for i_sign, (ineq, cos_t) in enumerate(zip("><", [cos_t_ipos, cos_t_ineg])):
        rs, rp = reflectance(cos_i, cos_t, eta)

        axes[0, i_sign].plot(theta_i, rs.real, label=r"$\Re r_s$")
        axes[0, i_sign].plot(theta_i, rs.imag, label=r"$\Im r_s$")
        axes[0, i_sign].plot(theta_i, rp.real, label=r"$\Re r_p$")
        axes[0, i_sign].plot(theta_i, rp.imag, label=r"$\Im r_p$")
        axes[0, i_sign].set_xlabel(r"$\theta_i$")
        axes[0, i_sign].set_ylabel("reflectance (complex amplitude)")
        axes[0, i_sign].set_title(f"$\\Im \\cos \\theta_t {ineq} 0$")
        axes[0, i_sign].legend()

        Rs = np.abs(rs*rs)
        Rp = np.abs(rp*rp)

        axes[1, i_sign].plot(theta_i, Rs, label=r"$R_s$")
        axes[1, i_sign].plot(theta_i, Rp, label=r"$R_p$")
        axes[1, i_sign].set_xlabel(r"$\theta_i$")
        axes[1, i_sign].set_ylabel("reflectance (power)")
        axes[1, i_sign].set_title(f"$\\Im \\cos \\theta_t {ineq} 0$")
        axes[1, i_sign].legend()

        delta = np.angle(rp/rs, deg=True)
        if delta[0] < -150 and delta[1] > 150:
            # For better visualization
            delta[0] += 360
        axes[2, i_sign].plot(theta_i, delta)
        axes[2, i_sign].set_xlabel(r"$\theta_i$")
        axes[2, i_sign].set_ylabel(r"$\delta_p - \delta_s$")
        axes[2, i_sign].set_title(f"$\\Im \\cos \\theta_t {ineq} 0$")

    plt.tight_layout()

eta = 1/1.5
show_delta(eta)

eta = 0.183 - 3.43j
show_delta(eta)

eta = 0.183 + 3.43j
show_delta(eta)

In my understanding, plots in Figures 12 and 14 in the Mitsuba documentation are based on considering phase as $\delta$ in $e^{+i\left(\omega t+\delta\right)}$. Consistently to my claims, my reproduced plots match to the plots in the Mitsuba documentation when taking imaginary parts of $\cos\theta_t$ to be negative. Note that taking the imaginary parts to be positive is relavent to considering phase as $\delta$ in $e^{i\left(kx-\omega t + \delta\right)}$, which yields plots with signed flipped values.

I finally found some time to play with this a bit, and indeed something seems to be wrong about the conductor case. Given that you raise two issues at once, let's start with the first one regarding the documentation text.

I'm not sure I completely follow this. Do you have a concrete suggestion of how to rephrase the description? And maybe a good source to cite in addition? My understanding here was the following:

A wave can be described by

E(t) = E_0 \, \exp(-i \, t \, \alpha \, n).

For real IOR $n$, this is simply going to be a non-absorptive wave, with real component

\Re[E(t)] = E_0 \, \cos(t \, \alpha \, n).

An imaginary component of $\eta$ can be introduced to model absorption. And given the " $exp(-i)$ " convention, this works out by choosing $n - k i$:

\exp(-i \, t \, \alpha \, (n - k \, i)) =
\exp(-i \, t \, \alpha \, n) \, \exp(i \, t \, \alpha \, k \, i) =
\exp(-i \, t \, \alpha \, n) \, \exp(-t \, \alpha \, k)

\Re[E(t)] = E_0 \, \cos(t \, \alpha \, n) \exp(-t \, \alpha \, k),

which amounts to a decaying wave when $k > 0$.

The problem what I want to point out comes from that the current Mitsuba 3 document explains formulae for waves only with time-dependency term, without spatial dependency. In general, a wave can be described by (now ignoring vector amplitude $\mathbf E_0$, for simplicity)

$$ E\left(\mathbf x, t\right) = E_0 \exp\left(\mp \boldsymbol\kappa \cdot \mathbf x \pm i \alpha t \right), \tag{2-1} $$

with the relation

$$ \eta^2 \alpha^2 - \boldsymbol\kappa \cdot \boldsymbol\kappa = 0. \tag{2-2} $$

Here $\alpha$ is a temporal frequency (written as $\omega$ in the my original post and other physics textbooks. However, it will be confused with direction vectors, so I guess it is the reason that Mitsuba 3 uses the symbol $\alpha$ rather than $\omega$), $\eta$ is a (real or complex) refractive index, and $\boldsymbol\kappa$ is a wave number ($2\pi$ over the wavelength). Eq. (2-2) is identical to Eq. (1) in my original post (substituting $\tilde \epsilon\mu \leftarrow \eta ^2$, $\omega \leftarrow \alpha$).

Note that Eq. (2-2) does not contain $\eta$ term. Assume that the notations $\alpha$ and $\boldsymbol\kappa$ in Eq. (2-2) denote the frequency and the wave number of an incident wave, which lives in the medium of $\eta = 1$. For the transmitted wave at a medium with $\eta \ne 1$, Eq.~\eqref{2-2} will be change to be contain $\eta$ dependency.

The first part I want to point out is that the $\eta$ dependency should appear in the wave number $\boldsymbol\kappa$ term, not the frequency $\alpha$ term. Concretely, not $\exp\left(\mp \boldsymbol\kappa \cdot \mathbf x \pm i \alpha t /\eta\right)$, we should use $\exp\left(\mp \boldsymbol\kappa_2 \cdot \mathbf x \pm i \alpha t \right)$ for some solution $\boldsymbol\kappa_2$ for $\boldsymbol\kappa_2 \cdot \boldsymbol\kappa_2=\eta^2 \boldsymbol\kappa \cdot \boldsymbol\kappa$ (since wave vectors are not scalar, it is more complicated than simply substituting $\boldsymbol\kappa_2 = \eta \boldsymbol\kappa$). Here, writing $i\alpha \eta t$ rather $i \alpha t / \eta$ in Mitsuba document seems to be another minor mistake. Whether the $\eta$ dependency appears in the spatial $\boldsymbol\kappa \cdot \mathbf x$ or the temporal $\alpha t$ term is a matter for physical meaning and boundary conditions of Maxwell's equations, while the discussion of whether "decay" or "diverge" is a matter for a later stage. While solving the exact boundary conditions will be more complicated, it can be summarized as "values of $E$ for incident, reflected, refracted waves should match at the interface of the two media with some constant multiplicative factor." Note that the fact that $\eta$ dependency appears in the spatial $\boldsymbol\kappa \cdot \mathbf x$ term indicates wavelengths will be changed, while frequencies will not be changed. It is a famous fact which often visualized as illustrations such as https://galileo-unbound.blog/2020/01/27/snells-law-the-five-fold-way/, where matching wavefronts at the interface implicitly implies the boundary conditions.

I also want to write a concrete suggestion for rephrasing the Mitsuba documentation, but I apologize I would be able to write that suggestion and more discussions at least next week due to my personal schedule.

I think I see where you're getting at. I found a slightly more gentle explanation in "Optics" by Hecht (5th Edition) Section 4.8 "Optical Properties of Metals" ... (Notation adjusted to avoid naming conflicts.)

Harmonic wave propagating along spatial dimension $x$ and time $t$: $$E(x,t) = E_0\exp\left(i (\omega t - \alpha x)\right)$$

with $\omega$ and $\alpha$ being the angular temporal and spatial frequencies respectively. I.e.

$$\omega = \frac{2\pi}{\lambda}, \alpha = \frac{2\pi}{\tau},$$

when parameterizing with temporal and spatial periods $\lambda$ and $\tau$ instead.

Spatial and temporal propagation are actually related by considering the propagation speed $v$ as $\alpha = \omega/v$. As the IOR is defined as a ratio of speed and the speed of light ($\eta = c / v$) we have $\alpha = (\omega \eta) / c$ and

$$ E(x,t) = E_0 \exp\left(i (\omega t - \frac{\omega \eta}{c} x)\right). $$

Now, choosing the convention to represent conductor IORs with a negative complex component $\eta = n - k i$ will actually result in exponential attenuation as the wave propagates into the medium.

$$E(x,t) = E_0 \exp(i \omega t) \exp(-i \frac{\omega \eta}{c} x)$$

$$= E_0 \exp(i \omega t) \exp(-i \frac{\omega (n - k i)}{c} x)$$

$$= E_0 \exp(i \omega t) \exp(-i \frac{\omega n}{c} x) \exp(i^2 \frac{\omega k}{c} x)$$

$$= E_0 \exp(i \omega t) \exp(-i \frac{\omega n}{c} x) \exp(-\frac{\omega k}{c} x)$$

$$= E_0 \exp(i (\omega t - \frac{\omega n}{c} x)) \exp(-\frac{\omega k}{c} x).$$

This is actually where we started, just with $\Re[\eta] = n$ describing the harmonic wave and $\Im[\eta] = k$ describing the additional attenuation term.

All that being said, I'm not convinced the discussion of both spatial and temporal dimension here adds a lot to the high-level insight behind the complex IOR. We can always be more precise and more correct but need to draw the line somewhere.

In particular, the temporal aspect actually doesn't play any role in these equations. So a compromise could be to change the documentation text referencing "increase in time" to "increase in distance into the medium" and replace the reasoning with the equations above---but just illustrating the spatial dimension. (I.e. remove " $\exp(i \omega t)$ " throughout.)

A reference to Hecht should then be enough for the more interested reader to dig deeper into the topic.

Sorry for late response

About your latest comment

Your understanding through [Hecht] looks right. (while denoting temporal and spatial periods by $\tau$ and $\lambda$ will be more consistent to other references, your equations in your notations are correct, except a minor errata $\Im\left[\eta\right]=k$ -> $\left|\Im\left[\eta\right]\right| = k$)

All that being said, I'm not convinced the discussion of both spatial and temporal dimension here adds a lot to the high-level insight behind the complex IOR. We can always be more precise and more correct but need to draw the line somewhere. ... A reference to Hecht should then be enough for the more interested reader to dig deeper into the topic.

I totally agree with you. One of the reasons I wrote such a long derivation is just clarification for verification. Since it is so confusing, I was concerned that I might be wrong, and wanted to clarify the details to avoid another error. I agree that writing brief explanation in the 1D wave case which you wrote above or simply stating "$n-ki$ leads to an exponential decay of the wave as it travels inside the conductor material. For the derivation see p. xx [Hecht]" in the documentation is enough.

While obliquely incident rays with vector spatial frequency rather than the scalar one will be more complicated, (unfortunately the direction of exponential decay and the direction of transmitted plane wave will be different in general), it will not be need to be explained in the documentation.

In particular, the temporal aspect actually doesn't play any role in these equations. So a compromise could be to change the documentation text referencing "increase in time" to "increase in distance into the medium" and replace the reasoning with the equations above---but just illustrating the spatial dimension. (I.e. remove " " throughout.)

I understand what you want say, and actually the temporal term itself seems not be matter to determine the sign of complex numbers. It will be confusing.

Note that for a general wave function formed $f\left(x - vt\right)$, its propagation direction is same as $-(\text{coefficient of }t)/(\text{coefficient of }x)$. You did test $E\left(x=\infty, t\right)=0$ to show the wave decays as it travels since its propagation direction has the $+$ sign: $-\left(i\omega\right)/\left(-i\frac{\omega n}c\right)>0$. If we had a wave with the formula $E\left(x,t\right) = E_0 \exp\left(-i\omega t - i\frac{\omega\eta }c x\right)$ which has the same sign for the coefficient of $x$ but the opposite sign for one of $t$, we should test $E\left(x=-\infty, t\right) = 0$ to show the wave decays since its propagation direction is $-\left(-i\omega\right)/\left(-i\frac{\omega n}c\right)<0$, $-x$ direction. Then it need the opposite imaginary sign of the complex IOR $\eta=n+ki$.

In summary, using $\omega >0$. Thus, for a harmonic wave with propagation direction $\propto \boldsymbol\alpha$, we have two choices of conventions: $\exp\left(i\left(\omega t - \boldsymbol\alpha \cdot\mathbf x\right)\right)$ and $\exp\left(i\left(\boldsymbol\alpha \cdot\mathbf x - \omega t \right)\right)$. (when we decide to make the symbol $\boldsymbol\alpha$ denote the direction of propagation, we cannot use $\exp\left(i\left(\omega t + \boldsymbol\alpha \cdot\mathbf x\right)\right)$ or $\exp\left(-i\left(\omega t+ \boldsymbol\alpha \cdot\mathbf x\right)\right)$ since these have the propagation direction $-\boldsymbol\alpha$) A way to concretely avoid confusions is calling the two conventions with exact dependencies of both temporal and spatial terms, but many previous literature including [Mueller] calls the two conventions as "$\exp\left(+i\omega t\right)$ and $\exp\left(-i\omega t\right)$ conventions", which implicitly determine spatial terms $\exp\left(-i\boldsymbol\alpha\cdot \mathbf x\right)$ and $\exp\left(+i\boldsymbol\alpha\cdot \mathbf x\right)$, resp. I recommend to use the same name (" $\exp \left( +i \omega t \right)$ and $\exp\left(-i\omega t\right)$ conventions") to refer to those conventions as previous literature.

Concrete suggestion for the documentation text

In conclusion, my suggestion is changing the following current phrases:

The same expressions also hold in the conductor case where the indices of refraction are complex valued, i.e. $\eta = n-ki$ with real and imaginary parts $n$ and $k$. The use of complex numbers here comes from the expression $\exp\left(-i\alpha \eta t\right)$ that describes a plane wave propagating with increasing time and a constant $\alpha$ related to the wave frequency. When introducing a value of $k>0$, this now leads to an exponential decay of the wave as it travels inside the conductor material: $\exp\left(-i n \alpha t \right) \cdot \exp\left(-\alpha k t\right)$. For that reason, $k$ is also known as the extinction coefficient.

There exist alternate conventions of this description where the direction of the time dependence is reversed ($\exp\left(+i\alpha\eta t\right)$) and as a consequence, the sign of the complex component is flipped to ($n+ki$) [Muller69]. We use the former convention with the negative sign ($n-ki$) in this document.

into the following new phrases:

The same expressions also hold in the conductor case where the indices of refraction are complex valued, i.e. $\eta = n-ki$ with real and imaginary parts $n$ and $k$. The use of complex numbers here comes from the expression $\exp\left(+i\alpha t\right)$ that describes a plane wave propagating with increasing time and a constant related to the wave frequency. When introducing a value of $k>0$, this now leads to an exponential decay of the wave as it travels inside the conductor material: $\exp\left(i\alpha\left( t - \frac{\eta}{c}x\right)\right)=\exp\left(i\alpha\left( t - \frac{n}{c}x\right)\right)\exp\left( - \frac{\alpha k}{c}x\right)$. We refer to Section 4.8 in [Hecht] for derivations. For that reason, $k$ is also known as the extinction coefficient.

There exist alternate conventions of this description where the direction of the time dependence is reversed ( $\exp\left(-i\alpha t\right)$ ) and as a consequence, the sign of the complex component is flipped to ($n+ki$) (Table 4 on page 17 of [Muller69]). We use the former convention with the negative sign ($n-ki$) in this document.

Remained issue: the imaginary sign for "cosine refracted angle"

One of the reasons I wrote the long derivation is that the confusing imaginary signs appear in both complex refractive index and complex refracted angle.

For complex refractive index, it is just fine to change the documentation text (mainly "comes from the expression $\exp\left(-i\alpha\eta t\right)$" to "comes from the expression $\exp\left(+i\alpha t\right)$"). The source code is currently uses $n-ki$ IOR values consistently.

Our discussion during this several weeks has mainly focused on the complex refractive index issue. However, for complex refracted angle for total internal reflection and reflection on conductors, what I found is that the values produced by the source code appear to be incorrect, and the code should be revised.

If you (or one of the people who manage Mitsuba 3) agree the problem in computing complex refracted angles in https://github.com/mitsuba-renderer/mitsuba3/blob/3013adb4f12a491f8dd37c32bcedf55c7998f9e8/include/mitsuba/render/fresnel.h#L227, I can make a PR for this issue.

except a minor errata $\Im\left[\eta\right]=k$ -> $\left|\Im\left[\eta\right]\right| = k$)

Good catch! And thanks for the additional explanations. So technically the time component is not needed, but it implicitly defines the sign of the spatial component. I think given all this discussion it's safe to say this maybe requires a bit more explanation in the documentation after all.

I suggest the following text:

The same expressions also hold in the conductor case where the indices of refraction are complex valued, i.e. $\eta = n - k i$ with real and imaginary parts $n$ and $k$. The use of complex values here is a mathematical trick to encode both the harmonic oscillation of a wave together with its decay when travelling into the conductive medium. As illustration, consider a wave travelling along spatial coordinate $x$ and time $t$:

$$ \exp\left(i\omega(t - \frac{\eta}{c}x)\right) = \exp\left( i\omega(t - \frac{n}{c}x)\right)\cdot\exp\left(-\omega\frac{k}{c}x\right) $$

$\omega$ denotes the angular spatial frequency and $c$ is the speed of light. A value $k > 0$ introduces an exponential falloff term---which is why $k$ is also known as the extinction coefficient. We refer to "Optical Properties of Metals" in [Hecht] (Section 4.8 in the 5th edition) for a derivation and further details.

There also exist alternate conventions of the above behaviour where the direction of time dependence is reversed, i.e. $\exp\left(i\omega(\frac{\eta}{c}x - t)\right)$. As a consequence, the sign of the imaginary part needs to flip accordingly to $n + k i$. These are sometimes referred to as the " $\exp(+i \omega t)$ and $\exp(-i \omega t)$ conventions " in the literature [Muller1969]. In this document we use the former, so $\eta = n - k i$.

Yes, please open a separate issue or PR for the actual code issue. Aplogogies for delaying that (arguably more important) discussion a bit. I just wanted to first clarify the confusion around the documentation.

Your suggestion for the text is good and I cannot find any problem, except a minor suggestion.

There also exist alternate conventions of the above behaviour where the direction of time dependence is reversed, i.e.

This sentence is not wrong, it's fine. But I'm a little concerned that the above phrase might be misread as some physical meaning in time. My suggestion is

There also exist alternate conventions of the above behaviour with a formula of the complex conjugate, i.e.

I think this is slightly better to clarify that the difference between those two formulae is nothing but a convention, without any change in physical meaning and behaviour.

This discussion has been too long so that it seems difficult to read the two issues (complex IOR, complex refracted angle) separately. From this issue, we can conclude that Mitsuba 3 uses the $\exp\left(+i\omega t\right)$ convention and $\eta = n-ik$ IOR. Based on this conclusion, I will open a new PR with reorganising only the content about complex refracted angles from this issue, also including my suggestion of code correction. The new PR may take a couple of days.

Sounds good. I'll prepare a PR for the documentation changes tomorrow.

I have created PR #1161 . The bug of the previous implementation of fresnel_polarized(Float, dr::Complex<Float>) is much obscure as I expected since it works for some range of copmlex IOR but produces wrong values for some other range. This PR includes such a discussion.

Sorry that this hasn't been getting too much attention. We're undergoing a big refactoring, we're nearing the end of it. I should then have the time to read through your thorough discussions. :pray:

No problem. Thanks for taking care of this issue.

mitsuba-renderer / mitsuba3