Effective electrical resistivity in a square array of oriented square inclusions

Square tessellation, i.e. a two-dimensional repetition of a square unit cell across the plane, is amongst the most straightforward engineering designs, naturally also found in nanofabrication. While hexagonal tiling is theoretically more efficient [1], and may exhibit marginally superior physical properties [e.g. 2], quadratic arrays with quadratic sub-features are often easier to manufacture and operate, as evidenced by their prevailing occurrence in micro-optoelectronic devices and electromagnetic metamaterial applications [3–11]. The broad range of target materials considered in this paper includes various types of biosensors, plasmonic materials, and optoelectronic devices (figure 1), all of which can be modeled as bi- or multi-composites as follows. The simplest imaginable bi-composite lies on the nanohole-nanofishnet continuum (figures 1(a)–(c)), where a uniform thin film (material 1) is perforated by square holes, constant both in size and spatial offsets with respect to one another (material 2 is void). A more sophisticated bi-composite results from filling that void with a different material, termed an inclusion (e.g. nanopatches in figure 1(d)). Planar multi-composites may be formed by the recursive nesting of inclusions within inclusions (e.g. single and double split-ring resonators in figures 1(e)–(f), respectively). Eventually, considerable complexity can be achieved by warping either the nanohole or the nanopatch arrays into the third dimension (e.g. nanowell/nanopyramid arrays in figures 1(g)–(h), and μLED array in figure 1(i), respectively). Since the electrical properties of all the abovementioned structures may be of key importance to their proper functioning, the premise of this paper is to make this broad range of materials and structures amenable to the theory of effective electrical transport, briefly outlined below.

Zoom In Zoom Out Reset image size

Most of the modern treatment of effective transport through composite media is due to Lord Rayleigh, who in 1892 considered the effects of circular or spherical inclusions of resistivity ρ₂, periodically embedded within a host medium of resistivity ρ₁ [12]. Serving the 1960s boom of microfabrication, further seminal theory was developed by Keller [13] and later Dykhne [14], who proved that the effective sheet resistance of an infinite (and possibly random) planar checkerboard bi-composite obeys ρ_eff = $\sqrt{{\rho }_{1}{\rho }_{2}}.$ The seminal mathematical analysis of effective transport in a bi-composite with oriented square inclusions [15] (i.e. inclusions aligned in parallel to the principal axes of the square unit cell), appears to be due to Mortola and Stefeé [16]. A decade and a half later, one of their conjectures was rigorously proven by Obnosov [17], who introduced the variable transformation:

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{1-{\rho }_{2}/{\rho }_{1}}{1+{\rho }_{2}/{\rho }_{1}}\end{eqnarray} \tag{ 1 }$$

and expressed the effective resistivity of a doubly-periodic bi-composite, whose inclusions occupy precisely c = 1/4 of the unit cell, as:

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{{\rm{e}}{\rm{f}}{\rm{f}}}}{{\rho }_{1}}=\sqrt{\displaystyle \frac{2-{\rm{\Delta }}}{2+{\rm{\Delta }}}}.\end{eqnarray} \tag{ 2 }$$

In order to appreciate equation (2), its key behavior and selected assumptions [17] are visualized in figure 2. The square unit cell (figure 2(a)) consists of a thin sheet, composed of two media with isotropic resistivities ρ₁ (host) and ρ₂ (inclusion), c being the fractional area of the inclusion in the unit cell. In the analytical solution for c = 1/4, the effective resistance (equation (1), figure 2(b)) varies between ρ₁/$\sqrt{3}$ to $\sqrt{3}$ ρ₁ for ρ₂ ≪ ρ₁ and ρ₁ ≪ ρ₂, respectively. Note, that the response is symmetric about ρ₂ = ρ₁ [13], and that the factor $\sqrt{3}$ is purely geometric [17]. Three representative scenarios in figure 2(b) are further shown in figures 2(c)–(e), visualizing the preferential current flow in the low-resistivity medium (black arrows, with size proportional to current density), alongside equipotential lines (grey contours, 0.05 V apart).

