diff --git a/docs/photonics/examples/detector.py b/docs/photonics/examples/detector.py
new file mode 100644
index 00000000..6903a6cd
--- /dev/null
+++ b/docs/photonics/examples/detector.py
@@ -0,0 +1,244 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.15.2
+#   kernelspec:
+#     display_name: femwell
+#     language: python
+#     name: python3
+# ---
+
+# # Photodetector field profiles
+#
+# ```{caution}
+# **It seems that in its current formulation, the detector basis can only represent up to 95% of the input mode (see plot below), indicating it is not complete. This does not occur if the germanium index is set purely real.**
+# ```
+#
+# ```{caution}
+# **This example ignores the reflection to cladding at the end of the detector. For long enough detectors where most of the light is absorbed, this should be a small effect.**
+# ```
+#
+# ```{caution}
+# **In this example, the scipy solver returns diverging high-order modes. The slepc solver does not seem to have this issue.**
+# ```
+#
+# Mode solving can be used to calculate the optical intensity profile along the length of an absorber. This is useful to calculate optical generation terms for semiconductor simulations.
+
+# +
+from collections import OrderedDict
+
+import matplotlib.pyplot as plt
+import numpy as np
+import shapely
+from skfem import Basis, ElementTriP0
+from skfem.io.meshio import from_meshio
+from tqdm import tqdm
+
+from femwell.maxwell.waveguide import Mode, calculate_hfield, compute_modes
+from femwell.mesh import mesh_from_OrderedDict
+from femwell.visualization import plot_domains
+
+# -
+
+# Consider a typical germanium-on-silicon vertical photodetector profile, potentially with some sidewalls. We create a single mesh that has both the silicon input width as well as the larger silicon width under the germanium to avoid mesh interpolation errors:
+
+# +
+silicon_core_thickness = 0.22
+germanium_thickness = 0.5
+germanium_sidewall_angle = 10
+clad_vertical_offset = 3
+
+silicon_input_width = 1.5
+silicon_detector_width = 3
+germanium_width = 1
+clad_horizontal_offset = 3
+
+silicon_input = shapely.geometry.box(
+    -silicon_input_width / 2, -silicon_core_thickness, silicon_input_width / 2, 0
+)
+silicon_detector = shapely.geometry.box(
+    -silicon_detector_width / 2, -silicon_core_thickness, silicon_detector_width / 2, 0
+)
+germanium = shapely.Polygon(
+    (
+        (-germanium_width / 2, 0),
+        (
+            -germanium_width / 2
+            + germanium_thickness / np.tan(np.rad2deg(germanium_sidewall_angle)),
+            germanium_thickness,
+        ),
+        (
+            germanium_width / 2
+            - germanium_thickness / np.tan(np.rad2deg(germanium_sidewall_angle)),
+            germanium_thickness,
+        ),
+        (germanium_width / 2, 0),
+    )
+)
+clad = shapely.geometry.box(
+    -silicon_detector_width / 2 - clad_horizontal_offset,
+    -silicon_core_thickness - clad_vertical_offset,
+    silicon_detector_width / 2 + clad_horizontal_offset,
+    germanium_thickness + clad_vertical_offset,
+)
+
+polygons = OrderedDict(
+    silicon_input=silicon_input, silicon_detector=silicon_detector, germanium=germanium, clad=clad
+)
+
+resolutions = dict(
+    silicon_input={"resolution": 0.025, "distance": 5},
+    silicon_detector={"resolution": 0.025, "distance": 5},
+    germanium={"resolution": 0.025, "distance": 5},
+)
+
+mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.5))
+mesh.draw().show()
+
+plot_domains(mesh)
+plt.show()
+# -
+
+# By tagging the materials appropriately, we can generate the cross-section of the input silicon:
+
+# +
+wavelength = 1.55
+
+# Using refractive indices at 1.55 um
+si_index = 3.45
+ge_index = 4.6530  # + 1j * 0.29800
+clad_index = 1.44
+
+# +
+basis0 = Basis(mesh, ElementTriP0(), intorder=4)
+
+epsilon_input = basis0.zeros(dtype=complex)
+
+for subdomain, n in {
+    "silicon_input": si_index,
+    "silicon_detector": clad_index,
+    "germanium": clad_index,
+    "clad": clad_index,
+}.items():
+    epsilon_input[basis0.get_dofs(elements=subdomain)] = n**2
+
+fig, axs = plt.subplots(1, 2)
+for ax in axs:
+    ax.set_aspect(1)
+axs[0].set_title(r"$\Re\epsilon$, input")
+basis0.plot(epsilon_input.real, colorbar=True, ax=axs[0])
+axs[1].set_title(r"$\Im\epsilon$, input")
+basis0.plot(epsilon_input.imag, shading="gouraud", colorbar=True, ax=axs[1])
+plt.show()
+# -
+
+# The fundamental mode of this geometry will be our input mode:
+
+input_modes = compute_modes(
+    basis0,
+    epsilon_input,
+    wavelength=wavelength,
+    num_modes=1,
+    order=1,
+    radius=np.inf,
+    solver="slepc",
+)
+
+# +
+input_mode = input_modes[0]
+
+input_mode.plot(input_mode.E.real, colorbar=True, direction="x")
+plt.show()
+
+# +
+basis0 = Basis(mesh, ElementTriP0(), intorder=4)
+
+epsilon_detector = basis0.zeros(dtype=complex)
+
+for subdomain, n in {
+    "silicon_input": si_index,
+    "silicon_detector": si_index,
+    "germanium": ge_index,
+    "clad": clad_index,
+}.items():
+    epsilon_detector[basis0.get_dofs(elements=subdomain)] = n**2
+
+fig, axs = plt.subplots(1, 2)
+for ax in axs:
+    ax.set_aspect(1)
+axs[0].set_title(r"$\Re\epsilon$, detector")
+basis0.plot(epsilon_detector.real, colorbar=True, ax=axs[0])
+axs[1].set_title(r"$\Im\epsilon$, detector")
+basis0.plot(epsilon_detector.imag, shading="gouraud", colorbar=True, ax=axs[1])
+plt.show()
+# -
+
+detector_modes = compute_modes(
+    basis0,
+    epsilon_detector,
+    wavelength=wavelength,
+    num_modes=40,
+    order=1,
+    radius=np.inf,
+    solver="slepc",
+)
+
+for i, mode in enumerate(detector_modes):
+    if not i % 5:
+        print(f"Mode index: {i}")
+        mode.plot(mode.E.real, colorbar=True, direction="x")
+        plt.show()
+
+#
+
+# +
+overlaps = []
+ks = []
+
+sum_modes = detector_modes[0].basis.zeros(dtype=complex)
+
+for mode in tqdm(detector_modes):
+    overlaps.append(mode.calculate_overlap(input_mode))
+    ks.append(mode.k)
+
+    sum_modes += mode.E * mode.calculate_overlap(input_mode)
+# -
+
+detector_modes[0].show(sum_modes.real, colorbar=True)
+
+# Here we evaluate the completeness of the detector basis, and see that it fails to capture about 5% of the incoming mode. We can also see the effect of modes with different symmetry and polarization to the input mode (plateaus in cumulative overlap).
+
+plt.plot(np.cumsum(np.abs(overlaps) ** 2))
+plt.ylim([0, 1.1])
+plt.axhline(y=1, color="k", linestyle="--")
+plt.xlabel("Mode index")
+plt.ylabel("Cumulative power overlap")
+
+# Following eigenmode expansion, the propagating field inside the detector cross-section given a pure input mode can be expressed as (ignoring reflections):
+#
+# $$ E(x,y,z) = \sum_k \braket{E_k|E_0} e^{-i \beta_k z} \ket{E_k} $$
+#
+# where $\braket{\bm{x}|E_k} = E_k(x,y)$ are the cross-sectional mode profiles.
+
+# +
+z = 0
+
+
+def field_profile_at_z(z):
+    field = detector_modes[0].basis.zeros(dtype=complex)
+    for i, (k, overlap, mode) in enumerate(zip(ks, overlaps, detector_modes)):
+        field += overlap * np.exp(-1j * (np.real(k) - 1j * np.abs(np.imag(k))) * z) * mode.E
+    return field
+
+
+# -
+
+for z in np.linspace(0, 2, 11):
+    field = field_profile_at_z(z)
+    detector_modes[0].show(field.real, colorbar=True)
+
+#
diff --git a/docs/photonics/examples/detector_perturbation_symmetry.py b/docs/photonics/examples/detector_perturbation_symmetry.py
new file mode 100644
index 00000000..90fc24f8
--- /dev/null
+++ b/docs/photonics/examples/detector_perturbation_symmetry.py
@@ -0,0 +1,272 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.15.2
+#   kernelspec:
+#     display_name: femwell
+#     language: python
+#     name: python3
+# ---
+
+# # Photodetector field profiles
+#
+# ```{caution}
+# **It seems that in its current formulation, the detector basis can only represent up to 95% of the input mode (see plot below), indicating it is not complete. This does not occur if the germanium index is set purely real.**
+# ```
+#
+# ```{caution}
+# **This example ignores the reflection to cladding at the end of the detector. For long enough detectors where most of the light is absorbed, this should be a small effect.**
+# ```
+#
+# ```{caution}
+# **In this example, the scipy solver returns diverging high-order modes. The slepc solver does not seem to have this issue.**
+# ```
+#
+# Mode solving can be used to calculate the optical intensity profile along the length of an absorber. This is useful to calculate optical generation terms for semiconductor simulations.
+
+# +
+from collections import OrderedDict
+
+import matplotlib.pyplot as plt
+import numpy as np
+import shapely
+from skfem import Basis, ElementTriP0
+from skfem.io.meshio import from_meshio
+from tqdm import tqdm
+
+from femwell.maxwell.waveguide import Mode, calculate_hfield, compute_modes
+from femwell.mesh import mesh_from_OrderedDict
+from femwell.visualization import plot_domains
+
+# -
+
+# Consider a typical germanium-on-silicon vertical photodetector profile, potentially with some sidewalls. We create a single mesh that has both the silicon input width as well as the larger silicon width under the germanium to avoid mesh interpolation errors:
+
+# +
+silicon_core_thickness = 0.22
+germanium_thickness = 0.5
+germanium_sidewall_angle = 10
+clad_vertical_offset = 3
+
+silicon_input_width = 2
+silicon_detector_width = 4
+germanium_width = 2
+clad_horizontal_offset = 3
+
+silicon_input = shapely.geometry.box(0, -silicon_core_thickness, silicon_input_width / 2, 0)
+silicon_detector = shapely.geometry.box(0, -silicon_core_thickness, silicon_detector_width / 2, 0)
+germanium = shapely.Polygon(
+    (
+        (0, 0),
+        (0, germanium_thickness),
+        (
+            germanium_width / 2
+            - germanium_thickness / np.tan(np.rad2deg(germanium_sidewall_angle)),
+            germanium_thickness,
+        ),
+        (germanium_width / 2, 0),
+    )
+)
+clad = shapely.geometry.box(
+    0,
+    -silicon_core_thickness - clad_vertical_offset,
+    silicon_detector_width / 2 + clad_horizontal_offset,
+    germanium_thickness + clad_vertical_offset,
+)
+
+polygons = OrderedDict(
+    silicon_input=silicon_input, silicon_detector=silicon_detector, germanium=germanium, clad=clad
+)
+
+resolutions = dict(
+    silicon_input={"resolution": 0.05, "distance": 5},
+    silicon_detector={"resolution": 0.05, "distance": 5},
+    germanium={"resolution": 0.05, "distance": 5},
+)
+
+mesh = from_meshio(
+    mesh_from_OrderedDict(
+        polygons, resolutions, default_resolution_max=0.5, filename="half_mesh.msh"
+    )
+)
+mesh.draw().show()
+
+plot_domains(mesh)
+plt.show()
+# -
+
+# By tagging the materials appropriately, we can generate the cross-section of the input silicon:
+
+# +
+wavelength = 1.55
+
+# Using refractive indices at 1.55 um
+si_index = 3.45
+ge_index = 4.6530 + 1j * 0.29800
+clad_index = 1.44
+
+# +
+# Add Neumann boundary on the left, and Dirichlet (condense) elsewhere
+manual_boundaries = {
+    "bottom": lambda xi: xi[1] == 0.0,
+    "top": lambda xi: xi[0] == germanium_thickness + clad_vertical_offset,
+    "right": lambda xi: xi[1] == silicon_detector_width / 2 + clad_horizontal_offset,
+}
+mesh = mesh.with_boundaries(manual_boundaries)
+
+basis0 = Basis(mesh, ElementTriP0(), intorder=4)
+
+epsilon_input = basis0.zeros(dtype=float)
+
+for subdomain, n in {
+    "silicon_input": si_index,
+    "silicon_detector": clad_index,
+    "germanium": clad_index,
+    "clad": clad_index,
+}.items():
+    epsilon_input[basis0.get_dofs(elements=subdomain)] = n**2
+
+fig, axs = plt.subplots(1, 2)
+for ax in axs:
+    ax.set_aspect(1)
+axs[0].set_title(r"$\Re\epsilon$, input")
+basis0.plot(epsilon_input.real, colorbar=True, ax=axs[0])
+axs[1].set_title(r"$\Im\epsilon$, input")
+basis0.plot(epsilon_input.imag, shading="gouraud", colorbar=True, ax=axs[1])
+plt.show()
+# -
+
+# The fundamental mode of this geometry will be our input mode:
+
+input_modes = compute_modes(
+    basis0,
+    epsilon_input,
+    wavelength=wavelength,
+    num_modes=1,
+    order=1,
+    radius=np.inf,
+    solver="scipy",
+    with_manual_boundaries=manual_boundaries,
+)
+
+# +
+input_mode = input_modes[0]
+
+input_mode.plot(input_mode.E.real, colorbar=True, direction="x")
+plt.show()
+
+# +
+basis0 = Basis(mesh, ElementTriP0(), intorder=4)
+
+epsilon_detector = basis0.zeros(dtype=float)
+
+for subdomain, n in {
+    "silicon_input": si_index,
+    "silicon_detector": si_index,
+    "germanium": np.real(ge_index),
+    "clad": clad_index,
+}.items():
+    epsilon_detector[basis0.get_dofs(elements=subdomain)] = n**2
+
+fig, axs = plt.subplots(1, 2)
+for ax in axs:
+    ax.set_aspect(1)
+axs[0].set_title(r"$\Re\epsilon$, detector")
+basis0.plot(epsilon_detector.real, colorbar=True, ax=axs[0])
+axs[1].set_title(r"$\Im\epsilon$, detector")
+basis0.plot(epsilon_detector.imag, shading="gouraud", colorbar=True, ax=axs[1])
+plt.show()
+# -
+
+detector_modes = compute_modes(
+    basis0,
+    epsilon_detector,
+    wavelength=wavelength,
+    num_modes=100,
+    order=1,
+    radius=np.inf,
+    solver="scipy",
+)
+
+for i, mode in enumerate(detector_modes):
+    if not i % 5:
+        print(f"Mode index: {i}")
+        mode.plot(mode.E.real, colorbar=True, direction="x")
+        plt.show()
+
+detector_modes[-1].k
+
+# Next, we "complexify" those modes by assuming the imaginary index of the germanium is an index perturbation:
+
+# +
+basis0 = Basis(mesh, ElementTriP0(), intorder=4)
+
+epsilon_detector_complex_perturbation = basis0.zeros(dtype=complex)
+
+for subdomain, n in {
+    "silicon_input": 0,
+    "silicon_detector": 0,
+    "germanium": np.imag(ge_index),
+    "clad": 0,
+}.items():
+    epsilon_detector_complex_perturbation[basis0.get_dofs(elements=subdomain)] = 1j * n**2
+
+fig, axs = plt.subplots(1, 2)
+for ax in axs:
+    ax.set_aspect(1)
+axs[0].set_title(r"$\Re\epsilon$, detector")
+basis0.plot(epsilon_detector_complex_perturbation.real, colorbar=True, ax=axs[0])
+axs[1].set_title(r"$\Im\epsilon$, detector")
+basis0.plot(epsilon_detector_complex_perturbation.imag, shading="gouraud", colorbar=True, ax=axs[1])
+plt.show()
+
+# +
+overlaps = []
+ks = []
+
+sum_modes = detector_modes[0].basis.zeros(dtype=complex)
+
+for mode in tqdm(detector_modes):
+    overlaps.append(mode.calculate_overlap(input_mode))
+    ks.append(mode.calculate_pertubated_k(epsilon_detector + epsilon_detector_complex_perturbation))
+    sum_modes += mode.E * mode.calculate_overlap(input_mode)
+# -
+
+detector_modes[0].show(sum_modes.real, colorbar=True)
+
+# Here we evaluate the completeness of the detector basis, and see that it fails to capture about 5% of the incoming mode. We can also see the effect of modes with different symmetry and polarization to the input mode (plateaus in cumulative overlap).
+
+plt.plot(np.cumsum(np.abs(overlaps) ** 2))
+plt.ylim([0, 1.1])
+plt.axhline(y=1, color="k", linestyle="--")
+plt.xlabel("Mode index")
+plt.ylabel("Cumulative power overlap")
+
+# Following eigenmode expansion, the propagating field inside the detector cross-section given a pure input mode can be expressed as (ignoring reflections):
+#
+# $$ E(x,y,z) = \sum_k \braket{E_k|E_0} e^{-i \beta_k z} \ket{E_k} $$
+#
+# where $\braket{\bm{x}|E_k} = E_k(x,y)$ are the cross-sectional mode profiles.
+
+# +
+z = 0
+
+
+def field_profile_at_z(z):
+    field = detector_modes[0].basis.zeros(dtype=complex)
+    for i, (k, overlap, mode) in enumerate(zip(ks, overlaps, detector_modes)):
+        field += overlap * np.exp(1j * k * z) * mode.E
+    return field
+
+
+# -
+
+for z in [500]:  # np.linspace(0, 50, 11):
+    field_basis = detector_modes[0].basis.zeros(dtype=complex)
+    field = field_profile_at_z(z)
+
+    detector_modes[0].show(np.real(field), colorbar=True)
diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index ae1782ce..fb309e85 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -241,6 +241,9 @@ def factor(w):
             )
         )
 
+    def calculate_pertubated_k(self, delta_epsilon):
+        return self.calculate_pertubated_neff(delta_epsilon) * 2 * np.pi / self.wavelength
+
     def calculate_pertubated_neff(self, delta_epsilon):
         return (
             self.n_eff
@@ -485,6 +488,7 @@ def compute_modes(
     num_modes=1,
     order=1,
     metallic_boundaries=False,
+    with_manual_boundaries=None,
     radius=np.inf,
     n_guess=None,
     solver="scipy",
@@ -544,6 +548,16 @@ def bform(e_t, e_z, v_t, v_z, w):
             ),
             solver=solver(k=num_modes, sigma=sigma),
         )
+    elif with_manual_boundaries is not None:
+        lams, xs = solve(
+            *condense(
+                -A,
+                -B,
+                D=basis.get_dofs(set(with_manual_boundaries.keys())),
+                x=basis.zeros(dtype=complex),
+            ),
+            solver=solver(k=num_modes, sigma=sigma),
+        )
     else:
         lams, xs = solve(
             -A,