diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c
index e182318..621bdcf 100644
--- a/qpms/scatsystem.c
+++ b/qpms/scatsystem.c
@@ -2128,7 +2128,7 @@ beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
     complex double omega_centre, double omega_rr, double omega_ri, 
     size_t contour_npoints, 
     double rank_tol, size_t rank_sel_min, double res_tol) {
-  qpms_ss_ensure_nonperiodic(ss);
+  qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsys_periodic_find_eigenmodes()");
   size_t n; // matrix dimension
   if (qpms_iri_is_valid(ss->sym, iri)) {
     n = ss->saecv_sizes[iri];
@@ -2149,3 +2149,47 @@ beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
   return result;
 }
 
+struct qpms_scatsys_periodic_eval_Beyn_ImTW_param {
+  const qpms_scatsys_t *ss;
+  const double *k; ///< Wavevector in cartesian coordinates.
+};
+
+/// Wrapper for Beyn algorithm (periodic system)
+static int qpms_scatsys_periodic_eval_Beyn_ImTW(complex double *target,
+    size_t m, complex double omega, void *params){
+  const struct qpms_scatsys_periodic_eval_Beyn_ImTW_param *p = params;
+  qpms_scatsys_at_omega_t *ssw = qpms_scatsys_at_omega(p->ss, omega);
+  QPMS_ENSURE(ssw != NULL, "qpms_scatsys_at_omega() returned NULL");
+  qpms_scatsys_at_omega_k_t sswk = {
+    .ssw = ssw,
+    .k = {p->k[0], p->k[1], p->k[2]}
+  };
+  QPMS_ASSERT(m == p->ss->fecv_size);
+  QPMS_ENSURE(NULL != 
+      qpms_scatsyswk_build_modeproblem_matrix_full(target, &sswk),
+      "qpms_scatsyswk_build_modeproblem_matrix_full() returned NULL");
+  qpms_scatsys_at_omega_free(ssw);
+  return QPMS_SUCCESS;
+}
+
+beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
+    const qpms_scatsys_t * const ss, const double k[3], 
+    complex double omega_centre, double omega_rr, double omega_ri, 
+    size_t contour_npoints, 
+    double rank_tol, size_t rank_sel_min, double res_tol) {
+  qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsys_finite_find_eigenmodes()");
+  size_t n = ss->fecv_size; // matrix dimension
+
+  beyn_contour_t *contour = beyn_contour_ellipse(omega_centre,
+      omega_rr, omega_ri, contour_npoints);
+
+  struct qpms_scatsys_periodic_eval_Beyn_ImTW_param p = {ss, k};
+  beyn_result_t *result = beyn_solve(n, n, 
+      qpms_scatsys_periodic_eval_Beyn_ImTW, NULL, (void *) &p,
+      contour, rank_tol, rank_sel_min, res_tol);
+  QPMS_ENSURE(result != NULL, "beyn_solve() returned NULL");
+
+  free(contour);
+  return result;
+}
+
diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h
index fc30bea..cf9bfab 100644
--- a/qpms/scatsystem.h
+++ b/qpms/scatsystem.h
@@ -519,6 +519,28 @@ qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(
 		const qpms_scatsys_at_omega_k_t *sswk
 		);
 
+/// Searches for periodic scattering system's eigenmodes using Beyn's algorithm.
+/**
+ * Currently, elliptical contour is used.
+ *
+ * TODO In the future, this will probably support irrep decomposition as well,
+ * but for the case of periodic / small systems,
+ * the bottleneck is the T-matrix and translation matrix evaluation
+ * rather than the linear algebra.
+ */
+struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
+		const qpms_scatsys_t *ss, 
+		/// Wavevector in cartesian coordinates (must lie in the lattice plane).
+		const double k[3],
+		complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
+		double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
+		double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
+		size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
+		double rank_tol, ///< (default: `1e-4`) TODO DOC.
+		size_t rank_sel_min, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
+		double res_tol ///< (default: `0.0`) TODO DOC.
+		);
+
 // ======================= Periodic system -only related stuff end =========================