Density-functional tight binding meets Maxwell: unraveling the mysteries of (strong) light–matter coupling efficiently

1 Introduction

Even a decade after the first experimental evidence showing that strong coupling of molecules in optical cavities can alter chemical properties, such as reaction rates [1], [2], [3], [4], [5], [6], [7], [8], energy [9], [10], [11] and charge [12], [13] transport, or self-assembly properties [14], [15], polaritonic chemistry remains a computationally and theoretically extremely challenging problem. While strong coupling arising from local plasmonic near-fields in single or few molecules has been measured and explained [16], [17], modified chemical reactivities in molecular ensembles in planar cavities are challenging to model. Indeed, we still lack a detailed mechanistic understanding and microscopic (local) insights on an ab initio level of theory into chemical ensembles under collective strong coupling [18], [19]. Collective insights beyond two level systems remain rare exceptions [20], [21], [22], [23], [24], [25], [26], [27], and many aspects remain open, e.g., the relevance of collective electron correlation effects that possess spin-glass features [25]. The reason for the computational complexity is that the collective strong coupling with an optical cavity introduces substantial feedback effects between different spatiotemporal scales, which are not present in traditional theoretical chemistry (see Figure 1). Consequently, established approximate methods from quantum chemistry, which have proven their reliability over decades, may no longer capture all the relevant physics. For this reason, different novel ab initio quantum electrodynamics (qed) methods have been developed such as quantum-electrodynamical density-functional theory (qedft) [28], [29], [30], QED-coupled cluster [31], [32], cavity Hartree–Fock, and many more [33], [34]. However, in practice, those accurate simulation methods are typically restricted to few-molecules strong-coupling simulations in the long-wavelength approximation. Consequently, the spatiotemporal modulations of the transverse electric fields are usually not resolved, and the physical complexity of realistic cavity mirrors is discarded. In particular, it is unclear how the interplay of these details affects the interplay between collective and local strong-coupling effects.

The Maxwell + time-dependent density-functional theory (tddft) methods that have been developed, both in local dipole approximation [30] and in full minimal coupling [35], are in principle able to consider all these interactions in the classical limit for the photons. However, the microscopic treatment of Maxwell’s equations in this formalism renders the application to molecular ensembles impractical due to their multiscale nature. In the past years, several quantum/classical methods have been developed to account for electromagnetic environments and strong coupling in nanoplasmonic arrangements [36], [37], [38], or STM junctions [39]; these methods, however, lack retardation effects that are needed for strong coupling with transverse modes. Similarly, hybrid quantum/classical treatments of near-field effects on molecular properties (from surface-enhanced Raman spectroscopy to strong coupling) have been proposed [40], [41], [42], which do not account for the radiative dynamics into the far field.

To overcome these spatiotemporal challenges in polaritonic chemistry, we have developed a computational framework that self-consistently couples Maxwell’s equations with dftb [43], a highly efficient and scalable electronic-structure method. This numerical framework allows to investigate nonperturbatively strong light–matter interactions of large polaritonic ensembles in realistic cavities and provides direct access to global and local properties of this multiscale problem. This method is a convenient fast approximation to the full solution of the Maxwell + tddft equations. Full minimal coupling Maxwell + tddft can be seen as the first principle reference method, which can account for beyond dipolar and magnetic field effects, e.g., giving access to Cherenkov emission [35]. In contrast, Maxwell + dftb is computationally more efficient, but simpler semi-empirical approximation, where local charges solely couple to the electric field vector [43].

On the ensemble scale, the method can directly compute two-dimensional optical spectra. This experimental technique has emerged as an useful tool for probing the complex dynamics of polaritonic systems, where the coherent hybridization of molecular excitations with confined electromagnetic modes creates new pathways for energy transfer, chemical reactivity, and nonlinear optical responses. Unlike conventional spectroscopic techniques, two-dimensional spectroscopy can directly map the coupling networks between polaritons, distinguish homogeneous from inhomogeneous broadening mechanisms [44], [45], and reveal quantum coherences that arise in cavity-modified chemistry [46], [47]. Such information is essential for the rational design of polaritonic devices, yet remains inaccessible to linear spectroscopy. The interpretation of such experiments, however, demands quantitative theoretical modeling that captures the full complexity of light–matter interaction beyond perturbative approximations [48], particularly when molecules are coupled to plasmonic metasurfaces [49] or nanostructured cavities, where field confinement, spatial inhomogeneity, and retardation effects fundamentally shape the polaritonic landscape. Current theoretical approaches, predominantly based on simplified cavity models or perturbative expansions, fail to capture the rich physics emerging at plasmonic interfaces, where near-field enhancement, radiative damping, and molecular response are intimately intertwined. The development of computational methods capable of simulating two-dimensional spectroscopic signals from first principles, accounting for the full electromagnetic complexity of metasurfaces and the quantum mechanical nature of molecular systems, thus represents a critical frontier in polaritonic chemistry – one that could bridge the gap between experimental observations and mechanistic understanding of cavity-controlled chemical processes.

