Floquet-Maxwell Simulations

• Floquet-Maxwell simulations are computational methods that model electromagnetic systems with periodic temporal or spatiotemporal modulations.
• They integrate Maxwell’s equations with Bloch or GSTC dynamics to analyze nonequilibrium quantum states, harmonic generation, and nonreciprocity.
• These methods enable efficient study of driven systems, with applications from quantum light–matter interactions to advanced photonic device design.

Floquet-Maxwell simulations constitute a class of computational and analytical methods for modeling electromagnetic systems subjected to periodic temporal and/or spatiotemporal modulations. Central to these approaches is the direct integration or rigorous Floquet analysis of Maxwell’s equations, often in combination with models for quantum or nonlinear matter, in settings where periodic driving leads to the emergence of nontrivial steady states and complex harmonic spectra. These frameworks are widely applied for investigating nonequilibrium phenomena in open quantum systems (as realized by Maxwell–Bloch models) and for analyzing harmonic generation, nonreciprocity, and wave-mixing in space–time modulated dispersive metasurfaces.

1. Governing Equations and Physical Models

The core of Floquet-Maxwell simulations is the coupling of electromagnetic field evolution (Maxwell’s equations) to periodically driven quantum or classical media, with dissipation and dispersion treated according to system requirements.

Maxwell–Bloch Model for Driven Quantum Systems:

In a non-magnetic dielectric without free charges or currents, Maxwell’s equations (SI units) simplify to

• Faraday’s law: ∇×E(r,t) = –∂ₜB(r,t)
• Ampère–Maxwell law: ∇×H(r,t) = ∂ₜD(r,t) with D(r,t)=ε₀E(r,t)+P(r,t), B(r,t)=μ₀H(r,t).

Elimination of H and B yields the driven wave equation:

$$\partial_t^2 E(r,t) - c^2 \nabla^2 E(r,t) = \mu_0 \partial_t^2 P(r,t)$$

For spatially uniform (“single-mode”) approximations, this simplifies to

$$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$

where $\gamma_\text{EM}$ accounts for cavity or propagation losses, $\Omega_\text{EM}$ the mode resonance.

The matter system (e.g., two-level quantum system) is governed by the Bloch equations with polarization coupling. Denoting the ground and excited state as $|g\rangle, |e\rangle$ with level splitting $\Delta$, the Hamiltonian in the electric dipole approximation is

$$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$

with $\mu$ the transition dipole. The dynamics of the density matrix are given by

$$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$

where $\Omega_R(t) = \mu E(t) / \hbar$, $$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$0 is the longitudinal relaxation time, and $$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$1 the decoherence (transverse relaxation) time (Sato et al., 2019).

Maxwell–Floquet Model for Space-Time Modulated Metasurfaces:

For zero-thickness Huygens’ metasurfaces with Lorentz-dispersive susceptibilities periodically modulated in space and time, the fields satisfy Maxwell’s curl equations with delta-function sheet polarization and magnetization at $$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$2: $$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$3

$$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$4

with Generalized Sheet Transition Conditions (GSTCs) at the interface, and polarization densities ($$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$5) following Lorentzian dynamics with space-time modulated parameters. All system coefficients (resonant frequencies, plasma frequencies, damping) are expanded in double Fourier series to enable systematic Floquet analysis (Tiukuvaara et al., 2020).

2. Floquet Ansatz and Harmonic Expansion

Periodic driving (either pure-time or space–time inhomogeneity) enables the application of Floquet’s theorem. For strictly periodic systems, field or matter solutions can be expanded in harmonics of the drive:

• Closed Quantum Systems: Solutions take the form $$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$6 with $$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$7 and $$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$8 the quasienergies. Observables and operator-valued quantities are expanded as Fourier series in $$\partial_t^2 E(t) + \gamma_\text{EM} \partial_t E(t) + \Omega_\text{EM}^2 E(t) = \mu_0 \partial_t^2 P(t)$$9.
• Metasurfaces: For systems periodic in both $\gamma_\text{EM}$0 and $\gamma_\text{EM}$1, the reflected and transmitted fields, as well as polarization densities, are expanded in double Floquet harmonic bases: $\gamma_\text{EM}$2 with

