FARR (Finite-difference time-domain ARRay) is a new high-performance 3D FDTD code designed for studying electromagnetic wave propagation in a magnetized collisional plasma. It is developed in C++ code optimized for studying computationally intensive problems. FARR is an open access project that is designed to be easily installed and used on a wide range of systems; from personal computers to supercomputers. To take advantage of highly parallel systems like next-generation supercomputers, FARR is domain decomposed and parallelized in three dimensions using Message Passing Interface (MPI). All I/O operations are handled using the ADIOS2 library (Godoy et al., 2020), and files are stored using either Binary Pack (BP) or HDF5 format. Figure 1 demonstrates strong scaling for FARR, running with up to 2048 processors on the Frontera system at the Texas Advanced Computing Center. Writing all vector quantities of E ⃗ and B ⃗ to disk can be time consuming, so custom I/O is available using the provided field monitor class, which defines a series of locations in the grid to save for later analysis.

FARR and its user manual can be easily accessed from its GitLab page at https://gitlab.com/longleywj/farr/. The fastest way to install FARR is using Conda to gather all requirements, which is detailed in the documentation. However, when using FARR for very large scale simulations, it can be beneficial to build against any prebuilt libraries/modules that are available. Post-processing of runs is handled through several easy and robust Python routines, including a lightweight class designed to interface with the output BP/HDF5 files. For complex 3D simulations, visualization through Paraview is available using either BP the files or a custom HDF5 loader.

FARR has a rich set of features that allow it to be configured to study a variety of different problems. The core FDTD algorithm that FARR uses is shown in Figure 2. In the rest of this section, we will provide a brief overview and accompanying examples of the key features that exist in FARR.

FARR currently has three different boundary conditions (BCs) implemented: reflecting, periodic, and a perfectly matched layer (PML) that absorbs all incident waves. We implement boundary conditions independently for each axis. For example, a user may wish to run with periodic BCs in the x and y directions and a PML in the z direction. This is easily set by flags in the input deck when running FARR.

Without explicit boundary conditions, the grid of an FDTD code is terminated by perfectly electrically conductive (PEC) boundary conditions that reflect all outwardly moving electromagnetic waves. This boundary condition is an unrealistic restriction in most cases and has pushed the development of more advanced absorbing boundary conditions which can dramatically lower the amount of energy reflected at the edges of the computational domain. The current state-of-the-art absorbing boundary condition is known as the perfectly matched layer (PML). The PML acts as a nonphysical lossy medium which exponentially decays outgoing waves.

In FARR we implement the convolutional, stretched-coordinate form of the PML, which is designed to be computationally efficient and independent of the material being simulated (Roden and Gedney, 2000). Figure 3 shows how the PML absorbs all incident waves from a point source at different times. This form of PML can handle inhomogenous, dispersive, anisotropic and nonlinear media. However, the perfectly matched layer has been shown to be unstable under certain conditions. When incident upon the PML, if the wave group and phase velocity contain anti parallel components, the wave will instead grow exponentially (Bécache et al., 2003). This situation can arise when simulating plane waves in a magnetized plasma (Chevalier et al., 2008), and it is partially for this reason that we choose to implement the Total-Field Scattered-Field formulation in FARR (Section 2.3).

In FDTD simulations it is often useful to make a distinction between the source wave field and any scattering, reflection, etc. that occurs as a result of the problem being simulated. The total-field scattered-field (TFSF) formulation is a boundary condition that is interior to and in addition to the BCs in Section 2.2. It is designed to separate the FDTD space into two distinct regions: the total-field and scattered-field regions, in which the normal FDTD algorithm is applied normally within each separate region (Merewether et al., 1980). At the boundary of these two regions, a correction term is included for the electric and magnetic field components, which will act as a source for the total field region. The wave source is analytically subtracted off at the TFSF boundary using the matched numerical dispersion technique of Guiffaut and Mahdjoubi (2000). The source wave will not enter the scattered field region unless it has been distorted/scattered in some way. In this way, we can implement arbitrary sources and distinguish them from the scattering processes that result from the simulation. For an in depth discussion of the total-field scattered formulation see Taflove and Hagnes (2005).

The TFSF is an integral part of FARR for several reasons, the most important of which is for introducing sources into the FDTD grid. Other methods for sources, such as hard sources, can cause spurious reflections and limit the practical uses of FDTD simulations (Taflove and Hagnes, 2005). FARR is designed to study ionospheric and radio science problems where we can assume the region of interest is in the far-field of any sources, and therefore we need to be able to model plane wave sources accurately. This creates an issue as plane wave sources are not directly compatible with the absorbing PML as the edges of the wavefront will drag through the PML and distort the wavefront or even cause a shift in the phase direction leading to the instability discussed in Section 2.2. The TFSF acts as an internal boundary condition that removes the source field before terminating the grid with the perfectly matched layer. Figure 4 shows how these two boundary conditions are applied in succession.