On the local scale, the numerical framework allows to investigate how the collective strong coupling affects microscopic details of molecules. Specifically, one can go beyond available results on collectively induced local strong-coupling effects, which commonly assume only a single, spatially homogeneous mode [20], [21], [22], [23], [24], [25], [26], and investigate the interplay between the spatial inhomogeneity of all the cavity modes and local effects. This allows crucial new insights into the question how collective strong coupling can change local chemical properties.

Another important aspect of the dftb + Maxwell framework is its computational simplicity, especially for one-dimensional cavities [43]. This allows to emphasize user-friendliness, an aspect that is often overlooked when developing novel simulation tools. Specifically, when interpreting experimental results, experimentalists favor simplified effective models over often complex ab initio light–matter simulation tools. However, to our knowledge, apart from very recent symmetry approaches [50], there is no effective model in polaritonic chemistry that can make reliable chemical predictions, beyond (collective) spectral properties (e.g., the collective Rabi splitting). For this reasons, we believe it is vital to develop and provide a user-friendly simulation framework that captures above features and provides reasonable ab initio predictions on a local computer within short simulation times. Providing such a tool will not only allow to quickly test and explore novel ideas, e.g., modified cavity geometries, losses and frequencies, or different solute–solvent configurations, before performing an experiment, but it will also help scientists to obtain a more intuitive understanding of the feedback between the different scales. We believe the dftb + Maxwell method will be capable to fill this gap.

In the following perspective article, we outline how the self-consistent combination of dftb with a fdtd Maxwell solver allows to address all these objectives. We illustrate its current capabilities with three compelling examples and eventually outline its versatility, focusing on different future developments that provide access to unexplored regimes of nonperturbative strong light–matter interaction. Overall, our goal is to introduce a simple but holistic simulation tool for the qualitative and quantitative design of polaritonic experiments and beyond.

Figure 1:

In polaritonic chemistry, the strong coupling of light and matter in optical cavities introduces feedback effects across different spatiotemporal scales, which are not captured by perturbative or standard computational chemistry methods (illustration taken from Ref. [51]).

2 Theory and computation

The coupling of the dftb Hamiltonian to an arbitrary external electric E follows Maxwell’s electromagnetic field equations. As a result, the time-evolution of the electric E and magnetic fields B can be determined by solving,

discretized on a real-space grid (see Figure 2). The vacuum speed of light is given by c ₀ and the vacuum permittivity by ϵ ₀. The coupling of the fields to the matter is due to the current density arising from the time-derivative of the polarization P .

Figure 2:

Illustration of the one-dimensional grid representation to solve the Maxwell equations for the propagating fields. Notice that the computational approach allows for a mixed representation of polarizable matter. For example, mirrors can be represented by the Drude model, whereas molecules can be described on a dftb level of theory in three dimensions.

In our implementation, the time-dependent polarization can be calculated on different levels of theory. For various metallo-dielectric interfaces such as mirrors, we have implemented several approaches. The simplest representation employs the continuum Drude model for free electrons, while more accurate treatments of material dispersion utilize Drude–Lorentz parametrizations that include interband transitions alongside the free-electron response [52]. We have also implemented modified Lorentz models that can account for additional resonances and damping mechanisms beyond the standard formulations [53]. For scenarios requiring the description of nonlinear optical phenomena, we employ a semiclassical hydrodynamics model for conduction electrons that provides quantitative treatment of nonlinear effects in metals while maintaining computational efficiency within the Maxwell solver framework [54].

Alternatively, for a more accurate microscopic description of matter, molecular polarizations can be determined by placing single molecules on the Maxwell grid, whose response to the transverse fields is determined from dftb calculations. Notice that in our implementation, the transverse fields are modified (locally) by the coupling to the matter, thus also affecting the coupling of matter to light at later time-steps. This allows to investigate microscopically, the feedback effects between light and matter.

The dftb Hamiltonian is derived from Kohn–Sham density functional theory by expansion of the energy functional around a reference density ρ ₀, which is truncated at a certain order (for this study, second order) around the density fluctuations [55], [56]. Then, the molecular orbitals ψ _i = ∑ _ν C _iν ϕ _ν are represented as a linear combination of atomic orbitals (LCAO), and a minimal basis set of pseudo-atomic orbitals {ϕ _ν } is used. After applying the two-center approximation on the Hamiltonian, and writing the density as a superposition of atomic contributions, the Hamiltonian matrix elements can be written as follows:

