Farfield and Purcell Factor for the usecase of a uLED

[chuanhong98]

Dear all,

thank you for the help and responses on my discussion of various sources to approximate a dipole last time, I am currently still working on it.

In the process of simulating and optimizing an uLED for my master's thesis I wish to analyze the purcell factor and the farfield pattern for the enhancement in spontaneous emission and the directivity of the resulting uLED. Since a Gaussian function is used to approximate a dipole in the paper, is it perhaps possible to extract the purcell factor and have you tried attempting said task?

Other than that, I have also noticed that the a function for the farfield has already been implemented in fmmax. I was wondering if the number of points taken in a Brillouin zone have an impact on the accuracy of the far field pattern and if you have tested it out before. If yes, from which grid size and above the far field pattern are reliable and accurate?

I thank you in advance and look forward to hearing from you again. :)

[mfschubert]

Yes, our paper computed a per-dipole emitted power, which should be proportional to the Purcell factor. See this figure:

With regards to convergence of the farfield, this may depend on the details of the structure and a convergence study is suggested.

[mfschubert]

The closest thing to an example would be the tests. There is one test which computes the farfield for a single dipole, and compares it to an analytical expression for the farfield.

You may be able to modify this to suit your needs.

[chuanhong98]

Yes, the emitted power was normalized---if I recall correctly, it was to the peak emitted power of the unoptimized device.

Regarding your questions about the gaussian source, you should regard these functions as merely providing a desired spatial profile, which you can scale as needed for your particular calculation. The peak value is always 1, so it can be multiplied by whatever prefactor you like.

The differences of factors of 2 should all cancel, and it can be verified that the profile FWHM matches the desired value, e.g. using the code provided here. Here's an example of a Gaussian where FWHM of 0.2 was specified.

By the way, for computing an accurate Purcell factor, another option is to divide by the computed free space power of a dipole, rather than the analytical free space power. Of course they should agree, but one might be easier.

Hi Martin! I am just revisiting this post and wondering about the analytical free space power. I did try to obtain my Purcell factor by dividing the emitted power in a structure with the emitted power of a dipole in free space

$\text{Purcell Factor} = \frac{\text{Emitted Power of a dipole(with structure)}}{\text{Computed Power of a dipole (free space)}}$

and received a good agreement with results from MEEP with a dummy test structure I have made.

After varying the number in the FWHM and the wavelength of the dipole, I realized that a change in the FWHM of the Gaussian does not affect the $\text{Computed Power of a dipole (free space)}$ but a change in the wavelength affects the $\text{Computed Power of a dipole (free space)}$ more significantly. To this end I have looked into the literatures and found 2 equations that might help: the power of a Gaussian beam and the total power of an electric dipole in a homogeneous medium.

The power of a Gaussian beam:
$P = \frac{1}{2}I_0\pi\omega_0^2 = \frac{1}{2}\frac{\mid E_0\mid^2}{2\eta}\pi\omega_0^2 = \frac{1}{2}\frac{\mid E_0\mid^2 n}{2\eta_0} \pi\omega_0^2$

The total power of an electric dipole in a homogeneous medium:
$P = \frac{\mu_0 \omega^4 p_0^2}{12 \pi c n}$

Taking into account the norms that FMMAX uses, that $\varepsilon_0, \mu_0, c$ equals 0, I still do not obtain the analytical emitted power of a dipole in free space that I desire. May I know which is the equation that you first referred to when writing the code?

Thanks a lot in advance and I look forward to hearing back from you again.

[mfschubert]

Hi @chuanhong98 , could you set up the calculation that is giving you an unexpected result, and then specify exactly what you are expecting? You can paste the code here.

I gather that you are computing the free space power of a dipole, and are expecting a wavelength dependence that is not actually observed. Having the code will help make this discussion efficient. Thanks!

[chuanhong98]

import functools
import time
from typing import Sequence, Tuple

import jax.numpy as jnp

from fmmax import basis, fields, fmm, scattering, sources


