The maximum voltage drop in an on-chip power distribution network: analysis of square, triangular and hexagonal power pad arrangements

A mathematical model of the voltage drop which arises in on-chip power distribution networks is used to compare the maximum voltage drop in the case of different geometric arrangements of the pads supplying power to the chip. These include the square or Manhattan power pad arrangement which currently predominates, as well as equilateral triangular and hexagonal arrangements. In agreement with findings in the literature and with physical and SPICE models, the equilateral power pad arrangement, independent of the underlying power mesh configuration, is found to minimize the maximum voltage drop. This headline finding is a consequence of relatively simple formulas for the voltage drop, with explicit error bounds, which are established using complex analysis techniques, and elliptic functions in particular.


Introduction
Control of the maximum voltage drop between power distribution pads is a factor of increasing importance in the design of the power distribution network (PDN) of modern IC computer chips. The voltage drop between power pads depends both on the current flowing in the power mesh between the pads and on the electrical resistance in the power mesh. The physical layout of a computer chip and the interaction between the chip and its power distribution network are described in detail by Shakeri and Meindl [7] in the context of both wire-bond and flip-chip PDN design. They focus on the dominant paradigm in which the power pads and the power mesh are arranged in a square grid, which is known as the Manhattan architecture, they derive the equations governing the voltage drop and provide the leading terms of the solution. The Y-architecture, in which pads are arranged in an equilateral lattice and the power mesh is also arranged in an equilateral grid, is considered by Chen et al. [3]. Analytical and simulation results are obtained which indicate a 5% reduction in the maximum voltage drop in the case of a single layer Y-architecture compared to the single layer Manhattan architecture.
Aquareles et al. [1] put the mathematical aspects of the work of Shakeri and Meindl on a firm footing. They obtain an asymptotic formula for the maximum voltage drop in terms of the size of the pads, including higher order terms that would seem to be beyond the techniques in [7]. The main mathematical tool they use is that of matched asymptotic expansions. In the present work, we use a complex analysis method to derive an expression for the maximum voltage drop in the case of the square pad arrangement. This method is simpler and more direct than the approach in [1] and covers, without additional effort, the case of pads arranged in an equilateral triangular array. With a little extra work, the method extends to treat the case of pads arranged in a hexagonal pattern.
The results we obtain suggest that the smaller maximum voltage drops observed by Chen et al. in [3] are due to the arrangement of the pads in an equilateral array and are Date: May 1, 2014. 1 independent of the configuration of the underlying power mesh. That is, for an equilateral disposition of the power pads superposed over a fine Manhattan power mesh there will be a similar voltage drop to that observed in the Y-architecture.
We also obtain formulas for the maximum voltage drop in each of these configurations (square, triangular and hexagonal). It is found than the hexagonal pad arrangement has the largest voltage drop of the three configurations considered. Even so, it may also be useful to have an explicit formula for the voltage drop in this case since, however important, control of the maximum voltage drop is but one of several constraints in the design of an on-chip PDN. Finally, the availability of explicit formulas makes it possible to accurately predict the maximum voltage drop at an early point in the circuit design stage, thereby obviating the need for costly re-design.