Here, the subindices I, J indicate different atoms, and Δq _K = q _K − Z _K , where Z _K indicates the charge of a neutral atom of species K and q _K are the time-dependent Mulliken charges, defined as q _K = Tr _K [ρS]. The reduced one-electron density matrix is indicated by ρ, S _μν ( R ) = ⟨ϕ _μ |ϕ _ν ⟩ are the overlap matrix elements and Tr _A [⋅] the partial trace over the orbitals centered on atom A. The γ _IK analytical function described the electrostatic interactions in DFTB, and it depends on the chemical hardness, which is computed from DFT calculations [57]. The non–self-consistent tight-binding Hamiltonian H ₀, defined as

are explicitly precomputed from DFT as a function of the ineratomic distance R . Together with the overlap matrix elements S _μν ( R ) these quantities are available in Slater–Koster tables, which are used and interpolated for the needed interatomic distance in the calculations. This allows for a rapid computational evaluation of the eigenvalue problem for the molecular orbitals

as well as for the update of the time-dependent quantities when doing a time-propagation (see below). Notice, however, that the atomic charge fluctuations Δq _I = q _I ( C ) − Z _I depend on the orbital coefficients C _μν . Consequently, Eq. (5) must be solved self-consistently.

Based on this formalism, a real-time time-dependent dftb (rt-TD-DFTB) formulation can be derived from the Runge–Gross theorem [58], [59]. When an additional external electric field E (for example, given by Maxwell equations) is added, the Hamiltonian matrix can be written as [59]

Here, the external potential is applied under the electric dipole approximation, where R _I is the position of atom I. When coupling to an external field calculated from a Maxwell simulation simultaneously, the value of the external field at the position of the molecule r _mol on the Maxwell grid is considered [43].

The time-propagation is given in terms of the one-electron density-matrix operator ρ μ ν = ∑ i ∈ o c c C ̃ i μ C ̃ i ν * , and by propagating the well-known Liouville–von Neumann equation for a nonorthogonal basis

Here, H ̃ μ ν = H μ ν − i D μ ν is a shorthand notation for an effective Hamiltonian that includes both electronic and nuclear dynamics contributions. The nonadiabatic coupling matrix elements defined as D μ ν = ⟨ ϕ μ | ϕ ̇ ν ⟩ are considered for Ehrenfest electron–ion dynamics (otherwise they are neglected). In practice, these terms are calculated as D μ ν = ∑ α ( ∂ α S μ ν ) R ̇ I α , ν ∈ I , α = { x , y , z } , where R ̇ I α are the components of the velocity of the nucleus of atom I [59]. The Ehrenfest dynamics at the level of dftb has been successfully applied to the study of light-driven nuclear dynamics in molecules [60], [61], metallic clusters [62], and nanoparticles [63]. Alternatively, Born–Oppenheimer molecular dynamic simulations can also be performed based on the solution of the stationary electronic-structure problem given in Eq. (5), while including the external field in a quasi-static fashion.

Both the Ehrenfest and Born–Oppenheimer dynamics are implemented in the DFTB + package [57], which we use as a library and coupled to the Maxwell equations by suitable routines of the DFTB + library that have been implemented for this purpose. Finally, the extended tight binding (xTB) methods are also supported in DFTB+ [57], [64], [65], giving access to simulations for molecules involving almost all elements of the periodic table.

Eventually, the time-dependent molecular polarizations are computed by

and plugged into Eq. (2). Here, N _M is a concentration term, which gives us the number of molecules represented by the simulated molecule in a single grid point and ⟨ μ ̂ ⟩ = Tr ρ ̂ μ ̂ . Notice, our current assignment of molecules to the Maxwell grid, as well as the scaling argument of Eq. (8) assumes no overlap of the electronic structure between different Maxwell-grid points. Thus, any intermolecular interactions remain entirely classical (no exchange and no correlation effects).

3 Illustrative applications

To demonstrate the capabilities of the dftb + Maxwell approach, we begin by briefly discussing three distinct applications. The first application focuses on global properties, which can be accessed by two-dimensional spectroscopy. Afterward, we illustrate local molecular changes in optical cavities as well as how to engineer cavity mirrors to reach a desired local effect. In all instances, the Maxwell equations are solved on a grid; however, different quantum descriptions for the matter are applied.

3.1 Two-dimensional spectroscopy

