The 2-d resonant cavity

Figure 83 shows a 2-d resonant cavity consisting of a hollow, rectangular, perfectly conducting channel of dimensions $L_x\times L_y$. Suppose that the walls of the channel are aligned along the $x$- and $y$-axes. We shall excite this cavity in a rather artificial manner by imposing a $z$-directed alternating current pattern of frequency $f$, which has the same spatial structure as the mode in which we are interested. Let us calculate the electric and magnetic field patterns excited within the cavity by such a current pattern.

Figure 83: A 2-d resonant cavity.

The electric and magnetic fields within the cavity can be written ${\bf E} = [0,0,$
$E_z(x,y,t)]$, and ${\bf B} = [B_x(x,y,t), B_y(x,y,t), 0]$, respectively. It follows from Maxwell's equations that

$$\frac{\partial H_x}{\partial t} + c\,\frac{\partial E_z}{\partial y} = 0,$$ (268)

$$\frac{\partial H_y}{\partial t} + c\,\frac{\partial E_z}{\partial x} = 0,$$ (269)

$$\frac{\partial E_z}{\partial t} + c\,\frac{\partial H_y}{\partial x} + c\,\frac{\partial H_x}{\partial y} = J_z,$$ (270)

where $c$ is the velocity of light, $H_x=c\,B_x$, $H_y=-c\,B_y$, and $J_z =-\mu_0\,c^2\,j_z$. Note that the above system of equations takes the form of three coupled advection equations with a source term. The boundary conditions are that the tangential electric field and the normal magnetic field must be zero at the conducting walls. It follows that

$$E_z=0$$ (271)

at all the walls (which are located at $x=0,L_x$ and $y=0,L_y$),

$$H_x = \frac{\partial H_y}{\partial x} = 0$$ (272)

at $x=0,L_x$, and

$$H_y = \frac{\partial H_x}{\partial y} = 0$$ (273)

at $y=0,L_y$. Finally, the normalized current pattern associated with the ($m$, $n$) mode takes the form

$$J_z(x,y,t) = J_0\,\sin(m\,\pi\,x/L_x)\,\sin(n\,\pi\,y/L_y)\,\sin(2\,\pi\,f\,t).$$ (274)

As usual, we discretize in time on the uniform grid $t_n=t_0+n\,\delta t$, for $n=0,1,2,\cdots$. Furthermore, in the $x$-direction, we discretize on the uniform grid $x_i = i\,\delta x$, for $i=0,I$, where $\delta x = L_x/I$. Finally, in the $y$-direction, we discretize on the uniform grid $y_j = j\,\delta y$, for $j=0,J$, where $\delta y = L_y/J$. Adopting a Crank-Nicholson temporal differencing scheme similar to that discussed in Sects. 7.4 and 7.6, Eqs. (268)-(270) yield

$$\frac{(H_x)_{i,j}^{n+1} -(H_x)_{i,j}^{n}}{\delta t} +\frac{c}{2}\... ...^{n+1}_{i,j} +\frac{c}{2}\left(\frac{\partial E_z}{\partial y}\right)^{n}_{i,j} = 0,$$ (275)

$$\frac{(H_y)_{i,j}^{n+1} -(H_y)_{i,j}^{n}}{\delta t} +\frac{c}{2}\... ...^{n+1}_{i,j} +\frac{c}{2}\left(\frac{\partial E_z}{\partial x}\right)^{n}_{i,j} = 0,$$ (276)

$$\frac{(E_z)_{i,j}^{n+1} -(E_z)_{i,j}^{n}}{\delta t} +\frac{c}{2}\... ...^{n+1}_{i,j} +\frac{c}{2}\left(\frac{\partial H_y}{\partial x}\right)^{n}_{i,j}$$

$$+\frac{c}{2}\left(\frac{\partial H_x}{\partial y}\right)^{n+1}_{i,j} +\frac{c}{2}\left(\frac{\partial H_x}{\partial y}\right)^{n}_{i,j} = (J_z)_{i,j}^n,$$ (277)

where $(H_x)_{i,j}^n\equiv H_x(x_i, y_j, t_n)$, etc.

Adopting a Fourier approach, we write

$$(E_z)_{i,j}^n = \sum_{i'=0,I}^{j'=0,J}\hat{E}_{i',j'}^n\,\sin(i\,i'\,\pi/I)\,\sin(j\,j'\,\pi/J),$$ (278)