$\gamma_\text{EM}$3

and $\gamma_\text{EM}$4 (Tiukuvaara et al., 2020).

In open quantum systems with dissipation, the nonequilibrium steady state $\gamma_\text{EM}$5 inherits the periodicity of the drive, and a Floquet fidelity can be defined based on the overlap of $\gamma_\text{EM}$6 orbitals with closed-system Floquet modes (Sato et al., 2019).

3. Dissipation, Decoherence, and Lindblad Formulation

Dissipation is incorporated via phenomenological relaxation times or Lindblad master equations. In Maxwell–Bloch formulations, relaxation is addressed with: $\gamma_\text{EM}$7 with

$\gamma_\text{EM}$8

$\gamma_\text{EM}$9

This is equivalent to the Lindblad form: $\Omega_\text{EM}$0 with $\Omega_\text{EM}$1, $\Omega_\text{EM}$2. Population relaxes via $\Omega_\text{EM}$3, while coherence is lost via $\Omega_\text{EM}$4, which directly impedes Floquet state formation. The explicit effect: Rabi-splitting remains visible if $\Omega_\text{EM}$5 ($\Omega_\text{EM}$6: Rabi period). Strong drive ($\Omega_\text{EM}$7 large, $\Omega_\text{EM}$8) can restore Floquet regimes even for $\Omega_\text{EM}$9 (Sato et al., 2019).

4. Computational Methods and Algorithms

Maxwell–Bloch Real-Time Integration:

• Initialize system parameters: $|g\rangle, |e\rangle$0, $|g\rangle, |e\rangle$1, drive frequency $|g\rangle, |e\rangle$2, amplitude $|g\rangle, |e\rangle$3, relaxation times $|g\rangle, |e\rangle$4, $|g\rangle, |e\rangle$5.
• Setup the Maxwell mode or use a prescribed drive $|g\rangle, |e\rangle$6.
• Initialize density matrix $|g\rangle, |e\rangle$7.
• Propagate $|g\rangle, |e\rangle$8 and $|g\rangle, |e\rangle$9 with a suitable ODE solver (e.g., RK4), using a time step $\Delta$0. For full coupling, update $\Delta$1 via finite-difference Maxwell each step.
• Continue until $\Delta$2 over $\Delta$3–$\Delta$4 driving cycles, indicating the steady state.
• Compute observables: Floquet-fidelity ($\Delta$5), spectral lines via probe coupling and Fourier analysis, quasienergy via matrix diagonalization in the closed-system limit.

Floquet–GSTC Matrix Approach for Metasurfaces:

• Expand all fields, polarizations, and Lorentz parameters as double Fourier series in space and time.
• Substitute expansions into Lorentz and GSTC equations, yielding a coupled linear system after equating harmonics.
• Truncate to $\Delta$6 harmonics for numerics; assemble unknown amplitudes $\Delta$7, $\Delta$8, $\Delta$9, $$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$0 into vectors.
• Solve the resulting linear matrix equation for the Floquet coefficients, e.g.,

$$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$1

with $$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$2 constructed from system parameters and Fourier coefficients (see Appendix A of (Tiukuvaara et al., 2020)).

• For arbitrary incident fields (e.g., Gaussian beams), decompose into plane waves via Fourier analysis in $$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$3 and superpose the solutions for each component.
• Adjust truncation orders $$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$4 for convergence.

Typical Parameter Choices (Maxwell–Bloch):

• $$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$5 (units), $$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$6; $$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$7 range $$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$8–$$H(t) = \frac{\Delta}{2} \sigma_z - \mu E(t) \sigma_x$$9
• $\mu$0; also $\mu$1
• Time-step $\mu$2–$\mu$3; check convergence as above.
• Calculation of $\mu$4 and mapping of quasienergy spectra as in (Sato et al., 2019).

5. Typical Results, Rules of Thumb, and Physical Insights

Quantum Floquet–Maxwell Regimes (Sato et al., 2019):