FARR is designed to study electromagnetic wave propagation in collisional plasmas like Earth’s ionosphere. To implement plasma effects, we must solve the electron momentum equation at each time step for the current J ⃗ . In finite-difference form, the momentum is given by Equation 1: J ⃗ n + 1 / 2 − J ⃗ n − 1 / 2 Δ t c + ν e J ⃗ n − 1 / 2 = ϵ 0 ω p e 2 E ⃗ − ω ⃗ c e × J ⃗ n + 1 / 2 + J ⃗ n − 1 / 2 2 (1) where Δ t c is the time step required for stability of the plasma routine given the local plasma ( ω p e ) , collision ( ν e ) , and cyclotron frequencies ( ω c e ) . This equation can be applied for each relevant species, however, over the short timescales of FARR simulations, we assume only electron motion is important. To solve for J n + 1 / 2 , we use the method described by Samimi and Simpson (2015), which is fully explicit and computationally efficient. However, when Δ t c is substantially lower than the simulation time step Δ t (required to satisfy the Courant-Friedrichs-Lewy condition), the current must to be updated multiple times before use in the main FDTD algorithm. This case typically occurs when the collision rate is high, and therefore the collisional routine from Pokhrel et al. (2018) is also implemented. This routine requires solving matrix equations, but can update the current density with only one step. All matrix updates are solved using the linear algebra library Armadillo (Sanderson and Curtin, 2019).

When running simulations in FARR, the proper algorithm and time step are selected automatically according to the global plasma, collision and cyclotron frequencies. As mentioned in Section 2.2, when simulating magnetized plasmas, the perfectly matched layer boundary condition can become unstable. In this case, to prevent instability, we recommend simulating with either periodic or reflecting boundary conditions, grading the plasma density to 0 at the boundaries, or containing the plasma within the total-field scattered-field domain (see Section 2.3).

There are several options for running a simulation with a background plasma in FARR. A global background magnetic field, plasma frequency and collision frequency can be specified in the input deck. FARR can also read spatially varying collision and plasma frequencies from an external HDF5 file that aligns with either the TFSF region or the entire simulation domain. Optionally, for studying time scales over which temporal variation of the plasma becomes important, the HDF5 file can contain a series of time steps to update the background plasma over the simulation. The best way to create an external HDF5 file for FARR is through the customizable Python adapter that reads in plasma densities from the HDF5 file and transforms the density to FARR’s grid. An example is provided in the source code to generate an HDF5 file from the output of an Electrostatic Parallel Particle-In-Cell (EPPIC) (Oppenheim and Dimant, 2004) simulation. This routine can be readily modified to read in plasma densities from other plasma simulators.

FARR includes a time domain near-to-far field transformation (NTFF) as part of the Python post-processing suite. The NTFF is based on the Surface Equivalence Principle, which holds that if the electric current, J ⃗ , and magnetic current, M ⃗ , are known everywhere on a closed surface, and there are no sources external to the surface, the complex amplitude anywhere outside the surface can be determined. Typical use cases for the near to far field transformation in FDTD include determining far field radiation patterns of antennas and measuring scattering cross sections.

In practice, this allows one to define a surface inside the FARR simulation and use it to compute the time varying electric and magnetic fields anywhere outside the computational domain. We use the method of Luebbers et al. (1991) to determine the time domain electromagnetic response at an arbitrary location in space. For the NTFF to be valid and work in FARR, all components of the electric and magnetic field on a closed surface must be output, and the surface should be defined to encompass all sources and scatterers. To perform the NTFF, we first define the closed surface to pass through E z nodes in the grid. We then linearly interpolate E ⃗ and H ⃗ to the surface, and determine the surface electric ( j s ⃗ = n ̂ × H s ⃗ ) and magnetic ( m s ⃗ = − n ̂ × E s ⃗ ) currents, where n ̂ is the surface outward normal on surface s . Integration and scaling of the surface currents (see Luebbers et al., 1991) yields the appropriate far field response. In Figure 5 we simulate a half-wave dipole antenna transmitting a 15 MHz sine wave. The electric and magnetic fields are post-processed using the NTFF routine to determine the far field radiation pattern 60 km outside of the computational domain, which is in good agreement with the theoretical pattern.

When the regions of interest for an FDTD simulation are in the near field, or when spatially varying waveforms are required, the total-field scattered-field can no longer be used as a source. Under these conditions we introduce two additional sources, individual antenna elements and antenna arrays. Elements and arrays are applied to the simulation domain with individual locations, phases and gains, giving the user complete control over the transmission pattern. Individual antenna elements are created within FARR using the thin-wire approximation (Umashankar et al., 1987), which allows for computationally efficient and accurate modeling of sub-cell features. This type of source is particularly useful when studying high frequency (HF) radio propagation in the bottomside ionosphere.

