Analytical solution for two-dimensional Laplace's equation in a shallow domain containing coplanar interdigitated boundaries
AAnalytical solution for two-dimensional Laplace’sequation in a shallow domain containing coplanarinterdigitated boundaries
Cristian F.
Guajardo Yévenes ID [email protected] Biological Engineering Program, Faculty of EngineeringPilot Plant Development and Training Institute
King Mongkut’s University of Technology Thonburi, Thailand
Werasak
Surareungchai ID [email protected] School of Bioresources and TechnologyNanoscience & Nanotechnology Graduate Program
King Mongkut’s University of Technology Thonburi, Thailand
January 7, 2021
Abstract
Laplace’s equation appears frequently in physical applications involving conservablequantities. Among these applications, miniaturized devices have been of interest, inparticular those using interdigitated arrays. Therefore, we solved the two-dimensionalLaplace’s equation for a shallow or finite domain consisting of interdigitated bound-aries. We achieved this by using Jacobian elliptic functions to conformally transformthe interdigitated domain into a parallel plates domain. The obtained expressionsfor potential distribution, flux density and flux allow for arbitrary domain height,different band widths and asymmetric potentials at the interdigitated array, besidesconsidering fringing effects at both ends of the array. All these expressions dependonly on relative dimensions, instead of absolute ones. With these results we showedthat the behavior in shallow or finite domains approaches that of a semi-infinite do-main, when its height is greater than the separation between the centers of consecutivebands. We also found that, for any desired but fixed flux, bands of equal width min-imize the total surface of the interdigitated array. Finally, we present approximateexpressions for the flux, based on elementary functions, which can be applied to easethe calculation of currents (faradaic or non-faradaic), capacitances and resistancesamong other possible applications. 1 a r X i v : . [ phy s i c s . c l a ss - ph ] J a n INTRODUCTION ρ Q j = − γ ∇ u uc mol m − ϕ = − D ∇ c mol m − s − – ρ H = cρ M T J m − q = − k ∇ T J m − s − T K ρ E = ∇ · D C m − j = − σ ∇ V C m − s − D = − (cid:15) ∇ V C m − Table 1:
Density of conservable quantity , flux density, related potentials or other aux-iliary variables, and their proportionality constants.
Generic density ρ Q , generic fluxdensity j , generic potential u and generic proportionality constant γ . Concentration c ,diffusion flux (density) ϕ , diffusion coefficient D . Heat density ρ H , heat flux (density) q , temperature T , specific heat capacity c , mass density ρ M and thermal conductivity k . Electric charge density ρ E , electric current density j , electric potential V , electricdisplacement D , electric conductivity σ and absolute permittivity (cid:15) . Several physical applications involve the density ρ Q of some conservable quantity Q (and a related potential), for example: concentration (which is also a potential in themathematical sense), heat density (temperature), or electric charge density (electricpotential). See Table 1 for details. If there is a flux density j modifying the amountof Q inside the domain of interest, then the continuity equation applies. In particular,if such flux density is proportional only to the gradient of some potential u , then thecontinuity equation is reduced to ∂ ρ Q ∂t = − ∇ · j and j = − γ ∇ u ⇒ ∂ ρ Q ∂t = γ ∇ u (1)In these cases, Laplace’s equation appears naturally when the density of the conserv-able quantity inside the domain of interest is zero ( ρ Q = 0 ) or, more commonly, whenthis density doesn’t change in time ( ∂ρ Q /∂t = 0 ) ∇ u = 0 (2)Among real applications, miniaturized (and more recently wearable) electronics,actuators and sensors have been of increasing interest. In particular, devices using interdigitated arrays (IDA) are a popular choice [1–3], due to ease of fabrication andadvantages in performance [4–6]. This motivates the use of Laplace’s equation in afinite or shallow domain with interdigitated boundaries.Solving this equation enables us to know the distribution of the potential u insteady state, which allows the computation of the flux density j and subsequently theflux i . An analytical expression for the latter is important in real applications, sincethe flux of a quantity can be measured experimentally when applying a difference ofpotentials at the IDA. Some examples of flux include: faradaic current (in electro-chemistry), rate of heat flow, and non-faradaic current (for electronic devices). Thelatter can be directly related to capacitantes and resistances.Analytically solving Laplace’s equation inside a finite or shallow domain with inter-digitated boundaries is challenging due to its mixed boundary conditions: Potentialsat the IDA bands (Dirichlet boundaries) alternate with insulations at the gaps (Neu-mann boundaries). Several attempts have been made to solve this problem, amongwhich we would like to highlight important results obtained for electrochemistry andresistive/capacite sensors.In electrochemistry, Aoki and colleagues [7] were probably the first in explaining,analytically and accurately, the faradaic current generated experimentally at IDA In vector calculus, the flux density j corresponds to the amount of quantity flowing per unitarea per unit time, while the flux i corresponds to the amount of quantity flowing per unit time. Intransport theory, in contrast, j is called flux while i is called flow rate . See Wikipedia’s entry Flux • . DESCRIPTION OF THE PROBLEM ± with respect to their ideal counterpart). Let an IDA be located at the bottom of a rectangular box of height H , whose wallsbehave as perfect insulators (Figure 1). The IDA consists of two arrays, A (black)and B (gray), of N A and N B bands each, where two consecutive bands are separatedby a distance W AB from each of their centers. The bands have width of w A and w B ,a common length L , and are assumed to have no elevation (or thickness) above thebottom of the domain.For the rectangular box we will consider two cases. In an ideal case, the IDA fitsperfectly inside the domain, such that all four vertical walls touch the borders of theIDA, and the two most exterior bands (one at each end) have only half width, whichresults in N A = N B . In a more practical configuration (Figure 1a), the vertical wallsdon’t touch the IDA, and the most exterior bands have full width. In this last case,we will consider N A +1 = N B in order to take advantage of the symmetry produced atboth ends of the IDA. However, the results can be extrapolated to the case N A = N B without much difficulty if desired. METHODS (a) Full view. w A w B W AB Hu A u B (b) Cross section. Figure 1:
Practical configuration of an interdigitated array (IDA) in a rectangularbox of finite height H . The IDA consists of two arrays of N A and N B bands each(shown N A + 1 = N B = 4 ) at which the potentials u A and u B are applied. Thetwo-dimensional domain can be approximated as an assembly of several interior cells (four shown) and two exterior cells at the ends of the IDA.The problem consists of finding the potential distribution u inside the rectangularbox, and the flux density j and flux i at the IDA, when we apply potentials u A and u B at the arrays A and B respectively. If the IDA fits exactly within an ideal domain, the rectangular box can be reduced totwo dimensions, due to the symmetry along the bands. This reduced domain can befurther modeled as an assembly of two-dimensional cells, due to the vertical symmetryboundaries at the center of each band.For a more practical case, the domain can be reduced to two dimensions (Figure1b) when the bands are sufficiently long, which ensures that fringing effects at theends of each band can be neglected. This can be further reduced to an assemblyof two-dimensional cells, if the number of bands is sufficiently large, which ensuresthat fringing effects at both ends of the IDA can also be neglected. In this case,the symmetry boundaries between cells remain vertical at the interior of the IDA.However, they will bend near the ends of the IDA, where our approximation of verticalsymmetry boundaries will produce some error. Choosing the last vertical symmetryboundary to be located one band before each end of the IDA helps to better accountfor bendings of this boundary at the last band. Regarding the bands, they are allowedto have a thickness or elevation different from zero, as long as it is sufficiently smallcompared with their widths.Thus the IDA in a practical domain consists of an assembly of several interiorcells and two exterior cells at the ends of the IDA (Figure 1b), and the solution ofLaplace’s equation will hold accurately. Whereas, the IDA in an ideal domain consistsonly of an assembly of interior cells , and the solution of Laplace’s equation will beexact.
In practice, directly solving Laplace’s equation for the interior and exterior cells iscomplicated. This is due to alternated boundaries of potential and insulation atthe bottom of the cells. However, by using domain transformations, it is possibleto rearrange these boundary conditions, so they can be placed at different walls ina transformed cell. Complex conformal transformations are useful, since they leaveLaplace’s equation (as well as potential and insulation/symmetry boundaries) invari-ant under domain changes [12, §5.7] [13, §6]. In particular,
Jacobian elliptic functions are of interest, since they conformally map a square domain into the upper half-plane [12, §2.5] [14, §22.18.ii]. Also the
Möbius functions are important, since theyare able to reorganize the upper half-plane, by mapping into itself [12, §2.3].
METHODS xz rr o r a r α r β r b r l r m r n WH w a w b v n v o v a v α v β v b v l v m − /k r − /k r vω a ω α ω β ω b ω l ω m ω n ω o −∞ /k ρ ωξζ ρρ α ρ β ρ b ρ l ρ m ρ n ρ o ρ a Figure 2: Complex transformation of the generic cell of IDA domain r = ( x, z ) into the conformal parallel-plates domain ρ = ( ξ, ζ ) , by using the auxiliary complexdomains v and ω . First, the transformation r → v maps the interior of the cell into theupper half-plane by using a Jacobian elliptic function. Later, the transformation v → ω reorganizes the structure of the upper half-plane through a Möbius transformation.Finally, the transformation ω → ρ maps the upper half-plane to the interior of theparallel-plates cell, by using a composition of squared root and an inverse Jacobianelliptic function.Let r = ( x, z ) = x + i z designate the coordinates of a generic cell (Figure 2) where i = − . This generic cell represents simultaneously interior and exterior cells ofan IDA, depending on how its parameters are chosen. In case of an interior cell , onechooses r o = r a = 0 , w a = w A / , w b = w B / , r b = W AB = r l = W (3a)whereas, in case of an exterior cell r o = r a = 0 , w a = w A / , w b = w B , r b = W AB + w B / < r l = W (3b)This generic domain can be conformally transformed by ρ = T ρr ( r ) = T ρω ◦ T ωv ◦ T vr ( r ) (4a)passing through the auxiliary domains v and ωv = T vr ( r ) = − cd (cid:18) K ( k r ) 2 r W , k r (cid:19) (4b) ω = T ωv ( v ) = ( v − v α )( v − v a ) ( v β − v a )( v β − v α ) (4c)into a parallel-plates domain of coordinates ρ = ( ξ, ζ ) = ξ + i ζ ρ = T ρω ( ω ) = 1 K ( k ρ ) arcsn( √ ω , k ρ ) (4d)where Laplace’s equation has a simple solution. Here, the moduli and their comple-ments are given by q r = Q ( k r ) = exp (cid:18) − π HW (cid:19) , q (cid:48) r = Q ( k (cid:48) r ) = exp (cid:18) − π W H (cid:19) (5a) k ρ = ( v b − v a )( v β − v α )( v b − v α )( v β − v a ) , k (cid:48) ρ = ( v b − v β )( v α − v a )( v b − v α )( v β − v a ) (5b) RESULTS r a , r α , r β and r b on the boundary of the IDA domain r are transformedto the auxiliary domain v as v a = − cd (cid:18) K ( k r ) 2 r a W , k r (cid:19) , v α = − cd (cid:18) K ( k r ) 2( r a + w a ) W , k r (cid:19) (6a) v b = − cd (cid:18) K ( k r ) 2 r b W , k r (cid:19) , v β = − cd (cid:18) K ( k r ) 2( r b − w b ) W , k r (cid:19) (6b)and finally to the parallel-plates domain ρ as ρ a = ξ a + i ζ a = 0 + i K (cid:48) ( k ρ ) K ( k ρ ) , ρ α = ξ α + i ζ α = 0 + i (7a) ρ b = ξ b + i ζ b = 1 + i K (cid:48) ( k ρ ) K ( k ρ ) , ρ β = ξ β + i ζ β = 1 + i (7b)For a detailed construction of this transformation see Supplementary information §S1 .Several special functions and definitions have been used for the conformal trans-formation T ρr : cd( z , k ) and arcsn( z , k ) correspond to Jacobian elliptic functions [14,Equations (22.2.8) and (22.15.12)] analogous to their circular counterparts cos() and arcsin() respectively, K ( k ) and K (cid:48) ( k ) correspond to the complete elliptic integral ofthe first kind and its associated function respectively [14, Equations (19.2.4), (19.2.8)and (19.2.9)], Q ( k ) corresponds to the elliptic nome function [14, Equation (22.2.1)],and k and k (cid:48) = √ − k correspond to the elliptic modulus and its complement, whichare used as parameters for elliptic functions. By using the transformation of the last section, Laplace’s equation and its boundaryconditions (potential and insulation/symmetry) can be written in the parallel-platesdomain ρ = ( ξ, ζ ) as ∂ u∂ξ ( ξ, ζ ) + ∂ u∂ζ ( ξ, ζ ) = 0 (8a) u (0 , ζ ) = u A , u (1 , ζ ) = u B , ∂ u∂ζ ( ξ,
0) = ∂ u∂ζ ( ξ, ζ a ) = 0 (8b)where ζ a = (cid:61) ρ a = K (cid:48) ( k ρ ) /K ( k ρ ) is the imaginary part of ρ a in Equations (7).This equation has a straightforward solution u ( ξ, ζ ) = u A + [ u B − u A ] ξ , whichcorresponds to a linear interpolation of the potentials at each plate. This result mustbe carried back to the generic IDA domain r = ( x, z ) to obtain the final solution u ( x, z ) = u A + [ u B − u A ] ξ ( x, z ) (9)where ξ ( x, z ) = (cid:60) T ρr ( x, z ) corresponds to the real part of the conformal transformationin Methods §3.2 . See Figure 3 for potential distributions at interior and exterior cells .Since the IDA bands are equipotential surfaces, the flux density at the bands mustbe perpendicular to them, and therefore it must possess only vertical component j ( x ) = − γ ∂ u∂z ( x,
0) = − γ [ u B − u A ] ∂ ξ∂z ( x,
0) = γ [ u B − u A ] (cid:61) ∂ ρ ∂ r ( x ) (10a)where the last equality is due to the Cauchy-Riemann equations [13, Theorem 3.2] ∂ ξ∂z = − ∂ ζ∂x = −(cid:61) ∂ ρ ∂x = −(cid:61) ∂ ρ ∂ r (10b) RESULTS x/W AB z / W A B . . . . . . ξ ( x, z ) 0.0 0.5 1.0 1.5 2.0 x/W AB z / W A B . . . . . . . ξ ( x, z ) (a) Tall domain ( H/W AB = 1 ). x/W AB z / W A B . . . . . . ξ ( x, z ) 0.0 0.5 1.0 1.5 2.0 x/W AB z / W A B . . . . . . . ξ ( x, z ) (b) Shallow domain ( H/W AB = 0 . ). Figure 3: Normalized potential distribution ξ ( x, z ) = [ u ( x, z ) − u A ] / [ u B − u A ] forequal band widths w A = w B = 0 . W AB and different domain heights H/W AB . Leftand right figures correspond to interior ( W = W AB ) and exterior ( W = 2 W AB ) cellsrespectively.After some algebraic calculations (see Supplementary information §S1.4 ), we obtainedthe derivative of the conformal transformation ρ ( r ) = T ρr ( r ) ∂ ρ ∂ r = i W K ( k r ) K ( k ρ ) ( v b − v α ) / ( v β − v a ) / ( v − v α ) / ( v − v β ) / (1 − v ) / (1 − k r v ) / ( v − v a ) / ( v b − v ) / (10c)where v is given in Equation (4b) and the parameters are given in Equations (3), (5)and (6). See Figure 4 for flux densities at the bottom of interior and exterior cells .Finally, we can obtain the flux per band by integrating the flux density on thesurface of a band E ∈ { A, B } i E = (cid:90) E j ( x ) L d x = − Lγ [ u B − u A ] (cid:90) E ∂ ξ∂z ( x,
0) d x (11a)which can be simplified using the Cauchy-Riemann identities in Equation (10b) − (cid:90) E ∂ ξ∂z ( x,
0) d x = (cid:61) (cid:90) E ∂ ρ ∂x ( x,
0) d x = (cid:61) (cid:90) E ∂ ρ ( r ) (11b)Due to symmetry, we can take twice the integral over half band for the interior cells (cid:61) (cid:90) A ∂ ρ ( r ) = 2 (cid:61) ρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) r α r a , (cid:61) (cid:90) B ∂ ρ ( r ) = 2 (cid:61) ρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) r b r β (11c)For the exterior cells , we just combine both cells to obtain a full band of A , while for B only one exterior cell suffices (cid:61) (cid:90) A ∂ ρ ( r ) = 2 (cid:61) ρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) r α r a , (cid:61) (cid:90) B ∂ ρ ( r ) = (cid:61) ρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) r b r β (11d) RESULTS x/W AB -4-2024 W = ρ ( x ) / r H/W AB = 2 . H/W AB = 1 . H/W AB = 0 . H/W AB = 0 . x/W AB -8-4048 W = ρ ( x ) / r H/W AB = 2 . H/W AB = 1 . H/W AB = 0 . H/W AB = 0 . Figure 4: Normalized flux density W (cid:61) ∂ ρ ( x ) /∂ r = W j ( x ) / ( γ [ u B − u A ]) for equalband widths w A = w B = 0 . W AB and different domain heights H/W AB . Left andright figures correspond to interior ( W = W AB ) and exterior ( W = 2 W AB ) cellsrespectively.Combining Equations (11) with (7) we have − i int A = i int B = L K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) int · γ [ u B − u A ] (12a)for bands in the interior cells , and − i ext A = 2 i ext B = L K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) ext · γ [ u B − u A ] (12b)for bands in the exterior cells . The value of k ρ can be obtained from Equations (3),(5) and (6) for interior and exterior cells . See Figure 5 for flux per band as a functionof the dimensions of interior and exterior cells . The flux depends on elliptic functions which are normally evaluated using specializedsoftware. However, expressions involving elementary functions are desirable for usingany general purpose software or calculator. For this reason, we simplified the expres-sion for the flux by approximating the ratio K (cid:48) ( k ρ ) /K ( k ρ ) . In particular, we chose tosimplify that of an interior cell , since this is the only kind of cell in an ideal domain ,and also the dominant one in case of a more practical domain (see Discussion §5.3 ).As shown in [7, Equation (29)] and [8, Equations (4) and (5)], this ratio can beapproximated for sufficiently small moduli k ρ or k (cid:48) ρ by K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) k ρ ≈ ≈ − π ln (cid:32) k ρ (cid:33) (13a) K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) k (cid:48) ρ ≈ ≈ − π ln (cid:32) k (cid:48) ρ (cid:33) − (13b)In our case, the moduli k ρ and k (cid:48) ρ in Equations (3a), (5b) and (6) also depend onelliptic functions k ρ = 2 (cid:34) cd (cid:18) K ( k r ) w A W AB , k r (cid:19) + cd (cid:18) K ( k r ) w B W AB , k r (cid:19)(cid:35)(cid:34) (cid:18) K ( k r ) w A W AB , k r (cid:19)(cid:35) (cid:34) (cid:18) K ( k r ) w B W AB , k r (cid:19)(cid:35) (14a) RESULTS w A /W AB = w B /W AB H / W A B . . . . . K ( k ρ ) /K ( k ρ ) 0.0 0.5 1.0 1.5 2.0 w A /W AB w B / W A B . . . . H/W AB = 0.3 0.0 0.5 1.0 1.5 2.0 w A /W AB w B / W A B . . . . H/W AB = 1.0 (a) Interior cell ( W/W AB = 1 ). w A /W AB = w B /W AB H / W A B . . . . . K ( k ρ ) /K ( k ρ ) 0.0 0.5 1.0 1.5 2.0 w A /W AB w B / W A B . . . . H/W AB = 0.3 0.0 0.5 1.0 1.5 2.0 w A /W AB w B / W A B . . . . H/W AB = 1.0 (b) Exterior cell ( W/W AB = 2 ). Figure 5: Selected two-dimensional slices of the normalized flux per band K (cid:48) ( k ρ ) /K ( k ρ ) = ( i A /L ) / ( γ [ u A − u B ]) as a function of the cells dimensions. Plotsconsidering exterior cells with ≤ W/W AB ≤ presented almost no variations (shownat Supplementary information §S2 ). k (cid:48) ρ = (cid:34) − cd (cid:18) K ( k r ) w A W AB , k r (cid:19)(cid:35)(cid:34) (cid:18) K ( k r ) w A W AB , k r (cid:19)(cid:35) (cid:34) − cd (cid:18) K ( k r ) w B W AB , k r (cid:19)(cid:35)(cid:34) (cid:18) K ( k r ) w B W AB , k r (cid:19)(cid:35) (14b)and therefore must be approximated too. To accomplish this, we found alternativerepresentations, which must be easier to approximate than their original counterparts.The one for k ρ depends on the gap g = W AB − w A = W AB − w B between bands ofequal width w A = w B k ρ = 4 sn (cid:18) K ( k r ) gW AB , k r (cid:19)(cid:34) (cid:18) K ( k r ) gW AB , k r (cid:19)(cid:35) (15a)due to quarter period translation cd( u + K ( k ) , k ) = − sn( u , k ) [14, Table 22.4.3] anddoble rotation of the argument sn( − u , k ) = sn( i u , k ) = − sn( u , k ) [14, Table 22.6.1].The one for k (cid:48) ρ can be expressed using different band widths k (cid:48) ρ = k (cid:48) r sd (cid:18) K ( k r )2 w A W AB , k r (cid:19) cn (cid:18) K ( k r )2 w A W AB , k r (cid:19) sd (cid:18) K ( k r )2 w B W AB , k r (cid:19) cn (cid:18) K ( k r )2 w B W AB , k r (cid:19) (15b)due to [1 − cd(2 u , k )] / [1 + cd(2 u , k )] = k (cid:48) sd( u , k ) / cn( u , k ) [15, Equations (1.10)and (4.1)]. Here the modulus k r is obtained from q r = Q ( k r ) = exp( − π H/W AB ) DISCUSSION sn , sd , cn and cd correspond to Jacobianelliptic functions analogous to their circular counterparts sin and cos [14, §22.2].Finally, we approximated k ρ and k (cid:48) ρ for the cases of tall and shallow domains. Incase of tall domains ( H is large) q r → + , k r → + and k (cid:48) r → − , which causes K ( k r ) → π/ , sn → sin , sd → sin and cd → cos [14, Equation (19.6.1) and Table22.5.3] or [16, Equations (10)] k ρ (cid:12)(cid:12)(cid:12)(cid:12) H → + ∞ g ≈ = 4 sin (cid:18) π gW AB (cid:19) (cid:34) (cid:18) π gW AB (cid:19)(cid:35) − ≈ (cid:18) π gW AB (cid:19) (16a) k (cid:48) ρ (cid:12)(cid:12)(cid:12)(cid:12) H → + ∞ w A ≈ w B ≈ = tan (cid:18) π w A W AB (cid:19) tan (cid:18) π w B W AB (cid:19) ≈ (cid:18) π w A W AB (cid:19) (cid:18) π w B W AB (cid:19) (16b)In case of shallow domains ( H is small) q (cid:48) r → + and k (cid:48) r → + , which causes sn → tanh , sd → sinh and cn → sech [14, 22.5.4]. However, to avoid the divergence in K ( k r ) → + ∞ [14, Equation (19.6.1)], it is more convenient to use the approximationsin [16, Equations (11)] as detailed in Supplementary information§S3 , which lead to k ρ (cid:12)(cid:12)(cid:12)(cid:12) H ≈ g ≈ ≈ (cid:18) π W AB H − ln √ (cid:19) tanh (cid:18) π gH (cid:19) (17a) k (cid:48) ρ (cid:12)(cid:12)(cid:12)(cid:12) H ≈ w A ≈ w B ≈ ≈
16 e − πW AB /H sinh (cid:18) π w A H (cid:19) sinh (cid:18) π w B H (cid:19) (17b)See also Table 2 at Discussion §5.5 to know the dimensions for which these approxi-mations hold.
The results for conformal transformation, potential distribution, flux density and fluxagree with experimentally validated expressions [7, 8, 11] and simulations [8, 17, 18]found in the literature. These expressions and simulations correspond to the casesof tall domain (with arbitrary band widths) and shallow domain (with equal bandwidths), which are particular cases of the results obtained in this report.In case of tall domain (large H ), the conformal transformation ( Methods §3.2 )agrees with that in [7, Equations (9), (10), (14), (15) and (19)]. The potential dis-tribution in Equation (9) corresponds to that in [7, Equations (20) and (21)] , whichdepends directly on the real part of the conformal transformation. The flux densityin Equations (10) agrees with that in [7, Equations (19) and (26)], which dependsdirectly on the imaginary part of the conformal transformation’s derivative. Thisderivative also agrees with that in [7, Equations (17) and (27)] by considering theidentity cos( α + β ) + cos( α − β ) = (1 + cos(2 α )) / (1 + cos(2 β )) / . The flux per bandin Equations (12) and their moduli in Equations (16) correspond to that in [7, Equa-tions (15) and (28)] and [8, Equations (2), (3) and (6)]. All these results agree asymp-totically when H → + ∞ , because k r → + , K ( k r ) → π/ and cd( · , k r ) → cos( · ) .In case of shallow domain, the flux per band in Equations (12) was evaluated andcompared with its simulated counterparts in [17, Figure 7a] [18, Figure 7a] and withthe experimentally validated expression in [11, Equation (13)]. All these results agreeas shown in the comparisons made in Figure 6. Note that [7, H in Equation (21)] has a typo, and the + sign, at the middle of the expression,should be replaced by a − sign. For facilitating the comparisons [7, c ∗ H in Equations (6) and (21)]and [7, nF D in Equation (26)] correspond to u B − u A and γ respectively. The parameter [7, p inEquation (15)] corresponds to k ρ , which is used instead of the modulus k ρ as the argument for allelliptic functions. DISCUSSION H/W AB K ( k ρ ) / K ( k ρ ) H/W AB K ( k ρ ) / K ( k ρ ) Figure 6: Normalized flux per band K (cid:48) ( k ρ ) /K ( k ρ ) = ( i A /L ) / ( γ [ u A − u B ]) for an interior cell with equal band widths w A /W AB = w B /W AB ∈ { . , . , . , . , . } (widths correspond to solid lines in order from bottom to top). Solid lines: Theoreticalexpresion in Equation (12a). Blue [18, Figure 7a] and orange [17, Figure 7a] dots:Literature simulations. Green dots: Theoretical expression [11, Equation (13)]. See Supplementary information §S4 for details on the data from the literature.
The conformal transformation T ρr in Methods §3.2 depends only on relative dimensionsand positions in the generic cell instead of absolute ones. The same occurs with itscomplex derivative in Equation (10c) when this is normalized as
W ∂ ρ /∂ r . Thissignifies that the potential u , normalized flux density W j and normalized flux i A /L in Results §4.1 will depend on relative dimensions of the IDA domain, and therefore,they will be shared among all IDA domains with proportional geometries. This isbecause these three quantities depend directly on the real/imaginary part of either T ρr or W ∂ ρ /∂ r .As general behavior, the flux in interior and exterior cells of an IDA domainincreases with both: the relative width of the bands ( w A /W AB and w B /W AB ) andthe relative height of the domain ( H/W AB ) (Figures 5 and 6). In case of an exteriorcell , the flux also increases with its relative width W/W AB , however any changeappears to be comparatively small for W/W AB (cid:38) . For a fixed domain height, theisolines of flux are symmetric with respect to the main diagonal of the plot (in caseof an interior cell ) and tend to collapse towards the anti-diagonal of the plot as thedomain height decreases (in case of interior and exterior cells ).In particular, for very shallow domains ( < H/W AB (cid:46) . ), the flux increasesslowly as the bands width increases (Figures 5 and 6). This is because the flux densitynear the center of each band is low (Figure 4), and therefore it cannot contribute withenough area for the flux, unless the bands are very wide. This low flux density nearthe center of the bands has its origin in the potential distribution (Figure 3), whichpresents almost no vertical gradient due to truncation by the roof of the domain.Most of the flux in this case comes from the flux density spikes at the edges of eachband (Figure 4), which originate mainly due to the horizontal gradient of potentialgenerated between neighboring bands (Figure 3). On the other hand, the flux forvery shallow domains increases quickly as the domain height increases (Figures 5 and6). This is because the flux density near the center of each band increases rapidly(Figure 4), boosted by a quick development of the vertical gradient of potential, asthe domain height increases (Figure 3).For shallow domains ( . (cid:46) H/W AB (cid:46) ), the flux increases more quickly as theband widths increases (Figures 5 and 6). This is because the flux density near thecenter of each band is higher (Figure 4), due to the vertical gradient of potential(Figure 3) that has been developed, and therefore it can contribute more area for theflux. On the other hand, the flux for shallow domains increases more slowly as thedomain height increases (Figures 5 and 6). The reason for this behavior is that theflux density near the center of each band also increases more slowly with the domainheight (Figure 4), since an important part of the vertical gradient of potential has DISCUSSION w A /W AB = w B /W AB H / W A B . . . . . . Relative error (a) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 0.3 (b) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 1.0 (c) Min: . . Max: . . Figure 7: Selected two-dimensional slices of the relative error i ext A /i int A − as a functionof the cells dimensions. Case of exterior cell with W/W AB = 2 . Plots considering exterior cells with ≤ W/W AB ≤ presented almost no variations (shown at Sup-plementary information §S3 ).been already developed (Figure 3).For tall domains (
H/W AB (cid:38) ), the flux (Figures 5 and 6) and flux density(Figure 4) approach their maximum limit and become independent of the domainheight (they depend only on the bands width). This is because the domain is so tallthat the potential distribution becomes uniform far from the surface of the IDA, anda fully developed potential gradient is present in the domain (Figure 3). In this case,the flux between neighboring bands is not affected by the roof of the domain andapproaches that of semi-infinite geometries.Similar results, in terms of geometry dependence, have also been found previouslythrough simulations for planar IDA configurations [17–21] as well as for elevated IDAconfigurations [19, 20, 22]. The total flux passing through an IDA is an important quantity to consider, since inmany cases it is possible to measure it experimentally.In case of an ideal domain , we have an assembly consisting only of interior cells with N A = N B , thus the total flux through each array is given by − N A i int A = N B i int B (18a)For a practical domain , we have an assembly of interior cells and two exterior cells such that N A + 1 = N B , thus the total flux through each array is given by − I A = I B (cid:117) − ( N A − i int A − i ext A = ( N B − i int B + 2 i ext B (18b)If we compare this total flux with the one of an ideal domain (both with the samenumber of bands for the array A ) I A N A i int A − (cid:117) N A (cid:32) i ext A i int A − (cid:33) (18c)we can see that the relative error of I A with respect to its ideal domain counterpart N A i int A decreases towards zero as the number of bands N A increases. This showsthat the flux through an interior cell is dominant compared to that of an exteriorcell , especially as the number of bands per array increases. This suggests that it ispossible to replace the more complex total flux of a practical domain by that of an ideal domain .In fact, for most practical cases, the relative error between the flux at exterior and interior cells is | i ext A /i int A − | ≤ . (Figure 7), and therefore the relative error DISCUSSION practical and ideal domains is | I A /N A i int A − | (cid:46) . /N A .This means that the number of digits in agreement for the total flux of both domainsrelates to the order of magnitude of the number of bands per array. Moreover, forthe cases where H/W AB ≤ w A /W AB = w B /W AB , we have | i ext A /i int A − | < . ,and then it is possible to obtain an additional digit in agreement for the flux of bothdomains.Considering the previous argument, we present some examples of total flux for thecase of an ideal domain . See Equation (12a).In case of an IDA of electric resistance R , an electric current I A is produced when applying a difference of electric potentials V A − V B to a medium of electric conductivity σ R = I A V A − V B = N A L K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) int σ (19a)In case of an IDA of capacitance C , an electric charge Q A (total flux of electricdisplacement) is accumulated when applying a difference of electric potentials V A − V B to a medium of absolute permittivity (cid:15)C = Q A V A − V B = N A L K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) int (cid:15) (19b)Note that if the IDA presents resitive and capacitive effects simultaneouly, then wehave RC = (cid:15)/σ (this relation is not limited to an IDA, but holds in general [10]). Incase of electrochemistry, a faradaic current I A is produced when there is a differenceof concentrations c A − c B of electrochemical species at the bands of the IDA I A = N A L K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) int F n e D [ c A − c B ] (19c)Here F corresponds to the Faraday’s constant, and the electrochemical species ex-change n e electrons in the faradaic process and have a diffusion coefficient D in themedium. The amount of material used for fabricating an IDA may be of importance in processesof additive manufacturing. This amount is proportional to m = w A /W AB + w B /W AB and can be minimized with ease in case of an ideal domain or a practical domain witha large amount of bands (Figure 5a).For fixed domain height and a desired flux, the minimum is obtained when w A = w B . This can be seen by selecting a desired isoline of flux, and intersecting it withthe anti-diagonal m = w A /W AB + w B /W AB . This produces two intersection points,which are symmetric with respect to the diagonal line w A /W AB = w B /W AB . Theminimum is obtained by moving the anti-diagonal towards the origin, which makesthe two intersection points on the isoline converge to a single point, thus obtaining w A = w B . The search for the regions of validity for Equations (13) (16) and (17) was donegraphically, and therefore, it was more convenient to consider bands of equal width,such that the resulting plots were two-dimensional.Approximations of K (cid:48) ( k ρ ) /K ( k ρ ) in Equations (13) and (16) for tall domains( H/W AB > ) were obtained first by Aoki and colleagues, for the case of wide bands( g ≈ ) [7, Eq. (32)], and later by Morf and colleagues, for the case of narrowbands [8, Eqs. (2), (6) and (7)]. The approximation for wide bands becomes accurate(error less than ± ) for bands of relative width > . . In the case of narrow bands,the approximation has an error less than ± for bands of relative width < . . CONCLUSION w A /W AB = w B /W AB H / W A B - . - . - . - . - . w A /W AB = w B /W AB H / W A B . . . . . . w A /W AB = w B /W AB H / W A B - . - . - . - . - . w A /W AB = w B /W AB H / W A B - . - . - . . . . . . . Figure 8: Relative error of the approximated ratio K (cid:48) ( k ρ ) /K ( k ρ ) for an ideal domain with respect to its exact counterpart. Top and bottom: Errors for tall and shallowdomains respectively. Left and right: Errors for narrow and wide bands respectively.In red: straight lines close to errors of ± .The regions of approximation for both sizes of bands overlap, thus covering all casesof H/W AB > with a relative error less than ± . See Table 2 for a summary.Approximations of K (cid:48) ( k ρ ) /K ( k ρ ) in Eqs. (13) and (17) for shallow domains( H/W AB ≤ ) are results that have not been published before. The approxima-tion for wide bands becomes accurate (error approximately less than ± ) for bandsof relative width approximately larger than − . H/W AB . In the case of narrowbands, the approximation has an error less than ± for bands of relative widthapproximately smaller than − . H/W AB . The regions of approximation for bothsizes of electrodes overlap, thus covering all cases of H/W AB ≤ with a relative errorless than ± . See Table 2 for a summary. It is possible to transform ideal and practical domains containing an IDA into aparallel-plates domain. In the latter, the solution of Laplace’s equation is simple,and corresponds to a linear interpolation of the potentials applied at both plates.This solution was transformed back into the IDA domain, leading to an analyticalexpression for the potential distribution (based on elliptic functions and integrals)from which flux density and flux were derived.All three expressions depend on relative dimensions of the IDA domain instead ofabsolute ones, meaning that the same results hold for proportionally similar IDAs.The flux density and flux at the IDA increase as the bands width and domain heightincrease, which correlates with an increase of gradients in the potential distribution.
ACKNOWLEDGEMENTS
Case of tall domains ( H/W AB > ) w E /W AB > . Equations (13a), (16a) w E /W AB < . Equations (13b), (16b)
Case of shallow domains ( H/W AB ≤ ) w E /W AB + 0 . H/W AB (cid:38) Equations (13a), (17a) w E /W AB + 0 . H/W AB (cid:46) Equations (13b), (17b)Table 2: Regions where the approximations hold with a relative error less than ± .Here w E := w A = w B . See Supplementary information §S5 for details.Their behavior approaches that of an IDA in a semi-infinite domain, once the do-main is sufficiently tall (approximately when the domain height is greater than theseparation between centers of consecutive bands).Among the three quantities, the flux plays an important role, since in many casesit is accessible for experimental measurement. Therefore, its analytical expressioncan be used during the design process of an IDA or as a benchmark of performancein measurements. In case of an ideal domain , its analytical expression holds exactly,while that of a practical domain approaches the latter as the number of bands perarray increases (the number of digits in agreement between both expressions correlateswith the order of magnitude of the number of bands per array). This is because theflux in interior cells of an IDA becomes dominant over that in exterior cells , makingfringing effects negligible.For an IDA behaving like one in an ideal domain , we found that bands of equalwidth minimize its total surface for any desired flux. Also, approximations for theexact flux were found. Elementary functions were used to approximate the cases oftall and shallow domains respectively, producing expressions which are accurate witha relative error smaller than ± with respect to their exact counterpart. Whenthese approximations are used in combination, they cover all possibilities of interestfor IDAs in confined domains. The authors deeply appreciate the aid and comments of Mithran Somasundrumand Sirimarn Ngamchana, which helped to improve the quality of this manuscript. cfgy gratefully acknowledges
Petchra Pra Jom Klao Ph. D. scholarship (Grant No.28/2558),
King Mongkut’s University of Technology Thonburi . The authors acknowl-edge the financial support provided by
King Mongkut’s University of TechnologyThonburi through the
KMUTT 55th Anniversary Commemorative Fund , and KingMongkut’s University and the
Research Network NANOTEC (RNN) program (grantno. P1851883) of the
National Nanotechnology Center (NANOTEC), NSTDA, Min-istry of Science and Technology, Thailand. ( pdf ) Details on calculations, additional figures and tables. ( zip ) Python scripts fornumerical calculations and plots [23]. [1] A. Mamishev, K. Sundara-Rajan, F. Yang, Y. Du, and M. Zahn, “Interdigitalsensors and transducers,”
Proceedings of the IEEE no. 5, (May, 2004)808–845, CiteSeerX:10.1.1.473.5989 .6[2] O. Atalay, “Textile-based, interdigital, capacitive, soft-strain sensor for wearableapplications,”
Materials no. 5, (May, 2018) 768. CC:BY/4.0.[3] S. Brosel-Oliu, N. Abramova, N. Uria, and A. Bratov, “Impedimetrictransducers based on interdigitated electrode arrays for bacterial detection – areview,” Analytica Chimica Acta (Dec., 2019) 1–19.[4] S. Szunerits and L. Thouin, “Microelectrode arrays,” in
Handbook ofElectrochemistry , C. G. Zoski, ed., ch. 10, pp. 391 – 428.[5] M. F. El-Kady and R. B. Kaner, “Scalable fabrication of high-power graphenemicro-supercapacitors for flexible and on-chip energy storage,”
NatureCommunications no. 1, (Feb., 2013) . OA.[6] S. Miserere and A. Merkoçi, “Microfluidic electrochemical biosensors:Fabrication and applications,” in Lab-on-a-chip devices and micro-total analysissystems , J. Castillo-León and W. E. Svendsen, eds., ch. 6, pp. 141–160.Springer International Publishing, 2015.[7] K. Aoki, M. Morita, O. Niwa, and H. Tabei, “Quantitative analysis of reversiblediffusion-controlled currents of redox soluble species at interdigitatedgitatedarray electrodes under steady-state conditions,”
Journal of ElectroanalyticalChemistry and Interfacial Electrochemistry no. 2, (Dec., 1988) 269–282. http://asura.apphy.u-fukui.ac.jp/~aoki/publication/PDF/JEC256_269.pdf .[8] W. E. Morf, M. Koudelka-Hep, and N. F. de Rooij, “Theoretical treatment andcomputer simulation of microelectrode arrays,”
Journal of ElectroanalyticalChemistry no. 1, (May, 2006) 47–56, rérodoc:16975 .[9] J. Wei, “Distributed capacitance of planar electrodes in optic and acousticsurface wave devices,”
IEEE Journal of Quantum Electronics no. 4, (Apr.,1977) 152–158.[10] W. Olthuis, W. Streekstra, and P. Bergveld, “Theoretical and experimentaldetermination of cell constants of planar-interdigitated electrolyte conductivitysensors,” Sensors and Actuators B: Chemical no. 1-3, (Mar., 1995) 252–256, CORE:11453584 . http://purl.utwente.nl/publications/15064 .[11] R. Igreja and C. J. Dias, “Analytical evaluation of the interdigital electrodescapacitance for a multi-layered structure,” Sensors and Actuators A: Physical no. 2-3, (May, 2004) 291–301.[12] T. A. Driscoll and L. N. Trefethen,
Schwarz-Christoffel mapping , vol. 8.Cambridge University Press, 2002.[13] P. J. Olver, “Complex analysis and conformal mapping,” June, 2020. .[14] F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F.Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A.McClain, eds.,
NIST Digital Library of Mathematical Functions . June, 2020. http://dlmf.nist.gov/ . Release 1.0.27 of 2020-06-15.[15] B. C. Carlson, “Symmetry in c, d, n of jacobian elliptic functions,”
Journal ofMathematical Analysis and Applications no. 1, (Nov., 2004) 242 – 253,
CORE:82142693 .[16] J. D. Fenton and R. S. Gardiner-Garden, “Rapidly-convergent methods forevaluating elliptic integrals and theta and elliptic functions,”
The ANZIAMJournal (July, 1982) 47–58. http://johndfenton.com/Papers/Fenton82b-Elliptic+Theta-functions.pdf .7[17] J. Strutwolf and D. E. Williams, “Electrochemical sensor design using coplanarand elevated interdigitated array electrodes. a computational study,” Electroanalysis no. 2, (Feb., 2005) 169–177.[18] C. Guajardo, S. Ngamchana, and W. Surareungchai, “Mathematical modelingof interdigitated electrode arrays in finite electrochemical cells,” Journal ofElectroanalytical Chemistry (Sept., 2013) 19–29, arXiv:1602.06428[physics.chem-ph] .[19] J.-I. Heo, Y. Lim, and H. Shin, “The effect of channel height and electrodeaspect ratio on redox cycling at carbon interdigitated array nanoelectrodesconfined in a microchannel,”
Analyst (2013) 6404–6411.[20] J.-I. Heo, “Development of electrochemical sensors based on carboninterdigitated array nanoelectrodes for high current amplification,” Master’sthesis, Ulsan National Institute of Science and Technology, June, 2014. .[21] Y. Kanno, T. Goto, K. Ino, K. Y. Inoue, Y. Takahashi, H. Shiku, andT. Matsue, “SU-8-based flexible amperometric device with ida electrodes toregenerate redox species in small spaces,”
Analytical Sciences no. 2, (2014)305–309, CORE:196707925 . https://kanazawa-u.repo.nii.ac.jp/?action=pages_view_main&active_action=repository_view_main_item_detail&item_id=9205&item_no=1&page_id=13&block_id=21 .[22] E. D. Goluch, B. Wolfrum, P. S. Singh, M. A. G. Zevenbergen, and S. G.Lemay, “Redox cycling in nanofluidic channels using interdigitated electrodes,” Analytical and Bioanalytical Chemistry no. 2, (May, 2009) 447–56,
CiteSeerX:10.1.1.611.2411 . http://resolver.tudelft.nl/uuid:15761e03-0f57-4763-8452-c02f04ca5ac8 .[23] C. F. Guajardo Yévenes, “Elektrodo,” Jan., 2021. https://doi.org/10.5281/zenodo.4422476 . upplementary information: Analytical solution for two-dimensional Laplace’sequation in a shallow domain containing coplanarinterdigitated boundaries
Cristian F.
Guajardo Yévenes ID [email protected] Werasak
Surareungchai ID [email protected] January 7, 2021
S1 Transformation of the generic cell of IDA domainto parallel plates
Alternated potential and insulation boundaries at the bottom of the cells make itdifficult to obtain an analytical solution for Laplace’s equation. Nevertheless, it ispossible to arrange these alternated boundary conditions, so they can be placed atdifferent walls in a transformed cell, by using conformal transformations.Complex conformal transformations leave Laplace’s equation, as well as potential(Dirichlet) and non-flux (Neumann) boundary conditions, invariant under domainchanges [1, §5.7]. In particular, Jacobian elliptic functions sn() and cd() are of interest,since they conformally transform a square domain into the upper half-plane [1, §2.5] [2,§22.18.ii]. Möbius functions are also important, since they reorganize the upper half-plane, by mapping it into itself [1, §2.3].Now we take a generic cell (Figure S1) representing simultaneously interior and exterior cells of an IDA, depending on how its parameters are chosen. In case of an interior cell , one chooses r o = r a = 0 , w a = w A / , w b = w B / , r b = r l = W AB (S1a)Whereas, in case of an exterior cell r o = r a = 0 , w a = w A / , w b = w B , r b = W AB + w B / < r l = W (S1b)Let r = x + i z designate the coordinates for the generic cell of the IDA domain,where i = − . We aim to find a complex function ρ = T ρr ( r ) that can transform thecell from the IDA domain into a parallel-plates domain of coordinates ρ = ξ + i ζ , asshown in Figure S1.The transformation T ρr will be decomposed in three stages T ρr = T ρω ◦ T ωv ◦ T vr (S2)First, T vr will transform the IDA domain r to the upper half-plane v by using theelliptic function cd() . Later, T ωv will reorganize the upper half-plane v into the upperhalf-plane ω by using a Möbius function. Finally, T ρω will transform the upper half-plane ω into the parallel plates domain ρ by using a combination of square root andinverse elliptic function arcsn() . S1 S2 xz rr o r a r α r β r b r l r m r n WH w a w b v n v o v a v α v β v b v l v m − /k r − /k r vω a ω α ω β ω b ω l ω m ω n ω o −∞ /k ρ ωξζ ρρ α ρ β ρ b ρ l ρ m ρ n ρ o ρ a Figure S1: Complex transformation ρ = T ρr ( r ) from the generalized IDA domain r = ( x, z ) into the conformal parallel-plates domain ρ = ( ξ, ζ ) , by using the auxiliarycomplex domains v and ω . S1.1 From r to v domain Considering the special values of the function arcsn() [2, Table 22.5.1] , k ) (S3a) ± K ( k ) = arcsn( ± , k ) (S3b) ± K ( k ) + i K (cid:48) ( k ) = arcsn( ± /k, k ) (S3c) i K (cid:48) ( k ) = arcsn( ∞ , k ) (S3d)one can construct a function ( T vr ) − that maps the upper half-plane of v into the IDAdomain r , as shown in Figure S1. The scale and translation of the function arcsn() must be chosen such that v o = − is mapped to r o = 0 and v l = 1 is mapped to r l = W , which leads to r = ( T vr ) − ( v ) = W (cid:20) K ( k r ) arcsn( v , k r ) (cid:21) (S4a) v = T vr ( r ) = sn (cid:18) K ( k r ) 2 r W − K ( k r ) , k r (cid:19) = − cd (cid:18) K ( k r ) 2 r W , k r (cid:19) (S4b)due to the quarter- and half-period properties [2, Table 22.4.3] sn( u − K ( k ) , k ) = sn( u + 3 K ( k ) , k ) = sn( u + K ( k ) + 2 K ( k ) , k ) = − cd( u , k ) (S5)where K ( k ) corresponds to one period.The appropriate modulus k r can be obtained by forcing v n = − /k r to be mappedto the upper-left corner r n = i H r n = ( T vr ) − ( v n ) ⇔ i H = i W K (cid:48) ( k r ) K ( k r ) This leads to a lattice parameter τ r [2, §20.1] τ r = i K (cid:48) ( k r ) K ( k r ) = i HW (S6)or, equivalently, to the elliptic nome and its associated counterpart q r = Q ( k r ) = exp (cid:18) − π HW (cid:19) (S7a) q (cid:48) r = Q ( k (cid:48) r ) = exp (cid:18) − π W H (cid:19) (S7b) S3by applying the nome function Q ( k ) , which is defined by the relations [3, §VI.3Equation (16)] [2, Equations (19.2.9) and (22.2.1)] ln Q ( k ) − π = − π ln Q ( k (cid:48) ) = K (cid:48) ( k ) K ( k ) (S8)The points v o = − , v l = 1 and v n = − /k r , define completely the function T vr in Equation (S4b), which transforms the IDA domain r into the upper half-plane v as shown in Figure S1. Thus all points of interest v a , v α , v β and v b are obtained byevaluating Equation (S4b) v a = − cd (cid:18) K ( k r ) 2 r a W , k r (cid:19) , v α = − cd (cid:18) K ( k r ) 2( r a + w a ) W , k r (cid:19) (S9a) v b = − cd (cid:18) K ( k r ) 2 r b W , k r (cid:19) , v β = − cd (cid:18) K ( k r ) 2( r b − w b ) W , k r (cid:19) (S9b) S1.2 From v to ω domain One can use a Möbius transformation T ωv for reorganizing the upper half-plane of v into the upper half-plane of ω . See Figure S1. This transformation is constructedsuch that: (i) it maps the band interval v ∈ [ v a , v α ] to the negative real axis of ω and the band interval v ∈ [ v β , v b ] to the real interval ω ∈ [1 , ω b ] , and (ii) it ensuresthe mapping of the upper half-plane of v into the upper half-plane of ω .Condition (i) can be satisfied by choosing the Möbius function as below ω = T ωv ( v ) = ( v − v α )( v − v a ) ( v β − v a )( v β − v α ) (S10a)since this maps v a , v α and v β into ω a = ∞ , ω α = 0 and ω β = 1 . Condition (ii) canbe ensured by rewritting T ωv as a composition of translations, rotations/scalings, andinversion ω = T ωv ( v ) = (cid:20) v a − v α v − v a (cid:21) (cid:34) v β − v a v β − v α (cid:35)(cid:124) (cid:123)(cid:122) (cid:125) p (S10b)The fact that v a − v α < and p > ensures that the upper half-plane of v is mappedto the upper half-plane of ω . See [4, Equation (5.7.3)], [3, §V.2 Equations (6), (7) and(10)] or [5, Examples 5.3, 5.4 and 5.7] for more details on decomposition and mappingof Möbius functions.Here most points of interest are already known from the definition of T ωv ω a = ∞ , ω α = 0 , ω β = 1 , (S11a)and the remaining point ω b can be obtained by direct evaluation ω b = ( v b − v α )( v b − v a ) ( v β − v a )( v β − v α ) (S11b) S1.3 From ω to ρ domain The last transformation T ρω is in charge of mapping the upper half-plane of ω intothe parallel-plates domain ρ . See Figure S1. This is achieved in two stages: (i) Theupper half-plane of ω is mapped to the first quadrant of √ ω , such that the negativereal axis of ω is mapped to the positive imaginary axis of √ ω and the real interval ω ∈ [1 , ω b ] is mapped to the real interval √ ω ∈ [1 , √ ω b ] . (ii) The first quadrant of √ ω is mapped to the parallel-plates domain ρ by using the special values of arcsn() in Eqs. (S3) or [2, Table 22.5.1]. S4The scaling and translation of the function arcsn() are chosen such that √ ω α = 0 is mapped to ρ α = 0 and √ ω β = 1 is mapped to ρ β = 1 , which leads to ρ = T ρω ( ω ) = 1 K ( k ρ ) arcsn( √ ω , k ρ ) (S12)The appropriate modulus k ρ is obtained by choosing √ ω b = 1 /k ρ , such that √ ω b bemapped to the upper right corner ρ b of the parallel-plates domain, which leads to thefollowing expressions for the modulus and complementary modulus k ρ = 1 ω b = ( v b − v a )( v b − v α ) ( v β − v α )( v β − v a ) (S13a) k (cid:48) ρ = 1 − k ρ = ( v b − v β )( v b − v α ) ( v α − v a )( v β − v a ) (S13b)Two points of interest, ρ α and ρ β , are already known from the definition of T ρω .The other two points, ρ a and ρ b , are defined when choosing the modulus k ρ and byusing the special values of arcsn() in Equations (S3) or [2, Table 22.5.1] ρ a = ξ a + i ζ a = 1 K ( k ρ ) arcsn( √ ω a , k ρ ) = 0 + i K (cid:48) ( k ρ ) K ( k ρ ) , ρ α = 0 (S14a) ρ b = ξ b + i ζ b = 1 K ( k ρ ) arcsn( √ ω b , k ρ ) = 1 + i K (cid:48) ( k ρ ) K ( k ρ ) , ρ β = 1 (S14b) S1.4 Derivative of the domain transformation
For convenience, we first summarize the domain transformation ρ = T ρr ( r ) = T ρω ◦ T ωv ◦ T vr ( r ) (S15a) ρ = T ρω ( ω ) = 1 K ( k ρ ) arcsn( √ ω , k ρ ) (S15b) ω = T ωv ( v ) = ( v − v α )( v − v a ) ( v β − v a )( v β − v α ) (S15c) v = T vr ( r ) = − cd (cid:18) K ( k r ) 2 r W , k r (cid:19) (S15d)To obtain its derivative, we decompose it by taking the derivatives for each of thetransformation steps ∂ ρ ∂ r = ∂ ρ ∂ ω ∂ ω ∂ v ∂ v ∂ r (S16a) ∂ ρ ∂ ω = 12 K ( k ρ ) 1 ω / (1 − ω ) / (1 − k ρ ω ) / (S16b) ∂ ω ∂ v = ( v α − v a )( v − v a ) ( v β − v a )( v β − v α ) (S16c) ∂ v ∂ r = 2 K ( k r ) W (1 − v ) / (1 − k r v ) / (S16d)where ∂ ρ /∂ ω is due to ∂ arcsn( u , k ) /∂ u = (1 − u ) − / (1 − k u ) − / in [2, Equation(22.15.12)], and ∂ v /∂ r is due to the combination of ∂ cd( u , k ) /∂ u = − k (cid:48) sd( u , k ) nd( u , k ) in [2, Table 22.13.1] and k (cid:48) sd( u , k ) = 1 − cd( u , k ) and k (cid:48) nd( u , k ) = 1 − k cd( u , k ) in [2, Equation (22.6.4)]. The value of each factor in the denominator of ∂ ρ /∂ ω is ω = ( v − v α )( v β − v a )( v − v a )( v β − v α ) (S17a) (1 − ω ) = − ( v − v β )( v α − v a )( v − v a )( v β − v α ) (S17b) (1 − k ρ ω ) = ( v b − v )( v α − v a )( v − v a )( v b − v α ) (S17c) S5since /k ρ = ω b = T ωv ( v b ) .Later, we combine the first two derivatives ∂ ρ /∂ ω and ∂ ω /∂ v ∂ ρ ∂ ω ∂ ω ∂ v = ± i K ( k ρ ) ( v b − v α ) / ( v β − v a ) / ( v − v α ) / ( v − v β ) / v − v a ) / ( v b − v ) / which leads to an expression that consists of two complex branches: ± i .By introducing the third derivative ∂ v /∂ r we arrive to the final expression ∂ ρ ∂ r = ± i W K ( k r ) K ( k ρ ) ( v b − v α ) / ( v β − v a ) / ( v − v α ) / ( v − v β ) / (1 − v ) / (1 − k r v ) / ( v − v a ) / ( v b − v ) / (S18)In case u B > u A , the flux density at each band B must be positive j ( x ) = γ [ u B − u A ] (cid:61) ∂ ρ ∂ r ( x ) ≥ , where the physical constant γ > (S19)In order to achieve this, the + i branch of ∂ ρ /∂ r in Equation (S18) must be chosenas the one carrying physical meaning. This is because /k r > ≥ v b > v > v β > v α > v a when v ∈ B , as it can be seen in Figure S1. S2 Flux per band for interior and exterior cells
The normalized flux per band for interior and exterior cells i int A /Lγ [ u A − u B ] = 2 K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) int , i ext A /Lγ [ u A − u B ] = 2 K (cid:48) ( k ρ ) K ( k ρ ) (cid:12)(cid:12)(cid:12)(cid:12) ext (S20)and the relative error of the flux for an exterior cell with respect to that of an interiorcell i ext A i int A − K (cid:48) ( k ρ ) /K ( k ρ ) | ext K (cid:48) ( k ρ ) /K ( k ρ ) | int − (S21)as functions of their geometrical parameters ( w A /W AB , w B /W AB , H/W AB and W/W AB )are difficult to visualize graphically, since they corresponds to three- and four-dimensionalscalar fields for the cases of interior and exterior cells respectively.Therefore, Figures S2 to S7 (as well as Figures 5 and 7 at the main text ) show two-dimensional slices and contour lines of K (cid:48) ( k ρ ) /K ( k ρ ) and i ext A /i int A − , as functionsof the relative dimensions of the cell. These slices were chosen such that they provideenough information to understand how K (cid:48) ( k ρ ) /K ( k ρ ) and i ext A /i int A − behave fordifferent domain dimensions. S3 Approximations for the interior cell
Calculating the flux per band at an interior cell requires the evaluation of K (cid:48) ( k ) /K ( k ) and its modulus; Equations (S7), (S9) and S13; which can be done by using specializedsoftware like Python , SageMath , Matlab and
Mathematica among others. However,an expression for the flux involving simpler elementary functions is highly desirable,since this would enable the use of any general purpose software or calculator.One convenient way to find such approximations is by using the nome function q = Q ( k ) , shown in Eq. (S8), as a way to compute the ratio K (cid:48) ( k ) /K ( k ) K (cid:48) ( k ) K ( k ) = ln Q ( k ) − π = − π ln Q ( k (cid:48) ) this ratio is closely related to the lattice parameter τ = i K (cid:48) ( k ) /K ( k ) and the nome q = exp( i πτ ) [2, §20.1, §22.1, and Equations (22.2.1) and (22.2.12)]. S6 w A /W AB = w B /W AB H / W A B . . . . . K ( k ρ ) /K ( k ρ ) (a) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 0.3 (b) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 1.0 (c) Min: . . Max: . . Figure S2: Normalized flux per band K (cid:48) ( k ρ ) /K ( k ρ ) = ( i A /L ) / ( γ [ u A − u B ]) as afunction of the cell dimensions. Case of exterior cell with W/W AB = 2 . w A /W AB = w B /W AB H / W A B . . . . . K ( k ρ ) /K ( k ρ ) (a) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 0.3 (b) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 1.0 (c) Min: . . Max: . . Figure S3: Normalized flux per band K (cid:48) ( k ρ ) /K ( k ρ ) = ( i A /L ) / ( γ [ u A − u B ]) as afunction of the cell dimensions. Case of exterior cell with W/W AB = 3 . w A /W AB = w B /W AB H / W A B . . . . . K ( k ρ ) /K ( k ρ ) (a) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 0.3 (b) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 1.0 (c) Min: . . Max: . . Figure S4: Normalized flux per band K (cid:48) ( k ρ ) /K ( k ρ ) = ( i A /L ) / ( γ [ u A − u B ]) as afunction of the cell dimensions. Case of exterior cell with W/W AB = 5 . S7 w A /W AB = w B /W AB H / W A B . . . . . . Relative error (a) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 0.3 (b) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 1.0 (c) Min: . . Max: . . Figure S5: Relative error of flux for an exterior cell with respect to that of an interiorcell i ext A /i int A − as a function of the cell dimensions. Case of exterior cell with W/W AB = 2 . w A /W AB = w B /W AB H / W A B . . . . . . Relative error (a) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 0.3 (b) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 1.0 (c) Min: . . Max: . . Figure S6: Relative error of flux for an exterior cell with respect to that of an interiorcell i ext A /i int A − as a function of the cell dimensions. Case of exterior cell with W/W AB = 3 . w A /W AB = w B /W AB H / W A B . . . . . . Relative error (a) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 0.3 (b) Min: . . Max: . . w A /W AB w B / W A B . . . . H/W AB = 1.0 (c) Min: . . Max: . . Figure S7: Relative error of flux for an exterior cell with respect to that of an interiorcell i ext A /i int A − as a function of the cell dimensions. Case of exterior cell with W/W AB = 5 . S8This is because the Taylor expansion of the nome function Q ( k ) is known and con-verges relatively fast [6, below Equation (12)] [2, Equation (19.5.5)] Q ( k ) = k
16 + 8 (cid:32) k (cid:33) + 84 (cid:32) k (cid:33) + 992 (cid:32) k (cid:33) + O ( k ) (S22)Therefore, for a sufficiently small modulus k or k (cid:48) , the ratio K (cid:48) ( k ) /K ( k ) can beapproximated with enough accuracy by using only the first term of the previousseries K (cid:48) ( k ) K ( k ) (cid:12)(cid:12)(cid:12)(cid:12) k ≈ = − π ln Q ( k ) (cid:12)(cid:12)(cid:12)(cid:12) k ≈ ≈ − π ln (cid:32) k (cid:33) (S23a) K (cid:48) ( k ) K ( k ) (cid:12)(cid:12)(cid:12)(cid:12) k (cid:48) ≈ = − π [ln Q ( k (cid:48) )] − (cid:12)(cid:12)(cid:12)(cid:12) k (cid:48) ≈ ≈ − π ln (cid:32) k (cid:48) (cid:33) − (S23b)as shown before in the literature [7, above Equation (32)] [8, Equation (5)].To aid in the process of approximation, it is convenient to rewrite the modulus k ρ and its complement k (cid:48) ρ . Since for an interior cell we have w a = w A , w b = w B , and r b = W AB = r l = W the points of interest of domain v in Equations (S9) reduce to v a = − , v α = − cd (cid:18) K ( k r ) w A W AB , k r (cid:19) (S24a) v b = 1 , v β = − cd (cid:18) K ( k r ) − w B W AB , k r (cid:19) = cd (cid:18) K ( k r ) w B W AB , k r (cid:19) (S24b)due to half period translation cd( u + 2 K ( k ) , k ) = − cd( u , k ) [2, Table 22.4.3] anddue to double rotation of the argument cd( − u , k ) = cd( i u , k ) = cd( u , k ) [2, Table22.6.1]. This also reduces q r and q (cid:48) r in Equations (S7) to q r = Q ( k r ) = exp (cid:18) − π HW AB (cid:19) , q (cid:48) r = Q ( k (cid:48) r ) = exp (cid:18) − π W AB H (cid:19) (S25)and the moduli k ρ and k (cid:48) ρ in Equations (S13) to k ρ = 2 (cid:34) cd (cid:18) K ( k r ) w A W AB , k r (cid:19) + cd (cid:18) K ( k r ) w B W AB , k r (cid:19)(cid:35)(cid:34) (cid:18) K ( k r ) w A W AB , k r (cid:19)(cid:35) (cid:34) (cid:18) K ( k r ) w B W AB , k r (cid:19)(cid:35) (S26a) k (cid:48) ρ = (cid:34) − cd (cid:18) K ( k r ) w A W AB , k r (cid:19)(cid:35)(cid:34) (cid:18) K ( k r ) w A W AB , k r (cid:19)(cid:35) (cid:34) − cd (cid:18) K ( k r ) w B W AB , k r (cid:19)(cid:35)(cid:34) (cid:18) K ( k r ) w B W AB , k r (cid:19)(cid:35) (S26b)This further leads to an alternative representation for k ρ , which depends on the gap g = W AB − w A = W AB − w B between bands of equal width w A = w B k ρ = 4 sn (cid:18) K ( k r ) gW AB , k r (cid:19)(cid:34) (cid:18) K ( k r ) gW AB , k r (cid:19)(cid:35) (S27a) S9due to quarter period translation cd( u + K ( k ) , k ) = − sn( u , k ) [2, Table 22.4.3] anddoble rotation of the argument sn( − u , k ) = sn( i u , k ) = − sn( u , k ) [2, Table 22.6.1].Also k (cid:48) ρ has an alternative representation, which can be expressed using different bandwidths k (cid:48) ρ = k (cid:48) r sd (cid:18) K ( k r )2 w A W AB , k r (cid:19) cn (cid:18) K ( k r )2 w A W AB , k r (cid:19) sd (cid:18) K ( k r )2 w B W AB , k r (cid:19) cn (cid:18) K ( k r )2 w B W AB , k r (cid:19) (S27b)due to [1 − cd(2 u , k )] / [1 + cd(2 u , k )] = k (cid:48) sd( u , k ) / cn( u , k ) [9, Equations (1.10)and (4.1)]. Here the special functions sn , sd , cn and cd correspond to Jacobian ellipticfunctions analogous to their circular counterparts sin and cos [2, §22.2].Finally, with this alternative representations for k ρ and k (cid:48) ρ , we can obtain approx-imations based on elementary functions for the cases of tall and shallow domains. S3.1 Approximations for tall domains
Consider the following approximations for sufficiently small nome q = Q ( k )lim q → + sn( u , k ) = lim q → + sd( u , k ) = sin (cid:18) π u K ( k ) (cid:19) (S28a) lim q → + cn( u , k ) = lim q → + cd( u , k ) = cos (cid:18) π u K ( k ) (cid:19) (S28b) lim q → + dn( u , k ) = lim q → + nd( u , k ) = 1 (S28c)[2, Equation (19.6.1) and Table 22.5.3] or [10, Equations (10)].In case of tall domains ( H is large) q r → + and k r → + , therefore we obtain k ρ (cid:12)(cid:12)(cid:12)(cid:12) H → + ∞ g ≈ = 4 sin (cid:18) π gW AB (cid:19)(cid:34) (cid:18) π gW AB (cid:19)(cid:35) ≈ (cid:18) π gW AB (cid:19) (S29a) k (cid:48) ρ (cid:12)(cid:12)(cid:12)(cid:12) H → + ∞ w A ≈ w B ≈ = tan (cid:18) π w A W AB (cid:19) tan (cid:18) π w B W AB (cid:19) ≈ (cid:18) π w A W AB (cid:19) (cid:18) π w B W AB (cid:19) (S29b) S3.2 Approximations for shallow domains
Consider the following approximations when the associated nome q (cid:48) = Q ( k (cid:48) ) is suffi-ciently small [10, Equations (11)] sn( u , k ) (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) ≈ ≈ ( k ) − / tanh (cid:18) π u K (cid:48) ( k ) (cid:19) (S30a) cn( u , k ) (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) ≈ ≈ (cid:32) k (cid:48) k q (cid:48) (cid:33) / sech (cid:18) π u K (cid:48) ( k ) (cid:19) (S30b) dn( u , k ) (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) ≈ ≈ (cid:32) k (cid:48) q (cid:48) (cid:33) / sech (cid:18) π u K (cid:48) ( k ) (cid:19) (S30c)also of the subsidiary function sd = sn / dnsd( u , k ) (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) ≈ ≈ (cid:18) q (cid:48) k k (cid:48) (cid:19) / sinh (cid:18) π u K (cid:48) ( k ) (cid:19) (S30d) S10and of the modulus ( k ) / [10, begining of §7] ( k ) / (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) ≈ ≈ − q (cid:48) q (cid:48) = tanh (cid:18) π K ( k ) K (cid:48) ( k ) −
12 ln 2 (cid:19) (S30e)where q (cid:48) = Q ( k (cid:48) ) = exp( − πK ( k ) /K ( k (cid:48) )) .In case of shallow domains ( H is small) q (cid:48) r → + and k (cid:48) r → + , which leads to k ρ (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) r ≈ g ≈ ≈ (cid:18) K ( k r ) gW AB , k r (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) r ≈ ≈ (cid:18) π K ( k r ) K (cid:48) ( k r ) gW AB (cid:19) tanh (cid:18) π K ( k r ) K (cid:48) ( k r ) − ln √ (cid:19) (S31a) k (cid:48) ρ (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) r ≈ ≈ q (cid:48) r sinh (cid:18) π K ( k r ) K (cid:48) ( k r ) w A W AB (cid:19) sinh (cid:18) π K ( k r ) K (cid:48) ( k r ) w B W AB (cid:19) (S31b)where the last equality was obtained through the identity u ) cosh( u ) = sinh(2 u ) ,as shown below for E ∈ { A, B } k (cid:48) r sd (cid:18) K ( k r )2 w E W AB , k r (cid:19) cn (cid:18) K ( k r )2 w E W AB , k r (cid:19) ≈ q (cid:48) r sinh (cid:18) π K ( k r ) K (cid:48) ( k r ) w E W AB (cid:19) sech (cid:18) π K ( k r ) K (cid:48) ( k r ) w E W AB (cid:19) = 4 q (cid:48) r sinh (cid:18) π K ( k r ) K (cid:48) ( k r ) w E W AB (cid:19) (S32)Here K (cid:48) ( k r ) /K ( k r ) = 2 H/W AB was obtained from q (cid:48) r = Q ( k (cid:48) r ) = exp( − πK ( k r ) /K ( k (cid:48) r )) =exp( − πW AB / H ) , due to Equations (S8) and (S25). S4 Comparison with previous results
In this section, simulations and theoretical expressions from the literature are writtenusing the notation of this work, in order to compare against the normalized flux perband in Equations (12) of the main text i A /Lγ [ u A − u B ] = 2 K (cid:48) ( k ρ ) K ( k ρ ) (S33)Note that all these results from the literature consider arrays with bands of equalwidth w E := w A = w B and an interior cell .Table S1 shows simulations obtained in [11, Fig. 7a] for the normalized faradaiccurrent in an interior cell . The expression [11, | G ss | corresponds to | G g | = | G c | in§2.2 and Equation (3) at steady state], which was rewritten using the notation of thistext | G ss | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N A i A LN A n e DF (cid:124) (cid:123)(cid:122) (cid:125) γ [ c b (cid:124)(cid:123)(cid:122)(cid:125) u A − (cid:124)(cid:123)(cid:122)(cid:125) u B ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 2 K (cid:48) ( k ρ ) K ( k ρ ) (S34)where n e is the number of electrons exchanged in the Faradaic process, D is the diffu-sion coefficient of the redox species, F is the Faraday’s constant and c b corresponds tothe bulk concentration of species in the electrochemical cell. These simulated resultscorrespond to the case of a process of limiting current (with counter electrode externalto the IDA), where the concentration at the bands of one array equals c b , while theconcentration at the bands of the other array equals zero. For plotting this data inFigure 6 at the main text , the horizontal variable [11, w e /h c in Figure 7a] was alsorewritten according to the notation of this text, that is w E /H in Table S1, and laterrelated to H/W by HW = (cid:34) w E H · (cid:18) gw E + w E w E (cid:19)(cid:35) − (S35) S11 w E /g = 0 . w E /g = 1 w E /H | G ss | w E /H | G ss | | G ss | = 2 K (cid:48) ( k ρ ) /K ( k ρ ) for an interior cell with bands of equal width w E = w A = w B [11, Figure 7a]. Data obtained by simu-lating for different values of w E /H and w E /g ∈ { . , } .where g is the gap between electrodes such that W AB = w E + g . Only the data forthe case w E = g ( w E /W AB = 0 . ) was plotted in Figure 6 at the main text .Table S2 shows exhaustive simulations obtained in [12, Fig. 7a] for the normalizedfaradaic current in an interior cell . The expression for the average flux density [12, ¯ φ lim λ ] was related to the current per array N A i A = N A Lw E F n e ¯ φ lim λ [12, Equation(12)], and then replaced in the expression of [12, Fig. 7a] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ¯ φ lim λ w E / π D ¯ c (cid:96), (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 1 π (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i A /LF n e D (cid:124) (cid:123)(cid:122) (cid:125) γ [2¯ c (cid:96), (cid:124)(cid:123)(cid:122)(cid:125) u A − (cid:124)(cid:123)(cid:122)(cid:125) u B ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 1 π · K (cid:48) ( k ρ ) K ( k ρ ) (S36)showing that the normalized currents in [12] and in this work are related by a factor of π ≈ . Here ¯ c (cid:96), is the smallest initial concentration between both redox species in w E /W AB H/W AB . /π . /π . /π . /π /π /π /π /π /π (2 /π ) K (cid:48) ( k ρ ) /K ( k ρ ) for bands of equal width w E = w A = w B [12, Figure 7a]. Data obtained by simulating exhaustively for differentvalues of H/W AB and w E /W AB . S12the electrochemical cell. These simulated results correspond to the case of a processof limiting current (with one of the arrays acting as counter electrode), where theconcentration at the bands of one array equals c (cid:96), , while the concentration at thebands of the other array equals zero [12, Figure 2]. See Figure 6 at the main text fora graphical representation of Table S2.The expression [13, C I in Equation (13)] shows the capacitance of what corre-sponds to half interior cell in this work. In [13] it is assumed that the potentialsapplied at bands of each array are symmetric ( u A = V and u B = − V ). Despite this,the difference of potentials inside half interior cell is only V , because this containsonly half band of one array. Also, since the charge (flux of electric displacement)at one band corresponds to i A , then in half interior cell we have a charge of i A / .Therefore, the normalized capacitance in half interior cell can be rewritten using thenotation of this text as C I (cid:15)L = 1 (cid:15)L i A / V = i A /L(cid:15) (cid:124)(cid:123)(cid:122)(cid:125) γ [ V (cid:124)(cid:123)(cid:122)(cid:125) u A − − V (cid:124)(cid:123)(cid:122)(cid:125) u B ] = 2 K (cid:48) ( k ρ ) K ( k ρ ) (S37)where (cid:15) = (cid:15) (cid:15) r is the absolute permittivity of the dielectric. See Figure 6 at the maintext for a comparison between [13, C I / ( (cid:15)L ) in Equation (13)] with K (cid:48) ( k ρ ) /K ( k ρ ) . S5 Regions of validity for approximated flux
The relative error between the approximated and exact fluxes ( N A ˜ i A and N A i A ) foran ideal domain corresponds to that between their normalized counterpartsrelative error = N A ˜ i A N A i A − (cid:34) K (cid:48) ( k ρ ) K ( k ρ ) (cid:35) app. (cid:34) K (cid:48) ( k ρ ) K ( k ρ ) (cid:35) − − (S38)This was used to find the combinations of band widths and domain height that pro-duce a relative error less than ± .
05 = ± , thus defining regions of parameterswhere the approximations are valid. The search for these regions of validity was donegraphically, and therefore, it was more convenient to consider bands of equal width w E := w A = w B , such that the resulting plots were two-dimensional.When expressing the ratio K (cid:48) ( k ) /K ( k ) in terms of the nome function Q ( k ) , orequivalently Q ( k (cid:48) ) due to Equation (S8), one can realize that its approximations aredue to a two-step approximation process (cid:34) K (cid:48) ( k ρ ) K ( k ρ ) (cid:35) app. := ln ˜ Q (˜ k ρ ) − π , for k ρ ≈ − π ln ˜ Q (˜ k (cid:48) ρ ) , for k (cid:48) ρ ≈ (S39)First, approximation of the nome function due to Equation (S22) by ˜ Q ( k ) := k / ,and later the approximations ˜ k ρ and ˜ k (cid:48) ρ of the moduli, which are given in Equations(S29) and (S31).Therefore, it is possible to analyze the propagation of the relative error, for k ρ ≈ ,by looking at the following total ratio N A ˜ i A N A i A = (cid:34) K (cid:48) ( k ρ ) K ( k ρ ) (cid:35) app. (cid:34) K (cid:48) ( k ρ ) K ( k ρ ) (cid:35) − = ln ˜ Q (˜ k ρ )ln Q ( k ρ ) = ln ˜ Q (˜ k ρ )ln ˜ Q ( k ρ ) ln ˜ Q ( k ρ )ln Q ( k ρ ) (S40)where a relative error of ± .
05 = ± corresponds to the total ratios . and . .The ratios ln ˜ Q ( k ρ ) / ln Q ( k ρ ) and ln ˜ Q (˜ k ρ ) / ln ˜ Q ( k ρ ) show the individual contributionsof the approximated nome ˜ Q and the approximated modulus ˜ k ρ to the total ratio ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) , and therefore, to the relative error. S13 w A /W AB = w B /W AB H / W A B . . . . . ln ˜ Q ( k ρ ) / ln Q ( k ρ ) 0.0 0.2 0.4 0.6 0.8 1.0 w A /W AB = w B /W AB H / W A B . . . . . . k ρ w A /W AB = w B /W AB H / W A B .
50 0 . . . . ln ˜ Q (˜ k ρ ) / ln ˜ Q ( k ρ ) 0.0 0.2 0.4 0.6 0.8 1.0 w A /W AB = w B /W AB H / W A B . . . . . . ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) Figure S8: Propagation of the ratio of approximation ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) for the caseof tall domains and wide bands. Contour lines of selected values are shown in black.Line w E /W AB = 0 . is shown in red.This is shown in Figure S8 for the case of tall domains and wide bands. The firstratio ln ˜ Q ( k ρ ) / ln Q ( k ρ ) < . (relative error of ) for band widths w E /W AB > . ( k ρ (cid:46) . ), which corresponds to very extreme cases of wide bands. This agreeswith the relative error of ≈ in [7, before Equation (32)] when approximatingonly the nome function ˜ Q ( k ρ ) for gap widths g/W AB < . ⇔ w E /W AB > . ( k ρ < . ). Fortunately, the contribution due to the approximation of the modulus ˜ k ρ in Equation (S29a) makes the second ratio ln ˜ Q (˜ k ρ ) / ln ˜ Q ( k ρ ) < in almost thewhole domain ( w E /W, H/W AB ) . This helps to extend the region where the total ratio ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) < . (relative error less than ), which corresponds to morereasonable band widths w E /W AB > . .A similar analysis of the propagation of the relative error can be done for k (cid:48) ρ ≈ when considering the following total ratio N A ˜ i A N A i A = (cid:34) K (cid:48) ( k ρ ) K ( k ρ ) (cid:35) app. (cid:34) K (cid:48) ( k ρ ) K ( k ρ ) (cid:35) − = ln Q ( k (cid:48) ρ )ln ˜ Q (˜ k (cid:48) ρ ) = ln ˜ Q ( k (cid:48) ρ )ln ˜ Q (˜ k (cid:48) ρ ) ln Q ( k (cid:48) ρ )ln ˜ Q ( k (cid:48) ρ ) (S41)where the ratios ln Q ( k (cid:48) ρ ) / ln ˜ Q ( k (cid:48) ρ ) and ln ˜ Q ( k (cid:48) ρ ) / ln ˜ Q (˜ k (cid:48) ρ ) show the individual contri-butions of the approximated nome ˜ Q and the approximated complementary modulus ˜ k (cid:48) ρ to the total ratio ln Q ( k (cid:48) ρ ) / ln ˜ Q (˜ k (cid:48) ρ ) , and therefore, to the relative error.This is shown in Figure S9 for the case of tall domains and narrow bands. Un-like the previous case, the first ratio ln Q ( k (cid:48) ρ ) / ln ˜ Q ( k (cid:48) ρ ) > . (relative error of − ) for widths w E /W AB < . ( k (cid:48) ρ (cid:46) . ), which corresponds to a largeportion of the parameter domain. This agrees with the statement in [8, beforeEquation (5)], which says that the approximation converges rapidly for gap widths S14 w A /W AB = w B /W AB H / W A B . . . . . ln Q ( k ρ ) / ln ˜ Q ( k ρ ) 0.0 0.2 0.4 0.6 0.8 1.0 w A /W AB = w B /W AB H / W A B . . . . . k ρ w A /W AB = w B /W AB H / W A B . . . . . ln ˜ Q ( k ρ ) / ln ˜ Q (˜ k ρ ) 0.0 0.2 0.4 0.6 0.8 1.0 w A /W AB = w B /W AB H / W A B . . . . . ln Q ( k ρ ) / ln ˜ Q (˜ k ρ ) Figure S9: Propapation of the ratio of approximation ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) for the caseof tall domains and narrow electrodes. Contour lines of selected values are shown inblack. Line w E /W AB = 0 . is shown in red. g > . w E ⇔ w E /W AB < . (no quantitative error is specified). Unfortunately,the contribution due to the approximation of the complementary modulus ˜ k (cid:48) ρ makesthe second ratio ln ˜ Q ( k (cid:48) ρ ) / ln ˜ Q (˜ k (cid:48) ρ ) < in most of the domain ( w E /W AB , H/W AB ) .This approximation helps to simplify the expression for the complementary modu-lus in Equation (S29b) at the expense of shrinking the domain where the total ratio ln Q ( k (cid:48) ρ ) / ln ˜ Q (˜ k (cid:48) ρ ) > . (relative error of approximation less than − ), whichcorresponds to electrode widths w E /W < . .The propagation of the relative error for cases of shallow cells have not beenreported previously in the literature and they are given below as a reference.Figure S10 shows the propagation of the total ratio, in Equation (S40), for thecase of shallow domains and wide bands. The first ratio ln ˜ Q ( k ρ ) / ln Q ( k ρ ) < . (relative error of ) for combinations of ( w E /W AB , H/W AB ) roughly to the rightof the line joining the points (1 , → (0 . , ( k ρ (cid:46) . ), which again correspondto very wide electrodes. The contribution due to the approximation of the modulus ˜ k ρ in Equation (S31a) makes the second ratio ln ˜ Q (˜ k ρ ) / ln ˜ Q ( k ρ ) < in most ofthe domain ( w E /W AB , H/W AB ) . This also extends the region where the total ratio . < ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) < . (relative error less than ± ), which correspondsto combinations of ( w E /W AB , H/W AB ) approximately at the right of the line joiningthe points (1 , → (0 . , .Figure S11 shows the propagation of the total ratio, in Equation (S41), for thecase of shallow domains and narrow bands. The first ratio ln Q ( k (cid:48) ρ ) / ln ˜ Q ( k (cid:48) ρ ) > . (relative error of − ) for combinations of ( w E /W AB , H/W AB ) roughly to the leftof the line joining the points (1 , → (0 . , ( k (cid:48) ρ (cid:46) . ), which correspondsto a large portion of the parameter domain. The approximation of the comple- S15 w A /W AB = w B /W AB H / W A B . . . . . ln ˜ Q ( k ρ ) / ln Q ( k ρ ) 0.0 0.2 0.4 0.6 0.8 1.0 w A /W AB = w B /W AB H / W A B . . . . . . k ρ w A /W AB = w B /W AB H / W A B . . . . . ln ˜ Q (˜ k ρ ) / ln ˜ Q ( k ρ ) 0.0 0.2 0.4 0.6 0.8 1.0 w A /W AB = w B /W AB H / W A B . . . . . . . . . ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) Figure S10: Propapation of the ratio of approximation ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) for the caseof shallow domains and wide bands. Contour lines of selected values are shown inblack. Line joining the points (1 , → (0 . , is shown in red.mentary modulus ˜ k (cid:48) ρ makes the second ratio ln ˜ Q ( k (cid:48) ρ ) / ln ˜ Q (˜ k (cid:48) ρ ) < in most of thedomain ( w E /W AB , H/W AB ) . This approximation helps to simplify the expressionfor the complementary modulus in Equation (S31b) at the expense of shrinking thedomain where the total ratio ln Q ( k (cid:48) ρ ) / ln ˜ Q (˜ k (cid:48) ρ ) > . (relative error of approxi-mation less than − ), which corresponds to combinations of ( w E /W AB , H/W AB ) approximately at the left of the line joining the points (1 , → (0 . , .In summary, the regions of validity for the approximated flux (within ± relativeerror) correspond to w E /W AB > . (wide bands) and w E /W AB < . (narrowbands) for the case of tall domains H/W AB > . In the case of shallow domains H/W AB ≤ , these regions of validity approximately correspond to the right of theline connecting the points (1 , → (0 . , (wide bands) and to the left of the lineconnecting connecting the points (1 , → (0 . , (narrow bands), where the linethat passes through the points (1 , and ( p, is given by w E W AB + (1 − p ) HW AB = 1 (S42)that is, w E + 0 . H (cid:38) W AB for wide bands and w E + 0 . H (cid:46) W AB for narrowbands. S6 References [1] T. A. Driscoll and L. N. Trefethen,
Schwarz-Christoffel mapping , vol. 8.Cambridge University Press, 2002.16 w A /W AB = w B /W AB H / W A B . . . . . ln Q ( k ρ ) / ln ˜ Q ( k ρ ) 0.0 0.2 0.4 0.6 0.8 1.0 w A /W AB = w B /W AB H / W A B . . . . . k ρ w A /W AB = w B /W AB H / W A B . . . . ln ˜ Q ( k ρ ) / ln ˜ Q (˜ k ρ ) 0.0 0.2 0.4 0.6 0.8 1.0 w A /W AB = w B /W AB H / W A B . . . . . ln Q ( k ρ ) / ln ˜ Q (˜ k ρ ) Figure S11: Propapation of the ratio of approximation ln ˜ Q (˜ k ρ ) / ln Q ( k ρ ) for the caseof shallow domains and narrow bands. Contour lines of selected values are shown inblack. Line joining the points (1 , → (0 . , is shown in red.[2] F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F.Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A.McClain, eds., NIST Digital Library of Mathematical Functions . June, 2020. http://dlmf.nist.gov/ . Release 1.0.27 of 2020-06-15.[3] Z. Nehari,
Conformal Mapping . International series in pure and appliedmathematics. McGraw-Hill, 1952.
GoogleBooks:9W24AAAAIAAJ .[4] M. J. Ablowitz and A. S. Fokas,
Complex variables: Introduction andapplications . Cambridge University Press, 2nd ed., Apr., 2003.[5] P. J. Olver, “Complex analysis and conformal mapping,” June, 2020. .[6] A. Kneser, “Neue untersuchung einer reihe aus der theorie der elliptischenfunktionen.”
Journal für die reine und angewandte Mathematik (1927)209–218. http://gdz.sub.uni-goettingen.de/dms/load/toc/?PID=PPN243919689_0158 .[7] K. Aoki, M. Morita, O. Niwa, and H. Tabei, “Quantitative analysis of reversiblediffusion-controlled currents of redox soluble species at interdigitatedgitatedarray electrodes under steady-state conditions,”
Journal of ElectroanalyticalChemistry and Interfacial Electrochemistry no. 2, (Dec., 1988) 269–282. http://asura.apphy.u-fukui.ac.jp/~aoki/publication/PDF/JEC256_269.pdf .17[8] W. E. Morf, M. Koudelka-Hep, and N. F. de Rooij, “Theoretical treatment andcomputer simulation of microelectrode arrays,”
Journal of ElectroanalyticalChemistry no. 1, (May, 2006) 47–56, rérodoc:16975 .[9] B. C. Carlson, “Symmetry in c, d, n of jacobian elliptic functions,”
Journal ofMathematical Analysis and Applications no. 1, (Nov., 2004) 242 – 253,
CORE:82142693 .[10] J. D. Fenton and R. S. Gardiner-Garden, “Rapidly-convergent methods forevaluating elliptic integrals and theta and elliptic functions,”
The ANZIAMJournal (July, 1982) 47–58. http://johndfenton.com/Papers/Fenton82b-Elliptic+Theta-functions.pdf .[11] J. Strutwolf and D. E. Williams, “Electrochemical sensor design using coplanarand elevated interdigitated array electrodes. a computational study,” Electroanalysis no. 2, (Feb., 2005) 169–177.[12] C. Guajardo, S. Ngamchana, and W. Surareungchai, “Mathematical modelingof interdigitated electrode arrays in finite electrochemical cells,” Journal ofElectroanalytical Chemistry (Sept., 2013) 19–29, arXiv:1602.06428[physics.chem-ph] .[13] R. Igreja and C. J. Dias, “Analytical evaluation of the interdigital electrodescapacitance for a multi-layered structure,”
Sensors and Actuators A: Physical112