Mathematical model of the voltage drop
In this section we describe the mathematical model of the power distribution network and the associated voltage drop as derived by Shakeri and Meindl [7].
The surface of the integrated circuit is modeled as an infinite complex plane in which the power pads of the power distribution network are modeled as circular disks. Power to the chip is supplied through these power pads and distributed through a fine grid of wires called the power mesh. The square and triangular arrangement of the pads are displayed below. The planar region consisting of the complex plane with these circular disks removed is denoted by Ω. Under the assumption of uniform current flow between pads, the voltage drop satisfies the equation ∆u = c as the power mesh (triangular or square) gets finer.
The constant c on the right hand side of this partial differential equation codes for the resistance properties of the wires of the mesh and the current drawn from the power network. In order to make a fair comparison between the voltage drop across different power distribution network configurations, the resistance properties of the underlying integrated circuits (IC) and the current drawn need to be the same, that is we need to use the same constant c in all cases. Moreover, since we measure the relative change in the maximum voltage drop across different arrangements of the power pads, and since the solution to ∆u = c is proportional to c, it suffices to take the common value c = 1 in the modeling equation. Next, the power distribution pads are held at constant voltage, which we may take to equal zero. Thus the governing partial differential equation for the voltage in the region Ω between the power pads is ∆u = 1 in Ω, u = 0 on ∂Ω.
The voltage between the pads will then be negative since u is subharmonic and the pads themselves are held at voltage 0, while the voltage drop relative to the pads will simply be −u. It is interesting to note that the solution of the partial differential equation ∆u = −2 in a domain D, also with zero Dirichlet boundary conditions, describes the expected exit time of Brownian motion from the domain. Thus, the problem of determining the maximum voltage drop is mathematically equivalent to determining the maximum expected lifetime of Brownian motion in the domain complementary to the power pads.
The partial differential equation (2.1) obeys a scaling law: If u(z) is the solution of ∆u = 1 in a domain D then v(w) = r 2 u(w/r) is the solution of ∆v = 1 in the domain rD. Thus, if the radius of the power pads and the spacing between their centres both change by a factor of r then the maximum voltage drop changes by a factor r 2 . If we know the voltage drop for all values of the radius of the power pads for some fixed spacing between their centres then we can scale this result to determine the voltage drop in the case of any power pad radius and any spacing between their centres. Next, in order to make a fair comparison between different geometric power pad configurations, the proportion of the area on the chip occupied by the power pads (let's call it p) should be the same in each case. Notice that p, the area of the power pads per unit area on the chip, does not change under the scaling z → rz discussed above, whereas the voltage drop changes by a factor r 2 . Thus, even for prescribed areal density p of the power pads, the voltage drop can be made as small as one wishes by taking smaller pads closer together. Thus, in order to make a fair comparison between different configurations, it is necessary not only to ensure that the aereal density of the power pads are the same in each configuration but also to specify the radius ε of each pad. The values of p and ǫ then determine the spacing between the pads. (Alternatively, one could instead specify the spacing between the pads rather than their radius, but this seems less natural.) Referring to Figure 1, each pad in the square arrangement configuration lies at the centre of a square of side d 1 which does not overlap with the corresponding square for any other pad. Thus the aereal density p = πε 2 /d 2 1 in this case. For an equilateral triangular arrangement, each pad lies at the centre of a diamond of area √ 3 d 2 2 /2 which does not overlap with the corresponding diamond for any other pad. Thus the aereal density of the pads in the equilateral configuration is p = 2πε 2 /( √ 3 d 2 2 ). For prescribed common radius ε of the pads, the aereal density of the pads will be the same in both configurations once Assuming, therefore, that in the square arrangement we have power pads of radius ε whose centres are unit distance apart, in the equilateral arrangement we should have power pads of radius ε whose centres are d 2 = √ 2/ 4 √ 3 ≃ 1.0745 apart. In the case of the hexagonal configuration, as shown in Figure 2, each pad lies at the centre of an equilateral triangle of sidelength √ 3 d 3 which does not overlap with the corresponding triangle for any other pad. The area of this triangle is , so that the aereal density for the hexagonal configuration is p = 4πε 2 /(3 √ 3 d 2 3 ). In order that this agrees with the aereal density p = πε 2 for the previous configurations, we need In this case, each hexagon has area 2.

Main numerical results
Analytic formulas for the voltage drop in the case of each of the arrangements of the pads considered above are established in Sections 4 and 5. These yield the following bounds for the maximum voltage drop. In terms of the radius ε of the pads, the maximum voltage drop V S max in the case of the square arrangement is The maximum voltage drop V T max in the case of the triangular configuration is In the case of the hexagonal configuration, the voltage drop at the centre of a hexagon is It is notable that, apart from the error term, the maximum voltage drop has the same dependence on the pad size in all three cases, the only difference being in the constant term. The conclusion is that the hexagonal pad arrangement has the worst voltage drop among the configurations which we consider, the best being the triangular lattice with the standard square lattice being in an intermediate position.
One intuitive explanation of this situation is that though in the hexagonal arrangement there are six disks around the origin they are, crucially, further separated from the origin than in the other configurations considered. It is possible to fit a bigger disk around the origin which does not meet the boundary of Ω and this allows the Brownian motion to increase its expected lifespan.
The analytical results to come in Sections 4 and 5 yield explicit error bounds in (3.1), (3.2) and (3.3), which are displayed graphically in Figure 3. The curves represent upper and lower bounds for the maximum voltage drops V S max (ε) and V T max (ε) which take account of the error terms. The graph in the hexagonal case shows the voltage drop V H (ε) at Equilateral configuration, upper bound Equilateral configuration, lower bound Square configuration, lower bound Square configuration, upper bound Hexagonal configuration, lower bound the centre of a hexagon. Presumably, this is the maximum voltage drop, that is the maximum voltage drop presumably occurs at the centre of a hexagon, but in any case the maximum voltage drop is at least this large. Thus, even at the limits of the error bounds, the equilateral arrangement outperforms the square and the hexagonal arrangements for all pad sizes. Note also that the error bounds are seen to be quite tight in both the square and the equilateral configurations, so that the formulas (3.1) and (3.2) are accurate. Note that the range of pad size ε (from 0.1 to 0.3) relative to the distance between the distance between the centres of the pads (d 1 , d 2 , d 3 , each of which is about unit size) is informed by industry norms (see [7, Table III]).
In order to test the robustness of these analytical results we assembled two boards, each with a rectangular mesh of resistances. A constant current sink was connected at each node. In one of the boards the voltage distribution was through a collection of pads in a square configuration and in the other the pads were in a triangular configuration. All pads were held at 5V. The maximum voltage drop was measured for each board. It was 1.91V in the triangular pad setting versus 2.03V in the square setting. SPICE simulations with the same configuration gave voltage drops of 1.94V in the triangular case and 2.05V in the square case. The difference between the on board measurements and the SPICE simulations may be due to less than perfect current sinks. In this section, analytic expressions for the voltage drop in both the square and the triangular pad arrangements are obtained. Both configurations correspond to lattices in the plane, permitting direct use of the standard theory of elliptic functions. We next set out those aspects of the theory that we will need, as well as the special results that pertain for the square and equilateral lattices, drawing on the classic text by Hille [6, Section 13.2] as a standard general reference.
4.1. The square and equilateral lattices. A lattice of points in the complex plane consists of all integer linear combinations 2w 1 m + 2w 3 n (m, n ∈ Z) of two given complex numbers 2w 1 and 2w 3 for which w 3 /w 1 has positive imaginary part. We immediately specialize to the case in which so that α is a q th -root of unity. In this case the lattice is described by The resulting lattice is invariant under multiplication by α precisely when there are integers k and j such that It is not difficult to see, for example by examining the resulting equations for the real and imaginary parts separately, that α will satisfy such an identity only in the cases q = 4 and q = 6. The case q = 4, with α = i, 2w 3 = id, α 2 = −1, corresponds to the square lattice. The case q = 6, with α = e πi/3 , 2w 3 = e πi/3 d, α 2 = α − 1, corresponds to the triangular lattice. The values of d in each case are governed by (2.2) which guarantee that the areal densities of the pads agree.
Much of the analysis in the next subsections is essentially unchanged whether we work with the square or with the triangular lattice. We will therefore retain the notation q, d, α with the understanding that in the case of the square lattice, in the case of the triangular lattice, (4.3) the advantage being that in this way we can treat both configurations simultaneously.

4.2.
The Weierstrass σ-function for a plane lattice. The Weierstrass elliptic Pfunction associated with a lattice is doubly periodic with periods 2w 1 and 2w 3 , and is analytic except for double poles at each of the lattice points. The corresponding σfunction is defined by where ′ denotes the product over all lattice points with zero omitted. The Weierstrass ζfunction is defined by where ′ denotes the sum over all lattice points with 0 omitted. 6 A quasi-periodicity property of the function σ plays a key role in our analysis. Set Then, [6, Identity 13.2.19], These two identities lead to the full quasi-periodicity property for any integers m and n. To proceed further, we need to compute η 1 and η 3 explicitly for the square and the triangular lattices, at which point the quasi-periodicity property (4.6) will become explicit in these cases. While these results are known, we give the explicit computations here for completeness.

4.3.
Computation of η 1 and η 3 for the square and triangular lattices. In the case of a general lattice, η 1 and η 3 are related by the identity 2w In the case of either of our lattices, this identity becomes (see (4.1)) The invariance of the lattice under multiplication by α = e 2πi/q and its powers, with q = 4 for the square lattice and q = 6 for the triangular lattice, leads to a second linear relationship between η 1 and η 3 as follows. By definition, Replacing λ by α k λ, k = 1, . . . , q − 1, gives a total of q expressions for η 1 . Adding these leads to we find that This procedure is repeated for η 3 , which is given by Replacing λ by α k λ, k = 1, . . . , q − 1, and adding all q expressions for η 3 , leads to Together, (4.9) and (4.10) yield Solving the simultaneous equations (4.7) and (4.11) gives .
(4.12) Lemma 1. In the case of the square lattice (d = 1, α = i) while in the case of the triangular lattice (d = √ 2/ 4 √ 3, α = e πi/3 ), These results are easily verified, in view of (4.1), by replacing α by i and d by 1 in (4.12) in the case of the square lattice to obtain (4.13). In the case of the triangular lattice, replace α by e πi/3 in (4.12) and use to obtain (4.14).

4.4.
Quasi-periodicity and true periodicity for the square and triangular lattices. These values for η 1 and η 3 lead to a simple form of the general quasi-periodicity relation (4.6) for the σ-function in the case of the square and the triangular lattices. Surprisingly, perhaps, this relation has the same form in both cases, thereby unifying the analysis required to derive an analytic expression for the IR-drop. With an eye to (4.6), recall that a general lattice point is λ m,n = 2w 1 m + 2w 3 n = md + nαd, where m and n are integers. Then, in the case of the square lattice with d = 1 and using the η-values given by (4.13), λ m,n = m + in and 2mη 1 + 2nη 3 = mπ − inπ = π λ m,n . 8 In the case of the triangular lattice with d = √ 2/ 4 √ 3, we use the η-values given by (4.14) to obtain The quasi-periodicity property (4.6) of the σ-function, in the case of either the square or the triangular lattice, therefore becomes σ z + λ m,n = (−1) m+n+mn exp (z + 1 2 λ m,n ) π λ m,n σ(z) = (−1) m+n+mn exp πλ m,n z + π 2 |λ m,n | 2 σ(z). This quasi-periodicity property of the σ-function leads to true periodicity of a related function.

Lemma 2. Set
In the case when either Λ is the square lattice or the triangular lattice, h is periodic in the sense that h(z + λ) = h(z), for z ∈ C \ Λ, λ ∈ Λ. Furthermore, the value of h doesn't change under reflection in any side of the relevant lattice.
Remark 1. The periodiciy of h in the case of square or triangular lattices also follows from the results in [4, Proposition 3.4] which builds upon work in [5]. Gröchenig and Lyubarskii have a more general periodicity result which is valid for all lattices and involves an explicit normalization factor in terms of η 1 and η 3 . The computation of η 1 and η 3 above shows that no normalization factor arises for triangular or square lattices.
Proof. Taking the logarithm of (4.15) with λ = λ m,n ∈ Λ leads to log σ(z + λ) = log |σ(z)| + Re πλz + π 2 |λ| 2 = log |σ(z)| + π 2 |λ| 2 + 2 Re λ z . But, This establishes the periodicity of h. To see that h is invariant under reflection in any side of the lattice, it suffices to show that h α 2k z = h(z) (4.17) where α = i and k = 0 or 1 in the case of the square lattice, while α = e πi/3 and k = 0, 1 or 2 in the case of the triangular lattice. In any of these cases, The invariance of the lattice under multiplication by α 2k , that is α 2k Λ = Λ, and then its invariance under complex conjugation, shows that On taking logarithms, the identity (4.17) follows.

4.5.
Analytic expressions for the IR-drop in the square and the triangular arrangements.
Let Ω ε = C \ λ∈Λ D(λ, ε) denote the region formed by removing from the plane a closed disk of radius ε about each lattice point. Our main result gives an analytic bound for the voltage drop in both the square and triangular arrangement of the pads. It continues to be possible to analyse both configurations simultaneously, which we do. After stating and proving the analytic bound, we derive the explicit numerical bounds (3.1) and (3.2) which prove, in particular, that the equilateral disposition outperforms the square arrangement. Before stating the main analytical result, Theorem 1, we need an estimate on the σ-function near the origin. Proof. Recall the expression (4.4) for the σ-function. By the symmetry of the lattice under multiplication by α k , we see that When these q expressions for σ(z) are multiplied together, one obtains where (4.8) leads to the elimination of the exponential terms, and the identity was used at the last step in (4.21). Taking the logarithm of (4.21) leads to log |σ(z)| = log |z| The power series expansion of the analytic function − log(1 − w) about 0 is Since |λ| ≥ 1 for λ ∈ Λ \ {0}, once |z| q ≤ 3 5 we can apply (4.23) with w = (z/λ) q to obtain 1 q Together with (4.22), this proves (4.18). The estimates (4.19) and (4.20) can be obtained numerically. may be written as

26)
and A q has the value given in the statement of Lemma 3.
Proof. Let h ε be the function which is harmonic on Ω ε and has boundary values By Lemma 2, these boundary values are periodic and therefore so too is h ε (that is, h ε (z + λ) = h ε (z) for z ∈ Ω ε and λ ∈ Λ). Define a function u ε by (4.25). Then ∆u ε = 1 in Ω ε , this because ∆ |z| 2 = 4 while log |σ(z)| is harmonic on Ω ε being the logarithm of the modulus of a non-vanishing analytic function there. Moreover, u ε vanishes on the boundary of Ω ε , so that u ε is the solution of (4.24).
Set D 0 to be the interior of the square with vertices 0, 1, 1 + i and i in the case of the square lattice and set D 0 to be the interior of the triangle with vertices 0, √ 2/ 4 √ 3 and √ 2e πi/3 / 4 √ 3 in the case of the triangular lattice. The bound (4.26) for h ε is obtained by applying the maximum principle to h ε on Ω ε ∩ D 0 . If h ε were to assume an extremal value on the closure of Ω ε ∩ D 0 at a point of Ω ε ∩ ∂D 0 then, by the symmetry of h ε in the sides of D 0 (see the final part of Lemma 2), h ε would have a local extremum there, contradicting the maximum principle. Thus h ε achieves its extremum values (over Ω ε or, equivalently, over Ω ε ∩ D 0 ) at a point of ∂Ω ε that is, again using the periodicity of h ε , at a point of C(0, ε). Taking account of the boundary values (4.27) and then Lemma 3 we see that, for |ζ| = ε, Thus, by the maximum principle, the harmonic function h ε satisfies the bound (4.26) throughout Ω ε .
Theorem 2. The maximum voltage drop V S max (ε), when the pads are arranged in a square lattice and with the parameters given in (4.3), satisfies where A 4 is given by (4.19) and correct to 12 decimal places. The maximum voltage drop V T max (ε), when the pads are arranged in an equilateral triangular lattice and with the parameters given in (4.3), satisfies where A 6 is given by (4.20) and correct to 12 decimal places.
Proof. In the case of the square arrangement of pads, the maximum voltage drop occurs at the point b s = (1 + i)/2 which lies at the centre of the square formed by the lattice points at 0, 1, 1 + i and i (see Section 6). The negative of the expression (4.25), evaluated at z = b s , is the maximum voltage drop. Since |b s | 2 = 1/2, Formulas 18.14.7 and 18.14.9 in Abramowitz and Stegun [2] give Scaling by t = 2 √ π/Γ 2 ( 1 4 ) so that w 1 = 1/2, and noting that the function σ also scales linearly, we find that Then, log |σ(b s )| = π 4 + log(2 √ 2π) − 2 log Γ( 1 4 ), so that (4.29) follows from (4.32), and then (4.28) follows from the bound (4.26) for h s .
In the case of the triangular pad arrangement, the maximum voltage drop occurs at the point b t = 3 −3/4 √ 2e πi/6 which lies at the centre of the equilateral triangle with vertices 0, 2w 1 = d, Formulas 18.13.15 and 18.13.28 in Abramowitz and Stegun [2] give the value of σ at the centre of the equilateral triangle as and w 3 = iw 1 . Then, , (4.31) follows from (4.33), and then (4.30) again follows from the bound (4.26) for h s .