$$(H_x)_{i,j}^n = \sum_{i'=0,I}^{j'=0,J} \hat{X}_{i',j'}^n\,\sin(i\,i'\,\pi/I)\,\cos(j\,j'\,\pi/J),$$ (279)

$$(H_y)_{i,j}^n = \sum_{i'=0,I}^{j'=0,J} \hat{Y}_{i',j'}^n\,\cos(i\,i'\,\pi/I)\,\sin(j\,j'\,\pi/J),$$ (280)

$$(J_z)_{i,j}^n = \sum_{i'=0,I}^{j'=0,J} \hat{J}_{i',j'}^n\,\sin(i\,i'\,\pi/I)\,\sin(j\,j'\,\pi/J)$$ (281)

which automatically satisfies the boundary conditions (271)-(273). Equations (275)-(277) yield

$$\hat{X}_{i,j}^{n+1}-\hat{X}_{i,j}^n + j\,D_y\,(\hat{E}_{i,j}^{n+1}+\hat{E}_{i,j}^{n}) = 0,$$ (282)

$$\hat{Y}_{i,j}^{n+1}-\hat{Y}_{i,j}^n + i\,D_x\,(\hat{E}_{i,j}^{n+1}+\hat{E}_{i,j}^{n}) = 0,$$ (283)

$$\hat{E}_{i,j}^{n+1}-\hat{E}_{i,j}^n - i\,D_x\,(\hat{Y}_{i,j}^{n+1}+\hat{Y}_{i,j}^{n}) - j\,D_j\,(\hat{X}_{i,j}^{n+1}+\hat{X}_{i,j}^{n}) = \delta t\,\hat{J}^{n}_{i,j},$$ (284)

for $i=0,I$ and $j=0,J$, where $D_x=\pi\,c\,\delta t/(2\,L_x)$ and $D_y=\pi\,c\,\delta t/(2\,L_y)$. It follows that

$$\hat{E}^{n+1}_{i,j} = \frac{(1-i^2\,D_x^2-j^2\,D_y^2)\,\hat{E}_{i,j}^n + 2\,j\,D_y\, \h... ...\,D_x\,\hat{Y}_{i,j}^n +\delta t \,\hat{J}^n_{i,j}} {1 + i^2\,D_x^2+j^2\,D_y^2}$$ (285)

$$\hat{X}_{i,j}^{n+1} = \hat{X}_{i,j}^n - j\,D_y\,(\hat{E}_{i,j}^{n+1}+\hat{E}_{i,j}^{n}),$$ (286)

$$\hat{Y}_{i,j}^{n+1} = \hat{Y}_{i,j}^n - i\,D_x\,(\hat{E}_{i,j}^{n+1}+\hat{E}_{i,j}^{n}).$$ (287)

The routine listed below solves the 2-d wave equation in a resonant cavity using the Crank-Nicholson scheme discussed above. The routine first Fourier transforms $H_x$, $H_y$, $E_z$, and $J_z$ in both the $x$- and $y$-directions, takes a time-step using Eqs. (285)-(287), and then reconstructs $H_x$, $H_y$, and $E_z$ via an double inverse Fourier transform.

// Wave2D.cpp

// Function to evolve 2-d wave equation:

//  d H_x / dt + c d E_z / dy = 0

//  d H_y / dt + c d E_z / dx = 0

//  d E_z / dt + c d H_y / dx + c d H_x / dy = J_z

// in region 0 < x < L_x and 0 < y < L_y

// Boundary conditions:

//  E_z(0, y) = E_z(L_x, y) = E_z(x, 0) = E_z(x, L_y) = 0

//  H_x(0, y) = H_x(L_x, y) = d H_y(0, y) / dx =  d H_y(L_x, y) / dx = 0

//  H_y(x, 0) = H_y(x, L_y) = d H_x(x, 0) / dy =  d H_x(x, L_y) / dy = 0

// Matrices Hx, Hy, Ez, Jz assumed to be of extent I+1, J+1.
// Now, (i,j)th elements of matrices correspond to

//  x_i = i * dx    i=0,I

//  y_j = j * dy    j=0,J

// Here, dx = L_x / I is grid spacing in x-direction,

// and dy = L_y / J is grid spacing in x-direction.

// Now, Dx = pi c dt / (2 L_x) and Dy = pi c dt / (2 L_y),

// where dt is time-step.

// Uses Crank-Nicholson scheme.

#include <blitz/array.h>

using namespace blitz;

void fft_forward_cos (Array<double,1> f, Array<double,1>& F);
void fft_backward_cos (Array<double,1> F, Array<double,1>& f);
void fft_forward_sin (Array<double,1> f, Array<double,1>& F);
void fft_backward_sin (Array<double,1> F, Array<double,1>& f);