The Maxwell integrator (see Eqs. (1) & (2)) of dftb plus fdtd is under active development, representing a significant advancement in computational electrodynamics for complex light–matter systems. Current implementations support a wide range of dispersive material models, including Drude, Lorentz, and Drude–Lorentz formulations, enabling accurate description of metallic, dielectric, and hybrid structures. The framework handles both linear optical phenomena and nonlinear optical processes including second- and third-harmonic generation, four-wave mixing, and intensity-dependent refractive index changes. These capabilities make the platform particularly suited for studying nanophotonic devices, metamaterials, and strongly coupled light–matter systems. One of the most exciting recent developments is a novel numerical scheme that enables the calculation of two-dimensional coherent spectra for arbitrary optical systems. Two-dimensional spectroscopy has emerged as a useful experimental technique for unraveling complex dynamics in condensed matter systems, providing information about energy transfer pathways, coherent couplings, and many-body interactions that remain hidden in conventional linear spectroscopy. In strongly coupled systems, where light and matter degrees of freedom hybridize to form new quasiparticles, two-dimensional spectroscopy can reveal quantum coherence, distinguish homogeneous from inhomogeneous broadening, and map energy relaxation cascades with unprecedented detail. Despite its experimental success, theoretical modeling of two-dimensional spectra has remained challenging. The majority of existing theoretical approaches rely on perturbative expansions, phenomenological models, or simplified system-bath interactions that may not capture the full complexity of realistic systems. More critically, truly self-consistent simulations that solve Maxwell’s equations simultaneously with quantum mechanical equations of motion have been conspicuously absent from the literature. This gap exists not due to lack of interest, but rather because of a fundamental technical challenge: in fdtd-based approaches where electromagnetic fields and material polarizations are evolved self-consistently, the time-domain electric and magnetic fields inherently contain contributions from all orders of the optical response – linear absorption, nonlinear mixing, scattered fields, and cavity feedback effects are all superimposed in the raw simulation data. This mixing makes it seemingly impossible to isolate the specific third-order nonlinear signals that constitute the two-dimensional spectroscopic response, as traditional postprocessing techniques cannot disentangle these contributions after they have been combined in the time-domain propagation.

The new numerical technique we are developing relies on the propagation of several coupled Maxwell–Schrödinger equations in parallel. Figure 3 illustrates the method using a numerical experiment with two pumps and a probe. The signal of interest is the third-order nonlinear process, which is obtained by propagating equations: (1) with all three pulses; (2) with pump 1 only; (3) with pump 2 only. The resulting signal is obtained by subtracting the electromagnetic fields obtained in (2) and (3) from the fields simulated in (1). All calculations are performed in the time domain. A representative example is shown in Figure 3(a), comprising two metal slabs sandwiching a molecular layer consisting of two-level emitters treated using standard Bloch equations. The parameters of the system are chosen such that the linear transmission exhibits two well-resolved polariton states shown in Figure 3(c). The results of the two-pump/probe sequence are presented in Figure 3(d)–(f). The two-dimensional spectrum reveals the rich nonlinear optical response of the strongly coupled cavity-polariton system through several distinct spectral features. The prominent diagonal peaks at (1.575, 1.575) eV and (1.817, 1.817) eV correspond to lower polariton (LP) and upper polariton (UP) self-interactions, respectively, confirming efficient excitation of both polariton branches with a Rabi splitting of 242 meV. A clear cross-peak at (LP, UP) demonstrates coherent coupling between LP and UP states, while the absence of the corresponding (UP, LP) peak indicates asymmetric relaxation favoring energy flow from upper to lower polariton branches. Most intriguingly, a narrow peak appears at (LP, 1.698) eV, where the signal frequency matches the bare emitter transition, suggesting the important role played by the dark states.

Figure 3:

Two-dimensional spectroscopy simulations: Panel (a) shows a schematic setup. Panel (b) illustrates the two pumps + probe sequence. Panel (c) shows the linear transmission spectrum indicating lower and upper polaritons. Panels (d)–(f) show the resulting two-dimensional signal as a function of the probe and signal frequencies. Absolute value of the signal is shown in panel (d), and panels (e) and (f) show its real and imaginary parts, correspondingly.

The computational framework demonstrated here for two-dimensional spectroscopy of strongly coupled cavity-emitter systems represents just one facet of the broader capabilities enabled by interfacing the Maxwell integrator with dftb electronic-structure methods. By combining the nonperturbative propagation of the full Maxwell-matter dynamics with dftb’s efficient treatment of molecular electronic structure, our approach surpasses the limitations of current methodologies that rely on perturbative expansions [66] – limitations that become particularly severe in the strong coupling regime where light–matter hybridization fundamentally alters both optical and molecular properties [67]. The self-consistent treatment of both the electromagnetic cavity modes and molecular degrees of freedom within a unified computational framework eliminates the need for artificial partitioning between “system” and “bath” that underlies most contemporary approaches, instead capturing the full complexity of polaritonic states and their nonlinear optical response from first principles. This capability to perform quantitative simulations of multidimensional spectroscopy without invoking rotating wave approximations, Markov assumptions, or perturbative field-matter coupling represents a significant advance toward predictive modeling of cavity-modified chemistry and opens new avenues for designing and optimizing polaritonic devices, where strong coupling effects dominate the system dynamics.

