Methods for Controlling Electrostatic Discharge and Electromagnetic Interference in Materials

Abstract

Methods for controlling electromagnetic fields in materials are presented that mitigate effects such as electrostatic discharge and electromagnetic/radio frequency interference. The first method determines the effective response of composite materials using a d-dimensional effective medium theory. The material consists of inhomogeneous two-layer inclusions with hyperspherical geometry. Non-integer dimensions represent fractal limits. The material medium is composed of a low hypervolume fraction of inclusions that are randomly distributed inside it. The effective response of the dielectric function is obtained using a virial expansion of the Maxwell–Garnett theory. The other method uses the transformation medium theory and involves the transformation of the material’s permittivity and permeability tensors so that the material exhibits a predefined effective response. By selecting appropriate transformations, a homogeneous material medium is transformed into an inhomogeneous version, forcing the electromagnetic fields to propagate along geodesic paths. These geodesics determine the behaviour of the fields inside the material. As a result, the material can be made to exhibit similar physical characteristics as those of a material composed of hyperspherical inclusions. The theoretical analysis presented is further studied and validated via the use of full-wave numerical simulations of Maxwell’s equations.

1. Introduction

Electrostatic discharge (ESD) and electromagnetic/radio frequency interference (EMI/RFI) are two critical phenomena that can have significant impact on various materials, particularly in the context of electronics, sensitive apparatus, and sensors. ESD is concerned with the sudden flow of electricity between two electrically charged objects caused by contact, an electrical short, or dielectric breakdown. EMI/RFI refers to the disturbance generated by an external source that affects an electrical circuit through electromagnetic induction, electrostatic coupling, or conduction. The work presented in this paper is motivated by the need to address these electromagnetic effects in space materials, but the results are also generic and apply to materials in general. Materials used in space platforms, in particular, have to be robust and resilient enough to cope with the harsh environment of space. Interaction with the surrounding plasma causes the accumulation of charged particles that induce an electric field in the materials. This results in such problems as electrostatic discharge (ESD) and electromagnetic/radio frequency interference (EMI/RFI) to occur. Electromagnetic discharge, in particular, can easily destroy the operation of electronic components within and around material structures [1,2]. In particular, ESD accounts for more than 50% of space platform anomalies and malfunctions [3]. Issues arising from EMI/RFI can come about via ESD, arcing, and other transient EM fields [4]. Studies into ESD and EMI/RFI have been conducted by many so that such processes are better understood and mitigated [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27].

Section 2 of the paper introduces some of the basic ideas that are explored in order to mitigate ESD and EM/RF interference. In Section 3, the effective medium theory (EMT) [28,29,30,31,32,33,34] is modified to a general d-dimensional form consisting of hyperinclusion geometries. In fact, a particular hypergeometry will be studied, which is that of hyperspheres because of their geometric symmetry [35]. The geometry of the hyperspheres is inhomogeneous because it is composed of two layers as depicted in Figure 1. It is shown, among other things, that parameters such as the dielectric function of a material can be controlled and mathematically determined via a virial expansion of the Maxwell–Garnett two-phase effective medium theory [36]. The EMT approach is written in d dimensions and to the first order in the hypervolume fraction of inclusions c.

In Section 4, an alternative method is used based on the transformation medium theory (TMT) [37,38,39,40,41,42,43,44]. The method allows the control of electromagnetic fields inside materials via a predefined form for the tensor parameters ${ϵ}^{ij}$ and ${\mu }^{ij}$ representing the permittivity and permeability, respectively. These describe how a material geometry can be designed with inhomogeneous regions inside it for controlling the path of fields. These field paths are in fact geodesics according to Fermat’s principle. The two alternative approaches complement each other, and they produce similar control of the electromagnetic fields so that problems such as ESD and other effects can be mitigated. Finally, Section 5 presents conclusions.

The theory developed and presented in the paper is examined via a full-wave numerical simulation of Maxwell’s equations using the Comsol Multiphysics [45] and FEKO [46] software.

2. Electrostatic Discharge and EM/RF Interference in Materials

One of the problems that has to be addressed in space material design is the control of electromagnetic (static) fields in order to mitigate electrostatic discharge (ESD) that causes arcing currents to form. By addressing ESD, the associated problem of electromagnetic and radio frequency interference (EMI/RFI) can be dealt with as well. In particular, arcing, caused by a charge buildup on space materials, can be understood using a relatively simple model. Suppose a negative charge builds up on a surface due to electrons. Consider another surface which now has a positive charge to counter the negative charge on the first surface. In simplistic terms, this is just a capacitor. If the region between the two surfaces is empty space, an electric field is induced between the surfaces, ${\mathbf{E}}_{app}$, due to the charges. The field direction is from the positive charges towards the negative charges. The empty space or vacuum between the surfaces has a dielectric constant or permittivity equal to that of free space, namely ${ϵ}_0$. The surfaces are separated by a distance d. The electric field is related to the charge density or charge per unit area $\sigma$ and the voltage V between the two surfaces as follows:

$$\begin{array}{c}\hfill E={\displaystyle \frac{\sigma }{{ϵ}_0}}={\displaystyle \frac{V}{d}}\end{array}$$

(1)

In addition, the force between any two charges is given by

$$\begin{array}{c}\hfill F=-{\displaystyle \frac{1}{4\pi {ϵ}_0}}{\displaystyle \frac{q_1{q}_2}{r^2}}\end{array}$$

(2)

since $q_1=-q_2$. Furthermore, the capacitance between the surfaces takes the form:

$$\begin{array}{c}\hfill C={\displaystyle \frac{{ϵ}_{0}A}{d}}\end{array}$$

(3)

The electric field, according to (1), builds up due to the charges, and if it passes a certain threshold it will cause arcing that may destroy electric/electronic components inside a structure. The requirement is to avoid this scenario or to radically reduce the electric field between the surfaces in order to eliminate or vastly reduce the risk of this happening. One approach that can facilitate this is to introduce a material between the two charged surfaces that has a very high dielectric constant $ϵ$. A material that has a high dielectric constant reduces the electric field since:

$$\begin{array}{c}\hfill {\mathbf{E}}_{eff}={\mathbf{E}}_{app}-{\mathbf{E}}_{pol}={\displaystyle \frac{\sigma }{ϵ{ϵ}_0}}\end{array}$$

(4)

Hence, the higher the value for $ϵ$, the smaller the effective electric field becomes. Such a material medium between two charged surfaces is made up of ‘particles’ which are polarised by the applied field ${\mathbf{E}}_{app}$ created between the two surfaces. The polarised particles then induce an opposing field ${\mathbf{E}}_{pol}$ so that the overall effective field is given by (4). As a result there is a reduction in the probability of arcing. In addition, it means that the force between the charges is reduced by a factor of $1/ϵ$ as well, i.e.,

$$\begin{array}{c}\hfill F=-{\displaystyle \frac{1}{4\pi ϵ{ϵ}_0}}{\displaystyle \frac{q_1{q}_2}{r^2}}\end{array}$$

(5)

Other material parameters are affected too. For example, the dielectric constant also depends on the temperature of a surface in an inverse way:

$$\begin{array}{c}\hfill ϵ\propto {\displaystyle \frac{1}{T}}\end{array}$$

(6)

where T is the temperature. When the temperature is high in a space material (facing the sun), the dielectric constant of the material goes down, and conversely, if the temperature is low (facing away from the sun), the dielectric of the material goes up. One problem that is addressed in this paper is the control and elimination of electromagnetic interference (EMI) and radio frequency interference (EFI). Transient fields in a material can interfere with other electronic components or instruments. It is necessary to mitigate these and many other effects via the manipulation of the permittivity/permeability parameters or effective material parameters such as the dielectric constant.

For example, to maximise the dielectric constant of a material in order to reduce or eliminate electrostatic discharge (ESD) effects, it is a matter of using a material that has high dielectric behaviour, but unfortunately, that is not so simple. Such a material may not be strong or robust enough to survive the space environment and other physical effects. An alternative would be required that has a high dielectric constant and at the same time has the ability to withstand the rigours of the space environment.

Two methods for achieving this are presented in this paper. The first is via the use of composite materials that contain very low concentrations of inhomogeneous hyperinclusion geometries. This means that the host material in which they are embedded in can be very strong and robust in order to cope with the space environment. A small concentration of the hyperinclusions is sufficient to prevent ESD amongst other things, while there is no compromise on the type of material that can be used for a particular structure. The d-dimensional framework is useful for examining the behaviour of fractal geometries which can be used in very narrow regions of a material structure such as corners where there is a larger concentration of charges. A $d=3$ hyperinclusion, i.e., a spherical inclusion, may not be a practical geometry for eliminating the electric field but a geometry that is in between a disc and a sphere might prove more useful, e.g., $d=2.5$. The second approach consists of a method that allows the design of materials such that their effective response is able to control the fields that exist inside them in some predefined way. This approach can be described as the transformation medium theory (TMT).

3. Effective Medium Theory

The effective medium theory (EMT) is a theoretical framework used in physics and materials science to describe the macroscopic properties of composite materials. These composite materials are heterogeneous materials composed of multiple phases or components, such as a mixture of different substances or materials. The primary idea behind the effective medium theory is to simplify the complex interactions between the individual components of a composite material by treating the material as a homogeneous medium with effective macroscopic properties. In other words, the EMT allows an approximation of the behaviour of a heterogeneous material as if it were a uniform medium with certain effective properties. The EMT is particularly useful when the microscopic details of the composite material are either unknown or too complicated to model accurately. By considering the average properties of the constituent materials and the geometry of the composite, the effective medium theory enables the prediction of various macroscopic physical properties, such as electrical conductivity, thermal conductivity, or optical properties of the composite material.

The effective medium theory is essentially the averaging of the interaction between each inclusion with the surrounding material they are embedded in. It takes into account all the interactions between the inclusions and the host material for a given concentration of inclusions c. Several models and approaches are employed in the EMT, each making different assumptions about the microstructural details of the composite material. Some widely used models include the Maxwell–Garnett model, the Bruggeman model, and the self-consistent approximation. These models provide valuable insights into the effective properties of composite materials, allowing researchers to make informed decisions about material design and optimisation for various applications. Before extending the Maxwell–Garnett (M-G) theory to a d-dimensional form for hyperspherical inclusions in a material, it is worth deriving the conventional M-G approach for a homogeneous dielectric sphere as a starting point.

Consider a homogeneous spherical inclusion inside a material experiencing an electric field $\mathbf{E}\left(\mathbf{r}\right)$. The electrostatic potential and the polarisability for such a configuration can be calculated using spherical polar coordinates $(r,\theta ,\varphi )$ because of the rotational symmetry. The electrostatic potential $\varphi$ is solved using the Laplace equation with no charges present:

$$\begin{array}{c}\hfill \Delta \varphi \left(\mathbf{r}\right)=0\end{array}$$

(7)

Using spherical coordinates means that Legendre polynomials $P_l(cos\theta )$ are required. The assumption is that there is no field excitation in the angular direction, so the $\varphi$ dependence is left out of the solutions. The overall potential is a superposition of the ‘excited’ potential and the potential of the spherical inclusion. The following ansatz for the potentials is made:

$$\begin{array}{c}\hfill {\varphi }_{in}(r,\theta )=\sum _{l=0}^{\infty }{A}_l{r}^l{P}_l(cos\theta )\end{array}$$

(8)

and

$$\begin{array}{c}\hfill {\varphi }_{out}(r,\theta )=-E_{0}r{P}_1(cos\theta )+\sum _{l=0}^{\infty }{\displaystyle \frac{C_l{P}_l(cos\theta )}{r^{l+1}}}\end{array}$$

(9)

Here, $A_l$ and $C_l$ are constants that have to be found by the boundary conditions at the surface of the sphere. Set $z=rcos\left(\theta \right)\equiv r{P}_1(cos\theta )$, and because $\mathbf{E}\left(\mathbf{r}\right)=E_0{\widehat{\mathbf{e}}}_z=-\nabla \varphi \left(\mathbf{r}\right)$, it leads to the form of the excited potential. The regions that are in and outside of the spherical inclusion are denoted ‘in’ and ‘out’, respectively. Since $\nabla \times \mathbf{E}\left(\mathbf{r}\right)=0$, the tangential part of the electric field is continuous, whereas $\nabla ·\mathbf{D}\left(\mathbf{r}\right)=0$ implies that the normal part of $\mathbf{D}\left(\mathbf{r}\right)$ is continuous. Here, $\mathbf{D}\left(\mathbf{r}\right)=ϵ\left(\mathbf{r}\right)\mathbf{E}\left(\mathbf{r}\right)$ is the constitutive relation which results in the continuity of $ϵ\left(\mathbf{r}\right)\mathbf{E}\left(\mathbf{r}\right)$. Using $\mathbf{E}\left(\mathbf{r}\right)=-\nabla \varphi \left(\mathbf{r}\right)$ gives:

$$\begin{array}{ccc}\hfill -{\displaystyle \frac{1}{r}}{\partial }_{\theta }{\varphi }_{in}(r=r_s,\theta )& =& -{\displaystyle \frac{1}{r}}{\partial }_{\theta }{\varphi }_{out}(r=r_s,\theta )\hfill \\ \hfill -{ϵ}_i{\partial }_r{\varphi }_{in}(r=r_s,\theta )& =& -{ϵ}_0{\partial }_r{\varphi }_{out}(r=r_s,\theta )\hfill \end{array}$$

(10)

where $r_s=a$ is the radius of the homogeneous sphere. Note also that ${\partial }_x\equiv \partial /\partial x$. Using the fact that

$$\begin{array}{c}\hfill P_l^1(cos\theta )={\displaystyle \frac{d}{d\theta }}{P}_l(cos\theta )\end{array}$$

(11)

gives

$$\begin{array}{ccc}\hfill {\partial }_{\theta }{\varphi }_{in}(r=r_s,\theta )& =& \sum _{l=1}^{\infty }{A}_l{r}_s^l{P}_l^1(cos\theta )\hfill \\ & =& -E_0{r}_s{P}_1^1(cos\theta )+\sum _{l=1}^{\infty }{\displaystyle \frac{C_l{P}_l^1(cos\theta )}{r_s^{l+1}}}\hfill \end{array}$$

(12)

The summation starts at $l=1$ because $P_0^1\left(x\right)=0$. From the orthogonality of the Legendre polynomials $P_l^m$, the following holds:

$$\begin{array}{ccc}\hfill A_1{r}_s& =& -E_0{r}_s+{\displaystyle \frac{C_1}{r_s^2}}\hfill \\ \hfill A_n{r}_s^{2n+1}& =& C_n\hfill \end{array}$$

(13)

for $n\ge 2$. Then, from the boundary conditions,

$$\begin{array}{ccc}\hfill {ϵ}_{in}{\partial }_r{\varphi }_{in}(r=r_s,\theta )& =& {ϵ}_{in}\sum _{l=1}^{\infty }l{A}_l{r}_s^{l-1}{P}_l(cos\theta )\hfill \\ & =& -{ϵ}_{out}\left[E_0{P}_1(cos\theta )+\sum _{l=0}^{\infty }{\displaystyle \frac{(l+1)C_l{P}_l(cos\theta )}{r_s^{l+2}}}\right]\hfill \end{array}$$

(14)

This implies that

$$\begin{array}{ccc}\hfill C_0& =& 0\hfill \\ \hfill {ϵ}_{in}{A}_1& =& -{ϵ}_{out}\left(E_0+2{\displaystyle \frac{C_1}{r_s^3}}\right)\hfill \\ \hfill {ϵ}_{in}n{A}_n{r}_s^{2n+1}& =& -{ϵ}_{out}(n+1)C_n\hfill \end{array}$$

(15)

for $n\ge 2$. The above system of equations can only hold if $A_n=C_n=0$∀$n\ne 1$. Finally, setting

$$\begin{array}{c}\hfill \left(\begin{array}{cc}{r}_s^3& -1\\ {\displaystyle \frac{{ϵ}_{in}}{{ϵ}_{out}}}{r}_s^3& 2\end{array}\right)\left(\begin{array}{c}{A}_1\\ C_1\end{array}\right)=\left(\begin{array}{c}-E_0{r}_s^3\\ -E_0{r}_s^3\end{array}\right)\end{array}$$

(16)

allows $A_1$ and $C_1$ to be solved so that

$$\begin{array}{c}\hfill A_1=-{\displaystyle \frac{3E_0}{2+\frac{{ϵ}_{in}}{{ϵ}_{out}}}}\end{array}$$

(17)

and

$$\begin{array}{c}\hfill C_1=r_s^3{E}_0{\displaystyle \frac{{ϵ}_{in}-{ϵ}_{out}}{{ϵ}_{in}+2{ϵ}_{out}}}\end{array}$$

(18)

The potential can now be written as:

$$\begin{array}{c}\hfill {\varphi }_{in}(r,\theta )=-{\displaystyle \frac{3E_0}{2+\frac{{ϵ}_{in}}{{ϵ}_{out}}}}r{P}_1(cos\theta )\end{array}$$

(19)

and

$$\begin{array}{c}\hfill {\varphi }_{out}(r,\theta )=-E_{0}r{P}_1(cos\theta )+E_0\left({\displaystyle \frac{{ϵ}_{in}-{ϵ}_{out}}{{ϵ}_{in}+2{ϵ}_{out}}}\right){\displaystyle \frac{r_s^3{P}_1(cos\theta )}{r^2}}\end{array}$$

(20)

The polarisability of the spherical inclusion $\alpha$ is related to the dipole moment $\mathbf{p}=\alpha \mathbf{D}$. The electrostatic potential of a dipole is given by

$$\begin{array}{c}\hfill {\varphi }_{dip}\left(\mathbf{r}\right)={\displaystyle \frac{1}{4\pi {ϵ}_0{ϵ}_{out}}}{\displaystyle \frac{\mathbf{p}·\mathbf{r}}{r^3}}\end{array}$$

(21)

Comparing the dipole potential to the potential of the sphere ${\varphi }_{out}(r,\theta )$, and since $\mathbf{p}·\mathbf{r}=pr{P}_1(cos{\theta }^{\prime })$ with the angle ${\theta }^{\prime }$ being the angle between $\mathbf{p}$ and $\mathbf{r}$,

$$\begin{array}{c}\hfill E_0\left({\displaystyle \frac{{ϵ}_{in}-{ϵ}_{out}}{{ϵ}_{in}+2{ϵ}_{out}}}\right){\displaystyle \frac{r_s^3{P}_1(cos\theta )}{r^2}}={\displaystyle \frac{1}{4\pi {ϵ}_0{ϵ}_{out}}}{\displaystyle \frac{pr{P}_1(cos{\theta }^{\prime })}{r^3}}\end{array}$$

(22)

if $\theta ={\theta }^{\prime }$ and $\mathbf{p}=p{\widehat{\mathbf{e}}}_z$, this implies that

$$\begin{array}{c}\hfill \mathbf{p}=4\pi r_s^3{ϵ}_0{ϵ}_{out}{\mathbf{E}}_0\left({\displaystyle \frac{{ϵ}_{in}-{ϵ}_{out}}{{ϵ}_{in}+2{ϵ}_{out}}}\right)\end{array}$$

(23)

where

$$\begin{array}{c}\hfill \alpha =4\pi r_s^3\left({\displaystyle \frac{{ϵ}_{in}-{ϵ}_{out}}{{ϵ}_{in}+2{ϵ}_{out}}}\right)\end{array}$$

(24)

is the polarisability of the spherical inclusion. Finally, from the Clausius–Mossotti equation that relates the effective medium permittivity or dielectric to the polarisability, namely,

$$\begin{array}{c}\hfill {\displaystyle \frac{ϵ-{ϵ}_{out}}{ϵ+2{ϵ}_{out}}}={\displaystyle \frac{\alpha }{3V}}\end{array}$$

(25)

where V is the volume of the medium (material) and using (24) gives:

$$\begin{array}{c}\hfill {\displaystyle \frac{ϵ-{ϵ}_{out}}{ϵ+2{ϵ}_{out}}}={\displaystyle \frac{V_s}{V}}\left({\displaystyle \frac{{ϵ}_{in}-{ϵ}_{out}}{{ϵ}_{in}+2{ϵ}_{out}}}\right)\end{array}$$

(26)

Let $c=V_s/V$ be the volume fraction of inclusions, ${ϵ}_{out}\equiv {ϵ}_0$ and ${ϵ}_{in}\equiv {ϵ}_1$, then the final form becomes:

$$\begin{array}{c}\hfill {\displaystyle \frac{ϵ-{ϵ}_0}{ϵ+2{ϵ}_0}}=c\left({\displaystyle \frac{{ϵ}_1-{ϵ}_0}{{ϵ}_1+2{ϵ}_0}}\right)\end{array}$$

(27)

Equation (27) is the Maxwell–Garnett effective medium model describing the interaction between spherical inclusions and the host material medium they are embedded in.

3.1. Effective Medium Theory of Inhomogeneous Hyperspheres

The EMT approach considered is based on extending the Maxwell–Garnett formulation to a d-dimensional framework for spherical inclusions; see (27). The non-integer dimensions of d can also be thought of as being fractal dimensions. The Maxwell–Garnett theory for a two-phase medium can now be written in d-dimensional form as:

$$\begin{array}{c}\hfill {\displaystyle \frac{ϵ-{ϵ}_0}{ϵ+(d-1){ϵ}_0}}=c\left({\displaystyle \frac{{ϵ}_1-{ϵ}_0}{{ϵ}_1+(d-1){ϵ}_0}}\right)\end{array}$$

(28)

Here, $ϵ$ is the dielectric function of the material that contains homogeneous hyperspheres of permittivity ${ϵ}_1$, and ${ϵ}_0$ is the permittivity of the surrounding material. See Figure 1 where the radii are equal, i.e., $a_1=a_2$ for the homogeneous case. Solving (28) for the case where the hypersphere hypervolume fraction of inclusions c is small yields:

$$\begin{array}{c}\hfill {\displaystyle \frac{ϵ}{{ϵ}_0}}=1+d\left({\displaystyle \frac{{ϵ}_1-{ϵ}_0}{{ϵ}_1+(d-1){ϵ}_0}}\right)c+d{\left({\displaystyle \frac{{ϵ}_1-{ϵ}_0}{{ϵ}_1+(d-1){ϵ}_0}}\right)}^2{c}^2+O\left(c^3\right)\end{array}$$

(29)

Effective medium theory or mean-field approximations are one-body interaction theories that account for the interaction of inclusions with the surrounding medium (material). This is represented by the linear term to order $O\left(c\right)$ quite accurately, provided that the random distribution of the inclusions inside the medium is sparse. This typically means values of $c=0.0-0.4$. As a consequence, any interactions between the hyperinclusions themselves is negligible and can be ignored. For greater values of c, the hyperinclusions are relatively close and higher-order interactions between them do matter. The coefficient of $O\left(c^2\right)$ represents two-body interactions, while $O\left(c^3\right)$ represents three-body interactions, and so on. The two-body interaction term in (29) for homogeneous hyperspherical inclusions is invalid, and more complicated theories that incorporate the correct interactions are required. However, for low concentrations of hyperspheres, mean-field approximations like the EMT are considerably accurate, and evaluating the linear term in (29) is sufficient for obtaining the material parameters such as the dielectric function, for example. According to (29), the following virial expansion for the dielectric function can be considered:

$$\begin{array}{c}\hfill ϵ={ϵ}_0(1+d{\gamma }_{0}c+O\left(c^2\right))\end{array}$$

(30)

For a random distribution of hyperspheres in a material medium, the total dipole moment induced by an electric field ${\mathbf{E}}_{\mathbf{0}}$ is given by:

$$\begin{array}{c}\hfill {\mathbf{p}}_{tot}={\gamma }_0{\mathbf{E}}_0\end{array}$$

(31)

where ${\gamma }_0$ is the effective polarisability that has the dimensional form:

$$\begin{array}{c}\hfill {\gamma }_0={\displaystyle \frac{{ϵ}_1-{ϵ}_0}{{ϵ}_1+(d-1){ϵ}_0}}\end{array}$$

(32)

For example, when $d=3$, the effective polarisability corresponding to spherical inclusions is given by:

$$\begin{array}{c}\hfill {\gamma }_0={\displaystyle \frac{{ϵ}_1-{ϵ}_0}{{ϵ}_1+2{ϵ}_0}}\end{array}$$

(33)

For the general case of inhomogeneous hyperspheres, the total d-dimensional dipole moment has a more complicated form compared to (32). It is derived to be:

$$\begin{array}{c}\hfill \gamma ={\displaystyle \frac{{\chi }_1+{\chi }_2{\beta }^d\left[1+(d-2){\chi }_1\right]}{1+(d-1){\chi }_1{\chi }_2{\beta }^d}}\end{array}$$

(34)

where the parameter $\beta$ is defined as the ratio of the inner over the outer radius of the inhomogeneous hyperinclusion, i.e., $\beta =a_1/a_2$. In addition, let ${\gamma }_1={ϵ}_2/{ϵ}_0$ and ${\gamma }_2={ϵ}_1/{ϵ}_2$; refer to Figure 1. The parameters ${\chi }_1$ and ${\chi }_2$ are then defined to be:

$$\begin{array}{c}\hfill {\chi }_1={\displaystyle \frac{{\gamma }_1-1}{{\gamma }_1+d-1}}\end{array}$$

(35)

and

$$\begin{array}{c}\hfill {\chi }_2={\displaystyle \frac{{\gamma }_2-1}{{\gamma }_2+d-1}}\end{array}$$

(36)

respectively. The dielectric function of the medium in now obtained via a virial expansion similar to (30):

