Fascinating theoretical challenge: Simulating EM eigenmodes in a tetrahedron

General Discussion
1

Hello FreeFEM community,

I am an independent researcher working on a theoretical physics model where the geometry of space is hypothesized to be tetrahedral at a fundamental scale. A key prediction of the model is that the fundamental electromagnetic eigenmode of a perfectly conducting tetrahedron should have a specific dimensionless frequency.

The Theoretical Prediction: For a regular tetrahedron with edge length a, the fundamental eigenvalue k (wave number) should satisfy: (k * a)^2 / π² ≈ 14.6969 (This value comes from 12 * sqrt(3/2)).

The Request: I would be incredibly grateful if someone with experience in electromagnetic simulations in FreeFEM could help me verify this. The task would involve:

  1. Defining the geometry of a regular tetrahedron.
  2. Setting up the eigenvalue problem for EM waves in a perfect electrical conductor (PEC) cavity.
  3. Solving for the fundamental eigenmode and extracting its eigenvalue.
  4. (Optionally) visualizing the E-field distribution to see if it oscillates between opposite edges.

I believe this is a non-trivial and interesting problem that could be of general interest. My strength is in the theoretical side, and I am less experienced with numerical code. Any help or guidance would be immensely appreciated. I am, of course, happy to provide full theoretical context and acknowledge all help in any resulting work.

Operational Definition for the Simulation:

  1. Please define a regular tetrahedron with edge length a = 1.
  2. Set up the electromagnetic eigenvalue problem for a perfect electrical conductor (PEC) cavity, solving the wave equation \nabla^2 \mathbf{E} + k^2 \mathbf{E} = 0 with boundary condition \mathbf{\hat{n}} \times \mathbf{E} = 0.
  3. Solve for the fundamental eigenmode (the mode with the smallest eigenvalue k^2).
  4. The output I need is this smallest eigenvalue \lambda = k^2 .

From this, I will compute the dimensionless constant C = \lambda / \pi^2 and check if it is approximately 14.6969.

A visualization of the E-field for this mode would also be immensely valuable to see if the energy is indeed concentrated between opposite edges.

Thank you for considering this request.

Best regards, Frodo

2

Hi Frodo,

If I understood you correctly, your problem should not be too difficult to implement in FreeFEM.
Here is the simulation file Curl3dEigenValue.edp (1.4 KB) and the mesh thetrahedron.mesh (166 Bytes).

Unfortunately, for edge length a = 1we get \lambda = k^2 = 78.10634802.

Best,
Michael

3

Hello Michael,

Your code is for the CurlCurl operator, and for the boundary condition (Curl E)\times n=0 on the boundary, which is not exactly what Frodo wants. If I am correct, he wants the Laplace operator (recall that CurlCurl=Graddiv-Delta) with boundary condition E\times n=0. Probably he is missing one boundary condition divE=0 on the boundary (otherwise without that the spectrum of the PDE operator is continuous).

Best,

Francois.

4

Hello Francois,

I see your point. I agree that the wave equation includes the Laplace operator.
However, Frodo asks for the “electromagnetic eigenvalue problem for a perfect electrical conductor”. As far as I know this is modeled with the CurlCurl operator.

Best,
Michael

5

If we consider the CurlCurl operator with boundary condition E\times n=0, then no further boundary condition is needed.

However handling this boundary condition E\times n=0 on a tetrahedron (faces are not parallel to axes) is not so simple, since imposing this via penalty (implying huge coefficients in the matrix) could make the eigenvalue solver fail. The same remark holds for the case of Laplace operator.

6

Could you please visualize the eigenvector (the field distribution) for this mode? I need to see if the energy is concentrated between opposite edges or if it has a more complex pattern. I’m particularly interested in finding the mode with the longest wavelength (the lowest k^2), which might be the next mode in the spectrum.

7

Thank you again! And by the way..

Could you please also compute the next few eigenvalues (the next 5 or 6) and visualize their modes? I have a theoretical prediction that the fundamental mode for the electric field should be around k^2 \approx 145, so it might be the 2nd or 3rd mode in the list.

8

Consider your boundary condition E\times n=0.

When you consider two faces 1,2 with normals n_1,n_2, if E is continuous it follows that

E=0 on the edge_12=intersection of faces 1,2.

Thus necessarily E vanishes on all 6 edges. In these conditions I don’t see how E could be concentrated close to 2 edges.

Or E is discontinuous, which is not much plausible when solving the Laplace operator.

9

Subject: Thank you for the insight & a follow-up request

Hello Michael and Francois