3.2 Collective strong coupling in optical cavities

Experimentally, collective strong coupling is manifested as a larger separation of the polaritonic peaks in the transmission spectrum (Rabi splitting) when the concentration of molecules inside the cavity is increased. To capture this phenomenon, we solve the Maxwell equations on a one-dimensional grid, while different numbers of molecules are described using dftb. The molecules are always resolved in three dimensions, providing a realistic representation of molecular structures. Allocating two sections of the one-dimensional grid to represent the mirrors, using the Drude model, it becomes possible to mimic a Fabry–Pérot cavity.

As a conceptual example, we work with nitrogen molecules (N₂), which possess a first electronic transition at 13.902 eV, according to dftb. Every molecule was placed on separate grid points spaced by 1 nm, forming molecular layers ranging from 2 nm to 80 nm in the center of the cavity. To construct a cavity that resonates with this transition, we set two Drude mirrors separated by 123 nm. The mirror thickness was chosen 20 nm, and the Drude parameters were set to Ω _p = 34 eV and γ = 0.181 eV. Under these conditions, the third mode of the cavity resonates with the electronic transition.^[1] Figure 4(a) shows the different transmission spectra obtained by explicitly varying the number of molecules simulated. As we can see, increasing the number of molecules leads to a larger Rabi splitting. The inset in Figure 4(a) shows the dependence of the Rabi splitting on the square root of the number of molecules N. It can be seen that the typical linear scaling with N breaks down for a certain amount of molecules. This is because the linear behavior is obtained assuming that all the molecules experience the same electric field, while in these simulations each molecule is subject to different field amplitudes depending on its position. For a small number of molecules located around the center of the cavity, the space occupied by them is less than half the wavelength of the excitation field, and therefore the field is almost constant and the linear dependence of the Rabi splitting with respect to N is reproduced. However, for a larger molecular ensemble, their spatial extent becomes comparable to the wavelength of the resonant cavity mode, breaking the N behavior. We thus see that the spatial distribution can have a profound impact on these relations.

Figure 4:

Ensemble of nitrogen molecules confined by a Drude mirror cavity: (a) Transmission spectra obtained after exciting the cavity and different number of simulated molecules. The gray dashed line indicates the first electronic excitation energy of N₂. Inset: Rabi splitting as a function of the square root of the number of molecules. (b) Fourier-transformed dipole moments of the simulated molecules as a function of their position in the cavity, for the system comprising 80 molecules. The color scale has been normalized with the maximum amplitude of the dipole moment.

To shed more light on the consequences of the spatial distribution of the cavity modes, we next analyze individual molecular properties. By Fourier transforming the time-dependent dipole moment, we can obtain information about the spectral response of each molecule. In Figure 4(b), the Fourier transform of the dipole moment for the case of 80 molecules is shown as a function of the position of the molecules. We can observe that the polaritonic peak intensities (at 13.7 eV and 14.06 eV) exhibit a spatial distribution that matches the expected profile of the cavity’s third mode, with nodes near one-third and two-thirds of the mirror separation distance (z ≈ ± 22 nm in the plot), evidencing local effects as a result of the realistic description of the molecular distribution in the cavity. In fact, there is experimental evidence of the position-dependent nature of the Rabi splitting when a J-aggregate layer is displaced inside the cavity, following a mode-like profile [68].

Figure 4(b) also reveals more features than are present in the transmission spectra (Figure 4(a)). The signals obtained above and below 13.9 eV correspond to dark states created by the electronic excitation of the molecules and to the off-resonant interaction between the molecules and higher modes, respectively. These results demonstrate the ability of our implementation to access local information that is lacking in the transmission spectrum.

3.3 Automatized cavity design

Inverse design is a powerful methodology for engineering optical devices that fulfill predefined characteristics. Among the many existing approaches, density-based topology optimization methods are the most flexible and versatile [69], [70]. We have included this tool in our implementation to design optical cavities in two dimensions, optimized for strong coupling with a selected molecule and to enhance the interaction between groups of the same molecule at different positions and scales.