$$\begin{array}{c}\hfill {\displaystyle \frac{ϵ}{{ϵ}_0}}=1+{\displaystyle \frac{d{\chi }_1+d{\chi }_2{\beta }^d\left[1+(d-2){\chi }_1\right]}{1+(d-1){\chi }_1{\chi }_2{\beta }^d}}c+O\left(c^2\right)\end{array}$$

(37)

The ability to control the behaviour of electromagnetic fields inside a material can be achieved via the randomly distributed inhomogeneous hyperspheres and their constitutive parameters. On the basis of (37), a general equation for predicting these parameters is given by:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{ij}={\lambda }_{ij}\left[1+(d-1){\chi }_1{\chi }_2{\beta }^d\right]-{\chi }_1\left[1+(d-2){\chi }_2{\beta }^d\right]-{\chi }_2{\beta }^d\end{array}$$

(38)

Equation (38) can be solved for particular limits as given by the pseudo-polarisation matrix ${\lambda }_{ij}$ with the condition that

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{ij}=0\end{array}$$

(39)

The limits of (38) are given by the matrix:

$$\begin{array}{c}\hfill {\lambda }_{ij}=\left[\begin{array}{cc}1& -{\displaystyle \frac{1}{d-1}}\\ \infty & -{\displaystyle \frac{1}{dc}}\\ 0& p\end{array}\right]\end{array}$$

(40)

Here, ${\lambda }_{11}=1$ is the perfect conducting limit, ${\lambda }_{12}=-{\displaystyle \frac{1}{d-1}}$ is the ‘holes’ limit, ${\lambda }_{21}=\infty$ maximises the dielectric function, also known as the Frohlich condition, ${\lambda }_{22}=-{\displaystyle \frac{1}{dc}}$ minimises the dielectric function, ${\lambda }_{31}=0$ eliminates the scattering of the hyperspheres or cloaks them from the fields, and ${\lambda }_{32}=p$ is any general value for the polarisation other than the limits that are already contained in the polarisation matrix. In what follows, three important limits are examined because of their importance in mitigating such things as electrostatic discharge (ESD), electromagnetic interference (EMI), and other effects in space materials.

3.2. The Limit ${\lambda }_{21}\to \infty$: Maximising the Dielectric Function

A material with hyperspherical inclusions can facilitate a large dielectric function which can help reduce electrostatic discharge (ESD) and arcing in materials. The dielectric of a material, also known as the relative permittivity, determines its ability to store electrical energy in an electric field. A material with a high dielectric function can effectively dissipate charges by creating an electric field that allows the charges to redistribute and neutralise. This helps prevent the buildup of static charges on the material’s surface, reducing the likelihood of ESD events. A dielectric function can confine and concentrate electric fields within a material. This confinement prevents field intensification at sharp edges or points, where arcing is more likely to occur. By confining the electric fields, the material can mitigate the risk of arcing and subsequent damage. A material with a high dielectric function can reduce the surface potential, which is the voltage difference between the material’s surface and its surroundings. Lower surface potentials minimise the likelihood of electrostatic discharge, as the potential difference required to initiate a discharge event is reduced. In addition, a material possessing a large dielectric function can act as an insulator, preventing the flow of electric current and reducing the risk of arcing. By effectively insulating conductive components or surfaces, the material can help maintain electrical isolation and prevent unintended electrical paths. When impedance matching is crucial, a material with a high dielectric function can facilitate the matching of impedance between different components or interfaces. This can help minimise reflections and ensure an efficient transfer of electrical energy, reducing the potential for ESD events.

On this basis, there is a need to obtain material structures that contain inclusions with a very high dielectric response. Using hyperspherical inclusions in a material, this effect can be studied by considering the limit ${\lambda }_{21}\to \infty$ given by (40). The idea is to increase the effective dielectric function or dielectric constant of the material via the manipulation of the polarisation, which is governed by the interactions between the hyperinclusions and the material they are embedded in. Taking the limit ${\lambda }_{21}\to \infty$ in (38) gives:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{21}\equiv 0\end{array}$$

(41)

Hence, (38) now becomes,

$$\begin{array}{c}\hfill -{\displaystyle \frac{1}{(d-1){\chi }_1{\chi }_2}}={\left({\displaystyle \frac{a_1}{a_2}}\right)}^d\end{array}$$

(42)

Setting the definitions

$$\begin{array}{c}\hfill \beta ={\displaystyle \frac{a_1}{a_2}},{\chi }_1={\displaystyle \frac{{\gamma }_1-1}{{\gamma }_1+d-1}},{\gamma }_1={\displaystyle \frac{{ϵ}_2}{{ϵ}_0}},{\chi }_2={\displaystyle \frac{{\gamma }_2-1}{{\gamma }_2+d-1}},{\gamma }_2={\displaystyle \frac{{ϵ}_1}{{ϵ}_2}}\end{array}$$

(43)

gives the following expression:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\gamma }_2-1}{{\gamma }_2+d-1}}=-{\displaystyle \frac{{\gamma }_1+d-1}{(d-1)({\gamma }_1-1){\beta }^d}}\end{array}$$

(44)

The requirement is to find a permittivity ${ϵ}_1$ or ${ϵ}_2$ (${\mu }_1$ or ${\mu }_2$) for a hyperinclusion as shown in Figure 1 that makes the polarisability, and hence the dielectric function as given by (37), approach infinity for a given hypervolume fraction of inclusions c. Suppose that ${ϵ}_1$ is the required constitutive parameter required in order to achieve this condition for any arbitrary value ${ϵ}_2$. Then, from (44):

$$\begin{array}{c}\hfill {ϵ}_1=(d-1){ϵ}_2\left[{\displaystyle \frac{({ϵ}_2-{ϵ}_0){\beta }^d-({ϵ}_2+(d-1){ϵ}_0)}{(d-1)({ϵ}_2-{ϵ}_0){\beta }^d+{ϵ}_2+(d-1){ϵ}_0}}\right]-\delta ϵ\end{array}$$

(45)

For any given value of ${ϵ}_2$ the corresponding constitutive parameter ${ϵ}_1$ that gives a polarisation and hence a dielectric function that tends to infinity is given by (45). Since the expression for ${ϵ}_1$ as given by (45) makes the dielectric function of the material go to infinity very quickly, it means that it is a mathematical singularity. By subtracting a small value of choice $\delta ϵ$ from the singular value, one can tune the value of ${ϵ}_1$ so that it can take large or small values. When $\delta ϵ=0$ the dielectric function of the material is infinite and for other values, for example $\delta ϵ=0.001$ the dielectric function of the material is very large but otherwise not infinite. For parameters as given in Figure 1, the variation in the constitutive permittivity ${ϵ}_1$ as a function of ${ϵ}_2$ is shown in Figure 2a for integer dimensional hyperinclusions that achieve very large values for the dielectric constant. Even larger values can be obtained when $\delta ϵ<<0.001$.

However, the values required for ${ϵ}_1$ to achieve very large effective dielectric constants for materials must be negative. This is shown in Figure 2b for the case where the inclusions are spheres $(d=3)$. As ${ϵ}_2$ increases in positive values, the corresponding ${ϵ}_1$ values are negative. Negative values imply the need to use a resonating material or structure that gives a negative permittivity ${ϵ}_1$. One way to produce negative values is through the use of metamaterials. In order to test the validity of the theory presented in this section, full-wave simulations were performed using the Comsol Multiphysics and FEKO software, which solve Maxwell’s equations numerically. The idea was to see what happened to a hyperinclusion inside a material experiencing an electromagnetic field. The effect on one inclusion can explain the effective response of the entire material consisting of multiple inclusions. The large effective dielectric response of the material was achieved via the elimination of the field inside the inclusions. An electromagnetic ‘hole’ was created within the inclusions where the field was zero or close to zero.

A composite material was modelled as shown in Figure 3. At a frequency of $f=0.8$ GHz, an electromagnetic wave was considered that propagated from left to right as shown in Figure 3a. In both Figure 3a,b an electromagnetic ‘hole’ was obtained where the field did not penetrate inside a spherical inclusion (2D cross section) but left the inside unaffected by the EM waves. This validates the notion that regions can be made to reduce large fields to zero or close to zero via the use of hyperinclusions inside the material. As a result, in such regions, the probability of electrostatic discharge is minimised, amongst other things.

3.3. The Limit ${\lambda }_{11}\to 1$: The Perfect Conducting Limit

Perfect conduction should be understood to mean in the limit as there is no material that possesses such a unique property. The properties of a material with perfect conducting hyperspherical inclusions are studied in this section. It is shown that such a material can reduce electrostatic discharge and arcing effects as well. The presence of perfectly conducting hyperspherical inclusions in a material can provide a path for the dissipation of electric charges. When an electrostatic charge builds up on the surface of the material, the conducting hyperinclusions can effectively distribute and dissipate the charge throughout the material, reducing the likelihood of a sudden discharge or arcing. By allowing the charges to spread out and reach equilibrium, the conducting inclusions help prevent the buildup of high potential differences that can lead to ESD or arcing. The conductivity of the hyperinclusions plays a crucial role in determining the effectiveness of the material in reducing ESD and arcing. For optimal performance, the inclusions should have a high electrical conductivity to efficiently distribute and dissipate the electric charges.

The distribution of conducting hyperinclusions within the composite material is important. A uniform dispersion throughout the material helps ensure a more consistent conductivity and charge dissipation in the material. Clustering or agglomeration of hyperinclusions should be minimised to maintain the desired electrical properties. This is another reason why in this paper, the hypervolume fraction of hyperinclusions c was taken to be $c=0.1$ or $10\%$ of the entire volume of a material. The material, which surrounds the conducting hyperinclusions, also plays a role in the overall performance of the material. The requirement is to affect its properties so that the material exhibits perfect conduction in order to mitigate ESD and other unwanted effects. This can be achieved by using hyperinclusions inside the material without needing a major alteration to its mechanical strength, chemical stability, and many other properties in order to maintain the integrity and functionality it was designed for.

In the limit ${\lambda }_{11}\to 1$, the hyperinclusions can be used inside a material medium to cancel the electric field. Using (38), the following equation needs to be solved by equating it to zero:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{11}={\lambda }_{11}\left[1+(d-1){\chi }_1{\chi }_2{\beta }^d\right]-{\chi }_1\left[1+(d-2){\chi }_2{\beta }^d\right]-{\chi }_2{\beta }^d\end{array}$$

(46)

Solving (46) gives the following expression:

$$\begin{array}{c}\hfill {\chi }_2={\beta }^{-d}\end{array}$$

(47)

From the definition of ${\chi }_2$, ${\gamma }_2$, and so on, the perfect conducting case is achieved when the constitutive parameter ${ϵ}_1$ is related to ${ϵ}_2$ (see Figure 1) via

$$\begin{array}{c}\hfill {ϵ}_1={ϵ}_2\left[{\displaystyle \frac{{\beta }^d+d-1}{{\beta }^d-1}}\right]\end{array}$$

(48)

for some ratio of the radii $\beta =a_1/a_2$. Thus, for a material to possess effective perfect conducting properties, the inhomogeneous ‘core’ of the hyperinclusions varies with $\beta$ for some fixed value of ${ϵ}_2$. An interesting limit is obtained when $a_1\to a_2$ or $\beta \to 1$. This corresponds to the case where the hyperinclusions become homogeneous and the permittivity ${ϵ}_1$ corresponds to the entire hyperinclusion. It also means that ${ϵ}_2\equiv {ϵ}_0$ of the surrounding material. Then, ${ϵ}_1$ tends to infinity:

$$\begin{array}{ccc}\hfill {ϵ}_1& =& \underset{\beta \to 1}{lim}{ϵ}_0\left[{\displaystyle \frac{{\beta }^d+d-1}{{\beta }^d-1}}\right]\hfill \\ & =& \infty \hfill \end{array}$$

(49)

For perfect conducting hyperinclusions, substitute (47) in (34), and after cancelling terms, it can be shown that the polarisation becomes $\gamma =1$. Substituting (47) in (37) also gives the effective dielectric constant for a perfect conducting material containing hyperinclusions as:

$$\begin{array}{c}\hfill ϵ={ϵ}_0+d{ϵ}_{0}c+O\left(c^2\right)\end{array}$$

(50)

For a material that has a perfect conducting response via its inclusions, the accumulation of charge in a region is suppressed and the electric field is zero. Numerical computations were performed to obtain the potential and electric field of a positive charge in a medium using

$$\begin{array}{c}\hfill \mathbf{E}=-\nabla \varphi \end{array}$$

(51)

In essence it is a matter of computing the equipotential hypersurfaces or in 2D, curves, which the electric field lines intersect perpendicularly. In curvilinear coordinates, the potential field $\varphi$ has to be obtained from the equation:

$$\begin{array}{ccc}\hfill {\nabla }^2\varphi & =& g^{ij}{\displaystyle \frac{{\partial }^2\varphi }{\partial {\xi }^i\partial {\xi }^j}}-g^{ij}{\mathrm{\Gamma }}_{ij}^k{\displaystyle \frac{\partial \varphi }{\partial {\xi }^k}}\hfill \\ & =& 0\hfill \end{array}$$