def simulate_uled(
    permittivity_epi: complex = (2.4869 + 0.0j) ** 2,                           ### GaN for epi @ 450nm
    thickness_layer_1: float = 25.0,
    thickness_layer_2: float = 25.0,
    thickness_layer_3: float = 25.0,
    thickness_layer_4: float = 25.0,
    pitch: float = 1802.0,
    wavelength: float = 450.0,
    dipole_fwhm: float = 20.0,
    approximate_num_terms: int = 3000,
    truncation: basis.Truncation = basis.Truncation.CIRCULAR,
    formulation: fmm.Formulation = fmm.Formulation.JONES_DIRECT_FOURIER,
    brillouin_grid_shape: Tuple[int, int] = (1, 1),
) -> Tuple[
    jnp.ndarray,
    jnp.ndarray,
    # jnp.ndarray,
    # jnp.ndarray,
    # Tuple[jnp.ndarray, jnp.ndarray, jnp.ndarray],
]:
 
    primitive_lattice_vectors = basis.LatticeVectors(
        u=basis.X * pitch, v=basis.Y * pitch
    )

    permittivities_before_source = [
        jnp.asarray(permittivity_epi)[jnp.newaxis, jnp.newaxis],
        jnp.asarray(permittivity_epi)[jnp.newaxis, jnp.newaxis],
    ]
    permittivities_after_source = [
        jnp.asarray(permittivity_epi)[jnp.newaxis, jnp.newaxis],
        jnp.asarray(permittivity_epi)[jnp.newaxis, jnp.newaxis],
    ]

    thicknesses_before_source = [
        jnp.asarray(thickness_layer_1),
        jnp.asarray(thickness_layer_2),
    ]
    thicknesses_after_source = [
        jnp.asarray(thickness_layer_3),
        jnp.asarray(thickness_layer_4),
    ]


    expansion = basis.generate_expansion(
        primitive_lattice_vectors=primitive_lattice_vectors,
        approximate_num_terms=approximate_num_terms,
        truncation=truncation,
    )

    in_plane_wavevector = basis.brillouin_zone_in_plane_wavevector(
        brillouin_grid_shape, primitive_lattice_vectors
    )
    assert in_plane_wavevector.shape[-1] == 2
    assert in_plane_wavevector.ndim == 3

    eigensolve = functools.partial(
        fmm.eigensolve_isotropic_media,
        wavelength=jnp.asarray(wavelength),
        in_plane_wavevector=in_plane_wavevector,
        primitive_lattice_vectors=primitive_lattice_vectors,
        expansion=expansion,
        formulation=formulation,
    )

    # Perform the layer eigensolve for each layer in the stack. Use the fact that
    # the first layer after the source is identical to the layer just before the
    # source in order to avoid one eigensolve.
    layer_solve_results_before_source = tuple(
        [eigensolve(permittivity=p) for p in permittivities_before_source]
    )
    layer_solve_results_after_source = tuple(
        [eigensolve(permittivity=p) for p in permittivities_after_source]
    )

    # Compute interior scattering matrices.
    s_matrices_interior_before_source = scattering.stack_s_matrices_interior(
        layer_solve_results_before_source, thicknesses_before_source
    )
    s_matrices_interior_after_source = scattering.stack_s_matrices_interior(
        layer_solve_results_after_source, thicknesses_after_source
    )
    s_matrix_before_source = s_matrices_interior_before_source[-1][0]
    s_matrix_after_source = s_matrices_interior_after_source[-1][0]

    dipole_locations = [
            (offset, 0)
            for offset in jnp.linspace(200,800,7)
        ]
    dipoles = sources.gaussian_source(
        fwhm=jnp.asarray(dipole_fwhm),
        location=jnp.asarray(dipole_locations),
        in_plane_wavevector=in_plane_wavevector,
        primitive_lattice_vectors=primitive_lattice_vectors,
        expansion=expansion,
    )

    zeros = jnp.zeros_like(dipoles)
    jx = jnp.concatenate([dipoles, zeros, zeros], axis=-1)
    jy = jnp.concatenate([zeros, dipoles, zeros], axis=-1)
    jz = jnp.concatenate([zeros, zeros, dipoles], axis=-1)

    (
        bwd_amplitude_ambient_end,
        fwd_amplitude_before_start,
        bwd_amplitude_before_end,
        fwd_amplitude_after_start,
        bwd_amplitude_after_end,
        fwd_amplitude_substrate_start,
    ) = sources.amplitudes_for_source(
        jx=jx,
        jy=jy,
        jz=jz,
        s_matrix_before_source=s_matrix_before_source,
        s_matrix_after_source=s_matrix_after_source,
    )

    # Compute the Poynting flux just before the source, after the source,
    # and in the ambient, to allow calculation of the extraction efficiency.
    bwd_amplitude_before_start = fields.propagate_amplitude(
        amplitude=bwd_amplitude_before_end,
        distance=s_matrix_before_source.end_layer_thickness,
        layer_solve_result=s_matrix_before_source.end_layer_solve_result,
    )
    fwd_flux_before_start, bwd_flux_before_start = fields.directional_poynting_flux(
        forward_amplitude=fwd_amplitude_before_start,
        backward_amplitude=bwd_amplitude_before_start,
        layer_solve_result=s_matrix_before_source.end_layer_solve_result,
    )

    fwd_amplitude_after_end = fields.propagate_amplitude(
        amplitude=fwd_amplitude_after_start,
        distance=s_matrix_after_source.start_layer_thickness,
        layer_solve_result=s_matrix_after_source.start_layer_solve_result,
    )
    fwd_flux_after_end, bwd_flux_after_end = fields.directional_poynting_flux(
        forward_amplitude=fwd_amplitude_after_end,
        backward_amplitude=bwd_amplitude_after_end,
        layer_solve_result=s_matrix_after_source.start_layer_solve_result,
    )

    _, bwd_flux_ambient_end = fields.directional_poynting_flux(
        forward_amplitude=jnp.zeros_like(bwd_amplitude_ambient_end),
        backward_amplitude=bwd_amplitude_ambient_end,
        layer_solve_result=s_matrix_before_source.start_layer_solve_result,
    )

    # Compute the total forward and backward flux resulting from the source. The
    # forward flux from the source is the difference between the forward flux just
    # after the source, and the forward flux just before the source. The backward
    # flux is defined analogously.
    fwd_flux_from_source = fwd_flux_after_end - fwd_flux_before_start
    bwd_flux_from_source = bwd_flux_before_start - bwd_flux_after_end

    # Compute total power emitted by the source by summing over the Brillouin
    # zone points and Fourier orders.
    assert fwd_flux_from_source.shape == brillouin_grid_shape + (
        2 * expansion.num_terms,
        jx.shape[-1],
    )
    forward_emitted_power = jnp.sum(fwd_flux_from_source, axis=(-4, -3, -2))
    backward_emitted_power = -jnp.sum(bwd_flux_from_source, axis=(-4, -3, -2))
    total_emitted_power = forward_emitted_power + backward_emitted_power

    total_extracted_power = -jnp.sum(bwd_flux_ambient_end, axis=(-4, -3, -2))
    assert total_extracted_power.shape == (jx.shape[-1],)

    # Calculate the per-dipole extraction efficiency.
    extraction_efficiency = total_extracted_power / total_emitted_power

    return dipole_locations, extraction_efficiency, total_emitted_power,