To demonstrate the use of antenna elements and phased arrays in FARR, we simulate arrays of half-wave dipole antennas operating at 10 MHz. The arrays have 2, 4 and 8 colinear elements, and are phased to have the main lobe at 90 ° (vertical) or 45 ° (oblique) elevation angles. The far field radiation patterns are determined using the near to far field transform (Section 2.5) and the results are shown in Figure 6. FARR includes Python routines to aid in setting the location and phasing of any number of antenna elements. When simulating arrays, all antenna elements use the same transmission with user specified phase offsets. For complete control over phase, location, transmission, etc., individual antenna elements should be used.

A key feature of FARR is the ability propagate waves through an inhomogeneous, magnetized plasma. To validate the plasma routines in 2.4, we run a set of simulations with a slowly varying plasma density and compare the results to the dispersion relations for O, X, R, and L mode waves. We perform these simulations in both collisionless and weakly collisional regimes to validate all routines. Table 1 shows the simulation parameters for FARR. The plasma density is set to increase linearly from 0 to n max ≈ 1.24 × 1 0 12 m − 3 over 2,000 m. The amplitude of an electromagnetic wave in a slowly varying plasma is described by the WKB approximation (Bellan, 2006). The WKB solution is a complex exponential, where the amplitude and phase change slowly as the wave propagates: E 1 x = 1 k x x exp i ∫ 0 x k x x ′ d x ′ (2)

In Equation 2, the amplitude and phase are determined by the wavenumber in the direction of the inhomogeneity; chosen as the x-direction here. For an O-mode wave, the wavenumber will vary spatially as the plasma frequency ω p e ( x ) varies, via the dispersion relation given in Equation 3: k x 2 x = ω 2 c 2 1 − ω p e 2 x ω 2 − k y 2 − k z 2 (3)

The simulations will propagate waves in the x-direction, so k y = k z = 0 . Therefore, as the wave propagates through the plasma, its amplitude at each position is a direct validation of the dispersion relation. We use the standard cold plasma dispersion relations for the R and L modes, given by Equations 4, 5 (Stix, 1992): R = k 2 c 2 ω 2 = 1 − ω p e 2 ω ω + Ω c e (4) L = k 2 c 2 ω 2 = 1 − ω p e 2 ω ω − Ω c e (5) where Ω c e carries a minus sign for electrons. The X mode dispersion relation is defined as Equation 6: k 2 c 2 ω 2 = 2 R L R + L (6)

The results of the WKB tests are shown in Figure 7. Because this is a time-domain code, truly monochromatic waves cannot be modeled. Instead, the input wave packet is a pulse of 100 wavelengths, with a Tukey window function applied. As the code models the full sinusoidal wave form, there is no information about the envelope of the wave. Therefore, the amplitude of the wave at each position is estimated by averaging the 500 highest values of | E | 2 over time. This estimated amplitude is then normalized by the initial amplitude of the wave at x = 0 , where n e = 0 . The results in Figure 7 are in general agreement with the WKB theory, therefore validating the dispersion relation for each propagation mode. The discrepancies between the simulations and the theory are due to the following:

1) Small oscillations in the simulation amplitudes will occur due to the finite bandwidth of the wave form. This is analogous to Gibbs phenomenon when trying to propagate a square wave in a fluid code.

As part of the validation of FARR, we compare simulations to the numerical ray tracing code PHaRLAP (Cervera and Harris, 2014) under a background ionosphere provided by IRI 2020 (Bilitza et al., 2022). Given the smoothly varying plasma density with only large scale gradients, the geometric optics approximation is applicable, and both FARR and PHaRLAP should provide similar results. The background plasma density from IRI begins at an altitude of approximately 65 km, with an foE = 1.15 MHz and hoE = 108 km. The altitude at which the plasma frequency matches 1 MHz is 96.75 km.

Using PHaRLAP we simulate a 1 MHz ray launched at a 60 ° elevation angle, and consider the height of reflection, time of flight, angle of arrival, and 1st hop distance. Here we use the 2D numerical ray tracing module in PHaRLAP which neglects the background magnetic field. We provide the same background plasma density to FARR and discretize the grid using 3,000 processors on the Stampede 3 supercomputer at the Texas Advanced Computing Center. The simulation parameters are provided in Table 2. We then launch a 1 MHz sine wave pulse using a 16 element linear array of half-wave dipole antennas, phased to 60 ° with a 7.16 ° half power beamwidth (HPBW). The source amplitude is modulated by a Gaussian envelope to provide a 100 kHz bandwidth. For comparison with PHaRLAP, we select the location in the grid with the maximum energy density at each timestep to be the center of the main lobe. Figure 8 compares the path taken by the main lobe in FARR with rays launched in PHaRLAP. With PHaRLAP, rays are launched spanning the HPBW of the antenna array. The location of the wave in FARR is in good agreement with the locations of PHaRLAP rays, with the exception of a small cluster of points at 100 km ground distance that are caused by a sidelobe of the antenna array. Table 3 summarizes the measured parameters, with the 1 hop distance for PHaRLAP being the ray launched at center of the main lobe. The general agreement between the simulations shows that FARR accurately reproduces important characteristics of ray-tracing when the geometric optics limit is valid.