(52)

The coordinates $(i,j)$ can take on values beyond $d=3$ dimensions. Thus, (52) can be written in the compact form:

$$\begin{array}{ccc}\hfill \Delta \varphi & =& {\displaystyle \frac{1}{\sqrt{\left|g\right|}}}{\displaystyle \frac{\partial }{\partial {\xi }^i}}\left(\sqrt{\left|g\right|}{g}^{ij}{\displaystyle \frac{\partial \varphi }{\partial {\xi }^j}}\right)\hfill \\ & =& 0\hfill \end{array}$$

(53)

Here, $g=det\left(g_{ij}\right)$, $g_{ij}$ is the Euclidean metric tensor relative to the new coordinates, and ${\mathrm{\Gamma }}_{ij}^k$ are its Christoffel symbols. Figure 4 and Figure 5 show $d=2$ homogeneous inclusions in a material with a charge. The inclusions are able to affect the electric field and potential induced by the charge. As the number of inclusions increases, a larger region is formed where both the potential and electric field are zero. Such configurations would be useful for reducing or eliminating ESD and other effects from regions inside a material.

3.4. The Limit ${\lambda }_{31}\to 0$: EMI/RFI Cancellation

In this section, the limit where scattering cancellation is achieved via the hyperspherical inclusions is considered. The aim is to shield against electromagnetic interference (EMI) and radio frequency interference (RFI) propagating inside a material. It is necessary to attenuate and absorb electromagnetic waves, reducing the impact of external electromagnetic fields on sensitive electronic components, e.g., magnetometers. This dual functionality of ESD protection and EMI/RFI shielding is critical for space materials. The limit corresponding to scattering cancellation is examined using hyperspherical inclusions in a material. In addition to the electromagnetic properties studied here, the same approach is valid for the manipulation of the thermal properties of materials. The idea is to enhance or eliminate the thermal conductivity of a material, facilitating the efficient control of heat distribution. This is particularly beneficial in applications where heat generation is a concern. Effective thermal management helps maintain the material’s performance and reliability, further enhancing its suitability for ESD and arcing prevention as well.

Using (38), the equation required for such a case is obtained from:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{31}={\lambda }_{31}\left[1+(d-1){\chi }_1{\chi }_2{\beta }^d\right]-{\chi }_1\left[1+(d-2){\chi }_2{\beta }^d\right]-{\chi }_2{\beta }^d\end{array}$$

(54)

for the limit ${\lambda }_{31}\to 0$ and ${\mathrm{\Omega }}_{31}=0$. From (54),

$$\begin{array}{c}\hfill {\chi }_2=-{\displaystyle \frac{{\chi }_1{\beta }^{-d}}{1+(d-2){\chi }_1}}\end{array}$$

(55)

Using the definition for ${\chi }_1$ and ${\chi }_2$, respectively, a relation between the permittivity ${ϵ}_1$ (or permeability ${\mu }_1$) of the inner layer of a hyperinclusion with respect to the permittivity ${ϵ}_2$ (or permeability ${\mu }_2$) of the outer layer can be determined that cancel transient fields. Using (55):

$$\begin{array}{c}\hfill {ϵ}_1={ϵ}_2\left({\displaystyle \frac{[{ϵ}_0+(d-1){ϵ}_2]{\beta }^d-(d-1)({ϵ}_2-{ϵ}_0)}{[{ϵ}_0+(d-1){ϵ}_2]{\beta }^d+{ϵ}_2-{ϵ}_0}}\right)\end{array}$$

(56)

A plot of (56) is shown in Figure 6 for different integer inclusion dimensions as a function of ${ϵ}_2$. Observe that for EM/RF interference cancellation, the constitutive parameters $({ϵ}_1,{ϵ}_2)$ are always positive. Figure 7 shows a full-wave numerical simulation using the FEKO software of a $d=3$ spherical inclusion surrounded by vacuum or air with relative permittivity of ${ϵ}_0=1$ and ${ϵ}_0\approx 1$, respectively. In vacuum or air, the value of ${ϵ}_1$ is negative; however, when the surrounding medium is ${ϵ}_0>1$, the values for ${ϵ}_1$ are always positive (see Figure 6). The perturbed transient field propagating from the left towards a homogeneous spherical inclusion in vacuum, shown in Figure 7a, is ‘smoothed-out’ when the inclusion is an inhomogeneous two-layer sphere, see Figure 7b. In fact, a more realistic scenario consists of transient fields inside a material medium. The fields can be eliminated by an array of two-layer inhomogeneous spherical inclusions as shown in Figure 8. Normal spherical inclusions fail to eliminate the transient fields, as they travel from left to right and contribute to their perturbation, see Figure 8a, but using inhomogeneous spherical inclusions with parameters as predicted via (56), the perturbed transient fields propagating to the right are eliminated; see Figure 8b.

4. Transformation Medium Theory

Previously, it was shown that the electromagnetic properties of a material could be manipulated by embedding inhomogeneous-layered hyperinclusions inside it. The effective properties of the material were then studied using the effective medium theory (EMT). In this section, a second approach is considered for controlling fields inside a material. It consists of ‘bending’ or ‘curving’ the electromagnetic fields in a region of interest so that they are made to behave in certain ways or even vanish. Instead of achieving this by inserting hyperinclusions in a material as before, the material itself is ‘designed’ to act as the control mechanism that changes electromagnetic behaviour inside it. This can only occur if the material contains regions that are inhomogeneous with respect to its permittivity and permeability. The mathematical properties of differential manifolds can be used to model these effects since the material can be thought of as being a three-dimensional manifold. In such materials (manifolds), the electromagnetic fields travel in paths that are not straight lines, but rather, the paths are curved due to the inhomogeneities inside the material. According to Fermat’s principle, these paths are the shortest distances between any two or more points in the material and are referred to as geodesics. In essence, material properties must be transformed from those of a homogeneous or flat material to those of an inhomogeneous one. This process can be referred to as the transformation medium theory (TMT) and can be used to calculate the physical properties of the material precisely.

The process involves the transformation of Maxwell’s equations from Cartesian coordinates, representing a flat or homogeneous material, to other coordinate systems of choice that represent the non-homogeneous version of the same material. Indeed, Maxwell’s equations are almost always invariant under coordinate transformations which means that the electromagnetic field parameters are generally the same, regardless of the coordinate system. However, two of the parameters appearing in Maxwell’s equations in materials do not remain invariant and must be transformed accordingly because of their strong dependence on the spatial variation in the material. These two parameters are the second-rank contravariant tensors for the permittivity ${ϵ}^{ij}$ and permeability ${\mu }^{ij}$, respectively. They are both directly connected to the refractive index of the material medium. These second-rank tensors, which describe how a material’s permittivity and permeability must transform in order to manipulate electromagnetic fields, rely on the metric tensor. The paths that fields travel along inside the inhomogeneous material must be optimal in the sense that they are the shortest paths or geodesics. The equations that describe the motion of the fields along these geodesics must be determined using Riemann differential geometry and tensor calculus. These equations can then be used to describe the bulk properties of the material. The mathematical equations that allow the permittivity/permeability tensors to be calculated, which control or manipulate electromagnetic fields inside materials, are now discussed and derived. In principle, one is then able to design and engineer a material to control EM fields as desired for practical applications.

4.1. Determining the Equations and Geodesics of Electromagnetic Fields in a Material

A process that determines the geodesics in a material medium is derived via the calculus of variations. From this, the equations of motion of the fields inside the material can also be obtained. These equations turn out to be complicated higher-order non-homogeneous differential equations which cannot be solved in closed form. Their solution usually requires a numerical computation. The equations of motion are derived first, and then, the geodesics that the fields propagate along inside the material are calculated. Consider the action A for which an extremum must be found with respect to the material coordinates ${\xi }^{m^{\prime }}=({\xi }^{1^{\prime }},{\xi }^{2^{\prime }},...,{\xi }^{M^{\prime }})=(x,y,z,...)$ and their gradients ${\dot{\xi }}^{m^{\prime }}$. Usually a ‘dot’ above the coordinate means differentiation with respect to time but here, it means differentiation as a function of an arbitrary parameter $\lambda$, i.e., ${\dot{\xi }}^{m^{\prime }}\equiv d{\xi }^{m^{\prime }}/d\lambda$, which could also be time. Then, the action that must be an extremum has the form:

$$\begin{array}{c}\hfill A={\int }_{{\lambda }_1}^{{\lambda }_2}L({\xi }^{m^{\prime }}\left(\lambda \right),{\dot{\xi }}^{m^{\prime }}\left(\lambda \right))d\lambda \end{array}$$

(57)

The coordinates in (57) represent many possible field paths that can be minimised, but there is always, in principle, one such path between the coordinate boundaries or beginning and end points of all paths such that ${\xi }^{m^{\prime }}\left({\lambda }_1\right)=0$ and ${\xi }^{m^{\prime }}\left({\lambda }_2\right)=0$. The extremum has to be the minimum if the shortest electromagnetic path between coordinates is required. To find this minimum path, the coordinates are perturbed from the non-optimised original path to that of an arbitrary one ${\eta }^m\left(\lambda \right)$ by an amount $ϵ$. That is, the action coordinates are perturbed by the transformations ${\xi }^{m^{\prime }}={\xi }^m+ϵ{\eta }^m$ and ${\dot{\xi }}^{m^{\prime }}={\dot{\xi }}^m+ϵ{\dot{\eta }}^m$, respectively, representing nearby paths. The minimum path is determined when the limit $ϵ\to 0$ is taken. The action then becomes:

$$\begin{array}{c}\hfill A={\int }_{{\lambda }_1}^{{\lambda }_2}L({\xi }^m+ϵ{\eta }^m,{\dot{\xi }}^m+ϵ{\dot{\eta }}^m)d\lambda \end{array}$$

(58)

where the dependence of ${\xi }^m$ and ${\dot{\xi }}^m$ on $\lambda$ has been dropped for brevity. In other words, the differentiation of the arguments of the integrand in (57) also depends on the perturbations (transformations) as shown in (58). Note that the coordinates ${\eta }^m$ and ${\dot{\eta }}^m$ are also functions of the arbitrary parameter $\lambda$. The perturbed paths also obey the conditions ${\eta }^m\left({\lambda }_1\right)=0$ and ${\eta }^m\left({\lambda }_2\right)=0$. Differentiating (57) with respect to $ϵ$ gives the extrema of the paths taken by electromagnetic fields inside a material medium,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dA}{dϵ}}& =& \underset{ϵ\to 0}{lim}{\displaystyle \frac{d}{dϵ}}{\int }_{{\lambda }_1}^{{\lambda }_2}L({\xi }^{m^{\prime }},{\dot{\xi }}^{m^{\prime }})d\lambda \hfill \\ & =& 0\hfill \end{array}$$

(59)

Thus, (59) becomes, as $ϵ\to 0$:

$$\begin{array}{ccc}\hfill \underset{ϵ\to 0}{lim}{\int }_{{\lambda }_1}^{{\lambda }_2}\left[{\displaystyle \frac{\partial L}{\partial {\xi }^{m^{\prime }}}}{\displaystyle \frac{\partial {\xi }^{m^{\prime }}}{\partial ϵ}}+{\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^{m^{\prime }}}}{\displaystyle \frac{\partial {\dot{\xi }}^{m^{\prime }}}{\partial ϵ}}\right]d\lambda & =& 0\to \hfill \\ \hfill \underset{ϵ\to 0}{lim}{\int }_{{\lambda }_1}^{{\lambda }_2}\left[{\displaystyle \frac{\partial L}{\partial {\xi }^{m^{\prime }}}}{\eta }^m+{\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^{m^{\prime }}}}{\dot{\eta }}^m\right]d\lambda & =& 0\hfill \end{array}$$

(60)

Equation (60) is called the weak form in the calculus of variations. Unfortunately, this form is not very useful for obtaining the minimum geodesic. However, the second term on the right can be integrated by parts to obtain:

$$\begin{array}{c}\hfill {\int }_{{\lambda }_1}^{{\lambda }_2}{\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^{m^{\prime }}}}{\dot{\eta }}^{m}d\lambda ={\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^{m^{\prime }}}}{\eta }^m\left(\lambda \right){|}_{{\lambda }_1}^{{\lambda }_2}-{\int }_{{\lambda }_1}^{{\lambda }_2}{\eta }^m{\displaystyle \frac{d}{d\lambda }}\left[{\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^{m^{\prime }}}}\right]d\lambda \end{array}$$

(61)

The first term on the right of (61) vanishes due to the boundary conditions ${\eta }^m\left({\lambda }_1\right)=0\equiv {\eta }^m\left({\lambda }_2\right)=0$ mentioned above. A substitution of (61) into () gives