The Hexagonal configuration
We estimate the voltage drop for the hexagonal lattice with the same aerial density of pads as in the case of the square and the equilateral power pad arrangements analysed in the previous section. The geometric setting is the following. We consider the domain where Λ is the set of vertices of the blue hexagonal grid shown in Figure 2.
It will be convenient to consider the set of centres Λ as the difference of two lattices: see Figure 2. The first lattice consists of the black and the red vertices in Figure 2, which we denote by BR, while the second lattice consists of the red vertices alone, which we denote by R. Thus Λ = BR \ R. Both BR and R are lattices that determine an equilateral grid. The main advantage of considering Λ as a difference of two lattices is that for any equilateral lattice we can construct an associated Weierstrass entire function with zeros on the lattice whose pseudo-periodicity properties were analysed in the previous section. Thus, instead of directly building an entire function with zeros on Λ we obtain more information by considering a quotient of two entire functions, one vanishing on BR and the other on R.
The maximum voltage drop corresponds to the minimum value of u, where u is the solution to ∆u = 1 in Ω ε and u = 0 in the boundary of Ω. The maximum voltage drop is, consequently, at least as big as −u(0) where 0 is at the centre of a hexagon.
Let us denote by σ(z) the Weierstrass σ-function associated with the equilateral triangular lattice with side length d 2 = √ 2/ 4 √ 3 as described in (4.3). The σ-function for the lattice BR, with sidelength d 3 = 2/ 4 √ 27, is then while the σ-function for the lattice R, with sidelength √ 3d 3 , is Clearly σ R vanishes on the vertices of R and σ BR vanishes on the vertices of BR.