Zoom In Zoom Out Reset image size

Micro four-point probe (M4PP) metrology is a well-established electrical technique, which enables measurement of key electrical material properties (e.g. sheet resistance or resistivity, carrier density and mobility, tunneling magnetoresistance, etc) on micrometer length scales [18, 19]. Given that the spatial resolution of the technique scales with the electrode pitch, sub-pitch resolution can be difficult to obtain [20]. For example, it has been shown that spatial variations in sheet resistance can be mapped reliably only when the spatial extent of such variations exceeds probe pitch by at least an order of magnitude [21]. Sub-pitch resolution may arguably be attained by deconvolution of data from a dense areal sampling [22]; however, this involves non-unique inversion of an underdetermined system, and thus is applicable only for simple (e.g. quasi-1D) patterns. Although recent theoretical progress has enabled accurate M4PP measurements of small pads with dimensions barely in excess of the M4PP probe footprint [23, 24], there is still a clear lack of M4PP methodology for smaller (i.e. sub-pitch) motif sizes. Motivated to close this gap, the current study devises a data reduction scheme for accurate M4PP metrology of electrical transport in doubly-periodic bi- or multi-composites with sub-pitch periodicity.

To verify that equation (2) remains valid also for M4PP measurements, where the electrical current propagates in a markedly different manner to what is modeled in figure 2, we undertook a full-scale M4PP simulation (figure 3) via finite element modeling (FEM) using COMSOL 5.5 [25]. Over a large 59 × 59 square matrix of bordering unit cells, a four-point resistance measurement is simulated via four collinear and equidistant circular contacts of radius r₀, where current is passed between two contacts and the voltage drop is registered by the remaining two (figure 3(a)). From the 4! = 24 combinations of electric terminal assignment to the four electrodes (I+, V+, V− and I−), we simulate only the three baseline configurations [18]; the resulting dual-configuration sheet resistances [18] are numerically indistinguishable, regardless of the chosen configuration pair.

Zoom In Zoom Out Reset image size

The numerically simulated ρ_eff/ρ₁ as a function of ρ₂/ρ₁ (symbols in figure 3(b)) match theory (solid line in figure 3(b)) to within a few percent (figure 3(c)) and thus confirm the relevance of equation (2) for M4PP array-metrology. The deviations in figure 3(c) arise not from numerical artefacts (which are <0.017% as will be shown in section 2.3), but mainly from the M4PP method itself (i.e. a finite electrical contact size r₀, and the position of each contact with respect to inclusions). Nevertheless, even for a probe whose pitch is only marginally larger (factor 1.24) than the unit cell size (circles in figures 3(a)–(c)), these deviations clearly remain second-order effects, overlain on top of an otherwise well-behaved theoretical trend (equation (2)). At this point, we can conclude that while equation (2) seems readily applicable to M4PP array-metrology, its inherent geometric constraint (c = 1/4) is specific, and is of limited practical use. However, anticipating that the symmetry of equation (2) may extend to other inclusion sizes, we take a more in-depth look at a single unit cell, as shown next.

For the rigorous mathematical derivation of equation (2), applicable only for c = 0.25, the reader is referred to [17]. Since from a mathematical standpoint, a closed-form analytical solution for c ≠ 0.25 seems improbable [26], here we resort to physical first principles, and make simplifying assumptions to obtain an approximate expression for ρ_eff for arbitrary values of c. To further generalize the treatment, in all our subsequent algebra we operate with isotropic sheet resistances R = ρ/d, i.e. resistivities ρ scaled by their corresponding materials' thicknesses d. This makes composites, whose components might have different thicknesses (all negligible in comparison to the M4PP probe pitch), amenable to the very same treatment.

