Periodic Reporting for period 3 - BET3D (March-on-in-Time: Boundary Element Time-Domain Domain Decomposition Methods)
Reporting period: 2024-09-01 to 2026-02-28
Summary of the context and overall objectives of the project
Modelling almost entirely replaced prototyping as a design methodology. In electromagnetics and optics, it has played a central role in the rapid development of communication and imaging, with applications in medicine, security, and energy. New technological requirements have always been met by breakthroughs in modelling.
Advances in TeraHertz technology have driven researchers to consider systems that are highly non-linear and strongly radiating. Modelling these systems requires methods that are based on time-stepping, and that can accurately describe unbounded regions. These requirements can only be met by time-domain boundary element methods (TD-BEMs). The use of other methods such as the finite element method, the method of moments, or the finite difference time-domain method – if at all applicable – results in unacceptable computational costs, and large errors in dispersion and radiation characteristics.
Unfortunately, TD-BEMs are not nearly as mature as other methods. They lack the ability to model all but the simplest systems. Attempts to use TD-BEMs to model realistic devices containing multiple materials, ports, coatings, or two-dimensional materials like graphene lead to instabilities, rendering the result of simulations completely useless. This state of affairs has persisted for over 50 years but has become an urgent problem now.
BET3D will aggressively accelerate the development of the time-domain boundary element method into a method that can model highly non-linear, strongly radiating systems and that unlocks the capability to model emerging technologies in THz communications and imaging, and beyond.
As part of the benchmarking and validation for this project, we will model the generation of THz radiation by an array of graphene resonant tunnelling diode oscillators, including feeds, ports, packaging, and antenna structure. The results of BET3D will enable modelling of non-linear, radiating systems, also in acoustics, elastodynamics, and fluid dynamics.
The objectives of this project are:
• Design a Boundary Element Time-Domain Domain-Decomposition Method: this allows modelling of systems of all shapes and compositions
• Stabilize the BET3D method in the late-time and multi-scale regimes: we want to make sure that errors do not overtake the real solution computed by our modelling tools
• Improve accuracy of the BET3D method by designing novel approximation and integration methods: the more accurate the description of the currents and the fields in the device we wish to study, the higher the quality of the solution, and the lower the chance that instabilities will occur
• Enable coupling to circuits (0D), transmission lines (1D), thin sheets (2D) and inhomogeneous materials (3D): modern systems contain two-dimensional materials such as graphene and are powered by circuits connected through gates. An integrated modelling approach should be able to take into account all these subcomponents.
• Design a flexible and robust parallel algorithm for the distributed execution of BET3D and validate on the modelling of an array of GRTD oscillators: to be able to model highly complicated real-life devices, the required computational costs should be minimized and distributed across an entire network of collaborating computers.
Advances in TeraHertz technology have driven researchers to consider systems that are highly non-linear and strongly radiating. Modelling these systems requires methods that are based on time-stepping, and that can accurately describe unbounded regions. These requirements can only be met by time-domain boundary element methods (TD-BEMs). The use of other methods such as the finite element method, the method of moments, or the finite difference time-domain method – if at all applicable – results in unacceptable computational costs, and large errors in dispersion and radiation characteristics.
Unfortunately, TD-BEMs are not nearly as mature as other methods. They lack the ability to model all but the simplest systems. Attempts to use TD-BEMs to model realistic devices containing multiple materials, ports, coatings, or two-dimensional materials like graphene lead to instabilities, rendering the result of simulations completely useless. This state of affairs has persisted for over 50 years but has become an urgent problem now.
BET3D will aggressively accelerate the development of the time-domain boundary element method into a method that can model highly non-linear, strongly radiating systems and that unlocks the capability to model emerging technologies in THz communications and imaging, and beyond.
As part of the benchmarking and validation for this project, we will model the generation of THz radiation by an array of graphene resonant tunnelling diode oscillators, including feeds, ports, packaging, and antenna structure. The results of BET3D will enable modelling of non-linear, radiating systems, also in acoustics, elastodynamics, and fluid dynamics.
The objectives of this project are:
• Design a Boundary Element Time-Domain Domain-Decomposition Method: this allows modelling of systems of all shapes and compositions
• Stabilize the BET3D method in the late-time and multi-scale regimes: we want to make sure that errors do not overtake the real solution computed by our modelling tools
• Improve accuracy of the BET3D method by designing novel approximation and integration methods: the more accurate the description of the currents and the fields in the device we wish to study, the higher the quality of the solution, and the lower the chance that instabilities will occur
• Enable coupling to circuits (0D), transmission lines (1D), thin sheets (2D) and inhomogeneous materials (3D): modern systems contain two-dimensional materials such as graphene and are powered by circuits connected through gates. An integrated modelling approach should be able to take into account all these subcomponents.
• Design a flexible and robust parallel algorithm for the distributed execution of BET3D and validate on the modelling of an array of GRTD oscillators: to be able to model highly complicated real-life devices, the required computational costs should be minimized and distributed across an entire network of collaborating computers.
Work performed from the beginning of the project to the end of the period covered by the report and main results achieved so far
BET3D has been focusing on the design of integral equation methods for the simulation of systems that can be made up from various materials that are arranged in any geometric configuration. The final goal of BET3D is to be able to perform these simulations in the time-domain (this means constructing the solution at subsequent times). This has been achieved for some of the simpler systems such as those composed entirely from perfect conductors or that are made up of a single material. For more complicated systems, a lot of preparatory work has been done in the frequency domain (where the impact of signals containing only a single frequency is investigated). In the frequency domain, modelling of fully general systems (in terms of geometry and in term of the materials from which it is built) has been covered.
Time domain simulations are notorious for being unstable: errors that enter the simulation process are liable to magnify as they are being used as inputs to compute the solution at subsequent time steps. This error magnification goes on until they completely dominate any useful information in the output of the simulation. For perfectly conducting systems and single body systems, we have constructed methods that are guaranteed to be free from such instabilities. What is more, these modelling methods can be used when both rapidly and extremely slowly varying inputs are present and when the geometry under investigation contains both large features and small details.
This work is central to enable the modelling of strongly non-linear and radiating systems such as encountered in the design of next-gen THz technology. Indeed, these systems are made up of various materials, and contain both large radiating structures and extremely small semi-conductors and two-dimensional components.
We have begun research in simulators that do not directly reconstruct the electric and magnetic fields in the devices of interest, but first reconstruct the so-called magnetic vector potential. This quantity is arguably more fundamental then the electric and magnetic fields. In quantum mechanics and quantum field theory, photons (the fundamental particles making up electromagnetic radiation) are described in terms of this vector potential and some famous experiments (such as the Aharonov-Bohm effect) provide evidence that the slightly more abstract vector potential is the true source of electromagnetic interactions. For modellers, expressing fields in terms of the potential has a number of computational advances: the behaviour of the vector potential can be studied more easily when multiple frequencies are present. We have been able to show that boundary integral equations for very general geometries can be constructed for the vector potential.
When modelling either the electromagnetic field or the vector potential in a discrete, computational setting, so-called basis functions that describe the behaviour of these fields in a small region in space are required. In general, a pair of dual sets of basis functions is required, such that one set can be easily expressed in terms of the other. Without such a pair it is not possible to build modelling tools that provide answers in a reasonable number of computations. When the unknown field lives on a surface, such a pair is known but unfortunately will introduce excess error on one of the fields that is being computed. We have built general pairs of basis functions that do not suffer from such a loss in accuracy for fields that live on surfaces, and pairs of basis functions that can be used for the modelling of fields that live in three-dimensional domains.
For the accurate modelling of fields in strongly non-linear, radiating systems, we do not only need appropriate basis functions, but also the ability to compute the fields radiated by a single basis function. This needs to be done with an extremely high accuracy to avoid that errors creep into the solution that can render that solution completely useless. We have designed exact integration routines for these so-called interaction integrals. These integrals combine contributions across an abstract four dimensional space. Our computation method adds these together by reducing the calculation to one across a three dimensional space and subsequently to a two-dimensional space for which solution methods are available. This allows to, for the first time, compute the result at the highest accuracy possible on modern computers. This in turn leads to simulation methods that enjoy the highest stability.
The strongly non-linear and radiating systems that are being encountered in research into next-gen THz devices often contain highly specialized components such as sheets of graphene. These special two-dimensional components are the source of the non-linearity and it is important that we model them with great accuracy. We have coupled traditional field solvers with solution methods for non-linear equations to arrive at a method that can enforce, in each point of the graphene sheet, the complicated relation between currents and fields that is the result of the sophisticated physics taking place at the sheet. It is important to use the correct set of basis functions to describe all relevant physical quantities: in this case the currents in the sheet and the fields at the sheet. We have demonstrated that this can be done for surfaces without boundaries and without sharp corners by, at each time step, solving a large non-linear system that has as many unknowns as there are geometric elements in our finite description of the sheet.
In many cases systems are made up from a small number of materials and each material occupies a separate geometric region. Sometimes, however, the materials can be mixed up with varying proportions, leading to a gradient in the material properties. An important example is organic tissue, where different materials can occupy the same regions or where they are interlaced in hard to separate and hard to model geometric patterns. In this case modelling happens by assigning different material parameters to the different three-dimensional pixels or voxels that make up the system. It is important that the description of fields in these regions can be coupled or connected to the description in regions that are occupied a by a single material only. We have designed modelling methods that can include both types of regions. We have cast the resulting equations in a form that can be quickly solved by modern computers. To minimise the size of these systems, we have designed a new class of modelling methods, which we call multi-trace single source methods. Single source methods, as opposed to classic or dual source methods, express the field dynamics in half the number of unknowns, resulting in a smaller number of computations required for the modelling of the system.
Realistic systems are quite complicated and to model them accurately a large number of variables are required. The number of variables can easily reach several millions. Traditional algorithms require a number of computations that scales as the square of the number of variables: this means that when the number of variables doubles, the number of computations required is multiplied by four. For complicated devices this leads to computational workloads that are simply not feasible. Fortunately, in 1999 the group of Michielssen and collaborators developed the plane wave time domain algorithm, a direct time-domain counterpart of the famous fast multipole method, to reduce the number of computations needed to almost proportional to the number of variables. The method revolves around the single shot computation of the field observed in one location in the device as radiated at another location over an extended time interval. The larger the separation of the points, the longer the intervals are allowed to be and the more efficient the resulting algorithm. We have revisited the plane wave algorithm and come up with the optimal lengths for each such time interval. The final goal is to redesign this algorithm for execution on a network of computers, even further optimising the simulation time.
Time domain simulations are notorious for being unstable: errors that enter the simulation process are liable to magnify as they are being used as inputs to compute the solution at subsequent time steps. This error magnification goes on until they completely dominate any useful information in the output of the simulation. For perfectly conducting systems and single body systems, we have constructed methods that are guaranteed to be free from such instabilities. What is more, these modelling methods can be used when both rapidly and extremely slowly varying inputs are present and when the geometry under investigation contains both large features and small details.
This work is central to enable the modelling of strongly non-linear and radiating systems such as encountered in the design of next-gen THz technology. Indeed, these systems are made up of various materials, and contain both large radiating structures and extremely small semi-conductors and two-dimensional components.
We have begun research in simulators that do not directly reconstruct the electric and magnetic fields in the devices of interest, but first reconstruct the so-called magnetic vector potential. This quantity is arguably more fundamental then the electric and magnetic fields. In quantum mechanics and quantum field theory, photons (the fundamental particles making up electromagnetic radiation) are described in terms of this vector potential and some famous experiments (such as the Aharonov-Bohm effect) provide evidence that the slightly more abstract vector potential is the true source of electromagnetic interactions. For modellers, expressing fields in terms of the potential has a number of computational advances: the behaviour of the vector potential can be studied more easily when multiple frequencies are present. We have been able to show that boundary integral equations for very general geometries can be constructed for the vector potential.
When modelling either the electromagnetic field or the vector potential in a discrete, computational setting, so-called basis functions that describe the behaviour of these fields in a small region in space are required. In general, a pair of dual sets of basis functions is required, such that one set can be easily expressed in terms of the other. Without such a pair it is not possible to build modelling tools that provide answers in a reasonable number of computations. When the unknown field lives on a surface, such a pair is known but unfortunately will introduce excess error on one of the fields that is being computed. We have built general pairs of basis functions that do not suffer from such a loss in accuracy for fields that live on surfaces, and pairs of basis functions that can be used for the modelling of fields that live in three-dimensional domains.
For the accurate modelling of fields in strongly non-linear, radiating systems, we do not only need appropriate basis functions, but also the ability to compute the fields radiated by a single basis function. This needs to be done with an extremely high accuracy to avoid that errors creep into the solution that can render that solution completely useless. We have designed exact integration routines for these so-called interaction integrals. These integrals combine contributions across an abstract four dimensional space. Our computation method adds these together by reducing the calculation to one across a three dimensional space and subsequently to a two-dimensional space for which solution methods are available. This allows to, for the first time, compute the result at the highest accuracy possible on modern computers. This in turn leads to simulation methods that enjoy the highest stability.
The strongly non-linear and radiating systems that are being encountered in research into next-gen THz devices often contain highly specialized components such as sheets of graphene. These special two-dimensional components are the source of the non-linearity and it is important that we model them with great accuracy. We have coupled traditional field solvers with solution methods for non-linear equations to arrive at a method that can enforce, in each point of the graphene sheet, the complicated relation between currents and fields that is the result of the sophisticated physics taking place at the sheet. It is important to use the correct set of basis functions to describe all relevant physical quantities: in this case the currents in the sheet and the fields at the sheet. We have demonstrated that this can be done for surfaces without boundaries and without sharp corners by, at each time step, solving a large non-linear system that has as many unknowns as there are geometric elements in our finite description of the sheet.
In many cases systems are made up from a small number of materials and each material occupies a separate geometric region. Sometimes, however, the materials can be mixed up with varying proportions, leading to a gradient in the material properties. An important example is organic tissue, where different materials can occupy the same regions or where they are interlaced in hard to separate and hard to model geometric patterns. In this case modelling happens by assigning different material parameters to the different three-dimensional pixels or voxels that make up the system. It is important that the description of fields in these regions can be coupled or connected to the description in regions that are occupied a by a single material only. We have designed modelling methods that can include both types of regions. We have cast the resulting equations in a form that can be quickly solved by modern computers. To minimise the size of these systems, we have designed a new class of modelling methods, which we call multi-trace single source methods. Single source methods, as opposed to classic or dual source methods, express the field dynamics in half the number of unknowns, resulting in a smaller number of computations required for the modelling of the system.
Realistic systems are quite complicated and to model them accurately a large number of variables are required. The number of variables can easily reach several millions. Traditional algorithms require a number of computations that scales as the square of the number of variables: this means that when the number of variables doubles, the number of computations required is multiplied by four. For complicated devices this leads to computational workloads that are simply not feasible. Fortunately, in 1999 the group of Michielssen and collaborators developed the plane wave time domain algorithm, a direct time-domain counterpart of the famous fast multipole method, to reduce the number of computations needed to almost proportional to the number of variables. The method revolves around the single shot computation of the field observed in one location in the device as radiated at another location over an extended time interval. The larger the separation of the points, the longer the intervals are allowed to be and the more efficient the resulting algorithm. We have revisited the plane wave algorithm and come up with the optimal lengths for each such time interval. The final goal is to redesign this algorithm for execution on a network of computers, even further optimising the simulation time.
Progress beyond the state of the art and expected potential impact (including the socio-economic impact and the wider societal implications of the project so far)
In the next reporting period, we will focus on those topics that are still required for the fully integrated modelling of strongly non-linear and radiating systems such as those encountered in the design of next-gen devices in THz generation and detection.
We will focus on building time-domain methods that can deal with almost completely general geometric configurations and that can be used when both large features and small geometric details are present in the device. The frequency-domain results achieved in the first period provide a firm foundation to base this work on. Our experience in the design of stable time-domain methods for perfectly conducting surfaces and single body penetrable systems have learnt us a great deal about how the temporal variation should be described.
Our modelling method of two-dimensional systems such as sheets of graphene currently only works well for surfaces without boundaries and without corners. This is not how these components typically appear in the design of novel technology. In particular, this does not leave room for the connection to a driving circuit that will provide the energy required to power these devices. In the next reporting period, we will focus on using so-called higher order basis functions that can describe the modelling of currents and fields with a much higher accuracy. To make this feasible, we will compute the interactions between different parts of the sheet using the so-called convolution quadrature method, which relies on computing the interactions first at a judiciously chosen set of frequencies, after which the full temporal interaction can be reconstructed by combining these contributions. For the computation at any given frequency, existing methods that are extremely flexible and accurate are available, so we expect to make rapid progress in this work package.
To enable the connection of the device of interest to a driving circuit, novel descriptions of the scattering by sheets connected to circuits will be developed. The challenge is to leave out the driving circuit, but not to create any spurious accumulations of current and charge at the ports where the system has been truncated. Recently, techniques that rely on a regularisation near the boundary has been presented in the context of stabilisation in the presence of small geometric details. We will use similar techniques to build realistic port models without introducing the need to include a full-wave simulation of the driving circuit.
Even though the plane wave time domain method with the resection strategy we have developed allows for the computation of time dependent fields using not more memory than strictly required, the finite amount of memory available on any given computer still puts bounds on the complexity of the devices that we can model. In the next reporting period, we will focus on the parallelisation of this algorithm on modern distributed hardware. To enable extremely large simulations, we will design an algorithm that can continue its computation, even if one or multiple computers fail or when during the simulation, additional computers are added to the network. This robust parallel algorithm will be designed to work with acoustic fields, electromagnetic fields, and the electromagnetic vector potential.
To demonstrate the effectiveness of the developed ecosystem of modelling tools, we will pursue the fully integrated modelling of a source of THz radiation that is based on the strongly non-linear surface dynamics in twisted bilayer graphene superlattices. This means we will take into account the effect of the various materials present in the device, and the effect of the ports that connect the device to the driving circuit. Because of its complex structure, the modelling of the device will require millions of variables.. The device contains many small geometric features and because of non-linear effects, the solution can exhibit features on different length scales. We will use our robust parallel plane wave time domain method, in conjunction with stabilisation methods in the presence of slowly varying fields and small geometric features, to compute a stable solution using the absolute minimal computational resources.
We will focus on building time-domain methods that can deal with almost completely general geometric configurations and that can be used when both large features and small geometric details are present in the device. The frequency-domain results achieved in the first period provide a firm foundation to base this work on. Our experience in the design of stable time-domain methods for perfectly conducting surfaces and single body penetrable systems have learnt us a great deal about how the temporal variation should be described.
Our modelling method of two-dimensional systems such as sheets of graphene currently only works well for surfaces without boundaries and without corners. This is not how these components typically appear in the design of novel technology. In particular, this does not leave room for the connection to a driving circuit that will provide the energy required to power these devices. In the next reporting period, we will focus on using so-called higher order basis functions that can describe the modelling of currents and fields with a much higher accuracy. To make this feasible, we will compute the interactions between different parts of the sheet using the so-called convolution quadrature method, which relies on computing the interactions first at a judiciously chosen set of frequencies, after which the full temporal interaction can be reconstructed by combining these contributions. For the computation at any given frequency, existing methods that are extremely flexible and accurate are available, so we expect to make rapid progress in this work package.
To enable the connection of the device of interest to a driving circuit, novel descriptions of the scattering by sheets connected to circuits will be developed. The challenge is to leave out the driving circuit, but not to create any spurious accumulations of current and charge at the ports where the system has been truncated. Recently, techniques that rely on a regularisation near the boundary has been presented in the context of stabilisation in the presence of small geometric details. We will use similar techniques to build realistic port models without introducing the need to include a full-wave simulation of the driving circuit.
Even though the plane wave time domain method with the resection strategy we have developed allows for the computation of time dependent fields using not more memory than strictly required, the finite amount of memory available on any given computer still puts bounds on the complexity of the devices that we can model. In the next reporting period, we will focus on the parallelisation of this algorithm on modern distributed hardware. To enable extremely large simulations, we will design an algorithm that can continue its computation, even if one or multiple computers fail or when during the simulation, additional computers are added to the network. This robust parallel algorithm will be designed to work with acoustic fields, electromagnetic fields, and the electromagnetic vector potential.
To demonstrate the effectiveness of the developed ecosystem of modelling tools, we will pursue the fully integrated modelling of a source of THz radiation that is based on the strongly non-linear surface dynamics in twisted bilayer graphene superlattices. This means we will take into account the effect of the various materials present in the device, and the effect of the ports that connect the device to the driving circuit. Because of its complex structure, the modelling of the device will require millions of variables.. The device contains many small geometric features and because of non-linear effects, the solution can exhibit features on different length scales. We will use our robust parallel plane wave time domain method, in conjunction with stabilisation methods in the presence of slowly varying fields and small geometric features, to compute a stable solution using the absolute minimal computational resources.