$$\begin{array}{c}\hfill \underset{ϵ\to 0}{lim}{\int }_{{\lambda }_1}^{{\lambda }_2}{\eta }^m\left(\lambda \right)\left[{\displaystyle \frac{\partial L}{\partial {\xi }^{m^{\prime }}}}-{\displaystyle \frac{d}{d\lambda }}\left[{\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^{m^{\prime }}}}\right]\right]d\lambda =0\end{array}$$

(62)

Here, the arbitrary function ${\eta }^m\left(\lambda \right)$ has been factored out. In fact, because ${\eta }^m\left(\lambda \right)$ is arbitrary, it can be chosen to be non-zero. Then, the only alternative that makes the integrand in (62) vanish to zero is

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial L}{\partial {\xi }^m}}-{\displaystyle \frac{d}{d\lambda }}\left[{\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^m}}\right]=0\end{array}$$

(63)

where the limit $ϵ\to 0$ has been taken in the coordinate expressions ${\xi }^{m^{\prime }}$ and ${\dot{\xi }}^{m^{\prime }}$. Equation (63) is the equation that must be used in order to obtain the optimal field path corresponding to a minimum separation between coordinates inside an inhomogeneous material medium. This equation is known as the Euler–Lagrange equation.

It is now possible to use the Euler–Lagrange equation to find the equations of motion and geodesics of the fields in the material. The geodesic distance S between coordinates inside the material is obtained on the basis of the covariant metric tensor $g_{ij}$ which is discussed more later. Then, S is given by

$$\begin{array}{c}\hfill S={\int }_{s_1}^{s_2}ds\end{array}$$

(64)

where the infinitesimal distance along the geodesic, $ds$, is obtained by

$$\begin{array}{c}\hfill d{s}^2=g_{ij}\left({\xi }^m\right)d{\xi }^{i}d{\xi }^j\end{array}$$

(65)

Here, repeated upper and lower indices are summed according to the Einstein convention. The distance element (65) can be parametrised using an arbitrary coordinate $\lambda$. Then,

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{ds}{d\lambda }}\right)}^2=g_{ij}\left({\xi }^m\right){\displaystyle \frac{d{\xi }^i}{d\lambda }}{\displaystyle \frac{d{\xi }^j}{d\lambda }}\end{array}$$

(66)

Substituting (66) into (64) results in the expression:

$$\begin{array}{c}\hfill S={\int }_{s_1}^{s_2}\sqrt{g_{ij}\left({\xi }^m\right){\displaystyle \frac{d{\xi }^i}{d\lambda }}{\displaystyle \frac{d{\xi }^j}{d\lambda }}}d\lambda \end{array}$$

(67)

Equation (67) gives the EM field geodesics inside a material. However, there are numerous such paths, and what is required is a specific path that optimises (shortens) the distance between coordinates inside the material. To find the minimum path, the integrand of (67) must be treated as an extremum. This requires the use of the Euler–Lagrange equation derived above with

$$\begin{array}{c}\hfill L({\xi }^m\left(\lambda \right),\dot{{\xi }^m}\left(\lambda \right))=\sqrt{g_{ij}\left({\xi }^m\right){\displaystyle \frac{d{\xi }^i}{d\lambda }}{\displaystyle \frac{d{\xi }^j}{d\lambda }}}\end{array}$$

(68)

Using (68) and the second term in brackets of the Euler–Lagrange equation (63) gives the following:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^k}}& =& {\displaystyle \frac{\partial }{\partial {\dot{\xi }}^k}}\sqrt{g_{ij}{\dot{\xi }}^i{\dot{\xi }}^j}\hfill \\ & =& {\displaystyle \frac{1}{2\sqrt{g_{ij}{\dot{\xi }}^i{\dot{\xi }}^j}}}{\displaystyle \frac{\partial }{\partial {\dot{\xi }}^k}}\left[g_{ij}{\dot{\xi }}^i{\dot{\xi }}^j\right]\hfill \\ & =& {\displaystyle \frac{1}{2L}}\left(g_{ij}{\dot{\xi }}^j{\displaystyle \frac{\partial {\dot{\xi }}^i}{\partial {\dot{\xi }}^k}}+g_{ij}{\dot{\xi }}^i{\displaystyle \frac{\partial {\dot{\xi }}^j}{\partial {\dot{\xi }}^k}}\right)\hfill \\ & =& {\displaystyle \frac{1}{2L}}\left(g_{ij}{\delta }_k^i{\dot{\xi }}^j+g_{ij}{\delta }_k^j{\dot{\xi }}^i\right)\hfill \end{array}$$

(69)

where the argument of the metric tensor $g_{ij}\left({\xi }^m\right)$ has been dropped for convenience. From the Kronecker delta tensors above, it can be seen that the non-zero entries occur only whenever $j=k$ and $i=k$, respectively, so that ${\delta }_k^j\to {\delta }_k^k$ and ${\delta }_k^i\to {\delta }_k^k$:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^k}}={\displaystyle \frac{1}{2L}}\left[g_{kj}{\dot{\xi }}^j+g_{ik}{\dot{\xi }}^i\right]\end{array}$$

(70)

Since the summed indices in (70) are dummy variables, setting $j=i$ in the first term gives the symmetric component form:

$$\begin{array}{c}\hfill g_{ki}{\dot{\xi }}^i={\displaystyle \frac{1}{2}}\left[g_{ki}+g_{ik}\right]{\dot{\xi }}^i\end{array}$$

(71)

then, this gives

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^k}}& =& {\displaystyle \frac{1}{2L}}\left[g_{ki}{\dot{\xi }}^i+g_{ik}{\dot{\xi }}^i\right]\hfill \\ & =& {\displaystyle \frac{1}{2L}}\left[2g_{ki}{\dot{\xi }}^i\right]\hfill \\ & =& {\displaystyle \frac{1}{L}}{g}_{ki}{\displaystyle \frac{d{\xi }^i}{d\lambda }}\hfill \end{array}$$

(72)

In order to eliminate L from (72), use is made of (66) where $ds/d\lambda =L$. More specifically, if a function is given as $f=f\left(s\right(\lambda \left)\right)$, then it holds that:

$$\begin{array}{c}\hfill {\displaystyle \frac{df}{d\lambda }}={\displaystyle \frac{df}{ds}}{\displaystyle \frac{ds}{d\lambda }}\equiv L{\displaystyle \frac{df}{ds}}\end{array}$$

(73)

In operator form, this can be written as:

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{d\lambda }}=L{\displaystyle \frac{d}{ds}}\end{array}$$

(74)

Substitution into (72) gives the final form for the second term in the Euler–Lagrange equation:

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{d\lambda }}\left[{\displaystyle \frac{\partial L}{\partial {\dot{\xi }}^k}}\right]={\displaystyle \frac{d}{d\lambda }}\left[g_{ki}{\displaystyle \frac{d{\xi }^i}{ds}}\right]\end{array}$$

(75)

The next step is to evaluate the first term in the Euler–Lagrange equation, which becomes

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial L}{\partial {\xi }^k}}& =& {\displaystyle \frac{\partial }{\partial {\xi }^k}}\sqrt{g_{ij}{\dot{\xi }}^i{\dot{\xi }}^j}\hfill \\ & =& {\displaystyle \frac{1}{2L}}\left[{\displaystyle \frac{\partial }{\partial {\xi }^k}}{g}_{ij}{\displaystyle \frac{d{\xi }^i}{d\lambda }}{\displaystyle \frac{d{\xi }^j}{d\lambda }}\right]\hfill \end{array}$$

(76)

Substituting both (75) and (76) into the Euler–Lagrange equation gives

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{d\lambda }}\left[g_{ki}{\displaystyle \frac{d{\xi }^i}{ds}}\right]-{\displaystyle \frac{1}{2L}}{\displaystyle \frac{\partial g_{ij}}{\partial {\xi }^k}}{\displaystyle \frac{d{\xi }^i}{d\lambda }}{\displaystyle \frac{d{\xi }^j}{d\lambda }}=0\end{array}$$

(77)

The arbitrary parametrisation $d\lambda$ and L can be eliminated by choosing $\lambda$ to be along the geodesic path S so that $\lambda =s$. Hence, from (74)

$$\begin{array}{c}\hfill L{\displaystyle \frac{d}{ds}}\left[g_{ki}{\displaystyle \frac{d{\xi }^i}{ds}}\right]-{\displaystyle \frac{L}{2}}{\displaystyle \frac{\partial g_{ij}}{\partial {\xi }^k}}{\displaystyle \frac{d{\xi }^i}{ds}}{\displaystyle \frac{d{\xi }^j}{ds}}=0\end{array}$$

(78)

where it is important to remember that the argument of the metric tensor is dependent on the coordinates ${\xi }^m$, i.e., $g_{ij}\left({\xi }^m\right)$. The geodesic Equation (78) gives the equations for the minimum path between coordinates in a material medium and the parameter s is a parametrisation along the path itself. Alternatively, (78) can be expanded as follows:

$$\begin{array}{c}\hfill g_{ki}{\displaystyle \frac{d^2{\xi }^i}{d{s}^2}}+{\displaystyle \frac{\partial g_{ki}}{\partial {\xi }^m}}{\displaystyle \frac{d{\xi }^m}{ds}}{\displaystyle \frac{d{\xi }^i}{ds}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial g_{ij}}{\partial {\xi }^k}}{\displaystyle \frac{d{\xi }^i}{ds}}{\displaystyle \frac{d{\xi }^j}{ds}}=0\end{array}$$

(79)

Equation (79) can be represented in a more compact form involving the Christoffel symbols of the second kind. To transform it, let the index in the metric tensor be $i=j$ and $m=i$ in the second term, since there is an implied summation for both these indices and they can be changed arbitrarily. At the same time, the partial derivative in the second term can now be written in terms of its symmetric component form

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial g_{kj}}{\partial {\xi }^i}}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial g_{kj}}{\partial {\xi }^i}}+{\displaystyle \frac{\partial g_{ik}}{\partial {\xi }^j}}\right)\end{array}$$

(80)

to obtain:

$$\begin{array}{c}\hfill g_{ki}{\displaystyle \frac{d^2{\xi }^i}{d{s}^2}}+{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{\partial g_{kj}}{\partial {\xi }^i}}+{\displaystyle \frac{\partial g_{ik}}{\partial {\xi }^j}}-{\displaystyle \frac{\partial g_{ij}}{\partial {\xi }^k}}\right]{\displaystyle \frac{d{\xi }^i}{ds}}{\displaystyle \frac{d{\xi }^j}{ds}}=0\end{array}$$

(81)

Multiplying both sides of (81) with the rank-two contravariant form $g^{nk}$, noting that

$$\begin{array}{c}\hfill g_{ki}{g}^{nk}={\delta }_i^n\end{array}$$

(82)

is the Kronecker-delta tensor, then setting the index $i=n$ in (81), it becomes:

$$\begin{array}{c}\hfill {\displaystyle \frac{d^2{\xi }^n}{d{s}^2}}+{\displaystyle \frac{1}{2}}{g}^{nk}\left[{\displaystyle \frac{\partial g_{kj}}{\partial {\xi }^i}}+{\displaystyle \frac{\partial g_{ik}}{\partial {\xi }^j}}-{\displaystyle \frac{\partial g_{ij}}{\partial {\xi }^k}}\right]{\displaystyle \frac{d{\xi }^i}{ds}}{\displaystyle \frac{d{\xi }^j}{ds}}=0\end{array}$$

(83)

Identifying the Christoffel symbol

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{ij}^n={\displaystyle \frac{1}{2}}{g}^{nk}\left[{\displaystyle \frac{\partial g_{kj}}{\partial {\xi }^i}}+{\displaystyle \frac{\partial g_{ik}}{\partial {\xi }^j}}-{\displaystyle \frac{\partial g_{ij}}{\partial {\xi }^k}}\right]\end{array}$$

(84)

(81) can be written in the final form

$$\begin{array}{c}\hfill {\displaystyle \frac{d^2{\xi }^n}{d{s}^2}}+{\mathrm{\Gamma }}_{ij}^n{\displaystyle \frac{d{\xi }^i}{ds}}{\displaystyle \frac{d{\xi }^j}{ds}}=0\end{array}$$

(85)

which is a more useful form for finding the minimum field paths, i.e., the geodesics via the equations of motion of the fields inside the material. As mentioned before, the equations that must be solved in order to obtain the geodesics are non-linear due to the requirement that the material must be inhomogeneous. However, when the material is homogeneous, it is considered to be Euclidean or flat. Mathematically, the geodesic is a straight line, as expected for such a material (manifold). This is because for a Euclidean material, the symmetric metric tensor $g_{ij}$ in (84) is nothing more than the Kronecker-delta tensor:

$$\begin{array}{ccc}\hfill g_{ij}& =& {\delta }_{ij}\to \hfill \\ \hfill {\delta }_{ij}& =& diag(1,1,1,...,n)\hfill \end{array}$$

(86)