Let the periodic array have a unit cell of size a × a, with a centered square inclusion of size $a\sqrt{c}$  × $a\sqrt{c},$ as before. When electrical current is forced across the unit cell from its left to its right edge (figures 2(a) and 4(a)), the vertical edges are potential lines, while the horizontal edges are flux lines. Due to symmetry, the bisecting vertical and horizontal lines through the unit cell are also potential and flux lines, respectively; thus, all further treatment can be reduced to a single quadrant of the unit cell (e.g. its lower right quarter in figure 4(a)). The host and inclusion materials have sheet resistances R₁ and R₂, respectively. We now recast the problem into a resistor network, under two complementary assumptions. In the first, all equipotential lines through the quadrant are taken to be approximately vertical, and thus the effective resistance can be calculated from a simple network with two parallel resistors in series with a third resistor (figure 4(b)):

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{e}}{\rm{f}}{\rm{f}},1}}{{R}_{1}}=\displaystyle \frac{\left[{\left(1-\sqrt{c}\right)}^{2}+\sqrt{c}\right]\tfrac{{R}_{2}}{{R}_{1}}+\sqrt{c}\left(1-\sqrt{c}\right)}{\sqrt{c}+\left(1-\sqrt{c}\right)\tfrac{{R}_{2}}{{R}_{1}}}.\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

Under the alternative assumption, all flux lines are taken to be longitudinal, and thus the effective resistance can be calculated by considering two series resistors in parallel with a third resistor (figure 4(c)):

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{e}}{\rm{f}}{\rm{f}},2}}{{R}_{1}}=\displaystyle \frac{\sqrt{c}\tfrac{{R}_{2}}{{R}_{1}}+\left(1-\sqrt{c}\right)}{\left[{\left(1-\sqrt{c}\right)}^{2}+\sqrt{c}\right]+\tfrac{{R}_{2}}{{R}_{1}}\sqrt{c}\left(1-\sqrt{c}\right)}.\end{eqnarray} \tag{ 4 }$$

Although neither of equations (3)–(4) is accurate, a more adequate description of the effective current distribution in the plane may be obtained if we imagine an infinite (and possibly random) checkerboard of quadrants, each of which obeys either equation (3) or (4). Using Dykhne's formalism [14], we take the geometric mean of equations (3)–(4) and obtain:

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{e}}{\rm{f}}{\rm{f}}}}{{R}_{1}}=\sqrt{\displaystyle \frac{{R}_{{\rm{e}}{\rm{f}}{\rm{f}},1}}{{R}_{1}}\cdot \displaystyle \frac{{R}_{{\rm{e}}{\rm{f}}{\rm{f}},2}}{{R}_{1}}}=\sqrt{\left(\displaystyle \frac{{m}_{1}\tfrac{{R}_{2}}{{R}_{1}}+{m}_{2}}{{m}_{1}+{m}_{2}\tfrac{{R}_{2}}{{R}_{1}}}\right)\left(\displaystyle \frac{{m}_{4}\tfrac{{R}_{2}}{{R}_{1}}+{m}_{3}}{{m}_{4}+{m}_{3}\tfrac{{R}_{2}}{{R}_{1}}}\right)},\end{eqnarray} \tag{ 5 }$$

where ${m}_{1}=\sqrt{c},\,{m}_{2}=1-\sqrt{c},\,{m}_{3}=\sqrt{c}-c,$ and ${m}_{4}\,=1-\sqrt{c}+c.$ Such spatial averaging is valid if one imagines an electrical current spreading in 2D through many unit cells. As will be shown next, equation (5) performs equally well for even a single unit cell. At this point we note that equation (5) qualifies as a purely algebraic extension of equation (2), since it reduces to the latter for c = 1/4. To test the validity of equation (5) for other values of c, we compare it to finite element model simulations, as described next.