The algorithm solves the Maxwell equations in the frequency domain on a two-dimensional grid at a fixed frequency corresponding to a chosen molecular transition. The optimization seeks to maximize E ( ω ) 2 at one or more target points, when a point-like current source is located at a specified position in the simulation box. To this end, we calculate the gradient d E ( ω ) 2 / d ε ( r , ω ) , where ɛ(r, ω) is the corresponding permittivity at the frequency ω and in the grid point with coordinate r. Using this gradient, the permittivity at each grid point within a selected region is updated continuously in each optimization step. A binary projection is applied so that the permittivity takes values corresponding to either the chosen design material or vacuum. The obtained designs are tested afterward in our dftb + Maxwell implementation including the molecules explicitly.^[2] For example, Figure 5(a) and (b) shows two cavities made of silicon and gold, respectively, designed to be in resonance with pentacene molecules (electronic transition at 3.97 eV), and to couple three groups of molecules separated by 300 nm. The designs aims to maximize E ( ω ) 2 at two target points located at (−300 nm, 0 nm) and (300 nm, 0 nm), when a point-like current source is placed at (0 nm, 0 nm). In the same figure, a profile of the resulting E ( r , ω ) 2 in the whole simulation box is also plotted, demonstrating the desired increase at the selected points as well as at the source point. To test each cavity, we conducted three simulations including one, two, and three groups of 37 pentacene molecules, respectively, located at one target point, both target points, and both target points along with the source point. The largest simulation, propagating three groups of molecules over 700 fs, took about 15 h using the resources mentioned above. In all cases, the system is excited with a short point-like pulse, positioned at the source point and tuned at ω = 4.0 eV. The signal transmitted out of the cavity is utilized to obtain the transmission spectra presented in Figure 5(c) and (d), for the silicon and gold cavities, respectively. In both cases, we observe that the presence of just one group of molecules already yields a Rabi splitting, indicating that we are in the strong-coupling regime. Notice that even though we fixed the designing frequency, the cavity mode does not exactly match the dressed pentacene transition frequency. In the case of the gold cavity, the difference between the cavity mode and the dressed molecular transition causes an off-resonant interaction, which explains the asymmetry observed in the Rabi splitting. As we increase the number of groups of molecules, the separation between the polariton signals increases. This demonstrates that, despite the large distance between the target points, the different groups of molecules are interacting through the cavity mode and form a polaritonic ensemble.

Figure 5:

Computational cavity design: (a) and (b) Present two examples of cavities designed in 2D using two different materials, silicon and gold, respectively. The designed cavities were extended into the third dimension for visualization purposes. The obtained designs aim to maximize E ( ω ) 2 at two target points located at (−300 nm, 0 nm) and (300 nm, 0 nm), when a point-like current source is placed at (0 nm, 0 nm). The design frequency was 3.97 eV, corresponding to the first electronic transition of the pentacene molecule. The resulting E ( ω ) 2 distribution in the simulation box, for the respective cavities, is also shown in these panels. (c) and (d) Display the transmission spectra obtained from DFTB + FDTD simulations after exciting the system with a short point-like pulse, positioned at (0 nm, 0 nm) and tuned to ω = 4.0 eV, within the silicon and gold cavities, respectively. The spectra compare the empty cavity (blue curve) with cavities containing different numbers of groups (37 molecules each), located at one target point (orange curve), two target points (green curve), and two target points and source point (red curve).

Even though the presented example focused on the optimization of E ( ω ) 2 through modification of the permittivity ɛ(r, ω), our implementation can be extended to the optimization of other field magnitudes, as well as those related with the magnetic field B(r, ω). Some of these cases will require the magnetic permeability (μ(r, ω)) as an additional optimization variable, which can be included straightforwardly.

4 Outlook

While the dftb + Maxwell framework is already able to address various different problems of polaritonic chemistry and strong light–matter coupling, there are several directions that are actively developed and will enhance the capabilities of this theoretical toolkit even further. In the following, we briefly sketch the main directions:

1. Thermal Equilibrium: In the context of thermal effects, the self-consistently coupled dftb + Maxwell framework represents a more general problem than usually encountered in standard ab initio Molecular Dynamics (MD) setups (e.g., Ehrenfest or Born–Oppenheimer MD). Here, not only the nuclear degrees of freedom can be coupled to a thermal bath (e.g., by Langevin equations), but also the field can be thermostated, for example, by explicitly accounting for the black-body spectrum in the transverse fields [71]. The inclusion of transverse degrees of freedom breaks the usual conservative-field assumption and hence introduces the possibilities of novel effects. Notice that under collective strong-coupling conditions, the definition and characteristics of thermal (non)-equilibrium becomes a nontrivial problem of fundamental theoretical interest [19], [25], [72], but also practical relevance. For example, modified thermodynamics would alter chemical reaction rates and thus corresponding theoretical considerations [73], [74], [75], [76], [77], [78], [79] would require adjustment. For example, from a statistical point of view, the strongly coupled molecular ensemble can no longer be partitioned into weakly interacting subentities that can be regarded statistically independent and form the basis for a thermal statistical description of large molecular ensembles [72]. Moreover, on resonance, theoretical models show that cavity-mediated intermolecular electron correlations, i.e., London-dispersion forces, are substantially modified, which can lead to spin-glass-like effects (see point 4 below) that indicate the violation of standard fluctuation–dissipation relations [25]. The tackle these fundamental theoretical challenges requires efficient and versatile computational approaches. We believe the dftb + Maxwell framework is ideal to investigate the above-mentioned issues.