void Wave2D (Array<double,2>& Hx, Array<double,2>& Hy, Array<double,2>& Ez,
             Array<double,2> Jz, double Dx, double Dy, double dt)
{
  // Find I and J. Declare local arrays
  int I = Hx.extent(0) - 1;
  int J = Hx.extent(1) - 1;
  Array<double,2> X(I+1, J+1), XX(I+1, J+1), XXX(I+1, J+1);
  Array<double,2> Y(I+1, J+1), YY(I+1, J+1), YYY(I+1, J+1); 
  Array<double,2> E(I+1, J+1), EE(I+1, J+1), EEE(I+1, J+1);
  Array<double,2> K(I+1, J+1), KK(I+1, J+1);

  // Fourier transform solution in x-direction
  for (int j = 0; j <= J; j++)
    {
      Array<double,1> In(I+1), Out(I+1);

      // Fourier transform Hx
      for (int i = 0; i <= I; i++) In(i) = Hx(i, j);
      fft_forward_sin (In, Out);
      for (int i = 0; i <= I; i++) X(i, j) = Out(i);

      // Fourier transform Hy
      for (int i = 0; i <= I; i++) In(i) = Hy(i, j);
      fft_forward_cos (In, Out);
      for (int i = 0; i <= I; i++) Y(i, j) = Out(i);

      // Fourier transform Ez
      for (int i = 0; i <= I; i++) In(i) = Ez(i, j);
      fft_forward_sin (In, Out);
      for (int i = 0; i <= I; i++) E(i, j) = Out(i);

      // Fourier transform Jz
      for (int i = 0; i <= I; i++) In(i) = dt * Jz(i, j);
      fft_forward_sin (In, Out);
      for (int i = 0; i <= I; i++) K(i, j) = Out(i);
    }

  // Fourier transform solution in y-direction
  for (int i = 0; i <= I; i++)
    {
      Array<double,1> In(J+1), Out(J+1);

      // Fourier transform Hx
      for (int j = 0; j <= J; j++) In(j) = X(i, j);
      fft_forward_cos (In, Out);
      for (int j = 0; j <= J; j++) XX(i, j) = Out(j);

      // Fourier transform Hy
      for (int j = 0; j <= J; j++) In(j) = Y(i, j);
      fft_forward_sin (In, Out);
      for (int j = 0; j <= J; j++) YY(i, j) = Out(j);

      // Fourier transform Ez
      for (int j = 0; j <= J; j++) In(j) = E(i, j);
      fft_forward_sin (In, Out);
      for (int j = 0; j <= J; j++) EE(i, j) = Out(j);

      // Fourier transform Jz
      for (int j = 0; j <= J; j++) In(j) = K(i, j);
      fft_forward_sin (In, Out);
      for (int j = 0; j <= J; j++) KK(i, j) = Out(j);
    }
  
  // Evolve XX, YY, and EE
  for (int i = 0; i <= I; i++)
    for (int j = 0; j <= J; j++)
      {
        double x = double (i) * Dx;
        double y = double (j) * Dy;
        double fp = 1. + x*x + y*y;
        double fm = 1. - x*x - y*y;

        EEE(i, j) = fm * EE(i, j) + 2. * y * XX(i, j) +
          2. * x * YY(i,j) + KK(i, j);
        EEE(i, j) /= fp;
        XXX(i, j) = XX(i, j) - y * (EEE(i, j) + EE(i, j));
        YYY(i, j) = YY(i, j) - x * (EEE(i, j) + EE(i, j));  
      }

  // Reconstruct solution via inverse Fourier transform in y-direction
  for (int i = 0; i <= I; i++)
    {
      Array<double,1> In(J+1), Out(J+1);

      // Reconstruct Hx
      for (int j = 0; j <= J; j++) In(j) = XXX(i, j);
      fft_backward_cos (In, Out);
      for (int j = 0; j <= J; j++) X(i, j) = Out(j);

      // Reconstruct Hy
      for (int j = 0; j <= J; j++) In(j) = YYY(i, j);
      fft_backward_sin (In, Out);
      for (int j = 0; j <= J; j++) Y(i, j) = Out(j);

      // Reconstruct Ez
      for (int j = 0; j <= J; j++) In(j) = EEE(i, j);
      fft_backward_sin (In, Out);
      for (int j = 0; j <= J; j++) E(i, j) = Out(j);
    } 

  // Reconstruct solution via inverse Fourier transform in x-direction
  for (int j = 0; j <= J; j++)
    {
      Array<double,1> In(I+1), Out(I+1);

      // Reconstruct Hx
      for (int i = 0; i <= I; i++) In(i) = X(i, j);
      fft_backward_sin (In, Out);
      for (int i = 0; i <= I; i++) Hx(i, j) = Out(i);

      // Reconstruct Hy
      for (int i = 0; i <= I; i++) In(i) = Y(i, j);
      fft_backward_cos (In, Out);
      for (int i = 0; i <= I; i++) Hy(i, j) = Out(i);

      // Reconstruct Ez
      for (int i = 0; i <= I; i++) In(i) = E(i, j);
      fft_backward_sin (In, Out);
      for (int i = 0; i <= I; i++) Ez(i, j) = Out(i);
    }
}

