Effective electrical resistivity in a square array of oriented square inclusions

The continuing miniaturization of optoelectronic devices, alongside the rise of electromagnetic metamaterials, poses an ongoing challenge to nanofabrication. With the increasing impracticality of quality control at a single-feature (-device) resolution, there is an increasing demand for array-based metrologies, where compliance to specifications can be monitored via signals arising from a multitude of features (devices). To this end, a square grid with quadratic sub-features is amongst the more common designs in nanotechnology (e.g. nanofishnets, nanoholes, nanopyramids, μLED arrays etc). The electrical resistivity of such a quadratic grid may be essential to its functionality; it can also be used to characterize the critical dimensions of the periodic features. While the problem of the effective electrical resistivity ρ eff of a thin sheet with resistivity ρ 1, hosting a doubly-periodic array of oriented square inclusions with resistivity ρ 2, has been treated before (Obnosov 1999 SIAM J. Appl. Math. 59 1267–87), a closed-form solution has been found for only one case, where the inclusion occupies c = 1/4 of the unit cell. Here we combine first-principle approximations, numerical modeling, and mathematical analysis to generalize ρ eff for an arbitrary inclusion size (0 < c < 1). We find that in the range 0.01 ≤ c ≤ 0.99, ρ eff may be approximated (to within <0.3% error with respect to finite element simulations) by: ρeff=ρ1α(c)−1−ρ2/ρ11+ρ2/ρ1α(c)+1−ρ2/ρ11+ρ2/ρ1c·α(c), α(c)=1+0.97070.9193+c1−c2.12610.4671. whereby at the limiting cases of c → 0 and c → 1, α approaches asymptotic values of α = 2.039 and α = 1/c − 1, respectively. The applicability of the approximation to considerably more complex structures, such as recursively-nested inclusions and/or nonplanar topologies, is demonstrated and discussed. While certainly not limited to, the theory is examined from within the scope of micro four-point probe (M4PP) metrology, which currently lacks data reduction schemes for periodic materials whose cell is smaller than the typical μm-scale M4PP footprint.

The continuing miniaturization of optoelectronic devices, alongside the rise of electromagnetic metamaterials, poses an ongoing challenge to nanofabrication. With the increasing impracticality of quality control at a single-feature (-device) resolution, there is an increasing demand for array-based metrologies, where compliance to specifications can be monitored via signals arising from a multitude of features (devices). To this end, a square grid with quadratic sub-features is amongst the more common designs in nanotechnology (e.g. nanofishnets, nanoholes, nanopyramids, μLED arrays etc). The electrical resistivity of such a quadratic grid may be essential to its functionality; it can also be used to characterize the critical dimensions of the periodic features. While the problem of the effective electrical resistivity ρ eff of a thin sheet with resistivity ρ 1 , hosting a doubly-periodic array of oriented square inclusions with resistivity ρ 2 , has been treated before (Obnosov 1999 SIAM J. Appl. Math. 59 1267-87), a closed-form solution has been found for only one case, where the inclusion occupies c=1/4 of the unit cell. Here we combine first-principle approximations, numerical modeling, and mathematical analysis to generalize ρ eff for an arbitrary inclusion size (0<c<1). We find that in the range 0.01c0.99, ρ eff may be approximated (to within<0.3% error with respect to finite element simulations) by: 1. Introduction

Target materials and their effective transport parameters
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][4][5][6][7][8][9][10][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 nanoholenanofishnet 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 multicomposites 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 Figure 1. Nanostructures amenable to the hereby presented theory include, among others: (a) nanohole arrays [3], (b)-(c) fishnets [4,5], (d) plasmonic nanoarrays [6], (e)-(f) single and double split-ring resonators [7,8], (g) biosensor nanowells [9], (h) nanopyramids [10], and (i) μLED arrays [11]. Image credits: μ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. 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 ρ 2 , periodically embedded within a host medium of resistivity ρ 1 [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 = r r . 1 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: and expressed the effective resistivity of a doubly-periodic bicomposite, whose inclusions occupy precisely c=1/4 of the unit cell, as: 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 ρ 1 (host) and ρ 2 (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 ρ 1 / 3 to 3 ρ 1 for ρ 2 = ρ 1 and ρ 1 = ρ 2 , respectively. Note, that the response is symmetric about ρ 2 =ρ 1 [13], and that the factor 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).

Four terminal sensing of a doubly-periodic array of square inclusions
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 fourpoint resistance measurement is simulated via four collinear and equidistant circular contacts of radius r 0 , 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.
The numerically simulated ρ eff /ρ 1 as a function of ρ 2 /ρ 1 (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 arraymetrology. 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 0 , 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.

First-principles approximation
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 c ×a 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 1 and R 2 , 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)): 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)): 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 .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.