Consider the function defined in Ω by
where c ε and d ε are to be chosen appropriately. Both functions v BR and v R have many symmetries. In particular they are symmetric across any line that extends any of the sides of the hexagon which form the original grid. Thus v has the same symmetry. Moreover ∆v = 1 in Ω, so that v is close to the desired solution u to the problem. In fact, they differ by a harmonic function, in that u = v + h. The desired value u(0) can be approximated by the value of v at the centre of the hexagon. The error that we make, that is h(0), can again be estimated by the maximum principle, in that |h(0)| ≤ sup ∂Ω |h| = sup ∂Ω |v|. The constants c ε and d ε will now be chosen so that both sup ∂Ω v BR and sup ∂Ω v R are small. The selection of c ε required to make v BR small on the boundary is the more straightforward. By the symmetries of v BR , sup ∂Ω v BR = sup ∂D(0,ε) v BR . Observe that although ∂D(0, ε) is not part of the boundary of Ω, all disks around the vertices of the combined black and red triangular grid are equal if we restrict our attention to v BR . On ∂D(0, ε) the value of v BR is close to a constant. In fact we see from Lemma 3 that v BR (z) = 3 8 Thus, with the choice of c ε = 1 2π log 1 ε + 3 8 ε 2 , we obtain that |v BR (z)| ≤ Cε 6 on ∂Ω ε . We consider now the values of v R on the boundary of Ω ε which consists of disks of radius ε centred at the baricentres of the red triangles. The function v R has the same behaviour at each. Let us denote one of the baricentres by A. Then, as established in Section 6, v R has a local minimun at A. We can actually prove that Thus, if we choose d ε = 1 8 |A| 2 − 1 2π log |σ R (A)| + 1 8 ε 2 , then |v R (z)| ≤ Cε 3 on ∂D(A, ε) and therefore on ∂Ω ε .
Finally we have proved that sup ∂Ωε |h| = sup ∂Ωε |v| ≤ Cε 3 . The voltage drop at the centre of a hexagon is −u(0) = −v(0) − h(0), and so Thus the voltage drop at the centre of a hexagon is which is (3.3). The conclusion is that the hexagonal grid has the worst voltage drop among the ones that we considered, with the best being the triangular lattice and the standard square lattice being in an intermediate position.