To test equation (5) against numerical results, FEM of electrical currents in a single unit cell were performed in COMSOL 5.5 [25]. Electrical currents of opposite signs flow in and out of the left and right edges of the unit cell, with voltage fixed to zero at its center. Solving for the electric field in steady-state, we divide the potential drop between the inlet and outlet current terminals by the current and sheet thickness to obtain the effective sheet resistance of the composite medium. The model was solved using a 'physics-controlled mesh' with 'extremely small' element size, corresponding to a Delaunay triangulation of ∼25000 elements, with an average element size ∼4.5 orders of magnitude smaller than the unit cell. To verify the universality of the results, we operated on non-dimensionalised variables, testing the response of R_eff/R₁ as a function of R₂/R₁ for a range of varying unit cell dimensions (size and thickness), absolute sheet resistances (R₁ and R₂), and current intensities (I), all of which had no numerical effect on the results. At the next step, the nondimensional variables of interest were varied from c = 0.01 to c = 0.99 (step of 0.01), and R₂/R₁ from 10⁻⁶ to 10⁶ (10 steps per decade).

Predictions of equation (5) for a few representative values of c (red curves in figure 4(d)) are overlain by corresponding FEM simulations (dots), showing good agreement over a wide range of c and R₂/R₁. To quantify the error, figure 4(e) shows the maximum deviation of equations (3)–(5) with reference to the numerical model for all numerically simulated values of c. The approximations of either equations (3)–(4) lead to relative errors of up to 16%. In contrast, their geometrical mean (equation (5)) keeps the relative error under 5% over the entire range of c (figure 4(e)). Remembering that equation (5) with c = 0.25 is an exact solution, we can quantify the accuracy of the FEM model itself. For c = 0.25, the maximum relative error of the simulated R_eff/R₁ ranged from 10⁻⁸ at R₂ = R₁, to a maximum of 0.017% at the scenario extremities R₁≫R₂ and R₁≪R₂. The latter figure (0.017%) is thus a representative value for error due to purely numerical artefacts.

Recognizing the power law in equation (2), we conjecture a purely algebraic generalization:

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{e}}{\rm{f}}{\rm{f}}}}{{R}_{1}}={\left(\displaystyle \frac{\alpha -{\rm{\Delta }}}{\alpha +{\rm{\Delta }}}\right)}^{\beta }\end{eqnarray} \tag{ 6 }$$

and hypothesize that the particular parameters for c = 1/4 (α = 2 and β = 1/2 in equation (2)) might be drawn from a broader continuum of values. As it turns out, the fit quality of equation (6) to the numerically simulated results is outstanding (figure 5(a)). The fact, that the maximum error of equation (6) is comparable to that of numerical artefacts (figure 5(b)), suggests that equation (6) may be generally valid also for c ≠ 0.25 (subject to proof, not attempted here). While the trends of α and β both appear smooth and continuous with respect to c (figure 5(c)), and are even strongly anticorrelated (Pearson's ρ = –0.84), the truly remarkable observation is that their ratio β/α is numerically indistinguishable from c (figure 5(d)). Utilizing the discovered relationship β = αc to get rid of β, we rewrite equation (6) as:

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{e}}{\rm{f}}{\rm{f}}}}{{R}_{1}}={\left(\displaystyle \frac{\alpha -{\rm{\Delta }}}{\alpha +{\rm{\Delta }}}\right)}^{\alpha c}\end{eqnarray} \tag{ 7 }$$

and proceed to investigate whether the relationship between α and c is also amenable to parameterization.

Zoom In Zoom Out Reset image size