The numerical calculations discussed below were performed using the above routine. The electromagnetic fields $H_x$, $H_y$, and $E_z$ were all initialized to zero everywhere at $t=0$. Figure 84 shows the maximum amplitude of $E_z$ versus the frequency, $f$, for an $m=1/n=1$ driving current distribution. It can be seen that there is a clear resonance at $f\simeq 0.7$.

Figure 84: Electromagnetic waves in a 2-d resonant cavity. The maximum amplitude of $E_z(x=0.5,y=0.5)$ between $t=0$ and $t=20$ versus the driving frequency, $f$. Numerical calculation performed using $m=1$, $n=1$, $L_x=1$, $L_y=1$, $c=1$, $I=J=32$, and $\delta t = 10^{-2}$.

Figures 85 and 86 illustrate the typical time variation of $E_z$, $H_x$, and $H_y$ for a non-resonant and a resonant case, respectively. For the non-resonant case, the traces take the form of interference patterns between the directly driven response, which oscillates at the driving frequency $f$, and the transient response, which oscillates at the natural frequency $f_0$ of the cavity. Note that the transients never decay, since there is no dissipation in the present problem. Incidentally, it is easily demonstrated that

$$f_0 = \frac{c}{2}\sqrt{\frac{1}{(n\,L_x)^2}+\frac{1}{(m\,L_y)^2}}.$$ (288)

Hence, it follows that $f_0=1/\sqrt{2} = 0.7071$ for $n=m=L_x=L_y=c=1$, which corresponds very well to the resonant frequency found in Fig. 84. For the resonant case, the traces take the form of waves of ever increasing amplitude which oscillate at the natural frequency $f_0$.

Figure 85: Electromagnetic waves in a 2-d resonant cavity. Time traces of $E_z(x=0.5,y=0.5)$ (solid curve), $H_x(x=0.5,y=0.0)$ (dashed curve), and $H_y(x=0.0,y=0.5)$ (dotted curve--obscured by dashed curve). Numerical calculation performed using $m=1$, $n=1$, $L_x=1$, $L_y=1$, $c=1$, $I=J=32$, $\delta t = 10^{-2}$, and $f=0.6$.

Figure 86: Electromagnetic waves in a 2-d resonant cavity. Time traces of $E_z(x=0.5,y=0.5)$ (solid curve), $H_x(x=0.5,y=0.0)$ (dashed curve), and $H_y(x=0.0,y=0.5)$ (dotted curve--obscured by dashed curve). Numerical calculation performed using $m=1$, $n=1$, $L_x=1$, $L_y=1$, $c=1$, $I=J=32$, $\delta t = 10^{-2}$, and $f=0.7071$.

Finally, Figs. 86 and 87 illustrate the spatial variation of the electromagnetic fields driven within the cavity when $m=1$ and $n=1$.

Figure 87: Electromagnetic waves in a 2-d resonant cavity. Spatial variation of $E_z$ (solid curve) and $H_y$ (dashed curve) in the $x$-direction at $t=20$ and $y=0.5$. Numerical calculation performed using $m=1$, $n=1$, $L_x=1$, $L_y=1$, $c=1$, $I=J=32$, $\delta t = 10^{-2}$, and $f=0.7071$.

Figure 88: Electromagnetic waves in a 2-d resonant cavity. Spatial variation of $E_z$ (solid curve) and $H_x$ (dashed curve) in the $y$-direction at $t=20$ and $x=0.5$. Numerical calculation performed using $m=1$, $n=1$, $L_x=1$, $L_y=1$, $c=1$, $I=J=32$, $\delta t = 10^{-2}$, and $f=0.7071$.