Where does the maximum voltage drop occur?
We examine now where the maximal voltage drop takes place in the square lattice configuration and in the equilateral setting. Heuristically, one expects the voltage drop to be maximal in the center of the squares and in the barycenter respectively. This has been taken for granted in the literature, but we will nevertheless give a rigorous proof of this intuitive fact. The case of the square is the easiest one.
Proof. Consider the solution v in the unbounded domain Ω to the mixed Dirichlet-Neumann problem as in the figure 4: We want to prove that it has a minimum at z = 0. We will prove that the function v y > 0 when ℑz > 0 and v y < 0 when ℑz < 0. Clearly Thus v y is harmonic. Moreover in the "straight" pieces of the boundary v y = 0. On the half circles v = 0, thus ∇v is perpendicular to the circles. Therefore v y = ∇v, In the case of the triangular lattice we consider the domain as in figure 5. The domain Ω is the equilateral triangle where we remove the three disks of equal radius centered at the corners of the triangle. Let p be the barycenter of the triangle and define the function u such that ∆u = 1 in the interior of Ω, u = 0 in the part of the boundary of Ω defined by the arcs of circle and ∂u/∂n = 0 in the part of the boundary of Ω defined by the sides of the triangle. The claim is the following: Claim. There is only a minimum value of u in Ω and it is attained at p.
Proof. We will make this argument by a variation on the radius of the disks. It will be convenient to denote the domains Ω t to the domain obtained removing the disks of radius t and u t the corresponding solution. We will denote by v the Green function of the flat torus whose fundamental domain is twice the equilateral triangle. It follows from the definition that the Green function of this torus is the function v(z) = 1 4 |z| 2 − 1 2π log |σ(z)| as we saw in Lemma 2. In a sense we will see that u t is very close to v as t → 0. We are interested in the critical points of u t . The corresponding critical points for v have been identified in [8] and the only ones appearing are the trivial ones that can be identified by symmetry considerations. There is a local minimum of v at p and three saddle points in the midpoints of the sides of the triangle.
We are going to prove that a very similar structure arises in the case of u t : There is a minimum at p and three saddle points at the midpoints of the sides of the triangle.
Along all this discussion we will restrict ourselves to the case 0 < t < t 0 where t 0 is the biggest radius such that the disks defining Ω t are disjoint since this is the only relevant case.
We start by observing that at the barycenter p there is a critical point for u t for symmetry reasons. Moreover since u t (e 2πi/3 (z − p)) = u t (z − p) the Hessian of u at p must be a constant times the identity matrix. Since ∆u(p) = 1 it follows that u xx (p) = u yy (p) = 1/2.
Let d 1 , d 2 , d 3 be vectors pointing from p to one of the vertex of the triangle as in Figure 5. By symmetry again the gradient of u t in any point of the median of the triangle is a multiple of d j .
Assume, for the moment being, that there is a δ such that for a given t < δ we have proven that u t d j (x) > 0 for any x in the median joining p with a vertex (excluding the barycenter), i.e., along the median the gradient is pointing towards the vertices.
Under this assumption we concentrate our attention on the yellow region in the picture consisting of one third of the original domain Ω t limited by two of the medians. We will prove that on the yellow region the function u o 1 which is the derivative of u in the direction o 1 := −d 3 if strictly positive. This is clear because the function u t o 1 is an harmonic function (∆u t o 1 = 0) and in the boundary of the shaded region it is positive: on the medians it is positive by assumption, on the sides of the triangle it is actually 0 by the definition of u t and on the arcs of circles the gradient of u t is pointing towards the center of the disks (u t ≡ 0 on the boundary of the disks and it is negative in Ω t ), thus u o 1 is positive on the arcs of circle that limit the shaded region. Now any point q belonging to the yellow region has the property that u(p) < u(q) since we can follow a path from p to q consisting of segment over the median followed by a segment in the direction of o 1 and in both segments u t will be increasing.
It remains to prove that u t d j (x) ≥ 0 on the corresponding median. Let us assume for the moment being that this is the case for all t ≤ δ. We will prove then that this is true for all t < t 0 .
Let us denote by t * the biggest t such that u t d j (x) ≥ 0 on all points of the median. We will see now that if t * < t 0 we reach a contradiction. By continuity u t * d j (x) ≥ 0 on the median. If we prove that actually u t * d j (x) > δ > 0 (6.1) on the median we would have reached a contradiction since t * would not be maximal. We cannot prove (6.1) directly since u d j (p) = 0, but in a neighborhood of p u d j (x) > u d j (p) since u t d j d j (p) = 1/2. Thus if t * is maximal it maybe only for two reasons. Either there is a point q in the interior of the median different from p such that u t * d j (q) = 0 or the same thing happens for the point q ′ that is in intersection of the median with the boundary of Ω t . Let us examine these two cases separately. In the first case u t * d j ≥ 0 along the median but it vanishes in some intermediate position. By symmetry it will happen in u t * d 1 and u t * d 2 simultaneously. Thus u t * o 1 is an harmonic function in the yellow region that it is positive in the boundary (and strictly positive on some points in the boundary, for instance near p). Thus, by the maximal principle, it is a strictly positive function in the interior of the yellow region. Thus u t * o 1 is positive in the median that bisects the yellow region. By symmetry again u t * o 3 is positive in the piece of the median denoted by o 3 in the picture. Therefore finally u t * d 1 ≥ 0 on the region delimitated by o 1 , o 2 and Ω t . Finally since u t * d 1 is harmonic it implies that it is strictly positive on the interior, i.e. on the median d 1 . Thus such q does not exist.
On the other hand u t * d 1 cannot vanish on the endpoint e where the median d 1 meets the circle because we are assuming that t * < t 0 and therefore the expected lifetime near the boundary of the disk can be estimated from below by the expected lifetime of a corona around the disk. This has an explicit expression that has positive derivative on the boundary. Thus u t * d 1 (e) > 0. We have reached a contradiction. It only remains to prove that we can start the argument, i.e. that that there is a δ such that for a given t < δ we have that u t d j (x) > 0 for any x in the median joining p with a vertex (excluding the barycenter). This is the case when t = 0. In this case we define u 0 = v, the Green function. In this case the gradient v d j > 0 along the median because by the results of [8] v has p as unique critical point in the interior of Ω 0 . For very small t the Green function v has values in the circles around the vertices of the triangle very close to a constant. Thus u t can be obtained by correcting u t with an harmonic function that in the circles almost coincides with a constant. One can check that u d j is close to v d j and thus it is positive if t is small enough.