After some trial and error, a plausible relationship emerges between the transformed parameters α' = α – 1 and c' = c/(1 − c) (figure 6(a)), wherein α' approaches a horizontal asymptote for c' → 0, and α'→1/c' for c'→∞. This suggests a hyperbolic relationship of the form α' ∼ 1/(1 + c'), possibly requiring further modifications to achieve a smooth transition between the two asymptotic behaviors. In particular, we find that:

$$\begin{eqnarray}&&\alpha (c)=1+\displaystyle \frac{{\gamma }_{1}}{{\left[{\gamma }_{2}+{\left(\tfrac{c}{1-c}\right)}^{{\gamma }_{3}}\right]}^{{\gamma }_{4}}}\end{eqnarray} \tag{ 8 }$$

yields excellent (R² = 0.9997) fits to α as a function of c, such as the one shown in figure 6(b), for the particular coefficients of γ₁ = 0.9499, γ₂ = 0.8063, γ₃ = 2.158, and γ₄ = 0.4576. Aiming to minimize the maximum error of equation (7) across the entire FEM dataset (rather than the misfit of equation (8) to the bestfit values of α), we obtain slightly different coefficients γ₁ through γ₄, which results in:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{R}_{{\rm{e}}{\rm{f}}{\rm{f}}}}{{R}_{1}}\,=\,{\left(\displaystyle \frac{\alpha (c)-\tfrac{1-{R}_{2}/{R}_{1}}{1+{R}_{2}/{R}_{1}}}{\alpha (c)+\tfrac{1-{R}_{2}/{R}_{1}}{1+{R}_{2}/{R}_{1}}}\right)}^{\alpha (c)\cdot c},\\ \alpha (c)\,=\,1+\displaystyle \frac{0.9707}{{\left[0.9193+{\left(\tfrac{c}{1-c}\right)}^{2.1261}\right]}^{0.4671}}\end{array}\end{eqnarray} \tag{ 9 }$$

which is the core result of this work. Evaluating the goodness-of-fit of equation (9) to the FEM dataset, we find that within the range 0.01 ≤ c ≤ 0.99, equation (8) represents the FEM results with a maximum relative error of 0.3% (figure 6(d)). Note that the presence of four minima in the error of equation (8) (red curve in figure 6(d)) at c values of 0.01, 0.56, 0.87 and 0.98 arise from the four fitting parameters γ₁ through γ₄ in equation (8).

Zoom In Zoom Out Reset image size

As with Kirchhoff's circuit law, which can be iteratively applied on subsets of a linear resistor network, we hypothesize that equation (9) may be similarly used to predict the behavior of a nested system of square inclusions. We expect that this should be possible via a recursive application of equation (9), further shorthanded as f(c, R₂/R₁), from the innermost to the outermost inclusions. For simplicity, we constrain ourselves to two materials R₁ and R₂ as before, progressively nesting them within one another (bottom right of figure 7). Alongside a simple square inclusion as already treated before (equation (10a )), we consider a single (equation (10b )) and a double (equation (10c )) square ring inclusion, writing the corresponding recursion as:

$$\begin{eqnarray}&&{R}_{{\rm{e}}{\rm{f}}{\rm{f}}}/{R}_{1}={f}\left(0.7,\tfrac{{R}_{2}}{{R}_{1}}\right),\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{R}_{{\rm{e}}{\rm{f}}{\rm{f}}}/{R}_{1}={f}\left(0.7,\tfrac{{R}_{2}}{{R}_{1}}{f}\left(0.7,\tfrac{{R}_{1}}{{R}_{2}}\right)\right),\end{eqnarray} \tag{ 10b }$$

$$\begin{eqnarray}&&{R}_{{\rm{e}}{\rm{f}}{\rm{f}}}/{R}_{1}=f\left(0.8,\tfrac{{R}_{2}}{{R}_{1}}f\left(0.8,\tfrac{{R}_{1}}{{R}_{2}}f\left(0.7,\tfrac{{R}_{2}}{{R}_{1}}f\left(0.7,\tfrac{{R}_{1}}{{R}_{2}}\right)\right)\right)\right).\end{eqnarray} \tag{ 10c }$$

Zoom In Zoom Out Reset image size