To see this, consider two arbitrary coordinates $(u,v)$ in a material M: ${\xi }^m=({\xi }^1,{\xi }^2)=(u,v)$. Then, from (85), after replacing the metric tensor $g_{ij}$ with the Euclidean metric ${\delta }_{ij}$, the derivatives with respect to this metric vanish in the Christoffel symbols, and the only term remaining is

$$\begin{array}{c}\hfill {\displaystyle \frac{d^2{\xi }^n}{d{s}^2}}=0\end{array}$$

(87)

Since there are two parameters or coordinates $u,v$, $n=1,2$. Set $n=1$ to obtain the equation of the first parameter and $n=2$ for the second. Following this process and using (87) gives two sets of equations,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d^2{\xi }^1}{d{s}^2}}& =& 0\hfill \\ \hfill {\displaystyle \frac{d^2{\xi }^2}{d{s}^2}}& =& 0\hfill \end{array}$$

(88)

These correspond to the two equations

$$\begin{array}{c}\hfill {\displaystyle \frac{d^{2}u}{d{s}^2}}=0\mathrm{and}{\displaystyle \frac{d^{2}v}{d{s}^2}}=0\end{array}$$

(89)

where ${\xi }^1=u$ and ${\xi }^2=v$. It is easy to show that the solutions of the differential equations in (89) are

$$\begin{array}{c}\hfill u\left(s\right)=c_{1}s+c_2\mathrm{and}v\left(s\right)=c_{3}s+c_4\end{array}$$

(90)

where $c_1,c_3$ are the gradients of straight lines with intercepts at $c_2,c_4$. Hence, the solutions (90) indicate that for a flat or Euclidean material, the field geodesic Equation (85) reduces to solutions for straight lines, i.e., the fields travel in straight paths. The field geodesic distance between two coordinates in a homogeneous material can be shown via (65) to be

$$\begin{array}{ccc}\hfill d{s}^2& =& {\delta }_{ij}d{\xi }^{i}d{\xi }^j\hfill \\ & =& {\delta }_{11}{\left(d{\xi }^1\right)}^2+2{\delta }_{12}d{\xi }^{1}d{\xi }^2+{\delta }_{22}{\left(d{\xi }^2\right)}^2\hfill \\ & =& d{u}^2+d{v}^2\hfill \end{array}$$

(91)

The geodesic distance becomes, after substituting $ds$ above into (64),

$$\begin{array}{c}\hfill S={\int }_{s_1}^{s_2}\sqrt{d{u}^2+d{v}^2}\end{array}$$

(92)

which can be re-written as a derivative in s:

$$\begin{array}{ccc}\hfill S& =& {\int }_{s_1}^{s_2}\sqrt{{\left({\displaystyle \frac{du}{ds}}\right)}^2+{\left({\displaystyle \frac{dv}{ds}}\right)}^2}ds\hfill \\ & =& {\int }_{s_1}^{s_2}\sqrt{c_1^2+c_3^2}ds\hfill \end{array}$$

(93)

Setting the constants and the square root to $m=\sqrt{c_1^2+c_3^2}$, the geodesic distance of a Euclidean material becomes

$$\begin{array}{c}\hfill S=m(s_2-s_1)\end{array}$$

(94)

as expected, where m can be considered as a kind of gradient that can be set to unity without loss of generality. Equation (94) indicates that in homogeneous (Euclidean or flat) materials, the shortest distance the EM fields travel along is a straight path and equal to the difference between the two points $s_1$ and $s_2$ inside the material. This is of course intuitive for a Euclidean material, but this is not the case for non-Euclidean geodesics corresponding to inhomogeneous materials. Inhomogeneity is required for manipulating the fields inside the material. Consider the same problem as above but this time for a non-Euclidean material. This requires the use of (85), where the coordinate or parameter space is the same as before, ${\xi }^m=({\xi }^1,{\xi }^2)=(u,v)$. As a result, the summation of the repeated indices goes from $i,j,n=1,2$. Setting $n=1$ means that a non-linear second-order differential equation corresponding to the first coordinate is obtained; the terms in (85) are summed, so the equation becomes:

$$\begin{array}{c}\hfill {\displaystyle \frac{d^{2}u}{d{s}^2}}+{\mathrm{\Gamma }}_{11}^1{\left({\displaystyle \frac{du}{ds}}\right)}^2+\left({\mathrm{\Gamma }}_{12}^1+{\mathrm{\Gamma }}_{21}^1\right){\displaystyle \frac{du}{ds}}{\displaystyle \frac{dv}{ds}}+{\mathrm{\Gamma }}_{22}^1{\left({\displaystyle \frac{dv}{ds}}\right)}^2=0\end{array}$$

(95)

The first Christoffel symbol in (95) is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{11}^1={\displaystyle \frac{1}{2}}{g}^{11}{\displaystyle \frac{\partial g_{11}}{\partial u}}+g^{12}{\displaystyle \frac{\partial g_{12}}{\partial u}}-{\displaystyle \frac{1}{2}}{g}^{12}{\displaystyle \frac{\partial g_{11}}{\partial v}}\end{array}$$

(96)

The next Christoffel symbol in (95) takes on the following form:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{12}^1={\displaystyle \frac{1}{2}}\left(g^{12}{\displaystyle \frac{\partial g_{22}}{\partial u}}+g^{11}{\displaystyle \frac{\partial g_{11}}{\partial v}}\right)\end{array}$$

(97)

The Christoffel symbols in the third term in (85) are symmetric, i.e., ${\mathrm{\Gamma }}_{12}^1\equiv {\mathrm{\Gamma }}_{21}^1$, so that:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{12}^1+{\mathrm{\Gamma }}_{21}^1=g^{12}{\displaystyle \frac{\partial g_{22}}{\partial u}}+g^{11}{\displaystyle \frac{\partial g_{11}}{\partial v}}\end{array}$$

(98)

The final Christoffel symbol turns out to be:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{22}^1={\displaystyle \frac{1}{2}}{g}^{12}{\displaystyle \frac{\partial g_{22}}{\partial v}}+g^{11}{\displaystyle \frac{\partial g_{12}}{\partial v}}-{\displaystyle \frac{1}{2}}{g}^{11}{\displaystyle \frac{\partial g_{22}}{\partial u}}\end{array}$$

(99)

Equations (95)–(99) can be substituted into (85) noting that $g_{ij}$ is the covariant form of the metric tensor, and its contravariant form is given by $g^{ij}$. Both the covariant and contravariant forms of the metric tensor are symmetric. Equation (95) is the non-linear higher-order partial differential equation for coordinate ${\xi }^1=u$. Similarly for coordinate ${\xi }^2=v$ and using $n=2$, the following equation is obtained:

$$\begin{array}{c}\hfill {\displaystyle \frac{d^{2}v}{d{s}^2}}+{\mathrm{\Gamma }}_{11}^2{\left({\displaystyle \frac{du}{ds}}\right)}^2+\left({\mathrm{\Gamma }}_{12}^2+{\mathrm{\Gamma }}_{21}^2\right){\displaystyle \frac{du}{ds}}{\displaystyle \frac{dv}{ds}}+{\mathrm{\Gamma }}_{22}^2{\left({\displaystyle \frac{dv}{ds}}\right)}^2=0\end{array}$$

(100)

As before, the Christoffel symbols can be obtained for (100) starting with:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{11}^2={\displaystyle \frac{1}{2}}{g}^{21}{\displaystyle \frac{\partial g_{11}}{\partial u}}+g^{22}{\displaystyle \frac{\partial g_{12}}{\partial u}}-{\displaystyle \frac{1}{2}}{g}^{22}{\displaystyle \frac{\partial g_{11}}{\partial v}}\end{array}$$

(101)

The third term gives the following Christoffel symbols:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{12}^2+{\mathrm{\Gamma }}_{21}^2=g^{21}{\displaystyle \frac{\partial g_{11}}{\partial v}}+g^{22}{\displaystyle \frac{\partial g_{22}}{\partial u}}\end{array}$$

(102)

The final Christoffel symbol is determined to be:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{22}^2={\displaystyle \frac{1}{2}}{g}^{22}{\displaystyle \frac{\partial g_{22}}{\partial v}}+g^{21}{\displaystyle \frac{\partial g_{12}}{\partial v}}-{\displaystyle \frac{1}{2}}{g}^{21}{\displaystyle \frac{\partial g_{22}}{\partial u}}\end{array}$$

(103)

Note that when the metric tensor is that for a material that is homogeneous and flat (Euclidean), namely (86), all the Christoffel symbols vanish and (95) and (100) are equivalent to (89) and hence (90), i.e., linear solutions as expected. Note also that two coordinates were considered above for a $2D$ material where $n=1,2$ and $i,j=1,2$. For a three-dimensional material, the number of coordinates increases by one, i.e., $n=1,2,3$ and $i,j=1,2,3$. The contravariant metric tensor can be obtained from the covariant tensor by the cofactor method:

$$\begin{array}{c}\hfill g^{ij}={\displaystyle \frac{{(-1)}^{i+j}cdet\left(g_{ji}\right)}{det\left(g_{ij}\right)}}\end{array}$$

(104)

where $det(·)$ is the normal determinant, and $cdet(·)$ is the cofactor determinant. Here, $g_{ji}$ is the transpose and the $jth$ row and $ith$ column are crossed out before obtaining the determinant of whatever rows and columns remain. Similarly, the covariant metric tensor can be obtained from the contravariant version by:

$$\begin{array}{c}\hfill g_{ij}={\displaystyle \frac{{(-1)}^{i+j}cdet\left(g^{ji}\right)}{det\left(g^{ij}\right)}}\end{array}$$

(105)

In three dimensions and for coordinates $(u,v,w)$, the geodesic distance between coordinates in a material is given by the integration of the element:

$$\begin{array}{c}\hfill d{s}^2=g_{11}d{u}^2+g_{22}d{v}^2+g_{33}d{w}^2+2g_{12}dudv+2g_{13}dudw+2g_{23}dvdw\end{array}$$

(106)

Equation (106) gives the geodesic path of electromagnetic fields through a $3D$ material and requires the explicit form for the geodesic tensor. This is defined and discussed in the next section.

4.2. Transformation of the Constitutive Parameters in Maxwell’s Equations

A material’s structure and properties can be mathematically equivalent to, and described by, a geometrical manifold. In addition, the properties of a material also obey Maxwell’s equations regardless of whether they are homogeneous or inhomogeneous. The electromagnetic equations are form-invariant and independent of the coordinate system they are used in to describe the material properties. Hence, Maxwell’s equations can be written as:

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\sqrt{g}}}{\left(\sqrt{g}{E}^i\right)}_{,i}={\displaystyle \frac{\rho }{{ϵ}_0}};{\displaystyle \frac{1}{\sqrt{g}}}{\left(\sqrt{g}{B}^i\right)}_{,i}=0\\ \hfill {ϵ}^{ijk}{E}_{k,j}=-{\displaystyle \frac{\partial B^i}{\partial t}};{ϵ}^{ijk}{B}_{k,j}={\displaystyle \frac{1}{c^2}}{\displaystyle \frac{\partial E^i}{\partial t}}+{\mu }_0{j}^i\end{array}$$

(107)

Here, g is the determinant of the metric tensor and is discussed in detail in the following. The alternating Levi-Civita tensor is given by ${ϵ}^{ijk}$ and $E_{k,j}\equiv \partial E_k/\partial {\xi }^j$. The contravariant tensors (107) can be written in terms of covariant tensors which resemble more closely Maxwell’s equations in dielectric media:

$$\begin{array}{c}\hfill {\left(\sqrt{g}{g}^{ij}{E}_j\right)}_{,i}={\displaystyle \frac{\sqrt{g}\rho }{{ϵ}_0}};{\left(\sqrt{g}{g}^{ij}{B}_j\right)}_{,i}=0\\ \hfill \left[ijk\right]E_{k,j}=-{\displaystyle \frac{\partial }{\partial t}}(\pm \sqrt{g}{g}^{ij}{B}_j);\left[ijk\right]B_{k,j}={\displaystyle \frac{1}{c^2}}{\displaystyle \frac{\partial }{\partial t}}(\pm \sqrt{g}{g}^{ij}{E}_j)+{\mu }_0\sqrt{g}{j}^i\end{array}$$

(108)

The Levi-Civita tensor is written as

$$\begin{array}{c}\hfill {ϵ}^{ijk}=\pm {\displaystyle \frac{\left[ijk\right]}{\sqrt{g}}}\end{array}$$

(109)

and the positive sign corresponds to the right-handed coordinate system, while the negative sign implies the left-handed system. The form of Maxwell’s equations given by (108) can be further defined as

$$\begin{array}{c}\hfill D_{,i}^i=\rho B_{,i}^i=0\\ \hfill \left[ijk\right]E_{k,j}=-{\displaystyle \frac{\partial B^i}{\partial t}};\left[ijk\right]H_{k,j}={\displaystyle \frac{\partial D^i}{\partial t}}+j^i\end{array}$$

(110)

If $B_i$ is replaced by $H_i/{\mu }_0$ by a rescaling of the charge and current densities, then the constitutive equations can be defined as:

$$\begin{array}{c}\hfill D^i={ϵ}_0{ϵ}^{ij}{E}_j;B^i={\mu }_0{\mu }^{ij}{H}_j\\ \hfill {ϵ}^{ij}={\mu }^{ij}=\pm \sqrt{g}{g}^{ij}\end{array}$$

(111)

As a consequence, the empty-space Maxwell’s equations in arbitrary coordinates and geometries turn out to be the same as Maxwell’s equations in a right-handed Cartesian coordinate system. Dielectric media are in fact the same as these arbitrary geometries, and the permittivity tensor ${ϵ}^{ij}$ is identical to the permeability tensor ${\mu }^{ij}$. If a medium (material) does not contain an inhomogeneous geometry, the Riemann tensor is zero, and the medium is flat or Euclidean. This means that fields are described by coordinate transformations in Cartesian space as shown below. The form invariance of Maxwell’s equations is thus preserved under such transformations except for the constitutive tensors ${ϵ}^{ij}$ and ${\mu }^{ij}$. It is these tensors, which describe the effect that a transformation has on a material medium, that govern how the electromagnetic fields propagate inside the medium in the real sense.

The two constitutive tensors are second-rank contravariant tensors, and they are both directly connected to the refractive index of the medium. They can be written as

$$\begin{array}{c}\hfill {ϵ}^{i^{\prime }{j}^{\prime }}={\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{\partial {\xi }^{i^{\prime }}}{\partial {\xi }^k}}{\displaystyle \frac{\partial {\xi }^{j^{\prime }}}{\partial {\xi }^l}}{g}^{kl}{ϵ}^{ij}\end{array}$$

(112)

and

$$\begin{array}{c}\hfill {\mu }^{i^{\prime }{j}^{\prime }}={\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{\partial {\xi }^{i^{\prime }}}{\partial {\xi }^k}}{\displaystyle \frac{\partial {\xi }^{j^{\prime }}}{\partial {\xi }^l}}{g}^{kl}{\mu }^{ij}\end{array}$$

(113)

respectively, where $g=det\left({ϵ}^{i^{\prime }{j}^{\prime }}\right)=det\left({\mu }^{i^{\prime }{j}^{\prime }}\right)$. The indices $(i,j)$ are summed over the three spatial dimensions. A summation is implied if the same index appears in covariant–contravariant form or vice versa as per the Einstein convention. The primes indicate the new coordinate system. It should be pointed out that these equations also hold for the case where the Riemann tensor vanishes and the medium is locally flat. For a medium (material) that has an inhomogeneity, the tensors become more involved. Using these tensors one can consider a coordinate transformation by a design that transforms the permittivity and permeability tensors (${ϵ}^{i^{\prime }{j}^{\prime }},{\mu }^{i^{\prime }{j}^{\prime }})$ in such a way that the fields propagate in a predefined way or completely vanish inside the material. From either (112) or (113), the metric tensor can be identified, which is:

$$\begin{array}{c}\hfill g^{i^{\prime }{j}^{\prime }}={\displaystyle \frac{\partial {\xi }^{i^{\prime }}}{\partial {\xi }^k}}{\displaystyle \frac{\partial {\xi }^{j^{\prime }}}{\partial {\xi }^l}}{g}^{kl}\end{array}$$

(114)

The contravariant permittivity and permeability tensors can be written in the compact form:

$$\begin{array}{c}\hfill {ϵ}^{i^{\prime }{j}^{\prime }}={\displaystyle \frac{1}{\sqrt{g}}}{g}^{i^{\prime }{j}^{\prime }}{ϵ}^{ij}\end{array}$$

(115)

and

$$\begin{array}{c}\hfill {\mu }^{i^{\prime }{j}^{\prime }}={\displaystyle \frac{1}{\sqrt{g}}}{g}^{i^{\prime }{j}^{\prime }}{\mu }^{ij}\end{array}$$

(116)

It should be understood that the coordinates chosen to transform the homogeneous Euclidean material to an inhomogeneous one can be arbitrary and according to the desired design specifications one desires. In what follows, some examples are considered in order to highlight the approach in $(x,y,z)$ Cartesian coordinates. However, in specific cases it may be easier to consider other coordinate systems for simplicity such as spherical polar coordinates, cylindrical coordinates, and so on. It is shown via the numerical solution of the field equations and geodesics that the material medium can be made to eliminate or change the fields inside so that electromagnetic effects such as ESD/arcing and EM/RF interference can be minimised or even eliminated, for example.

4.3. Controlling EM Fields in Different Regions of a Material

In this section, equations and their solutions are obtained that transform a material medium so that electromagnetic fields are manipulated in some pre-required manner. Let the three-dimensional coordinates take values $m=1,2,3$ inside the material structure. The coordinates then become ${\xi }^{m^{\prime }}=({\xi }^{1^{\prime }},{\xi }^{2^{\prime }},{\xi }^{3^{\prime }})\equiv (x^{\prime },y^{\prime },z^{\prime })$ and ${\xi }^m=({\xi }^1,{\xi }^2,{\xi }^3)\equiv (x,y,z)$. Suppose that a design feature for the material is to weaken or eliminate the field in a given region so that the ESD/arcing probability is reduced dramatically or even eliminated. One way to achieve this is to ‘compress’ the internal structure of the material from a homogeneous or flat structure to an inhomogeneous one, as shown in Figure 9. Consider the following coordinate transformations:

$$\begin{array}{ccc}\hfill x^{\prime }& =& x\hfill \\ \hfill y^{\prime }& =& y\left[1-e^{-\left({\displaystyle \frac{x^2}{{\sigma }_x^2}}+{\displaystyle \frac{y^2}{{\sigma }_y^2}}\right)}\right]\hfill \\ \hfill z^{\prime }& =& z\hfill \end{array}$$

(117)

Figure 9a shows a homogeneous material. In the language of differential geometry, this is merely a flat or Euclidean manifold. Conversely, Figure 9b shows the same medium but this time, it is ‘compressed’ so that it now exhibits an inhomogeneous behaviour from the point of view of the fields in the material. The mathematical transformations that achieve this are given by (117). On the other hand, it is possible to ‘expand’ the material medium via the transformations given by (118):

$$\begin{array}{ccc}\hfill x^{\prime }& =& x\hfill \\ \hfill y^{\prime }& =& y\left[1-sin\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)\right]\hfill \\ \hfill z^{\prime }& =& z\hfill \end{array}$$

(118)

Figure 10a,b show the effect of such a transformation for the parameters given.

The transformations (117) and (118) can be written as a contravariant rank = 1 tensor, i.e., a vector. For example, (118) becomes:

$$\begin{array}{ccc}\hfill {\xi }^{m^{\prime }}& =& ({\xi }^{1^{\prime }}({\xi }^1,{\xi }^2,{\xi }^3),{\xi }^{2^{\prime }}({\xi }^1,{\xi }^2,{\xi }^3),{\xi }^{3^{\prime }}({\xi }^1,{\xi }^2,{\xi }^3))\hfill \\ & =& (x^{\prime }(x,y,z),y^{\prime }(x,y,z),z^{\prime }(x,y,z))\hfill \\ & =& (x,y-ysin\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right),z)\hfill \end{array}$$

(119)

In ESD/arcing effects, for example, the electric component of the fields is dominant and the magnetic component can be ignored because it is a quasi-static limit. Either way, the mathematical procedure considered is the same for the permeability tensor as well, which describes the magnetic field. The permittivity tensor is now obtained using (112) for the transformations given by (118) (the same process can be used in the case of (117)). The components of the permittivity tensor are determined by:

$$\begin{array}{c}\hfill {ϵ}^{1^{\prime }{1}^{\prime }}={\displaystyle \frac{1}{\sqrt{g}}}{g}^{1^{\prime }{1}^{\prime }}{ϵ}^{11}\end{array}$$

(120)

and

$$\begin{array}{c}\hfill {ϵ}^{1^{\prime }{2}^{\prime }}={\displaystyle \frac{1}{\sqrt{g}}}{g}^{1^{\prime }{2}^{\prime }}{ϵ}^{12}\end{array}$$

(121)

and so on. Here, $1^{\prime }{1}^{\prime }\equiv x^{\prime }{x}^{\prime }$, $1^{\prime }{2}^{\prime }\equiv x^{\prime }{y}^{\prime }$, etc., represent the Cartesian coordinates. Therefore, it is paramount to calculate the metric tensor using (114), which enables the full derivation of the permittivity tensor components. Here, the unprimed contravariant metric tensor $g^{ij}$ represents the homogeneous medium, since the medium is Euclidean or flat to begin with before applying the transformations. This means that $g^{ij}={\delta }^{ij}$, with the latter being the Kronecker tensor. Thus,

$$\begin{array}{ccc}\hfill g^{1^{\prime }{1}^{\prime }}& =& {\displaystyle \frac{\partial {\xi }^{1^{\prime }}}{\partial {\xi }^k}}{\displaystyle \frac{\partial {\xi }^{1^{\prime }}}{\partial {\xi }^l}}{\delta }^{kl}\hfill \\ & =& {\left({\displaystyle \frac{\partial x^{\prime }}{\partial x}}\right)}^2+{\left({\displaystyle \frac{\partial x^{\prime }}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial x^{\prime }}{\partial z}}\right)}^2\hfill \end{array}$$

(122)

Using the transformation (118), the first $x^{\prime }{x}^{\prime }$-component of the metric tensor becomes:

$$\begin{array}{c}\hfill g^{x^{\prime }{x}^{\prime }}=1\end{array}$$

(123)

Similarly,

$$\begin{array}{ccc}\hfill g^{1^{\prime }{2}^{\prime }}& =& {\displaystyle \frac{\partial {\xi }^{1^{\prime }}}{\partial {\xi }^k}}{\displaystyle \frac{\partial {\xi }^{2^{\prime }}}{\partial {\xi }^l}}{\delta }^{kl}\hfill \\ & =& {\displaystyle \frac{\partial x^{\prime }}{\partial x}}{\displaystyle \frac{\partial y^{\prime }}{\partial x}}+{\displaystyle \frac{\partial x^{\prime }}{\partial y}}{\displaystyle \frac{\partial y^{\prime }}{\partial y}}+{\displaystyle \frac{\partial x^{\prime }}{\partial z}}{\displaystyle \frac{\partial y^{\prime }}{\partial z}}\hfill \end{array}$$

(124)

The second component of the geometric tensor becomes:

$$\begin{array}{c}\hfill g^{x^{\prime }{y}^{\prime }}=-{\displaystyle \frac{2xy}{a}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)\end{array}$$

(125)

It can be shown by symmetry that the following holds for the geometric tensor components:

$$\begin{array}{c}\hfill g^{x^{\prime }{y}^{\prime }}=g^{y^{\prime }{x}^{\prime }};g^{x^{\prime }{z}^{\prime }}=g^{z^{\prime }{x}^{\prime }};g^{y^{\prime }{z}^{\prime }}=g^{z^{\prime }{y}^{\prime }};g^{z^{\prime }{z}^{\prime }}=g^{x^{\prime }{x}^{\prime }};\end{array}$$

(126)

The only non-trivial component is $g^{y^{\prime }{y}^{\prime }}$, which is obtained from:

$$\begin{array}{ccc}\hfill g^{2^{\prime }{2}^{\prime }}& =& {\displaystyle \frac{\partial {\xi }^{2^{\prime }}}{\partial {\xi }^k}}{\displaystyle \frac{\partial {\xi }^{2^{\prime }}}{\partial {\xi }^l}}{\delta }^{kl}\hfill \\ & =& {\left({\displaystyle \frac{\partial y^{\prime }}{\partial x}}\right)}^2+{\left({\displaystyle \frac{\partial y^{\prime }}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial y^{\prime }}{\partial z}}\right)}^2\hfill \end{array}$$

(127)

Hence,

$$\begin{array}{cc}{g}^{y^{\prime }{y}^{\prime }}={\displaystyle \frac{4x^2{y}^2}{a^2}}{cos}^2\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)+\hfill & \\ & \hfill {\left[1-sin\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)+{\displaystyle \frac{2y^2}{b}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)\right]}^2\end{array}$$

(128)

The contravariant metric tensor takes the final form:

$$\begin{array}{c}\hfill g^{i^{\prime }{j}^{\prime }}=\left[\begin{array}{ccc}1& -{\displaystyle \frac{2xy}{a}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)& 0\\ -{\displaystyle \frac{2xy}{a}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)& {\omega }^2+{\displaystyle \frac{4x^2{y}^2}{a^2}}{cos}^2\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)& 0\\ 0& 0& 1\end{array}\right]\end{array}$$

(129)

where

$$\begin{array}{c}\hfill \omega =1+{\displaystyle \frac{2y^2}{b}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)-sin\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)\end{array}$$

(130)

