Fiber Bragg Gratings: Coupled Mode Theory and Transfer Matrix Method

The spectral response of a Fiber Bragg Grating (FBG) relies on the coupling between forward and backward propagating core modes, governed by the Bragg condition (λ_B = 2 * n_eff * Λ). In MATLAB, the Transfer Matrix Method (TMM) is the standard approach for calculating reflection and transmission spectra.

Theoretical Background

In TMM, the grating is discretized into small uniform segments. Each segment is represented by a 2×2 matrix, and the overall system response is found by multiplying these matrices. For a uniform segment, the matrix is defined as:

M = [ cosh(γL) - i(δ/γ)sinh(γL) , -i(κ/γ)sinh(γL) ; i(κ/γ)sinh(γL) , cosh(γL) + i(δ/γ)sinh(γL) ]

Where γ = sqrt(κ² - δ²). Here, κ is the coupling coefficient, δ = β - π/Λ is the detuning factor, and L is the segment length.

MATLAB Implementation: Uniform FBG

The following script calculates and plots the reflection and transmission spectra for a uniform FBG using TMM logic.

% Initialization of physical parameters
neff = 1.45;           % Effective refractive index
period = 530e-9;       % Grating period in meters
gratingLength = 10e-3; % Total length of the grating
deltaIndex = 1e-4;     % Index modulation depth
wavelengths = linspace(1500e-9, 1600e-9, 1000); % Scan range

% Calculate coupling and detuning
couplingK = pi * deltaIndex ./ wavelengths;
detuningDelta = 2 * pi * neff * (1 ./ wavelengths - 1 / (2 * period));

% TMM Calculation
gammaVal = sqrt(couplingK.^2 - detuningDelta.^2);
coshTerm = cosh(gammaVal * gratingLength);
sinhTerm = sinh(gammaVal * gratingLength);

% Spectral Power
reflectance = (sinhTerm.^2) ./ (coshTerm.^2 - (detuningDelta ./ couplingK).^2);
transmittance = 1 ./ (coshTerm.^2 - (detuningDelta ./ couplingK).^2);

% Plotting Results
figure;
subplot(2,1,1);
plot(wavelengths * 1e9, reflectance, 'LineWidth', 2);
xlabel('Wavelength (nm)');
ylabel('Reflectance');
title('Uniform FBG Reflection Spectrum');
grid on;

subplot(2,1,2);
plot(wavelengths * 1e9, transmittance, 'LineWidth', 2);
xlabel('Wavelength (nm)');
ylabel('Transmittance');
title('Uniform FBG Transmission Spectrum');
grid on;

Extension: Chirped FBG Simulation

For a Chirped FBG where the period varies along the length, the grating is sliced into sections. Each section uses a local period Λ_i for its matrix calculation.

% Chirped FBG Parameters
p0 = 530e-9;           % Starting period
chirpRate = 0.1e-9;    % Period change per meter
totalLen = 10e-3;      % Length
pos = linspace(0, totalLen, 100); % Position vector
periods = p0 + chirpRate * pos;   % Period distribution

segLen = totalLen / 100;
M_accum = eye(2);      % Accumulated matrix

for idx = 1:100
    p_local = periods(idx);
    k_local = pi * deltaIndex / wavelengths; 
    d_local = 2 * pi * neff * (1 ./ wavelengths - 1 / (2 * p_local));
    g_local = sqrt(k_local.^2 - d_local.^2);
    
    % Local transfer matrix
    M_seg = [cosh(g_local*segLen) - 1i*d_local/g_local*sinh(g_local*segLen), ...
             -1i*k_local/g_local*sinh(g_local*segLen); ...
             1i*k_local/g_local*sinh(g_local*segLen), ...
             cosh(g_local*segLen) + 1i*d_local/g_local*sinh(g_local*segLen)];
    
    M_accum = M_seg * M_accum; % Matrix multiplication
end

% Extract coefficients
r_coeff = M_accum(2,1) / M_accum(1,1);
t_coeff = 1 / M_accum(1,1);
R_final = abs(r_coeff)^2;
T_final = abs(t_coeff)^2;

Planar Gratings: Fraunhofer Diffraction

Planar grating spectra are determined by diffraction from a periodic structure. The far-field intensity for an amplitude grating is described by the Fraunhofer diffraciton equation.

I(θ) = I₀ [ sin(Nπd sinθ/λ) / N sin(πd sinθ/λ) ]² [ sin(πa sinθ/λ) / sin(πd sinθ/λ) ]²

Where N is the number of slits, d is the period, a is the slit width, and θ is the diffraction angle.

MATLAB Implementation: Planar Grating

% Parameters for Planar Grating
lambda0 = 632.8e-9;    % He-Ne Laser wavelength
gratPeriod = 1e-6;     % Grating period
slitWidth = 0.5e-6;    % Slit width
numSlits = 100;        % Number of slits
angles = linspace(-pi/6, pi/6, 1000); % Angular range

% Intensity Calculation
I0 = 1;
interference = (sin(numSlits * pi * gratPeriod * sin(angles) / lambda0) ./ ...
                (numSlits * sin(pi * gratPeriod * sin(angles) / lambda0))).^2;
diffraction = (sin(pi * slitWidth * sin(angles) / lambda0) ./ ...
               sin(pi * gratPeriod * sin(angles) / lambda0)).^2;
intensity = I0 * interference .* diffraction;

% Visualization
figure;
plot(angles * 180/pi, intensity, 'LineWidth', 2);
xlabel('Diffraction Angle (degrees)');
ylabel('Relative Intensity');
title('Far-Field Diffraction Pattern');
grid on;

Volume Holographic Gratings (VHG): Kogelnik Theory

VHG performance is based on volume Bragg diffraction. Kogelnik's coupled-wave theory provides an analytical solution for diffraction efficiency in thick media.

For a reflection VHG, the efficiency is:

η = sinh²(ψ) / [cosh²(ψ) - (ξ/ψ)²]

Where ψ = π Δn d / (λ cosθ) and ξ is the dephasing parameter.

MATLAB Implementation: VHG Efficiency

% VHG Parameters
wl = 532e-9;           % Wavelength
n_media = 1.5;         % Substrate index
modIndex = 0.01;       % Index modulation
thickness = 10e-6;     % Grating thickness
periodVHG = 1e-6;      % Grating period
theta_in = 30;         % Incident angle in degrees

% Convert to radians
theta_rad = deg2rad(theta_in);

% Wave vector calculations
K_vec = 2 * pi / periodVHG;
k0 = 2 * pi / wl;
tau = K_vec * cos(0 - theta_rad) - 2 * pi * n_media / wl; 

% Kogelnik parameters
xi_param = thickness * tau / (2 * cos(theta_rad));
psi_param = pi * modIndex * thickness / (wl * cos(theta_rad));

% Diffraction Efficiency
eff = (sinh(psi_param)^2) / (cosh(psi_param)^2 - (xi_param / psi_param)^2);
fprintf('VHG Diffraction Efficiency: %.2f%%\n', eff * 100);

Comparison of Grating Simulation Techniques

Implementation Notes

• Numerical Stability: When the coupling coefficient is much larger than the detuning (κ >> δ), γ may become imaginary. Ensure your code handles complex numbers correct using hyperbolic functions.
• Polarization: For planar and sub-wavelength gratings, simulate both TE and TM polarizations separately as the response is polarization-dependent.
• Visualization: Always label the Bragg wavelength, 3dB bandwidth, and peak reflectance/transmittance in you're plots for clarity.