The predictions of equation (9a )–(9c ) are plotted as continuous curves in figure 7. After confirming that the theoretical behavior of nested inclusions discernibly diverges from that of a single-inclusion system, we proceeded to examine our theoretical predictions against a numerical simulation (symbols in figure 7). The match between theory and simulation in figure 7 confirms our hypothesis that equation (9) can be used recursively to predict the behavior of arbitrarily complex hierarchies of nested square inclusions, granted that they are all oriented and centered within each other. At this point, it is also important to note that since the approximation errors of the nested instances of equation (9) may be multiplicative, caution must be taken with increasing recusion depth.

Although the ability to model the effective resistivity of planar square grids is important, an increasing number of thin film devices involve out-of-plane nanostructuring in order to attain certain functionalities (e.g. flexibility in particular). An example of a doubly-periodic non-planar array of hollow frustums (lacking their base faces, but interconnected through their edges), whose top and side faces may have different electric properties, is shown in figure 8(a), alongside four M4PP electrodes probing that surface (to approximate scale; cf. figure 3(a)). Since current can flow only at the surface of the target structure, it is reasonable to assume that the problem of effective sheet resistance in figure 8(a) can be reduced to an equivalent 2D problem via conformal mapping [27]. Since the latter is angle-preserving, all the core mathematical description of the electrical field will be retained (i.e. derivatives and gradients) on an equivalent and simpler topology.

Zoom In Zoom Out Reset image size

To test this idea, a series of square truncated pyramids (frustums), with aspect ratios (height/base) ranging from 0.1 to 1, were constructed (figure 8(b)) and densely meshed via Delaunay triangulation (not shown). Each frustum is hollow, lacks its base face, and is a composite of medium R₁ at its sidewalls (colored), and R₂ at its top face (grey); the relative area of R₂ with respect to the total surface area is quoted below each shape in figure 8(b). To avoid a rigorous mathematical treatment at this proof-of-concept stage, we resorted to a quasi-conformal mapping technique prevalent in image analysis, known as Teichmüller mapping, and used an openly-available algorithm [28] to project the dense frustum triangulations onto a planar unit square (figure 8(c)), whereon the relative area of R₂ noticeably scales down (quoted below each shape in figure 8(c)). After verifying that the obtained angle distortions are small (0° ± 2.5°), we approximate c as the fractional area of R₂ in the Teichmüller map, and use equation (9) to predict the effective resistivity of the composite non-planar bi-composite as a function of R₂/R₁ (continuous curves in figure 8(d)). The original 3D frustum is then subjected to a FEM resistivity simulation on a unit cell (as in section 2, only on 3D topology), and the results (square symbols in figure 8(d)) are overlain on the theoretical expectation.

The error of the theoretical expectation relative to actual FEM simulations is <12% for the shallowest frustum considered, and <0.3% for the tallest. A reasonable explanation for the progressive departure of analytics from numerics for the shallowest frustums is that their originally strictly square top faces become mildly rounded on a Teichmüller map, thereby ceasing to be 'square'. For the majority of the taller frustums though, the area of R₂ in the Teichmüller map is so small that essentially it is indistinguishable from an equivalent square. However, a larger and non-square shape occupying the majority of the unit cell (e.g. two leftmost frustums in figures 8(b)–(c)) could appreciably affect the electric flow lines (see figures 2(c)–(e)), bringing about the observed departure of theory from simulations as observed in figure 8(d). While several other factors may be in play and should be considered (the accuracy of the Teichmüller mapping, and of the electric field simulation in 3D), a fractional error of <12% for a broad range of non-planar unit cell geometries appears quite promising, especially if it replaces a computationally-demanding numerical simulation with a robust semi-analytical approach. As a final remark, we mention that full-scale M4PP simulations, similar to those in section 2 yet on arrays of interconnected frustums, retain the same behavior as a single unit cell (figures 8(b), (d)), and thus strengthen the relevance of this approach for M4PP metrology as well.