• Rabi Splitting: Double-peak quasienergy splitting of size $\mu$5 remains robust if $\mu$6.
• Drive Strength vs Decoherence: For $\mu$7, Floquet features emerge when $\mu$8 is large enough, i.e., $\mu$9 (typically $$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$0–$$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$1).
• Resonant Drive ($$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$2): Weak fields cause heating and dissipative destruction of Floquet fidelity ($$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$3). Strong driving reinstates $$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$4 via dynamical stabilization.
• Off-Resonant Regime ($$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$5): Heating is suppressed; dissipation stays inactive and $$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$6 at low fields. Nonlinear excitation at higher fields triggers dissipation, temporarily reducing $$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$7 before strong dressing again stabilizes Floquet structure at large drive.
• Heating Suppression: The crossing of net heating rate ($$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$8) demarcates the regime of coherent Floquet state formation.

Metasurface Floquet–Maxwell Phenomena (Tiukuvaara et al., 2020):

• Pure Spatial Modulation ($$\begin{aligned} \partial_t \rho_{ee} &= -i[\Omega_R(t) (\rho_{ge} - \rho_{eg})] - (\rho_{ee} - \rho_{eq})/T_1 \ \partial_t \rho_{ge} &= -i \Delta \rho_{ge} - i \Omega_R(t)(\rho_{ee} - \rho_{gg}) - \rho_{ge}/T_2 \end{aligned}$$9): Harmonic diffraction patterns (cosine, sawtooth profiles) accurately reproduced by Floquet-GSTC method; strong spatial asymmetry via non-cosine profiles.
• Pure Temporal Modulation ($\Omega_R(t) = \mu E(t) / \hbar$0): Generation of temporal sidebands; strong modulation yields negative-frequency components recoverable via FFT.
• Space–Time Modulation: Standing-wave modulations preserve reciprocity (Onsager–Casimir holds), while traveling-wave modulations induce nonreciprocity, with up- and down-conversion between ports at different frequencies/angles. Output for Gaussian beams or complex excitations obtained by superposing solved plane-wave components.

Scenario	Parameter Regime / Key Effect	Reference
Rabi Splitting Survival	$\Omega_R(t) = \mu E(t) / \hbar$1	(Sato et al., 2019)
Floquet Recovery (Strong Drive)	$\Omega_R(t) = \mu E(t) / \hbar$2	(Sato et al., 2019)
Reciprocal Metasurface	Standing wave, $\Omega_R(t) = \mu E(t) / \hbar$3	(Tiukuvaara et al., 2020)
Nonreciprocal Metasurface	Traveling wave, $\Omega_R(t) = \mu E(t) / \hbar$4 breaks Onsager	(Tiukuvaara et al., 2020)

A plausible implication is that periodic suppression of heating, whether by detuning or field engineering, is as significant as minimizing material damage for preserving Floquet coherence in both quantum and classical driven systems.

6. Extensions and Applicability

Floquet-Maxwell methods are extensible to broader classes of systems:

• Multilevel Quantum Systems: Generalize by expanding $\Omega_R(t) = \mu E(t) / \hbar$5 to $\Omega_R(t) = \mu E(t) / \hbar$6 and using the corresponding operators, with the identical relaxation-time or more sophisticated microscopic Lindblad kernels (Sato et al., 2019).
• Multiresonator Surfaces: Floquet–GSTC schemes extend naturally to sums of Lorentz poles, enabling accurate modeling of advanced dispersive and nonlinear metasurfaces (Tiukuvaara et al., 2020).
• Arbitrary Periodic Modulation: Any analytical or synthesized periodic profile in space and/or time is admitted, provided the relevant harmonics are retained for convergence.
• Oblique or Structured Excitations: Fourier decomposition and superposition enable handling of Gaussian beams or complex incident waveforms in metasurface analysis.

The high accuracy and computational efficiency of Floquet-expansion-based approaches compared to brute-force time-domain solvers make them especially suitable for steady-state, high-dimensional, and highly resolved harmonic response studies. Versatile applications include quantum light–matter dynamics, harmonic generation, nonreciprocity, and custom beam engineering in advanced photonic platforms (Sato et al., 2019, Tiukuvaara et al., 2020).