2. Beyond Gaseous Phase: In our compact theoretical introduction of dftb + Maxwell, we have focused on gas-phase molecular setups. To account for liquid condition, i.e., longitudinal intermolecular interactions, it is straightforward to replace the dftb molecules placed on the Maxwell grid by explicit molecular ensembles on a dftb level of theory. For example, this would allow to explicitly account for the solvation shells surrounding the solute molecules. Moreover, using the basics of qedft, which makes dftb a limiting case of qedft, we can extend density-functional embedding methods to account also for longitudinal and transverse interactions between the molecular ensembles at different (macroscopic) Maxwell grid points.

   Using such an embedding strategy, it becomes possible to also include local (within one macrosopic grid point) transverse field effects in the dftb + Maxwell setup. However, these single-grid-point local-field effects might not be overly important. The fully quantized Pauli–Fierz theory reveals that transverse local-field effects (from the dipole self-energy term) are expected to be significant mainly on intermediate to large length scales, exceeding typical Coulomb correlation lengths [25]. Thus, the explicit modeling of these interactions via the macroscopic Maxwell equations as done in the Maxwell + dftb might be accurate enough.

3. Collective Quantum/Correlation Effects: Choosing the classical Maxwell equations to propagate the macroscopic transverse and longitudinal fields in time, a priori removes any physical effect that may emerge on this scale due to the quantization of the electro-magnetic field, e.g., spontaneous emission or Casimir forces. In particular, the restructuring of the continuum of vacuum modes in optical cavities cannot be resolved. However, this is considered a crucial ingredient for going from weak to strong coupling conditions. For example, in weak coupling, the vacuum-field fluctuations are responsible for the Casimir forces and the Purcell enhancement of spontaneous emission. Under (collective) strong coupling, the relevance of vacuum fluctuations is still under debate. In particular, it remains unclear to what extend and under which conditions the quantum nature of the electromagnetic field is key and surpasses standard thermal fluctuations. Experimental evidence suggest that in particular meeting collective-resonance conditions and (modified) London dispersion forces are crucial ingredients to explain changes of chemical properties under collective strong coupling. For example, modified bulk properties such as self-aggregation or changed conformational equilibrium constants were observed. In addition, rare events (tunneling and chemical reaction rates) could also be altered. Those chemical effects are believed to be a consequence of the quantized nature of the transverse electric fields. In more detail, as recently shown in Ref. [25], those effect could be triggered by the spin-glass-like nature of transverse electron correlations in optical cavities. Moreover, the explicit theoretical mapping onto a spin glass does not only support aforementioned experimental facts theoretically, but it also predicts phase-transitions, heating and hysteresis effects under resonant collective strong coupling [25]. All of which have been observed experimentally in various situations [15], [80], [81], [82], [83], [84], [85], [86], [87], [88], [89].

   In order to include these quantum effects in the dftb + fdtd scheme, we can again rely on the fact that dftb + Maxwell is the mean-field limit of qedft. Thus, by including transverse exchange–correlation functional approximations [90], [91], [92], [93] (on top of the standard longitudinal exchange–correlation functional approximations, which describe the Coulomb interaction, and are already included in dftb + Maxwell), we can efficiently and accurately simulate these novel quantum effects also in this framework.

4. Photon-field Fluctuations by Semiclassical Methods: Another research direction that can be explored with dftb + Maxwell methods is the possibility to recover quantum effects of light such as vacuum fluctuations and correlations via multitrajectory Ehrenfest methods. These approaches are based on propagating a swarm of Ehrenfest trajectories, where the classical photon field coordinates are sampled from a Wigner distribution, and the observables are ensemble averages of the individual trajectories [94], [95]. Among the different semiclassical methods that have been bench-marked, the multitrajectory Ehrenfest method qualitatively describes the dynamics of electron and photon populations and captures bound photonic states. While these calculations have been done considering a mode-resolved version of the Maxwell equations, they can be adapted to the case of spatially resolved fields. One of the advantages of dftb simulations for this resides in the large amount of trajectories that are needed, on the order of 10⁵, to achieve converged results. While this number becomes prohibitively expensive for full ab initio simulations of realistic systems, they become feasible by the use of dftb.

