"1d" Green's functions for near2far

[stevengj]

For a "1d" system described by a 2d (or 3d) simulation with Bloch-periodic boundary conditions in the x (or xy) directions, the Green's function is simply a planewave. It would be nice to support this directly in the near2far feature.

(We support a finite periodicity by a lattice sum, but in the limit of zero periodicity this is absurdly inefficient … the lattice sum effectively goes to a quadrature rule. We should just use the planewave formula.)

This is potentially useful for things like computing the Green's function (dipole emission pattern) for a multilayer film, which corresponds to simply integrating this "1d" near2far field over kx and ky (inside the light code). Of course, there are many specialized methods to compute multilayer Green's functions that will be vastly more efficient than Meep for this purpose (for free/open-source code, see e.g. https://github.com/modelics/strata and https://github.com/odmiller/MultilayerEM), but it shouldn't be too expensive in Meep.

[stevengj]

The first step is to write out all of the Green's functions in "1d".

For example, consider an electric current sheet $\vec{J} = \delta(z) \hat{x} e^{i(k_x x + k_y y) - i \omega t}$. This produces two planewaves, one for $z > 0$ propagating in the $+z$ direction, and one in $z < 0$ propagating in $-z$. These planewaves are proportional to $\sim \exp(i(k_x x + k_y y \pm k_z z))$ where $k_z = \sqrt{(\omega/c)^2 \varepsilon \mu - k_x^2 - k_y^2}$. (For $k_y=0$, by symmetry it must be "p" polarized with electric field in the $xz$ plane, but for $k_y \ne 0$ it's more messy.) One then just matches boundary conditions (there is a jump in the tangential H field across an electric current sheet) to work out the amplitudes of these two planewave and hence the field everywhere in space.

Similarly for an electric current polarized in the y direction, and similarly (swapping E and H) for magnetic currents.

Probably someone has written this all out somewhere if you can find it, but it is also not too hard to re-derive.

[oskooi]

The equation for planewave propagation through a 1D multilayer stack using the transfer-matrix method is derived in Section 4.1 "Complex amplitudes for reflection and transmission" of the notes "Multilayer optical calculations" by Steven J. Byrnes in arXiv:1603.02720:

[stevengj]

That's not what we're talking about. We're talking about (a) a homogeneous medium and (b) off-axis propagation.

For on-axis propagation, this is just the $e^{i\delta_n}$ term in the expressions you quoted (which are on-axis). We don't have multiple layers in near2far.

[stevengj]

There should be papers somewhere on the 4x4 off-axis transfer matrices that work all this out (e.g. I know Yariv did it in the cylindrical case, and the off-axis planar case might be in his book someplace). You still have to relate it to the current sources, but that's not too hard. (The whole thing is not too hard in principle, just tedious.)

[stevengj]

(Once you get the $k_y=0$ solution, you can just rotate them around the $z$ axis to get the other directions.)

[stevengj]

The underlying principle here is Brillouin zone integration, also called the "array scanning method" because you can think of it as a phased-array antenna that you integrate over the phases. On review article is Capolino et al. (2007). If you put in a "dipole" with Bloch-periodic boundary conditions (really a phased array of dipoles) and you integrate k over the Brillouin zone, you get the response of an isolated dipole.

In the case of continuous translational symmetry (zero period), the Brillouin zone is infinite, so you are just integrating over all k. That's the case for multilayer films.

[oskooi]

Probably someone has written this all out somewhere if you can find it, but it is also not too hard to re-derive.

Looks like the solution is derived in Section 1.5 "General Plane Wave Solutions" of Microwave Engineering (4th edition) by Prof. David M. Pozar. In the MKS units used in the book, the intrinsic impedance of free space is $\eta_0 = \sqrt{\mu_0 / \varepsilon_0} = 377$ Ω (page 17 in Section 1.4 "Planewaves in a Losless Medium") . $k_0 = |\vec{k}|$ is the wavevector in free space (equation 1.67 on page 21).

[oskooi]

Here is a summary of the calculation of the 1D Green's function.

[stevengj]

Again, you don't need to calculate $k_y \ne 0$ separately — it is equivalent to the $k_y = 0$ solution under a rotation.

To compute the far field, where you can neglect evanescent-wave contributions, you can just integrate $k_x, k_y$ over the ball $k_\Vert \le \omega n / c$ (for an ambient medium with index $n$).

[oskooi]

1. Set up the 1D simulation with $(k_x, k_y)$ Bloch-periodic boundaries. Time-step the fields and compute the DFT fields ($E_x, E_y, H_x, H_y$) at a single point at the monitor location. Meep calculations are done for a single fixed dipole (fixed source component).
2. Using the Principle of Equivalence, convert the DFT fields into equivalent electric and magnetic currents $J_x, J_y$ and $K_x, K_y$.
3. Rotate the coordinate system around the $z$ axis such that $k_y = 0$. This includes rotating the far-field points (i.e., hemisphere) and rotating the equivalent currents (i.e., multiply the currents by the rotation matrix). In the new coordinate system, we have the $\mathcal{S}$ and $\mathcal{P}$ decomposition.
4. Use the analytic formulas to obtain the amplitudes of the $\mathcal{S}$ and $\mathcal{P}$ planewaves from the equivalent currents.
5. Evaluate the planewaves (with the correct amplitudes) at the rotated far-field points to obtain $\vec{E}$ and $\vec{H}$ in the far field. If the original points are on a hemisphere then the rotated points are also on a hemisphere because we are rotating around the $z$ axis.
6. Rotate $\vec{E}$ and $\vec{H}$ back into the original location (i.e. multiply by the inverse of the rotation matrix from step 3).
7. Add the fields from step 6 to the running tally of the total fields at each point.

Repeat steps 1 - 7 by integrate over all $(k_x, k_y)$ in the light cone. Use cylindrical coordinates to integrate over $(k_r, k_\phi)$ using equally spaced points for $k_\phi$ (because it is periodic) and Gaussian quadrature points for $k_r$. (For periodic systems, we would be use equally spaced grid points but there is no periodicity here.)

[oskooi]

As a test, verify that the results from the 1D calculation for vacuum match the analytical Green's functions (green3d).

[oskooi]

Here is a first version of the calculation.

"""Radiation profile of an Ex dipole in free space."""

from enum import Enum
from typing import Tuple
import math
import warnings

import meep as mp
import numpy as np


Current = Enum("Current", "Jx Jy Kx Ky")

RESOLUTION_UM = 40
WAVELENGTH_UM = 1.0
FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM
DEBUG_OUTPUT = True

frequency = 1 / WAVELENGTH_UM


def dipole_in_vacuum(
        current_component: Current, kx: float, ky: float
) -> Tuple[complex, complex, complex, complex]:
    pml_um = 1.0
    air_um = 10.0
    size_z_um = pml_um + air_um + pml_um
    cell_size = mp.Vector3(0, 0, size_z_um)

    print(f"dipole_in_vacuum: kx = {kx}, ky = {ky}")

    pml_layers = [mp.Absorber(thickness=pml_um)]

    if current_component.name == "Jx":
        src_cmpt = mp.Ex
    elif current_component.name == "Jy":
        src_cmpt = mp.Ey
    elif current_component.name == "Kx":
        src_cmpt = mp.Hx
    else:
        src_cmpt = mp.Hy

    sources = [
        mp.Source(
            src=mp.GaussianSource(frequency, fwidth=0.1 * frequency),
            component=src_cmpt,
            center=mp.Vector3()
        )
    ]

    sim = mp.Simulation(
        resolution=RESOLUTION_UM,
        force_complex_fields=True,
        cell_size=cell_size,
        sources=sources,
        boundary_layers=pml_layers,
        k_point=mp.Vector3(kx, ky, 0)
    )

    dft_fields = sim.add_dft_fields(
        [mp.Ex, mp.Ey, mp.Hx, mp.Hy],
        frequency,
        0,
        1,
        center=mp.Vector3(0, 0, 0.5 * size_z_um - pml_um),
        size=mp.Vector3()
    )

    sim.run(
        until_after_sources=mp.stop_when_fields_decayed(
            10, src_cmpt, mp.Vector3(0, 0, 0.5423), 1e-6
        )
    )

    ex_dft = sim.get_dft_array(dft_fields, mp.Ex, 0)
    ey_dft = sim.get_dft_array(dft_fields, mp.Ey, 0)
    hx_dft = sim.get_dft_array(dft_fields, mp.Hx, 0)
    hy_dft = sim.get_dft_array(dft_fields, mp.Hy, 0)

    if DEBUG_OUTPUT:
        print(f"ex_dft:, {ex_dft.real:.4f}{ex_dft.imag:+.4f}*1j")
        print(f"ey_dft:, {ey_dft.real:.4f}{ey_dft.imag:+.4f}*1j")
        print(f"hx_dft:, {hx_dft.real:.4f}{hx_dft.imag:+.4f}*1j")
        print(f"hy_dft:, {hy_dft.real:.4f}{hy_dft.imag:+.4f}*1j")

    return ex_dft, ey_dft, hx_dft, hy_dft


def equivalent_currents(
        ex_dft: complex, ey_dft: complex, hx_dft: complex, hy_dft: complex
) -> Tuple[Tuple[complex, complex], Tuple[complex, complex]]:
    """Returns the equivalent electric and magnetic currents as a 2-tuple."""

    electric_current = (-hy_dft, hx_dft)
    magnetic_current = (ey_dft, -ex_dft)

    return electric_current, magnetic_current


def far_fields(
        kz: float,
        rz: float,
        current_amplitude: complex,
        current_component: Current,
) -> Tuple[complex, complex]:
    """Computes the S- or P-polarized far fields from a sheet current."""

    if current_component.name == "Jx":
        # Jx --> (Ex, Hy) [P pol.]                                                                               
        e_far = -0.5 * current_amplitude
        h_far = -0.5 * current_amplitude
    elif current_component.name == "Jy":
        # Jy --> (Hx, Ey) [S pol.]                                                                               
        e_far = -0.5 * current_amplitude
        h_far = 0.5 * current_amplitude
    elif current_component.name == "Kx":
        # Kx --> (Hx, Ey) [S pol.]                                                                               
        e_far = 0.5 * current_amplitude
        h_far = -0.5 * current_amplitude
    else:
        # Ky --> (Ex, Hy) [P pol.]                                                                               
        e_far = 0.5 * current_amplitude
        h_far = 0.5 * current_amplitude

    phase = np.exp(1j * kz * rz)
    e_far *= phase
    h_far *= phase

    return e_far, h_far


def spherical_to_cartesian(
        polar_rad,
        azimuthal_rad
) -> Tuple[float, float, float]:
    """Converts a point in spherical to Cartesian coordinates."""

    rx = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.cos(azimuthal_rad)
    ry = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.sin(azimuthal_rad)
    rz = FARFIELD_RADIUS_UM * np.cos(polar_rad)

    return rx, ry, rz

if __name__ == "__main__":
    # Jx --> (Ex, Hy) [P pol.]                                                                                   
    dipole_component = Current.Jx

    num_k = 21
    kxs = np.linspace(-frequency, frequency, num_k)
    kys = np.linspace(-frequency, frequency, num_k)

    # Far fields are defined on the surface of a hemisphere.                                                     
    num_polar = 20
    num_azimuthal = 30
    polar_rad = np.linspace(0, 0.5 * np.pi, num_polar)
    azimuthal_rad = np.linspace(0, 2 * np.pi, num_azimuthal)

    current_components = [Current.Jx, Current.Jy, Current.Kx, Current.Ky]
    farfield_components = ['Ex', 'Ey', 'Hx', 'Hy']

    total_farfields = {}
    for component in farfield_components:
        total_farfields[component] = np.zeros((num_polar, num_azimuthal), dtype=np.complex128)

    for kx in kxs:
        for ky in kys:

            # Skip wavevectors which are outside the light cone.                                                 
            if np.sqrt(kx**2 + ky**2) >= frequency:
                continue

            ex_dft, ey_dft, hx_dft, hy_dft = dipole_in_vacuum(dipole_component, kx, ky)

            # Rotation angle around z axis to force ky = 0. 0 is +x.                                             
            if kx:
                rotation_rad = -math.atan(ky / kx)
            else:
                if ky:
                    rotation_rad = -0.5 * np.pi
                else:
                    rotation_rad = 0

            if DEBUG_OUTPUT:
                print(f"rotation angle:, {math.degrees(rotation_rad):.2f}°")

            rotation_matrix = np.array(
                [
                    [np.cos(rotation_rad), -np.sin(rotation_rad)],
                    [np.sin(rotation_rad), np.cos(rotation_rad)]
                ]
            )

            electric_current, magnetic_current = equivalent_currents(
                ex_dft,
                ey_dft,
                hx_dft,
                hy_dft
            )

            if DEBUG_OUTPUT:
                print(f"electric_current:, {electric_current}")
                print(f"magnetic_current:, {magnetic_current}")

            electric_current_rotated = (
                rotation_matrix @ np.transpose(
                    np.array(
                        [
                            electric_current[0],
                            electric_current[1]
                        ]
                    )
                )
            )

            magnetic_current_rotated = (
                rotation_matrix @ np.transpose(
                    np.array(
                        [
                            magnetic_current[0],
                            magnetic_current[1]
                        ]
                    )
                )
            )

            if DEBUG_OUTPUT:
                print(f"electric_current_rotated:, {electric_current_rotated}")
                print(f"magnetic_current_rotated:, {magnetic_current_rotated}")

            kz = (frequency**2 - kx**2 - ky**2)**0.5

            if DEBUG_OUTPUT:
                print(f"(kx, ky, kz) = ({kx:.6f}, {ky:.6f}, {kz:.6f})")

            # Verify that ky of rotated wavevector is 0.                                                         
            k_rotated = rotation_matrix @ np.transpose(np.array([kx, ky]))
            print(f"k_rotated = ({k_rotated[0]:.2f}, {k_rotated[1]:.2f})")
            if k_rotated[1] == 0:
                print("rotated ky is zero.")

            current_amplitudes = [
                electric_current_rotated[0],
                electric_current_rotated[1],
                magnetic_current_rotated[0],
                magnetic_current_rotated[1],
            ]

            farfields = {}
            for component in farfield_components:
                farfields[component] = np.zeros((num_polar, num_azimuthal), dtype=np.complex128)

            for i, polar in enumerate(polar_rad):
                for j, azimuthal in enumerate(azimuthal_rad):
                    rx, ry, rz = spherical_to_cartesian(polar, azimuthal)

                    for current_component, current_amplitude in zip(current_components, current_amplitudes):
                        if abs(current_amplitude) == 0:
                            continue

                        e_far, h_far = far_fields(
                            kz,
                            rz,
                            current_amplitude,
                            current_component,
                        )

                        if current_component.name == "Jx" or current_component.name == "Ky":
                            # P polarization                                                                     
                            farfields['Ex'][i, j] += e_far
                            farfields['Hy'][i, j] += h_far
                        elif current_component.name == "Jy" or current_component.name == "Kx":
                            # S polarization                                                                     
                            farfields['Ey'][i, j] += e_far
                            farfields['Hx'][i, j] += h_far

            # apply inverse rotation matrix to total far fields                                                  
            inverse_rotation_matrix = np.linalg.inv(rotation_matrix)

            for i in range(num_polar):
                for j in range(num_azimuthal):
                    total_farfields['Ex'][i, j] = (
                        inverse_rotation_matrix[0, 0] * farfields['Ex'][i, j] +
                        inverse_rotation_matrix[0, 1] * farfields['Ey'][i, j]
                    )
                    total_farfields['Ey'][i, j] = (
                        inverse_rotation_matrix[1, 0] * farfields['Ex'][i, j] +
                        inverse_rotation_matrix[1, 1] * farfields['Ey'][i, j]
                    )
                    total_farfields['Hx'][i, j] = (
                        inverse_rotation_matrix[0, 0] * farfields['Hx'][i, j] +
                        inverse_rotation_matrix[0, 1] * farfields['Hy'][i, j]
                    )
                    total_farfields['Hy'][i, j] = (
                        inverse_rotation_matrix[1, 0] * farfields['Hx'][i, j] +
                        inverse_rotation_matrix[1, 1] * farfields['Hy'][i, j]
                    )

        if mp.am_master():
            np.savez(
                "dipole_farfields.npz",
                **total_farfields,
                RESOLUTION_UM=RESOLUTION_UM,
                WAVELENGTH_UM=WAVELENGTH_UM,
                FARFIELD_RADIUS_UM=FARFIELD_RADIUS_UM,
                num_polar=num_polar,
                num_azimuthal=num_azimuthal,
                polar_rad=polar_rad,
                azimuthal_rad=azimuthal_rad,
                kxs=kxs,
                kys=kys
            )

[stevengj]

A quick check: use exactly the method you are using now, but replace the integral over $k_x, k_y$ with just the $k_x = k_y = 0$ term. The result should be a planewave traveling vertically.

Similarly, you could keep just some other $k_x, k_y$ term and you should get a planewave traveling in another direction.

(You could look at the power, but even better just evaluate the far-field points on a grid and plot the planewave.)

[oskooi]

A quick check: use exactly the method you are using now, but replace the integral over $k_x, k_y$ with just the $k_x = k_y = 0$ term. The result should be a planewave traveling vertically.

I verified that a single ($k_x$, $k_y$) term of the Brillouin-zone integration in the setup produces a planewave with the correct fields and angle.

As an example, for a planewave with $\vec{k} = \left(\sin(45^\circ), 0, \sqrt{\omega^2 -\sin^2(45^\circ)}\right)$, a $\mathcal{J_x}$ current source produces nonzero $\mathcal{P}$-polarization fields for ($E_x$, $H_y$, $E_z$) where $\Re(E_x)$ in the far field in the $xz$ plane is shown below.