Maybe you can elaborate again why you need complex numbers in time domain??

the problem is that wheres for L-periodic structures in x-direction we would encounter a phase shift of the fields in frequency domain: u(x=0) = u(x=L) * exp(I * kx * L), for the pulses used in FDTD we would have a time shift in time-domain, which would sometimes involve fields at later times to do the updates and which would be furthermore frequency dependent.

The PBC I am trying to implement solves the problem of frequency-dependence by injecting a plane wave pulse into the simulation volume, which posseses a frequency-independent parallel wave-vector.

For a periodicity of Lx in x-direction this means in time domain: u(x=0, y, t+1) = A*sin(kx*Lx) * u(x=L, y, t) + B *cos(kx*Lx) * u(x=L,y, t). But the "weights" A and B of the sine and the cosine are not known apriori.

While writing I realize, that maybe it would be doable what you suggested earlier... to keep the real fields known to the main engine and multiply them by cos(kx * L) and save new values inside the extension, which would involve a multiplication by sin(kx * L).

See for example the book "Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology" page 84... the key-word is "Bloch-Waves" :

https://books.google.de/books?id=ynRM33 ... td&f=false]

At the moment I think one may initialize the real part as usual and once we reach one side of the computational grid, we obtain an imaginary part which is cos(kx * Lx) the real-part and and from there on each further timestep includes "mixing" between imaginary and real parts by means of multiplications by sin(kx * Lx) / cos(kx * Lx)

They are never set as all voltages (and currents) are initialized as zero and the engine just never is able to change the values and additonally the operator for the boundaries is zero too.

Let me elaborate on this to see if I understood correctly:

Since for example Ux(t+1,nx,ny,nz) only depends on Ux(t,nx,ny,nz) and Ix(t, nx, ny-1, nz-1) it needs: Ux(t,nx,ny,nz) = Ix(t, nx, ny-1, nz-1) = 0 in order for Ux(t,nx,ny,nz) = 0 for all timesteps.

What you say is, that for example Ux(., 0,ny,nz), U(.,nx,0,nz) and U(.,nx,ny,0) as well as those components at the other faces of the simulation volume are just never updated?

If I set the BC in -x to PMC, then all currents with index nx=0 are zero and a time-update (regardless of the boundary condition) only updates "inner" voltages,U(.,nx>0<Nx,ny>0<Ny,nz<Nz),

where the indices Nx, Ny, Nz indicate the greatest indices in the respective direction?