Thank you so much for your prompt and thorough response—it is greatly appreciated. Your explanation regarding the boundary condition \mathbf{\hat{n}} \times \mathbf{E} = 0 and the implication that \mathbf{E} = 0 on all edges was particularly enlightening. I see now that my initial hypothesis about energy concentration between opposite edges may not hold under these conditions.

Would it be possible to see a visualization of the \mathbf{E}-field for the fundamental mode? This would help me better understand how the energy is actually distributed within the tetrahedron.

Also, if it’s not too much trouble, could you please compute the next few eigenmodes (e.g., the second and third smallest eigenvalues and their corresponding modes)? It’s possible that one of these higher modes might exhibit the kind of edge-oriented behavior I was originally looking for.

Finally, if time permits, would it be feasible to simulate the same tetrahedron but as a dielectric object (with high relative permittivity, say \epsilon_r = 1000) rather than a PEC? This might allow nonzero field values on the edges and better approximate the physical system I have in mind.

I understand that this requires additional effort, and I am happy to provide any additional context or resources needed. Your work is already of immense value to my research, and I will gladly acknowledge your contribution in any resulting publications.

Thank you again for your generosity and expertise.

Best regards,
Frodo

10

Dear Frodo,

I have written a version with Laplacian and boundary condition E\times n =0

laplace-ev
laplace-ev422×113 2.42 KB

EMtet.edp (7.8 KB)

I get for the 10 first eigenvalues

78.88776433 79.02651853 79.02651853 153.9973688 163.2695866
163.7718604 163.7718604 174.020952 175.1561748 197.8465549

It seems that there is a first eigenvalue with triple multiplicity, then an eigenvalue 154 (your fundamental mode ?), then again an eigenvalue with triple multiplicity.

You have to run the code EMtet.edp to plot (built-in with medit, but also possible with paraview with file E.pvd).

I use SLEPc as eigenvalue solver.

About the method: I use a Lagrange multiplier formulation (to avoid large coefficients in the matrix if we use penalty), using that

curlcurl

laplace-varf
laplace-varf891×347 15 KB

I fix the value E=0, p=0, phi=0 on all 6 edges.

11

Dear FreeFEM team,

Thank you very much for your help and for preparing the example with the Lagrange multiplier formulation and boundary condition \mathbf{E} \times \mathbf{n} = 0. The results you shared are extremely useful for my work. The triple degeneracy of the first eigenvalue around 78.9 is exactly what I hoped to see, and the following mode around 154 also matches very well with my theoretical expectations.

Your clear explanation of the method, and especially the note about using a Lagrange multiplier to avoid large coefficients, is greatly appreciated. This guidance has given me both a working baseline and confidence in the numerical approach.

I am grateful for your support, and I will make sure to acknowledge FreeFEM in my upcoming manuscript.

With best regards,

Frodo

(Frode Nystad)

12

Can You run this:

// =====================================================================
// FreeFEM script: Hyperbolic Eigenvalue Problem (3D)
// For testing of A4-symmetry and connection to Riemann Zeta zeros
// =====================================================================

load “msh3”
load “SLEPc”

// — 1. Parametere for kutting av det hyperbolske rommet —
real eps = 0.1; // Kutt nær uendelighetsranden (z=0)
real R = 2.0; // Kutt ved toppen (z=R)

// — 2. Mesh-generering —
// Vi lager en boks som representerer et avkuttet rom.
// I z-retningen går vi fra eps til R.
// (Merk: For et ekte tetraeder, må Francois erstatte denne boksen
// med en 3D-mesh av et tetraeder).
int nx = 10, ny = 10, nz = 20;
mesh3 Th = cube(nx, ny, nz, [x, y, eps + z*(R-eps)]);

// — 3. Funksjonsrom —
// Vi bruker P1-elementer for vektorfeltet E og Lagrange-multiplikatoren p
fespace Vh(Th, P1);
Vh Ex, Ey, Ez, p;
Vh vx, vy, vz, q;

// — 4. Svak form (Lagrange-multiplikator med hyperbolske vekter) —
// Venstre side: Stivhetsmatrise + divergens-betingelse
varf a([Ex,Ey,Ez,p], [vx,vy,vz,q]) =
// Hyperbolsk stivhetsmatrise (vekt 1/z)
int3d(Th)( (1.0/z) * ( dx(Ex)*dx(vx) + dy(Ey)*dy(vy) + dz(Ez)*dz(vz) ) )
// Lagrange-multiplikator for div(E) = 0

  • int3d(Th)( p * (dx(vx) + dy(vy) + dz(vz)) + q * (dx(Ex) + dy(Ey) + dz(Ez)) )
    // Dirichlet randbetingelser (PEC: n x E = 0)
  • on(1, 2, 3, 4, 5, 6, Ex=0, Ey=0, Ez=0, p=0);