n_terms = jnp.array([])
time_array = jnp.array([])
extraction_efficiency_array = jnp.empty((0, 21))
power_array = jnp.empty((0, 21))

if __name__ == "__main__":

    for approximate_num_terms in (200, 500, 1000, 1200, 1400, 1600, 1800, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000):
        t0 = time.time()
        dipole_locations, extraction_efficiency, total_emitted_power = simulate_uled(
            approximate_num_terms=approximate_num_terms, brillouin_grid_shape=(1,1)
        )
        print(
            f"Results with approximate_num_terms = {approximate_num_terms} (time = {time.time() - t0:.2f}s):\n"
            f"    Extraction efficiency = {jnp.squeeze(extraction_efficiency * 100)}%\n"
            f"      Total emitted power = {jnp.squeeze(total_emitted_power)}"
        )
        n_terms = jnp.append(n_terms, jnp.array([approximate_num_terms]))
        time_array = jnp.append(time_array, jnp.array([time.time() - t0]))
        extraction_efficiency_array = jnp.append(extraction_efficiency_array, jnp.array([jnp.squeeze(extraction_efficiency * 100)]), axis=0)
        power_array = jnp.append(power_array, jnp.array([jnp.squeeze(total_emitted_power)]), axis=0)
        jnp.savez("GaN_uniform_FWHM_20nm_wavelength_450nm.npz", dipole_locations=dipole_locations, n_terms=n_terms, time_array=time_array, extraction_efficiency_array=extraction_efficiency_array, power_array=power_array)

Hello Martin, attached above is the code that I modified to calculate the emitted power of a dipole in a homogeneous medium. The 3 factors that affect this emitted power are: pitch size (crucial in determination of primitive lattice vector), FWHM of Gaussian and wavelength, which are intuitive. And I was mistaken before, the FWHM does affect the emitted power but not in the way I expected it to.

As the FWHM is increased, the emitted power decreases.

As the wavelength decreases, the emitted power increases.

I was just wondering which analytical formula does it corresponds to (whether the power of a dipole / power of a Gaussian beam) and how to predict this behaviour. Thank you!

[mfschubert]

OK, I I created a colab notebook which does the required calculation. It uses BZ integration in order to better replicate the situation of a dipole in free space, and I do see roughly the expected behavior.

To get this result, I had to account for some wavelength dependence in the dipole moment. As it stands, it seems the dipole moments of dipole sources produced by the gaussian_source function will have some wavelength dependence (since the function does not take wavelength as an argument).

Perhaps this correction was omitted in your code?