Numerical model
To test equation (5) against numerical results, FEM of electrical currents in a single unit cell were performed in COM-SOL 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  (2)). 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 1 as a function of R 2 /R 1 for a range of varying unit cell dimensions (size and thickness), absolute sheet resistances (R 1 and R 2 ), 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 2 /R 1 from 10 −6 to 10 6 (10 steps per decade).
Predictions of 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 1 ranged from 10 −8 at R 2 =R 1 , to a maximum of 0.017% at the scenario extremities R 1 ?R 2 and R 1 =R 2 . The latter figure (0.017%) is thus a representative value for error due to purely numerical artefacts.

Analytical approximation
Recognizing the power law in equation (2), we conjecture a purely algebraic generalization: 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: c eff 1 and proceed to investigate whether the relationship between α and c is also amenable to parameterization.
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: yields excellent (R 2 =0.9997) fits to α as a function of c, such as the one shown in figure 6(b), for the particular coefficients of γ 1 =0.9499, γ 2 =0.8063, γ 3 =2.158, and γ 4 =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 γ 1 through γ 4 , which results in: which is the core result of this work. Evaluating the goodnessof-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)

Nested structures
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 2 /R 1 ), from the innermost to the outermost inclusions. For simplicity, we constrain ourselves to two materials R 1 and R 2 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: 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.

Three-dimensional structures
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.
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 1 at its sidewalls (colored), and R 2 at its top face (grey); the relative area of R 2 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 2 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 2 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 2 /R 1 (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 2 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.

Conclusions
The determination of the effective properties of composites continues to be an active and exciting area of research [15,29,30], particularly when single-feature (or singledevice) metrology at the nanoscale becomes physicallyimprobable [31]. Our present work extends a former treatment of effective resistivity in square arrays of quadratic inclusions [17] to arbitrarily-sized inclusions. Apart from equation (9), which is a surprisingly compact and accurate generalization of [17], we have also devised calculation schemes for recursively-nested inclusions, and non-planar topologies. The maximum error of our obtained approximation (0.3%) is comparable to the current best-case reproducibility of M4PP measurements on advanced thin film material stacks [32]. Should the need for higher accuracy arise, the interpolation of tabulated values of α(c) (appendix) reduces the maximum error of equation (9) to 0.03%. As evident from our analysis, other variables related to the M4PP method may superpose a few additional percent deviation from the theoretical predictions [33], thus making equation (9) satisfactory for any basic data reduction in M4PP array-metrology. While an analogous treatment for hexagonal arrays (e.g. graphene nanomesh devices [34]) is easily conceivable, this remains outside of the scope of the present study.
While any theoretical inference from numerical simulations should be treated with caution, the empirically discovered relationship in equation (7) is unequivocal: the effective resistivity in doubly-periodic array of oriented square inclusions follows a power law, whose exponent is proportional to the inclusion's fractional area. Equation (7) remains a conjecture, which we hope could turn useful for a more rigorous analysis, as has already occurred in this field before [15]. As electronic devices and interconnects continue to shrink, the extraction of accurate and statistically robust information from a given array of devices is of extreme importance in nanofabrication [31]. We are thus hopeful that our current work extends the M4PP application portfolio to a broader set of materials and structures characterized by double periodicity, where the unit cell size may be (considerably) smaller than the M4PP footprint.

Acknowledgments
This work has been supported by grants 8054-00020B (Denmark's Innovation Fund) and 8048-00088B Sapere Aude (Independent Research Fund Denmark). BG is grateful to Yurii Obnosov for discussions on extending his seminal 1999 treatment, and to Frederik Østerberg, Lior Shiv, Kristoffer Kalhauge, Rong Lin, and Yakov Guralnik for feedback on work in progress.

Comment added in proof
In a L'esprit de l'escalier of sorts, we have noticed that by defining Δ eff =(1 -R eff /R 1 )/(1+R eff /R 1 ) (see equation (2)), a strikingly linear relationship Δ eff @ cΔ is observed. Rearranging the latter, we obtain R eff /R 1 =(1 -cΔ)/(1+cΔ), which does not only surpass equation (5) by brevity, but which is also marginally more accurate (<4.8% error for any c), thus probably qualifying as the best back-of-the-envelope calculation scheme to date.  (9) with c estimated from the quasi-conformal mapping (continuous curves); FEM and theory diverge by not more than 12% (relative).