According to (115) or (116), the tensor product of the metric tensor has to be taken with ${ϵ}^{ij}$. Define the latter as the orthogonal tensor in the Euclidean homogeneous material and set it to ${ϵ}^{ij}={ϵ}_{0}diag(1,1,1)={ϵ}_0{\delta }^{ij}$. Let ${ϵ}_0=1$ for free space, $a=1$ and $b=5$, then the permittivity tensor takes the form,

$$\begin{array}{c}\hfill {ϵ}^{i^{\prime }{j}^{\prime }}=\left[\begin{array}{ccc}{\displaystyle \frac{1}{g}}& -{\displaystyle \frac{2xy}{g}}cos\left(x^2-{\displaystyle \frac{y^2}{5}}\right)& 0\\ -{\displaystyle \frac{2xy}{g}}cos\left(x^2-{\displaystyle \frac{y^2}{5}}\right)& g+{\displaystyle \frac{4x^2{y}^2}{g}}{cos}^2\left(x^2-{\displaystyle \frac{y^2}{5}}\right)& 0\\ 0& 0& {\displaystyle \frac{1}{g}}\end{array}\right]\end{array}$$

(131)

The permeability tensor has the same form as the contravariant permittivity tensor, i.e., ${ϵ}^{i^{\prime }{j}^{\prime }}\equiv {\mu }^{i^{\prime }{j}^{\prime }}$ and the parameter $g=det\left(g^{i^{\prime }{j}^{\prime }}\right)\equiv {\omega }^2$ is the determinant of the metric tensor which is given by

$$\begin{array}{c}\hfill g=1-sin\left(x^2-{\displaystyle \frac{y^2}{5}}\right)+{\displaystyle \frac{2y^2}{5}}cos\left(x^2-{\displaystyle \frac{y^2}{5}}\right)\end{array}$$

(132)

The permittivity tensor as given by (131) is for an ‘expanded’ material. Figure 11 shows how the components of this tensor vary inside the material.

Using the same mathematical analysis for the expansion of the material medium given above, the permittivity or permeability tensors for a material experiencing compression can be derived as well. For brevity reasons, the mathematical analysis and final form of the tensor is not given; however, Figure 12 shows a plot of its components.

4.4. Propagation of the Fields along Their Geodesics in Inhomogeneous Materials

Consider the propagation of an EM field in the transformed medium. The fields are changed by the inhomogeneous medium to $E^{i^{\prime }}$ from the fields in the homogeneous medium via the Lorentz transform:

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_i^{i^{\prime }}={\displaystyle \frac{\partial {\xi }^{i^{\prime }}}{\partial {\xi }^i}}\end{array}$$

(133)

where $i=1,2,3$, such that:

$$\begin{array}{ccc}\hfill E^{i^{\prime }}& =& {\mathrm{\Lambda }}_i^{i^{\prime }}{E}^i\hfill \\ & =& {\displaystyle \frac{\partial {\xi }^{i^{\prime }}}{\partial {\xi }^i}}{E}^i\hfill \end{array}$$

(134)

The components of the transformed field can thus be obtained starting with $i=1=x$

$$\begin{array}{c}\hfill E^{x^{\prime }}={\displaystyle \frac{\partial x^{\prime }}{\partial x}}{E}^x+{\displaystyle \frac{\partial x^{\prime }}{\partial y}}{E}^y+{\displaystyle \frac{\partial x^{\prime }}{\partial z}}{E}^z\end{array}$$

(135)

Similarly, for $i=2=y$, the fields become

$$\begin{array}{c}\hfill E^{y^{\prime }}={\displaystyle \frac{\partial y^{\prime }}{\partial x}}{E}^x+{\displaystyle \frac{\partial y^{\prime }}{\partial y}}{E}^y+{\displaystyle \frac{\partial y^{\prime }}{\partial z}}{E}^z\end{array}$$

(136)

and for $i=3=z$,

$$\begin{array}{c}\hfill E^{z^{\prime }}={\displaystyle \frac{\partial z^{\prime }}{\partial x}}{E}^x+{\displaystyle \frac{\partial z^{\prime }}{\partial y}}{E}^y+{\displaystyle \frac{\partial z^{\prime }}{\partial z}}{E}^z\end{array}$$

(137)

Assuming a $TE$-polarised EM field in the x-direction of a lossless medium, namely, $E^y=exp(-ikx)$ with $k=2\pi /\lambda$ where $\lambda$ is the wavelength, means that all other fields are zero: $E^x=E^z=0$. That is, the only transformed field component that is non-zero becomes:

$$\begin{array}{ccc}\hfill E^{y^{\prime }}& =& {\displaystyle \frac{\partial y^{\prime }}{\partial y}}{E}^y\hfill \\ & =& {\displaystyle \frac{\partial y^{\prime }}{\partial y}}{e}^{-ikx}\hfill \end{array}$$

(138)

Figure 13 shows the propagation of the field in the expanded medium showing how the field is close to or equal to zero in the central regions of the medium. On the other hand, Figure 14 shows the field propagating in the compressed medium and it exhibits an electromagnetic ‘hole’ analogous to Figure 3 via the effective medium approach. Furthermore, Figure 15 shows that the size of the electromagnetic ‘hole’ in the compressed material can be manipulated so that a larger region with a zero EM field can be created. Changing the transformations and parameters to

$$\begin{array}{ccc}\hfill x^{\prime }& =& x\hfill \\ \hfill y^{\prime }& =& y^4(1-exp(-4(x^2+y^2)))\hfill \\ \hfill z^{\prime }& =& z\hfill \end{array}$$

(139)

induces the ‘splitting’ of the entire electromagnetic field as shown in Figure 16. Using (105), the covariant metric tensor is obtained from (129), and it becomes:

$$\begin{array}{c}\hfill g_{i^{\prime }{j}^{\prime }}=\left[\begin{array}{ccc}1+{\displaystyle \frac{4x^2{y}^2}{a^2{w}^2}}{cos}^2\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)& {\displaystyle \frac{2xy}{a{w}^2}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)& 0\\ {\displaystyle \frac{2xy}{a{w}^2}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)& {\displaystyle \frac{1}{w^2}}& 0\\ 0& 0& 1\end{array}\right]\end{array}$$

(140)

where

$$\begin{array}{c}\hfill \omega =1+{\displaystyle \frac{2y^2}{b}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)-sin\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)\end{array}$$

(141)

The electromagnetic fields travel along the geodesic in the material, and the path is given by:

$$\begin{array}{cc}d{s}^2=\left[1+{\displaystyle \frac{4x^2{y}^2}{a^2{w}^2}}{cos}^2\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)\right]d{x}^2+{\displaystyle \frac{4xy}{a{w}^2}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)dxdy+\hfill & \\ & \hfill {\displaystyle \frac{1}{w^2}}d{y}^2+d{z}^2\end{array}$$

(142)

which must be integrated in order to obtain the entire path the fields travel along in the material medium, namely,

$$\begin{array}{c}\hfill S=\int ds\end{array}$$

(143)

For example, the geodesic of the fields along the x-coordinate is determined from

$$\begin{array}{c}S={\int }_{x_1}^{x_2}[1+{\displaystyle \frac{4x^2{y}^2}{a^2{w}^2}}{cos}^2\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right)+{\displaystyle \frac{4xy}{a{w}^2}}cos\left({\displaystyle \frac{x^2}{a}}-{\displaystyle \frac{y^2}{b}}\right){\displaystyle \frac{dy}{dx}}+ \hfill \\ \hfill {\displaystyle \frac{1}{w^2}}{\left({\displaystyle \frac{dy}{dx}}\right)}^2+{\left({\displaystyle \frac{dz}{dx}}\right)}^2{]}^{{\displaystyle \frac{1}{2}}}dx\end{array}$$

(144)

and similarly for other coordinates. In order to obtain the final geodesic for the fields along the x-direction, the derivatives appearing in (144) must be solved. This is achieved by using (85), where this time, the indices are summed over three dimensions inside the material medium, i.e., $n,i,j=1,2,3$. Thus, if ${\xi }^m=({\xi }^1,{\xi }^2,{\xi }^3)=(x,y,z)$, then for $n=1$ the equation that corresponds to the x-direction becomes

$$\begin{array}{c}1+{\mathrm{\Gamma }}_{11}^1+{\mathrm{\Gamma }}_{22}^1{\left({\displaystyle \frac{dy}{dx}}\right)}^2+{\mathrm{\Gamma }}_{33}^1{\left({\displaystyle \frac{dz}{dx}}\right)}^2+\left({\mathrm{\Gamma }}_{12}^1+{\mathrm{\Gamma }}_{21}^1\right){\displaystyle \frac{dy}{dx}})+ \hfill \\ \hfill \left({\mathrm{\Gamma }}_{13}^1+{\mathrm{\Gamma }}_{31}^1\right){\displaystyle \frac{dz}{dx}})+\left({\mathrm{\Gamma }}_{23}^1+{\mathrm{\Gamma }}_{32}^1\right){\displaystyle \frac{dy}{dx}}{\displaystyle \frac{dz}{dx}})=0\end{array}$$

The same can be performed for y and z, i.e., for $n=2$ and $n=3$, respectively. The set of non-linear differential equations obtained must be solved if the path of the field in the material is required. The solutions can then be substituted into the geodesic expression that minimises the path of the field as a function of the coordinates such as the x-direction, which is given by (144). The Christoffel symbols of the second kind ${\mathrm{\Gamma }}_{ij}^n$ are determined via (84). Figure 17 shows a combination of compression and expansion transformations yielding a ‘warped’ material medium that concentrates the field in a central region but eliminates it in the surrounding regions of the material using the following transformations:

$$\begin{array}{ccc}\hfill x^{\prime }& =& x\hfill \\ \hfill y^{\prime }& =& {sinc}^2(x/2){sinc}^2(3y/2)\hfill \\ \hfill z^{\prime }& =& z\hfill \end{array}$$

(145)

This is of use to sensors or electric components in material structures that require sensing of external fields without the risk from such things as transient fields or ESD/arcing effects and so on. Similarly, as shown in Figure 13, Figure 14, Figure 15 and Figure 16, mathematical transformations can be used that manipulate the fields inside materials for similar purposes. The mathematical approach discussed in this paper for designing such materials can be engineered in the practical sense using $3D$-printing techniques for example. The printing of thin layers or surfaces that form the material can be designed according to the field transformation approach presented. Another manufacturing approach might utilise a combination of both hyperinclusions and material transformation design for obtaining a material that mitigates various unwanted electromagnetic field effects.

5. Conclusions

Two methods were presented for the manipulation of electromagnetic effects in space materials. In particular, the problems of electrostatic discharge/arcing and EM/RF interference were studied via the dielectric function of a material. The first method consisted in an effective medium theory derived in a multidimensional framework. The effective behaviour of a material with hyperspherical inclusions was studied, with non-integer dimensions corresponding to fractal inclusion geometries. The second method consisted in ‘designing’ the material medium itself via the use of differential calculus. The fields travelled along geodesic paths inside the material, and the constitutive parameters were derived via transformations using the metric tensor and other electromagnetic field components. The theory presented was studied and validated using full-wave numerical simulations of Maxwell’s equations using the Comsol Multiphysics and FEKO software.

Embedding hyperspherical inclusions in host materials can achieve the desired effective response using values for the hyperinclusion material permittivities/permeabilities that are easily found in nature. In cases requiring negative permittivity or permeability, metamaterials or similar constructs might be necessary. In the EMT section, the paper demonstrated how any material could be modified to exhibit a high dielectric function, useful for ESD cancellation. Naturally occurring materials with such large dielectric values are rare. Additionally, natural materials may be too heavy, weak, or lossy. Instead, using strong and light materials that meet the intended purposes is crucial. The required high dielectric response can be achieved by incorporating inclusions with the properties described in this paper, enhancing the material’s effective performance without altering its bulk properties. For designing material media using the transformation medium theory, 3D printing can create layers with graded permittivity or permeability as mathematically described in the paper. The 3D printer can be programmed to follow algorithms based on this paper’s mathematical results, printing out the layer properties and thickness accordingly.

This paper focused on electromagnetic properties, but the mathematical analysis of the dielectric function is interconnected with properties like electrical conductivity, thermal conductivity, and mechanical strength. The EMT approach can determine thermal conductivity using the Maxwell–Eucken model, which estimates the effective thermal conductivity of a composite material based on the constituent phases’ thermal conductivities and (hyper)volume fractions. It can also determine the effective elastic moduli (e.g., Young’s modulus, bulk modulus, and shear modulus) of composite materials using the Voigt, Reuss, and Hill models. The Voigt and Reuss bounds provide upper and lower limits for the effective moduli. The EMT approach can also be used to obtain estimates of the effective acoustic impedance of composite materials, which is useful in materials engineering and acoustic metamaterials design among other things.

Experimental fabrication and studies are required to validate the theoretical results and numerical simulations presented in the paper. This includes exploring and optimising material compositions and structures to enhance electromagnetic field control methods. Experiments with various inclusions, geometric configurations (fractal or conventional), and host material properties will further this research.