5. Nuclear Quantum Dynamics: So far our implementation of dftb + Maxwell only allows for classical trajectories of the nuclei. However, for molecules in a dilute gas phase, quantum dynamical effects can start to play a role. For example, the rotational dynamics of apolar diatomic gases can be well approximated by the quantum rigid rotors. In that sense, it would be very interesting to investigate the feedback between classical Maxwell fields and the quantized dynamics of nuclei. At least for dilute gases, our framework should be extendable to incorporate nuclear quantum dynamics for example by replacing the dftb part by multiconfiguration time-dependent Hartree propagations [96], [97] or hierarchical equations of motion [98], [99]. Clearly, the time-propagation would be restricted to small molecules and short simulation times. Still, collective feedback effects on the quantum dynamics should become accessible for many molecules thanks to the dilute gas limit. In that case, the different molecules can be considered interacting only via classical Maxwell fields. This effectively suppresses the scaling of the quantum dynamical complexity with respect to the number of molecules.

6. Simplicity and User-Friendliness: The gap between theory and experiments in polaritonic chemistry is enhanced due to the vastly different interpretational tools. Since currently most numerical ab initio QED tools are not easy to use without much prior experience, experimentalist favor simplified models to interpret the observed effects. We, therefore, think that an efficient and user-friendly numerical toolkit will be central to bridge this gap. That is, we aim to develop a framework that can be easily employed to interpret experimental results on the fly and even allow to optimize experimental details in situ. Thanks to the computational speed of the one-dimensional Maxwell solver in combination with semi-empirical performance of dftb, it becomes possible to aim for cloud platforms, which allow users to execute tasks directly from their web browser. This reduces the entry barrier for setting up a computational research environment [100]. In that regard, we plan the future integration of dftb + Maxwell into AiiDAlab and make our code open source. AiiDAlab is a web platform that enables computational scientists to package scientific workflows and computational environments and share them with their collaborators and peers. It allows for user-friendly graphical interfaces of key functionalities and is made available through the AiiDAlab app store [100]. Still more complex functionalities, that are only applicable on large high-performance computing facilities, can be developed concurrently.

7. Three-Dimensional Maxwell Grid: Our Maxwell solver incorporates advanced parallelization capabilities that enable simulations of large ensembles of quantum emitters in complex electromagnetic environments. The implementation employs a novel molecular-mapping technique built on top of the conventional three-dimensional domain decomposition methodology, which efficiently balances computational load across parallel processors and allows direct access to collective dynamics by drastically increasing the number of molecules that can be simulated simultaneously [101]. This parallelization framework is being merged with electronic-structure codes, enabling us to perform large-scale simulations of many molecules interacting with complex three-dimensional metasurfaces. The combination of efficient parallel computing with first-principles molecular calculations opens new possibilities for studying collective optical phenomena in realistic nanostructured environments, where both the electromagnetic field distribution and the quantum mechanical response of individual emitters are treated with high fidelity [102].

5 Summary and conclusion

This perspective article demonstrates the numerous capabilities and opportunities of a computationally efficient and versatile simulation framework for (strong) light–matter interaction in realistic configurations. The dftb + Maxwell framework bridges macroscopic classical field propagation with microscopic molecular resolution on a quantized level of theory, including their respective feedback effects in space and time. Capturing this delicate interplay offers nonperturbative, quantitative insights into (multidimensional) spectroscopic observables (see, e.g., Section 3.1) with detailed access to local chemical changes, for example, induced by the hybridization of light and matter inside optical cavities (see, e.g., Section 3.2). Eventually, this computational framework can be employed to specifically design the photonic environment on demand (see, e.g., Section 3.3). Based on these fundamental capabilities, we have highlighted different opportunities to enhance the predictive power of this toolkit, including temperature or collectively correlated quantum effects for coupled light and matter in a nonperturbative fashion. However, all of these effects can lead to an intricate feedback between different scales, and hence common simplified interpretational tools can fail to capture these important details. Therefore, an essential ingredient for the success of this computational framework is its computational simplicity and speed that extends to collective effects of strong light–matter interaction. Fortunately, for many practical cases, the one-dimensional representation of the Maxwell fields is sufficient together with a three-dimensional representation of molecular ensembles. In this case, the dftb + Maxwell method allows for an almost real-time exploration of chemical or physical parameters on a desktop machine. At this point, we would like to stress that a 3D Maxwell-field representation is possible using the here-presented scheme. Its implementation will be presented elsewhere, and it will be made accessible open source as a part of the octopus simulation suite [103]. Overall, we believe our approach provides a great opportunity to be developed into a user-friendly tool suitable for in situ interpretation of experimental setups and the design of novel industrial applications, without having to rely on oversimplified effective model systems or dedicated high-performance computing infrastructure.