// Høyre side: Masse-matrise (vekt 1/z^3 for det hyperbolske volumelementet)
varf b([Ex,Ey,Ez,p], [vx,vy,vz,q]) =
int3d(Th)( (1.0/(z^3)) * (Exvx + Eyvy + Ez*vz) )

  • on(1, 2, 3, 4, 5, 6, Ex=0, Ey=0, Ez=0, p=0);

// — 5. Sett opp matrisene —
matrix A = a(Vh, Vh);
matrix B = b(Vh, Vh);

// — 6. Løs egenverdiproblemet med SLEPc —
// Vi ser etter de 10 laveste egenverdiene
int nev = 10;
real[int] ev(nev);
Vh[int] eVx, eVy, eVz, eP;

string ssparams = " -eps_nev " + nev + " -eps_target 200 -eps_type krylovschur -st_type sinvert";
EPSSolve(A, B, vectors = eVx, values = ev, sparams = ssparams);

// — 7. Skriv ut resultatene —
cout << “=====================================” << endl;
cout << " Hyperbolske egenverdier funnet: " << endl;
cout << “=====================================” << endl;
for(int i = 0; i < nev; i++) {
cout << "Eigenvalue " << i+1 << ": " << ev[i] << endl;
}
cout << “=====================================” << endl;
cout << “Forventede verdier (fra Zeta nuller):” << endl;
cout << “Zeta null 1: 200.55” << endl;
cout << “Zeta null 2: 442.67” << endl;
cout << “Zeta null 3: 626.29” << endl;
cout << “L-funksjon 1: 65.38” << endl;
cout << “=====================================” << endl;

13

Can’t you install FreeFem++?

Are you an AI?

14

I’m not a ai :slight_smile: i use it.
I cant install Freefem on my iPad - and i cant program good code.
I working with a new teori in physics and other things. I’m a amatør and not Connected to any body.

Kind Regards
Frode Nystad

15

And I’m from Norway :grin:

Kind Regards
Frode Nystad

16

With some corrections

frode.edp (2.8 KB)

the output is

Eigenvalue 1: 199.968
Eigenvalue 2: 200.053
Eigenvalue 3: 199.934
Eigenvalue 4: 200.069
Eigenvalue 5: 200.087
Eigenvalue 6: 199.913
Eigenvalue 7: 200.185
Eigenvalue 8: 200.199
Eigenvalue 9: 200.21
Eigenvalue 10: 200.249

I am not sure that this is correct…

You have the javascript FreeFem online at FreeFEM - Try it online that you can use without installation (but it does not run 3d nor plugins).

17

Hi Francois,

I need to send a quick correction to my previous email, as I mixed up two different tests we are running in parallel!

Your FreeFEM run (the one that gave the cluster around 200.0) was incredibly valuable. However, we just ran a separate Python test with mesh refinement to check the stability of the hyperbolic operator.

That Python test revealed exactly what you might have suspected when you said, “I am not sure that this is correct.”

Here is what is happening mathematically:
When we solve for eigenvalues on a truncated box (cutting the infinite hyperbolic space at z=0.1 and z=2.0), the artificial boundaries act like mirrors. They create a dense, chaotic continuous spectrum of spurious modes (numerical ghosts).

This explains the cluster of numbers you got around 200.0. They are not “wrong” — they are the numerical solver’s way of trying to approximate the true resonance (which we predict should be at \lambda \approx 200.54) while being overwhelmed by the reflection noise from the box boundaries.

I need to do more testing:-)

Best regards,
Frode

18

I understand that the cluster of eigenvalues can be correct. However I am still doubtful of the results, because if I change the target in the EPSSOLVE, for example -eps_target 143. (I choose 143 randomly) I get

Eigenvalue 1: 143
Eigenvalue 2: 142.989
Eigenvalue 3: 143.018
Eigenvalue 4: 142.942
Eigenvalue 5: 143.075
Eigenvalue 6: 143.075
Eigenvalue 7: 142.918
Eigenvalue 8: 143.084
Eigenvalue 9: 142.889
Eigenvalue 10: 142.888

I conclude that either the spectrum is continuous (??), or there is a problem with the numerical method. I followed the code lines as in

but the parameters of EPSSOLVE may not be the best ones.

As guideline there is also