On the solution of Laplace's equation in the vicinity of triple-junctions

TL;DR

This paper analyzes the behavior of solutions to Laplace boundary integral equations near triple junctions in composite media, providing a theoretical characterization and developing efficient numerical discretization methods.

Contribution

It characterizes solution behavior near triple junctions and introduces an efficient discretization approach for boundary integral equations in such regions.

Findings

Solutions near triple junctions are well-approximated by functions of the form t^β.

The powers β depend only on material properties and boundary angles.

The proposed discretization method demonstrates high accuracy and efficiency.

Abstract

In this paper we characterize the behavior of solutions to systems of boundary integral equations associated with Laplace transmission problems in composite media consisting of regions with polygonal boundaries. In particular we consider triple junctions, i.e. points at which three distinct media meet. We show that, under suitable conditions, solutions to the boundary integral equations in the vicinity of a triple junction are well-approximated by linear combinations of functions of the form $t^\beta,$ where $t$ is the distance of the point from the junction and the powers $\beta$ depend only on the material properties of the media and the angles at which their boundaries meet. Moreover, we use this analysis to design efficient discretizations of boundary integral equations for Laplace transmission problems in regions with triple junctions and demonstrate the accuracy and efficiency of this algorithm with a number of examples.

Equations

$$\Delta u_i$$

$${\mu }_{\ell (i)}{u}_{\ell (i)}-{\mu }_{r(i)}{u}_{r(i)}$$

$${\nu }_{\ell (i)}\frac{\partial u_{\ell (i)}}{\partial n_i}-{\nu }_{r(i)}\frac{\partial u_{r(i)}}{\partial n_i}$$

$$\underset{|r|\to \infty }{\lim }\left(r\log (r)u_0^{\prime }(r)-u_0(r)\right)$$

$${\gamma }_i(t)={\mathbf{v}}_{\ell }+t\frac{{\mathbf{v}}_r-{\mathbf{v}}_{\ell }}{\Vert {\mathbf{v}}_r-{\mathbf{v}}_{\ell }\Vert }.$$

$${\mathbf{n}}_i=\frac{({\mathbf{v}}_r-{\mathbf{v}}_{\ell }{)}^{\perp }}{\Vert {\mathbf{v}}_r-{\mathbf{v}}_{\ell }\Vert },$$

$${\mathcal{S}}_{{\Gamma }_i}[\sigma ](\mathbf{y})=-\frac{1}{2\pi }{\int }_{{\Gamma }_i}\log \Vert \mathbf{x}-\mathbf{y}\Vert \sigma (\mathbf{x})d{S}_{\mathbf{x}}$$

$${\mathcal{D}}_{{\Gamma }_i}[\sigma ](\mathbf{y})=\frac{1}{2\pi }{\int }_{{\Gamma }_i}\frac{\mathbf{n}(\mathbf{x})\cdot (\mathbf{y}-\mathbf{x})}{\Vert \mathbf{x}-\mathbf{y}{\Vert }^2}\sigma (\mathbf{x})d{S}_{\mathbf{x}},$$

$${\mathcal{D}}_{{\Gamma }_i}^{\ast }[\sigma ](\mathbf{y})=\frac{1}{2\pi }{\int }_{{\Gamma }_i}\frac{\mathbf{n}(\mathbf{y})\cdot (\mathbf{x}-\mathbf{y})}{\Vert \mathbf{x}-\mathbf{y}{\Vert }^2}\sigma (\mathbf{x})d{S}_{\mathbf{x}},$$

$$K(\mathbf{x},\mathbf{y})=\frac{1}{2\pi }\frac{\mathbf{n}(\mathbf{x})\cdot (\mathbf{y}-\mathbf{x})}{\Vert \mathbf{x}-\mathbf{y}{\Vert }^2}.$$

$$-1+\alpha <\frac{\mathbf{x}-{\mathbf{x}}_0}{\Vert \mathbf{x}-{\mathbf{x}}_0\Vert }\cdot {\gamma }_i^{\prime }(t_0)<1-\alpha$$

$$\underset{\mathbf{x}\to {\mathbf{x}}_0^{\pm }}{\lim }{\mathcal{S}}_{{\Gamma }_i}[\sigma ](\mathbf{x})={\mathcal{S}}_{{\Gamma }_i}[\sigma ]({\mathbf{x}}_0)$$

$$\underset{\mathbf{x}\to {\mathbf{x}}_0^{\pm }}{\lim }{\mathcal{D}}_{{\Gamma }_i}[\rho ](\mathbf{x})=p.v.{\mathcal{D}}_{{\Gamma }_i}[\rho ]({\mathbf{x}}_0)\mp \frac{\rho ({\mathbf{x}}_0)}{2}$$

$$\underset{\mathbf{x}\to {\mathbf{x}}_0^{\pm }}{\lim }{\mathbf{n}}_i\cdot \nabla {\mathcal{S}}_{{\Gamma }_i}[\rho ](\mathbf{x})=p.v.{\mathcal{D}}_{{\Gamma }_i}^{\ast }[\rho ]({\mathbf{x}}_0)\pm \frac{\rho ({\mathbf{x}}_0)}{2},$$

$$\underset{\mathbf{x}\to {\mathbf{x}}_0^{\pm }}{\lim }{\mathbf{n}}_i\cdot \nabla {\mathcal{D}}_{{\Gamma }_i}[\rho ](\mathbf{x}),$$

$$u_i(\mathbf{x})=\frac{1}{{\mu }_i}\sum _{j=1}^m{\mathcal{S}}_{{\Gamma }_j}[{\rho }_j](\mathbf{x})+\frac{1}{{\nu }_i}\sum _{j=1}^m{\mathcal{D}}_{{\Gamma }_j}[{\sigma }_j](\mathbf{x})\mathbf{x}\in {\Omega }_i.$$

$$-\frac{1}{2}{\sigma }_i+\frac{{\mu }_{r(i)}{\nu }_{\ell (i)}-{\mu }_{\ell (i)}{\nu }_{r(i)}}{{\mu }_{r(i)}{\nu }_{\ell (i)}+{\mu }_{\ell (i)}{\nu }_{r(i)}}\sum _{\ell =1}^m{\mathcal{D}}_{{\Gamma }_{\ell }}[{\sigma }_{\ell }]$$

$$-\frac{1}{2}{\rho }_i+\frac{{\mu }_{r(i)}{\nu }_{\ell (i)}-{\mu }_{\ell (i)}{\nu }_{r(i)}}{{\mu }_{r(i)}{\nu }_{\ell (i)}+{\mu }_{\ell (i)}{\nu }_{r(i)}}\sum _{\ell =1}^m{\mathcal{D}}_{{\Gamma }_{\ell }}^{\ast }[{\rho }_{\ell }]$$

$${\int }_{\Gamma \setminus B_r({\mathbf{x}}_p)}K(\mathbf{x},\mathbf{y})\sigma (\mathbf{x})d{S}_{\mathbf{x}},{\int }_{\Gamma \setminus B_r({\mathbf{x}}_p)}K(\mathbf{y},\mathbf{x})\rho (\mathbf{x})d{S}_{\mathbf{x}},$$

$${\mathcal{D}}_{(\ell ,m);({\ell }^{\prime },m^{\prime })}[\sigma ](t)=p.v.{{\mathcal{D}}_{{\Gamma }_{(\ell ,m)}}[\sigma ]|}_{{\Gamma }_{({\ell }^{\prime },m^{\prime })}}$$

$${\mathcal{D}}_{(\ell ,m);({\ell }^{\prime },m^{\prime })}^{\ast }[\rho ](t)=p.v.{{\mathcal{D}}_{{\Gamma }_{(\ell ,m)}}^{\ast }[\rho ]|}_{{\Gamma }_{({\ell }^{\prime },m^{\prime })}}$$

$$d_{(1,2)}$$

$$d_{(2,3)}$$

$$d_{(3,1)}$$

$${\mathcal{K}}_{\text{dir}}[\mathbf{h}{]}_{(i,j)}=-\frac{1}{2}{h}_{(i,j)}+d_{(i,j)}\sum _{(\ell ,m)\in \mathcal{J}}{\mathcal{D}}_{(\ell ,m);(i,j)}[h_{(\ell ,m)}],$$

$${\mathcal{A}}_{\text{dir}}(a,b,\beta )=\left(\begin{array}{ccc}\sin (\pi \beta )& b\sin \beta (\pi -{\theta }_2)& -b\sin \beta (\pi -{\theta }_1)\\ (a+b)/(1+ab)\sin \beta (\pi -{\theta }_2)& \sin (\pi \beta )& -(a+b)/(1+ab)\sin \beta (\pi -{\theta }_3)\\ a\sin \pi \beta (1-{\theta }_1)& -a\sin \pi \beta (1-{\theta }_3)& \sin (\pi \beta )\end{array}\right).$$

$$h_{(i,j)}(t)$$

$$\sigma =\left[\begin{array}{c}{\sigma }_{1,2}(t)\\ {\sigma }_{2,3}(t)\\ {\sigma }_{3,1}(t)\end{array}\right]=\sum _{i=0}^N\sum _{j=0}^2{p}_{i,j}{\mathbf{v}}_{i,j}{t}^{{\beta }_{i,j}},$$

$$\underset{(i,j)\in \mathcal{J}}{\max }\left|h_{(i,j)}-{\mathcal{K}}_{\text{dir}}[\sigma {]}_{(i,j)}\right|\le C{t}^{N+1},$$

$${\mathcal{A}}_{\text{neu}}(a,b,\beta )=\left(\begin{array}{ccc}\sin (\pi \beta )& -b\sin \beta (\pi -{\theta }_2)& b\sin \beta (\pi -{\theta }_1)\\ -(a+b)/(1+ab)\sin \beta (\pi -{\theta }_2)& \sin (\pi \beta )& (a+b)/(1+ab)\sin \beta (\pi -{\theta }_3)\\ -a\sin \beta (\pi -{\theta }_1)& a\sin \pi \beta (1-{\theta }_3)& \sin (\pi \beta )\end{array}\right).$$

$$h_{(i,j)}(t)$$

$$\sigma =\left[\begin{array}{c}{\sigma }_{1,2}(t)\\ {\sigma }_{2,3}(t)\\ {\sigma }_{3,1}(t)\end{array}\right]=\sum _{i=1}^{N+1}\sum _{j=0}^2{p}_{i,j}{\mathbf{w}}_{i,j}{t}^{{\beta }_{i,j}-1},$$

$$\underset{(i,j)\in \mathcal{J}}{\max }\left|h_{(i,j)}-{\mathcal{K}}_{\text{neu}}[\sigma {]}_{(i,j)}\right|\le C{t}^{N+1},$$

$$Y=\{({\theta }_1,{\theta }_2,a,b):0<{\theta }_1,{\theta }_2<2\pi ,{\theta }_1+{\theta }_2<2\pi ,-1<a,b<1\}.$$

Full text

1 Introduction

Composite media, i.e. media consisting of multiple materials in close proximity or contact, are both ubiquitous in nature and fascinating in applications since their macroscopic properties can be substantially different than those of their components. One property of particular interest is the electrostatic response of composite media, typically the electric potential in the medium which is produced by an externally-applied time-independent electric field. In such situations one often assumes that the associated electric potential satisfies Laplace’s equation in the interior of each medium and that along each edge where two media meet one prescribes the jump in the normal derivative of the potential. Typically the potentials in these jump relations appear multiplied by coefficients depending on the electric permittivity. This leads to a collection of coupled partial differential equations (PDEs). In addition to classical electrostatics problems, the same equations also arise in, among other things, percolation theory, homogenization theory, and the study field enhancements in vacuum insulators (see, for example,  [1, 2, 3, 4, 5, 6, 7]).

Using classical potential theory this set of partial differential equations (PDEs) can be reduced to a system of second-kind boundary integral equations (BIEs). In particular, the solution to the PDE in each region is represented as a linear combination of a single-layer and a double-layer potential on the boundary of each subregion. If the edges of the media are smooth then the corresponding kernels in the integral equation are as well. Near corners, however, the solutions to both the differential equations and the integral equations can develop singularities.

Analytically, the behavior of solutions to both the PDEs and BIEs have been the subject of extensive analysis (see, for example  [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]). In particular, the existence and uniqueness of solutions in an $\mathbb{L}^{2}$-sense is well-known, under certain natural assumptions on the material properties [18, 19]. Moreover, the asymptotic form of the singularities in the vicinity of a junction has been determined for the solutions of both the PDE and its corresponding BIE [11, 12, 8, 20, 21].

Computationally the singular nature of the solutions poses significant challenges for many existing numerical methods for solving both the PDEs and BIEs. Typical approaches involve introducing many additional degrees of freedom near the junctions which can impede the speed of the solver and impose prohibitive limits on the size and complexity of geometries which can be considered. Recursive compressed inverse preconditioning (RCIP) is one way of circumventing the difficulty introduced by the presence of junctions in the BIE formulation [22]. In this approach, the extra degrees of freedom introduced by the refinement near the junctions are eliminated from the linear system. Moreover, the compression and refinement are performed concomitantly for multiple junctions in parallel. This approach gives an algorithm which scales linearly in the number of degrees of freedom added to resolve the singularities near the junction. The resulting linear system has essentially the same number of degrees of freedom as it would if the junctions were absent.

In this paper we restrict our attention to the case of triple junctions, extending the existing analysis by showing that under suitable restrictions the solution to the BIEs can be well-approximated in the vicinity of a triple junction by a linear combination of $t^{\beta_{j}}$, where $t$ is the distance from the triple junction and the $\beta_{j}$’s are a countable collection of real numbers defined implicitly by an equation depending only on the angles at which the interfaces meet and the material properties of the corresponding media. This analysis enables the construction of an efficient computational algorithm for solving Laplace’s equation in regions with multiple junctions. In particular, using this representation we construct an accurate and efficient quadrature scheme for the BIE which requires no refinement near the junction. The properties of this discretization are illustrated with a number of numerical examples.

This paper is organized as follows. In section 2 we state the boundary value problem for the Laplace triple junction transmission problem, summarize relevant properties of layer potentials, and describe the reduction of the boundary value problem to a system of boundary integral equations. In section 3 we present the main theoretical results of this work, the proofs of which are given in appendices A and B. In section 4 we discuss two conjectures extending the results of section 3 based on extensive numerical evidence. In section 5, we describe a Nyström discretization which exploits explicit knowledge of the structure of solutions to the integral equations in the vicinity of triple junctions, and in section 6 we demonstrate its effectiveness of numerical solvers. Finally, in section 7 we summarize the results and outline directions for future research.

2 Boundary value problem

Consider a composite medium consisting of a set of $n$ polygonal domains $\Omega_{1},\dots,\Omega_{n}$ (see Figure 1) with boundaries consisting of $m$ edges $\Gamma_{1},\dots,\Gamma_{m}$ and $k$ vertices ${\boldsymbol{v}}_{1},\dots,{\boldsymbol{v}}_{k}$. For a given edge $\Gamma_{i}$ let $L_{i}$ denote its length, $\boldsymbol{n}_{i}$ its normal, $\ell(i),r(i)$ the polygons to the left and right, respectively, and $\gamma_{i}$ be an arclength parameterization of $\Gamma_{i}$. Finally we denote the union of the regions $\Omega_{1},\dots,\Omega_{n}$ by $\Omega$ and denote the complement of $\Omega$ by $\Omega_{0}.$

Given positive constants $\mu_{1},\dots,\mu_{n}$ and $\nu_{1},\dots,\nu_{n}$ we consider the following boundary value problem

[TABLE]

where $f_{i}$ and $g_{i}$ are analytic functions on $\Gamma_{i},$ $i=1,\dots,m$, and $\ell(i),r(i)$ denote the regions on the left and right with respect to the normal of edge $\Gamma_{i}$.

Remark 2.1**.**

In this work we assume that all the normals $\boldsymbol{n}_{1},\dots,\boldsymbol{n}_{m}$ to $\Gamma_{1},\dots,\Gamma_{m}$ are positively oriented with respect to the parameterization $\gamma_{i}(t)$ of the edge $\Gamma_{i}$. Specifically, if $\Gamma_{i}$ is a line segment between vertices ${\boldsymbol{v}}_{\ell}$, ${\boldsymbol{v}}_{r}$, and $\gamma_{i}(t):[0,L_{i}]\to\Gamma_{i}$ is a parameterization of $\Gamma_{i}$, given by

[TABLE]

Then the normal on edge $\Gamma_{i}$, is given by

[TABLE]

where for a point $\boldsymbol{x}=(x_{1},x_{2})\in\mathbb{R}^{2}$, $\boldsymbol{x}^{\perp}=(x_{2},-x_{1})$.

Remark 2.2**.**

The existence and uniqueness of solutions to (1) is a classical result [19].

Remark 2.3**.**

In this paper we assume that no more than three edges meet at each vertex. Similar analysis holds for domains with higher-order junctions and will be published at a later date.

Remark 2.4**.**

Here we assume that $\mu_{1},\dots,\mu_{n}$, and $\nu_{1},\dots,\nu_{n}$ are positive constants. In principle the analysis presented here extends to the case where the constants are negative or complex provided the the constants $(\mu_{j}\nu_{i}+\mu_{i}\nu_{j})/(\mu_{j}\nu_{i}-\mu_{i}\nu_{j})$ across each edge are outside the closure of the essential spectrum of the double layer potential defined on the boundary, and the underlying differential equation admits a unique solution. Note that for non-negative coefficients this is always true, since these constants are all in magnitude greater than $1$, and the spectral radius of the double layer potential is bounded by $1$.

2.1 Layer potentials

Before reducing the boundary value problem (1) to a boundary integral equation we first introduce the layer potential operators and summarize their relevant properties.

Definition 2.1**.**

Given a density $\sigma$ defined on $\Gamma_{i},$ $i=1,\dots,m,$ the single-layer potential is defined by

[TABLE]

and the double-layer potential is defined via the formula

[TABLE]

Remark 2.5**.**

In light of the previous definition, evidently the adjoint of the double-layer potential is given by the formula

[TABLE]

Definition 2.2**.**

For $\boldsymbol{x}\in\Gamma$ we define the kernel $K(\boldsymbol{x},\boldsymbol{y})$ by

[TABLE]

The following theorems describe the limiting values of the single and double layer potential on the boundary $\Gamma_{i}$.

Theorem 2.1**.**

Suppose that $\boldsymbol{x}_{0}$ is a point in the interior of the segment $\Gamma_{i}$. Suppose the point $\boldsymbol{x}$ approaches a point $\boldsymbol{x}_{0}$ along a path such that

[TABLE]

for some $\alpha>0$. If $(\boldsymbol{x}-\boldsymbol{x}_{0})\cdot\boldsymbol{n}_{i}<0$, we will refer to this limit as $\boldsymbol{x}\to\boldsymbol{x}_{0}^{-}$, and if $(\boldsymbol{x}-\boldsymbol{x}_{0})\cdot\boldsymbol{n}_{i}>0$, we will refer to this limit as $\boldsymbol{x}\to\boldsymbol{x}_{0}^{+}$.

Then

[TABLE]

where p.v. refers to the fact that the principal value of the integral should be taken.

Moreover, both the limits

[TABLE]

exist and are equal.

Remark 2.6**.**

In the following we will suppress the ${\rm p.v.}$ from expressions involving layer potentials evaluated at a point on the boundary. Unless otherwise stated, in such cases the principal value should always be taken.

2.2 Integral representation

In classical potential theory the boundary value problem (1) is reduced to a boundary integral equation for a new collection of unknowns $\rho_{i},\sigma_{i}\in\mathbb{L}^{2}(\Gamma_{i})$, $i=1,\dots,m$ related to $u_{i}:\Omega_{i}\to\mathbb{R}$, $i=1,\dots,n$ in the following manner

[TABLE]

We note that by construction $u_{i}$ is harmonic in $\Omega_{i}$, $i=0,1,\dots,n$. Enforcing the jump conditions across the edges and applying theorem 2.1 yields the following system of integral equations for the unknown densities $\rho_{i}$ and $\sigma_{i}$

[TABLE]

for $i=1,\dots,m.$

We note that the preceding representation has several advantages. Firstly, the kernels of integral equations (14) and (15) are smooth except at the vertices. In particular, the weakly-singular and hypersingular terms arising from the single-layer potential and the derivative of the double-layer potential, respectively, are absent. Secondly, the equations for the single-layer density $\rho$ and the double-layer density $\sigma$ are completely decoupled and can be analyzed separately. Moreover, (15) is the adjoint of (14) and hence the structure of solutions to (15) can be inferred from the behavior of solutions to (14).

Remark 2.7**.**

The above representation also appears in [21] and is related to the work in [16]. It has been shown in [18] that the boundary integral equations (14),(15) are well-posed for $f_{i},g_{i}\in\mathbb{L}^{2}[\Gamma_{i}]$.

2.3 The single-vertex problem

The following lemma reduces the problem of analyzing the behavior of the densities $\rho$ and $\sigma$ in the vicinity of a triple junction with locally-analytic data to the analysis of an integral equation on a set of three intersecting line segments.

Lemma 2.1**.**

Let $\sigma,\rho$ satisfy the boundary integral equation (14) and (15), respectively. Consider three edges $\Gamma_{i},$ $\Gamma_{j},$ and $\Gamma_{k}$ meeting at a vertex $v_{p}$. If $\boldsymbol{x}_{p}$ denotes the coordinates of the vertex $v_{m}$ then there exists an $r>0$ such that

[TABLE]

are analytic functions of $\boldsymbol{y}$ for all $\boldsymbol{y}\in B_{r}(\boldsymbol{x}_{p})$. Here $B_{r}(\boldsymbol{x}_{p})$ denotes the ball of radius $r$ centered at $\boldsymbol{x}_{p}$.

Remark 2.8**.**

We note that by choosing $r$ sufficiently small we can assume that the intersection of all three-edges with $B_{r}(\boldsymbol{x}_{p})$ are of length $r$. Moreover, since Laplace’s equation is invariant under scalings the subproblem associated with the corner can be mapped to an integral equation on three intersecting edges of unit length.

In light of the preceding remark, in the remainder of this paper we restrict our attention to the geometry shown in Figure 2.

The following notation will be used in our analysis of triple junctions.

Remark 2.9**.**

Suppose that $\Gamma_{(\ell,m)}$ and $\Gamma_{(\ell^{\prime},m^{\prime})}$ are two (possibly identical) edges of a triple junction in which all edges are of length one. For $(\ell,m)$ and $(\ell^{\prime},m^{\prime})$ in $\{(1,2),(2,3),(3,1)\}$ and $t\in(0,1)$ let

[TABLE]

and

[TABLE]

for any $\sigma,\rho\in\mathbb{L}^{2}(\Gamma_{(3,1)}\cup\Gamma_{(1,2)}\cup\Gamma_{(2,3)}).$ Note that if $(\ell,m)=(\ell^{\prime},m^{\prime})$ then both quantities are identically zero for any $\sigma$ and $\rho.$ If $(\ell,m)\neq(\ell^{\prime},m^{\prime})$ then the principal value is not required.

Finally, in the following we will also denote the restrictions of $\sigma$ and $\rho$ to an edge $\Gamma_{(\ell,m)}$ by $\sigma_{(\ell,m)}$ and $\rho_{(\ell,m)},$ respectively.

3 Main results

In this section we state several theorems which characterize the behavior of the solutions $\sigma,\rho$ to  eqs. 14 and 15 for the single-vertex problem with piecewise smooth boundary data $f$ and $g.$ Before doing so we first introduce some convenient notation. To that end, let $\Gamma_{(1,2)},\Gamma_{(2,3)}$ and $\Gamma_{(3,1)}$ be three edges of unit length meeting at a vertex as in Figure 2. Let $\theta_{1},$ $\theta_{2},$ and $\theta_{3}$ be the angles at which they meet and suppose that $0<\theta_{1},\theta_{2},\theta_{3}<2\pi$ are real numbers summing to $2\pi$. Let $\Omega_{1}$ denote the region bordered by $\Gamma_{(3,1)}$ and $\Gamma_{(1,2)},$ $\Omega_{2}$ the region bordered by $\Gamma_{(1,2)}$ and $\Gamma_{(2,3)},$ and $\Omega_{3}$ denote the region bordered by $\Gamma_{(2,3)}$ and $\Gamma_{(3,1)}$. Finally, let $\mu_{i}$ and $\nu_{i}$ be the parameters corresponding to $\Omega_{i},$ $i=1,2,3$ and define the constants $d_{(1,2)},d_{(2,3)}$ and $d_{(3,1)}$ by

[TABLE]

Remark 3.1**.**

We note the following properties of $d_{(3,1)},d_{(1,2)},d_{(2,3)}$ which, for notational convenience, we will denote by $a,$ $b,$ and $c,$ respectively. Firstly, since $\mu_{i},\nu_{i}$ are positive real numbers, it follows that $a,b,c\in(-1,1)$. Secondly, a simple calculation shows that $c=-(a+b)/(1+ab)$. Thus, at each triple junction, there are two parameters $(a,b)$ which encapsulate the relevant information regarding material properties at that junction. For the rest of the paper, in a slight abuse of notation, we will refer to $(a,b)$ as the material parameters.

Next we define several quantities which will be used in the statement of the main results. Let $\mathcal{J}$ denote the set of indices $\{(1,2),(2,3),(3,1)\}$ and $X=\mathbb{L}^{2}(\Gamma_{(1,2)})\otimes\mathbb{L}^{2}(\Gamma_{(2,3)})\otimes\mathbb{L}^{2}(\Gamma_{(3,1)})$. Let $\mathcal{K}_{\text{dir}}:X\to X$, and $\mathcal{K}_{\text{neu}}:X\to X$ denote the bounded operators in eq. 14 and eq. 15 respectively. For any operator $A:X\to X$, $\boldsymbol{h}\in X$, and $(i,j)\in\mathcal{J}$, we denote the restriction of $A[\boldsymbol{h}]$ to the edge $\Gamma_{(i,j)}$ by $A[\boldsymbol{h}]_{(i,j)}$. For example, given $\boldsymbol{h}(t)=[h_{(1,2)}(t),h_{(2,3)}(t),h_{(3,1)}(t)]^{T}\in X$, and $(i,j)\in\mathcal{J}$,

[TABLE]

where the operators $\mathcal{D}_{(\ell,m);(i,j)}$ are defined in eq. 17.

We are interested in the following two problems:

for what collection of $\boldsymbol{h}\in X$ are $\mathcal{K}_{\text{dir}}[\boldsymbol{h}]$, and $\mathcal{K}_{\text{neu}}[\boldsymbol{h}]$ piecewise smooth functions on each of the edges $\Gamma_{(i,j)}$, $(i,j)\in\mathcal{J}$; and 2. 2.

given $h_{(i,j)}\in\mathcal{P}_{N}$, a polynomial of degree at most $N$, construct an explicit basis for $\mathcal{K}_{\text{dir}}^{-1}[\boldsymbol{h}]$ and $\mathcal{K}_{\text{neu}}^{-1}[\boldsymbol{h}]$   .

In section 3.1, we address these questions for $\mathcal{K}_{\text{dir}}$, while in section 3.2 we present analogous results for $\mathcal{K}_{\text{neu}}$.

3.1 Analysis of $\mathcal{K}_{\text{dir}}$

Suppose that $\boldsymbol{h}(t)=[h_{(1,2)}(t),h_{(2,3)}(t),h_{(3,1)}(t)]^{T}={\boldsymbol{v}}t^{\beta}$, where $t$ denotes the distance along the edge $\Gamma_{(i,j)}$ from the triple junction, and ${\boldsymbol{v}}\in\mathbb{R}^{3}$ and $\beta\in\mathbb{R}$ are constants.

In the following theorem, we derive necessary conditions on $\beta,{\boldsymbol{v}}$ such that $\mathcal{K}_{\text{dir}}[\boldsymbol{h}]_{(i,j)}$ is a smooth function on each edge $\Gamma_{(i,j)}$, $(i,j)\in\mathcal{J}$.

Theorem 3.1**.**

Let $\mathcal{A}_{\text{dir}}(a,b,\beta)\in\mathbb{R}^{3\times 3}$ denote the matrix given by

[TABLE]

Suppose that $\beta$ is a positive real number such that $\det{\mathcal{A}_{\text{dir}}(d_{(3,1)},d_{(1,2)},\beta)}=0$ and that ${\boldsymbol{v}}$ is a null-vector of $\mathcal{A}_{\text{dir}}(d_{(3,1)},d_{(1,2)},\beta)$. Let $\boldsymbol{h}(t)={\boldsymbol{v}}t^{\beta}$, $0<t<1$. Then $\mathcal{K}_{\text{dir}}[\boldsymbol{h}]_{(i,j)}$ is an analytic function of $t$, for $0<t<1$, on each of the edges $\Gamma_{(i,j)}$, $(i,j)\in\mathcal{J}$.

The above theorem guarantees that for appropriately chosen densities $\boldsymbol{h}\in X$, the potential $\mathcal{K}_{\text{dir}}[\boldsymbol{h}]$ is an analytic function on each of the edges.

We now consider the construction of a basis for $\mathcal{K}_{\text{dir}}^{-1}[\boldsymbol{h}]$, when $h_{(i,j)}\in\mathcal{P}_{N}$, $(i,j)\in\mathcal{J}$ for some $N>0$.

In order to prove this result, we require a collection of $\beta,{\boldsymbol{v}}$ satisfying the conditions of theorem 3.1. The following lemma states the existence of a countable collection of $\beta,{\boldsymbol{v}}$ which are analytic on a subset of $(-1,1)^{2}$.

Lemma 3.1**.**

Suppose that $\theta_{1},\theta_{2},\theta_{3}$ are irrational numbers summing to $2\pi$, and $(a,b)\in(-1,1)^{2}$. Then there exists a countable collection of open subsets of $(-1,1)^{2},$ denoted by $S_{i,j},$ as well as a corresponding set of functions $\beta_{i,j}:S_{i,j}\to\mathbb{R}$, $i=0,1,2,\ldots$, $j=0,1,2$ such that $\det{\mathcal{A}_{\text{dir}}}(a,b,\beta_{i,j})=0$ for all $(a,b)\in S_{i,j}$. The corresponding null-vectors ${\boldsymbol{v}}_{i,j}:S_{i,j}\to\mathbb{R}^{3}$ of $\mathcal{A}_{\text{dir}}(a,b,\beta_{i,j})$ are also analytic functions. Finally, for any $N>0$, $|\cap_{i=0}^{N}\cap_{j=0}^{2}S_{i,j}|>0$.

In the following theorem, we present the main result of this section which gives a basis for $\mathcal{K}_{\text{dir}}^{-1}[\boldsymbol{h}]$.

Theorem 3.2**.**

Consider the same geometry as in fig. 2, where $\theta_{1},\theta_{2}$, and $\theta_{3}$ sum to $2\pi$ and $\theta_{1}/\pi,$ $\theta_{2}/\pi,$ and $\theta_{3}/\pi$ are irrational. Let $\beta_{i,j},{\boldsymbol{v}}_{i,j},S_{i,j}$, $i=0,1,2,\ldots$, $j=0,1,2$ be as defined in lemma 3.1, and for any positive integer $N$, let $S_{N}$ denote the region of common analyticity of $\beta_{i,j},{\boldsymbol{v}}_{i,j}$, i.e., $S_{N}=\cap_{i=0}^{N}\cap_{j=0}^{2}S_{i,j}$. Finally, suppose that $h_{(i,j)}^{k}$, $(i,j)\in\mathcal{J},k=0,1,2,\ldots N$ are real constants, and define $h_{(i,j)}$ by

[TABLE]

$0<t<$1. Then there exists an open region $\tilde{S}_{N}\subset S_{N}\subset(-1,1)^{2}$ with $|\tilde{S}_{N}|>0$ such that the following holds. For all $(a,b)\in\tilde{S}_{N}$, there exist constants $p_{i,j}$, $i=0,1,\ldots N$, $j=0,1,2$, such that

[TABLE]

satisfies

[TABLE]

for $0<t<1$, where $C$ is a constant.

3.2 Analysis of $\mathcal{K}_{\text{neu}}$

Suppose that $\boldsymbol{h}(t)=[h_{(1,2)}(t),h_{(2,3)}(t),h_{(3,1)}(t)]^{T}={\boldsymbol{w}}t^{\beta-1}$, where $t$ denotes the distance on the edge $\Gamma_{(i,j)}$ from the triple junction, and ${\boldsymbol{w}}\in\mathbb{R}^{3}$ and $\beta$ are constants. In the following theorem, we discuss necessary conditions on $\beta,{\boldsymbol{w}}$ guaranteeing that $\mathcal{K}_{\text{neu}}[\boldsymbol{h}]_{(i,j)}$ is a smooth function on each edge $\Gamma_{(i,j)}$, $(i,j)\in\mathcal{J}$.

Theorem 3.3**.**

Let $\mathcal{A}_{\text{neu}}(a,b,\beta)\in\mathbb{R}^{3\times 3}$ denote the matrix given by

[TABLE]

Suppose that $\beta$ is a positive real number such that $\det{\mathcal{A}_{\text{neu}}(d_{(3,1)},d_{(1,2)},\beta)}=0$ and let ${\boldsymbol{w}}$ denote a corresponding null-vector of $\mathcal{A}_{\text{neu}}(d_{(3,1)},d_{(1,2)},\beta)$. Let $h={\boldsymbol{w}}t^{\beta-1}$, $0<t<1$. Then $\mathcal{K}_{\text{neu}}[h]_{(i,j)}$ is an analytic function of $t$, for $0<t<1$, on each of the edges $\Gamma_{(i,j)}$, $(i,j)\in\mathcal{J}$.

Before proceeding a few remarks are in order.

Remark 3.2**.**

We note that $\det{\mathcal{A}_{\text{dir}}(a,b,\beta)}=\det{\mathcal{A}_{\text{neu}}(a,b,\beta)}$. Thus, the existence of $\beta,{\boldsymbol{w}}$, which satisfy the conditions of theorem 3.3 is guaranteed by lemma 3.1.

Remark 3.3**.**

For a given $\beta$, if there exists a ${\boldsymbol{v}}\in\mathbb{R}^{3}$ such that $\mathcal{K}_{\text{dir}}[{\boldsymbol{v}}t^{\beta}]$ is piecewise smooth then there also exists a vector ${\boldsymbol{w}}\in\mathbb{R}^{3}$ such that $\mathcal{K}_{\text{neu}}[{\boldsymbol{w}}t^{\beta-1}]$ is also a smooth function. However, the requirement that ${\boldsymbol{w}}t^{\beta-1}\in X$ implies that for $\mathcal{K}_{\text{neu}}$, only $\beta$’s which satisfy $\beta>1/2$ are admissible.

For $\mathcal{K}_{\text{dir}}$, note that $\beta_{0,j}=0$ for $j=0,1,2$ (see proof of lemma 3.1 contained in section A.1). These densities are essential for the proof of theorem 3.2, since these are the only basis functions for which the projection of their image under $\mathcal{K}_{\text{dir}}$ onto the constant functions are non-zero.

However, since $\beta_{0,j}\not>1/2$, the densities ${\boldsymbol{w}}_{0,j}t^{\beta_{0,j}-1}$ are excluded from the representation for the solution to the equation $\mathcal{K}_{\text{neu}}[\sigma]=\boldsymbol{h}$. Note that, unlike $\mathcal{K}_{\text{dir}}[{\boldsymbol{v}}_{i,j}t^{\beta_{i,j}}]$, $\mathcal{K}_{\text{neu}}[{\boldsymbol{w}}_{i,j}t^{\beta_{i,j}-1}]$, $i=1,2,\ldots$, $j=0,1,2$, have a non-zero projection onto the constants (see lemma B.2).

The following theorem is a converse of theorem 3.3 under suitable restrictions.

Theorem 3.4**.**

Consider the same geometry as in fig. 2, where $\theta_{1},\theta_{2}$, and $\theta_{3}$ are irrational numbers summing to $2.$ Let $\beta_{i,j},{\boldsymbol{w}}_{i,j},S_{i,j}$, $i=0,1,2,\ldots$, $j=0,1,2$ be as defined in lemma 3.1. Let $T_{i,j}$ denote the open subset of $(-1,1)^{2}$ on which $\beta_{i,j}$ and ${\boldsymbol{w}}_{i,j}$ are analytic and $\beta_{i,j}>1/2.$ For any positive integer $N$, let $S^{\text{neu}}_{N}$ denote the region of common analyticity of $\beta_{i,j},{\boldsymbol{w}}_{i,j}$, i.e., $S^{\text{neu}}_{N}=\cap_{i=1}^{N+1}\cap_{j=0}^{2}T_{i,j}$. Finally, suppose that $h_{(i,j)}^{k}$, $(i,j)\in\mathcal{J},k=0,1,2,\ldots N$ are real constants, and define $h_{(i,j)}$ by

[TABLE]

$0<t<$1.

Then there exists an open region $\tilde{S}^{\text{neu}}_{N}\subset S^{\text{neu}}_{N}\subset(-1,1)^{2}$ with $|\tilde{S}^{\text{neu}}_{N}|>0$ such that the following holds. For all $(a,b)\in\tilde{S}^{\text{neu}}_{N}$, there exist constants $p_{i,j}$, $i=1,2,\ldots N+1$, $j=0,1,2$, such that

[TABLE]

satisfies

[TABLE]

for $0<t<1$, where $C$ is a constant.

4 Conjectures

There are four independent parameters that completely describe the triple junction problem, any two out of the three angles $\{\theta_{1},\theta_{2},\theta_{3}\},$ and any two of the parameters $\{d_{(1,2)},d_{(2,3)},d_{(3,1)}\}=\{b,c,a\}.$ Let $Y\subset\mathbb{R}^{4}$, denote the subset of $\mathbb{R}^{4}$ associated with the four free parameters that completely describe any triple junction given by

[TABLE]

When $\theta_{1},\theta_{2}$, are irrational multiples of $\pi$, and $(a,b)$ are in the neighborhoods of $a=0$, $b=0$, and $c=0$, the results theorems 3.2 and 3.4 construct an explicit basis of non-smooth functions for the solutions of $\mathcal{K}_{\text{dir}}[\sigma]=\boldsymbol{h}$, and $\mathcal{K}_{\text{neu}}[\sigma]=\boldsymbol{h}$, and show that this basis maps onto the space of boundary data given by piecewise polynomials on each of the edges meeting at the triple junction. However, extensive numerical studies suggest that both of these results can be improved significantly. In particular, we believe that this analysis extends to all $(\theta_{1},\theta_{2},a,b)\in Y$, except for a set of measure zero. Moreover, on the measure zero set where this basis is not sufficient, we expect the solution to have additional logarithmic singularities; including functions of the form $t^{\beta}\log{(t)}{\boldsymbol{v}}$ should be sufficient to fix the deficiency of the basis. We expect the analysis to be similar in spirit to the analysis carried out for the solution of Dirichlet and Neumann problems for Laplace’s equations on vicinity of corners (see [25, 26]).

In this section, we present a few open questions for further extending the results theorems 3.2 and 3.4, and present numerical evidence to support these conjectures.

4.1 Existence of $\beta_{i,j}$

The solutions $\beta_{i,j}$, $i=0,1,2,\ldots$, $j=0,1,2$, are constructed as the implicit solutions of $\det{\mathcal{A}_{\text{dir}}(a,b,\beta)}=0$ (recall that $\det{\mathcal{A}_{\text{dir}}(a,b,\beta)}=\det{\mathcal{A}_{\text{neu}}(a,b,\beta)}$). Note that $\det{\mathcal{A}_{\text{dir}}(a,b,\beta)}=\sin{(\pi\beta)}\cdot\alpha(a,b,c;\beta)$ where $\alpha$ is as defined in eq. 52. From this, it follows that $\beta_{i,0}=i$ always satisfies $\det{\mathcal{A}_{\text{dir}}(a,b,\beta)}=0$ for all $\theta_{1},\theta_{2}$, and that $\beta_{0,j}=0$ results in three linearly independent basis functions of the form $t^{\beta}{\boldsymbol{v}}$ since $\mathcal{A}_{\text{dir}}(a,b,0)=0$.

The remaining $\beta_{i,j}$, $i=1,2,\ldots$, $j=1,2$, are constructed in the following manner. $\alpha(a,b,c;\beta)$ simplifies significantly along $a=0$, $b=0$, and $c=0$, and the existence of $\beta_{i,j}$ which satisfy $\det{\mathcal{A}_{\text{dir}}(a,b,\beta)}=0$ is guaranteed based on the explicit construction detailed in [28]. The construction then uses the implicit function theorem to extend the existence of $\beta_{i,j}$ to a subset of $(a,b)\in(-1,1)^{2}$. The implicit function theorem is a local result and only guarantees existence in local neighborhoods of the initial points. However, extensive numerical evidence suggests that the $\beta_{i,j}$ are well-defined and analytic for all $(a,b)\in(-1,1)^{2}$ and all $\theta_{1},\theta_{2}$. In fig. 3, we plot a few of these functions to illustrate this result.

Conjecture 4.1**.**

There exists a countable collection of $\beta_{i,j}$, $i=1,2,\ldots,$, $j=1,2$ which satisfy $\alpha(a,b,c;\beta_{i,j})=0$. Moreover, these $\beta_{i,j}$ are analytic functions of $\theta_{1},\theta_{2},a$, and $b$, for all $(\theta_{1},\theta_{2},a,b)\in Y$.

An alternate strategy for proving this result is by making the following observation. For fixed $\theta_{1},\theta_{2}$, consider the curve $\gamma_{m}:(m,m+1)\to\mathbb{R}^{3}$ defined by

[TABLE]

where $m$ is an integer. This defines a curve in $\mathbb{R}^{3}$ for which $|\gamma_{m}|\to\infty$ for each $m$. Then consider the family of hyperboloids parameterized by $(a,b)$ given by

[TABLE]

It follows immediately that the solutions to $\alpha(a,b,c;\beta)=0$ can be characterized geometrically as points in the intersection of the hyperboloid $H(x,y,z;a,b)$ with the curve $\gamma_{m}$.

4.2 Completeness of the singular basis

Having identified the $\beta_{i,j}$, and the corresponding null vectors ${\boldsymbol{v}}_{i,j}$ for $\mathcal{A}_{\text{dir}}$, and ${\boldsymbol{w}}_{i,j}$ for $\mathcal{A}_{\text{neu}}$, the second part of the proof shows that every set of boundary data which is a polynomial of degree less than or equal to $N$ on each of the edges, has a solution to the integral equations eqs. 14 and 15 in the ${\boldsymbol{v}}_{i,j}t^{\beta_{i,j}}$ basis for $\mathcal{K}_{\text{dir}}$ and ${\boldsymbol{w}}_{i,j}t^{\beta_{i,j}-1}$ for $\mathcal{K}_{\text{neu}}$ which agrees with the boundary data with error $O(t^{N+1})$.

This part of the proof relies on constructing an explicit mapping from the coefficients of the density $\sigma$ in the ${\boldsymbol{v}}_{i,j}t^{\beta_{i,j}}$ to the coefficients of Taylor expansions for $\mathcal{K}_{\text{dir}}[\sigma].$ Then, along $a=0$, $b=0$, or $c=0$, based on the results in [28], we show that this mapping is invertible along these edges. It then follows from the continuity of determinants that the mapping is invertible for open neighborhoods of the line segments $a=0$, $b=0$, $c=0$. This implies that in the basis ${\boldsymbol{v}}_{i,j}t^{\beta_{i,j}}$ there exists a $\sigma$ such that $|\mathcal{K}_{\text{dir}}[\sigma]-\boldsymbol{h}|\leq O(t^{N+1})$, for all boundary data $f$ in the space of polynomials with degree less than or equal to $N$.

While we prove this result for an open neighborhood $(a,b)$ of the line segments $a=0$, $b=0$, $c=0$, when the angles $\theta_{1},\theta_{2}$ are irrational multiples of $\pi$, we expect the bases to have this property for all $(\theta_{1},\theta_{2},a,b)\in Y$ except for a measure zero set. Moreover, this measure zero set is the set of $(\theta_{1},\theta_{2},a,b)$ for which the multiplicity of $\beta_{i,j}$ as a repeated root of $\det{\mathcal{A}_{\text{dir}}(a,b,\beta_{i,j})}=0$ is not the same as the dimension of the null space of $\mathcal{A}_{\text{dir}}(a,b,\beta_{i,j})$.

Conjecture 4.2**.**

Suppose that 4.1 holds, i.e. $\beta_{i,j}:Y\to\mathbb{R}$ are analytic functions. Suppose further that $h_{(i,j)}^{k}$, $(i,j)\in\mathcal{J},k=0,1,2,\ldots N$ are real constants, and suppose that

[TABLE]

$0<t<$1. Then there exists a measure zero set $S$ such that for all $(\theta_{1},\theta_{2},a,b)\in Y\setminus S$ the following result holds. There exist constants $p_{i,j}$, $i=0,1,\ldots N$, $j=0,1,2$, such that

[TABLE]

satisfies

[TABLE]

for $0<t<1$, where $C$ is a constant.

In fig. 3, we plot sections of the zero measure set on which 4.2 does not hold.

5 Discretization of eqs. 14 and 15

In this section we discuss a numerical method for solving eqs. 14 and 15 for the unknown densities $\sigma,\rho$ which exploits the analysis of their behavior in the vicinity of triple junctions. There are two general approaches for discretizing these integral equations: Galerkin methods, in which the densities $\rho$ and $\sigma$ is represented directly in terms of appropriate basis functions, and Nyström methods, where the solution is represented in terms of its values at specially chosen discretization nodes. In this paper, we use a Nyström discretization for solving eq. 14, though we note that the expansions in theorems 3.2 and 3.4 can also be used to construct efficient Galerkin discretizations.

In [23], the authors developed a Nyström discretization for resolving the singular behavior of solutions to integral equations in the vicinity of corners. In this approach, the authors obtain a basis of solutions to the integral equation in the vicinity of the corner by solving a small number of local problems. Based on these families of solutions, discretization nodes capable of interpolating the span of these solutions, coupled with quadratures for handling far-field interactions (inner products of the basis of solutions with smooth functions), and special quadratures for handling near interactions (for resolving the near singular behavior of the kernel in the vicinity of the corner) are developed. This approach was later specialized for the solution of Laplace’s equation on polygonal domains to obtain universal discretization nodes, and quadrature rules [24].

Recent advances in the analysis of integral equations for Laplace’s equation have provided analytic representations of solutions to integral equations in the vicinity of the corners [25, 26], obviating the need for obtaining the span of solutions in the vicinity of corners through numerical means. Based on the approach above, these analytical results have been exploited to construct universal discretization and quadratures for solutions in vicinity of corners [27]. We briefly discuss the construction of the Nyström discretization in [27] below. Let $\mathcal{F}$ denote the family of functions

[TABLE]

Then there exist $t_{j}\in[0,1]$, $w_{j}>0$, an orthogonal basis $\phi_{j}(t)$, $j=1,2,\ldots\boldsymbol{k}_{\textrm{AB}}=36$, and a $\boldsymbol{k}_{\textrm{AB}}\times\boldsymbol{k}_{\textrm{AB}}$ matrix $V$ whose condition number is $O(1)$, with the following features. For any $f\in\mathcal{F}$, there exists $c_{j}$, such that

[TABLE]

Let $f_{j}=f(t_{j})\sqrt{w_{j}}$ denote the samples of the function at the discretization nodes scaled by the square root of the quadrature weights. The matrix $V$ maps $f_{j}$ to its coefficients $c_{j}$ in the $\phi_{j}$ basis. Finally the weights $w_{j}$ are such that

[TABLE]

Specialized quadrature rules for handling the near-singular interaction between corner panels which meet at the same vertex are also constructed. The Dirichlet problem for Laplace’s equation can then be discretized using panels with scaled Gauss-Legendre nodes for panels which are away from corners, and using scaled nodes $t_{j}$ for panels at corners.

In the vicinity of triple junctions, the behavior of the solution $\sigma$ of eq. 14 can be represented to high-order as a linear combination of functions in $\mathcal{F}$. Thus the discretization for the Dirichlet problem discussed above can be used to obtain a Nyström discretization for eq. 14. Unfortunately, the same is not true when solving eq. 15, since the singular behavior of $\rho$ is not contained in the span of $\mathcal{F}$. In particular, the leading order singularity in $\rho$ is of the form $t^{\beta}$ where $\beta\in(-\frac{1}{2},0)$. The nature of the singularity of $\rho$ is similar to the singular behavior of solutions to integral equations corresponding to the Neumann problem on polygonal domains.

Recall that  eq. 15 is the adjoint of eq. 14. Thus, formally, one could use the transpose of the Nyström discretization of eq. 14 to solve eq. 15. Specifically, if $\overline{\rho}=\{\rho_{j}\}_{j=1}^{N}$ are the unknown values of $\rho$ at the discretization nodes, and $\overline{g}=\{g_{j}\}_{j=1}^{N}$ denote the samples of the boundary data for eq. 15 at the discretization nodes, then we solve the linear system

[TABLE]

where $M$ is the matrix corresponding to Nyström discretization of eq. 14. The solution $\overline{\rho}$ is a high-order accurate weak solution for the density $\rho$ which can be used to evaluate the solution to eq. 15 accurately away from the corner panels of the boundary $\Gamma$. This weak solution can be further refined to obtain accurate approximations of the potentials in the vicinity of corner panels through solving a sequence small linear systems for updating the solution $\rho_{j}$ in the vicinity of the corner panels. This procedure is discussed in detail in [32].

6 Numerical examples

We illustrate the performance of the algorithm with several numerical examples. In each of the problems let $\Omega_{0}$ denote the exterior domain and $\Omega_{i},i=1,2,\ldots N_{r}$ denote the interior regions. Let $c_{j,k}$, $k=1,2,\ldots 10$, denote points outside of the region $\Omega_{j}$ for $j=1,2,\ldots N_{r}$. The results in sections 6.2 and 6.3 have been computed using dense linear algebra routines, while the results in sections 6.1 and 6.4 have been computed using GMRES where the matrix vector product computation has been accelerated using fast multipole methods [33].

6.1 Accuracy

In order to demonstrate the accuracy of our method we solve the PDE with boundary data corresponding to known harmonic functions using our discretization of the integral equation formulation. We set $u_{j}(x)=\sum_{k=1}^{10}\log|x-c_{j,k}|$ and set $u\equiv 0$ for $x\in\Omega_{0}$. We then compute the boundary data

[TABLE]

and solve for $\sigma,\rho$. Given the discrete solution for $\sigma,\rho$, we compare the computed solution and plot the error in the computed at targets in the interior of each of the regions. In figs. 4 and 5, we demonstrate the results for two sample geometries.

Remark 6.1**.**

Note that we do not use special quadratures for handling near boundary targets which is responsible for the loss of accuracy close to the boundary. For panels away from the corner, the potential at near boundary targets can be computed accurately using several standard methods such as Quadrature by expansion, or product quadrature (see [29, 30, 31]). In order to evaluate the solution at points lying close to a corner panel a different approach is required. A detailed description of a computationally efficient algorithm for evaluating the solution accurately arbitrarily close to a corner is presented in [32].

6.2 Condition number dependence on $\mu,\nu$

In this section, we discuss the dependence of the condition number of the discretized linear systems as a function of the material parameters of the regions. Recall that the condition number of a linear system $A$, which we denote by $\kappa(A)$, is the ratio of the largest singular value $s_{\text{max}}$ to the the smallest singular value $s_{\text{min}}$, i.e. $\kappa(A)=s_{\text{max}}/s_{\text{min}}$. As discussed in section section 3, for fixed angles the integral equation and the analytical behavior of integral equations eqs. 14 and 15 are solely a function of $d_{(1,2)},d_{(2,3)},d_{(3,1)}$ defined in eq. 19. Furthermore, $d_{(1,2)}$ can be expressed in terms of $d_{(3,1)},d_{(2,3)}$ which are contained in the interval $(-1,1)$. As before, let $a=d_{(3,1)}$ and $b=d_{(2,3)}$. Since the discrete linear system corresponding to eq. 15, is the adjoint of the linear system corresponding to eq. 14, it suffices to study the condition number for either linear system.

In fig. 6, we plot the condition number of the discretization of eq. 14 as we vary $(a,b)\in(-1,1)^{2}$, by holding the values of $\mu$ in each of the regions to be fixed. In particular, we set $\mu_{1}=0.37$, $\mu_{2}=0.81$, $\mu_{3}=1$, and $\nu_{3}=0.77$. $\nu_{1}$, $\nu_{2}$ can then be defined in terms of $(a,b)$ as

[TABLE]

We note that the problem is well-behaved for almost all values of $(a,b)$ and becomes ill-conditioned as we approach the line $b=-1$ and $a=1$. This behavior is expected since the underlying physical problem also has rank-deficiency along these limits since these values of the parameters correspond to interior Neumann problems in regions $1$ and $2$ respectively.

6.3 Condition number dependence on angles at the triple junction

In this section we discuss the the dependence of the condition number of the discretized linear systems as a function of the angles at the triple junction. Let $\theta_{1},\theta_{2},\theta_{3}$, denote the angles at the triple junction, then $\theta_{1}+\theta_{2}+\theta_{3}=2\pi$. The three angles at any triple junction can be parameterized by $\theta_{1},\theta_{2}$ in the simplex $\{(\theta_{1},\theta_{2}):\theta_{1}>0,\theta_{2}>0\,,\theta_{1}+\theta_{2}<2\pi\}$. Suppose that we split this simplex into $4$ regions as shown in fig. 7. By symmetry it suffices to vary the angles $(\theta_{1},\theta_{2})\in(0,\pi)^{2}$.

The physical problem as either of the angles approach [math] or $2\pi$ becomes increasingly ill-conditioned due to close-to-touching interactions on the entire edge (not just near the corner). In order to avoid these issues and to automate geometry generation as we vary the angles $\theta_{1},\theta_{2}$, we use two different types of geometries for regions I and IV which are shown in fig. 7.

Resolving the close-to-touching interactions has numerical consequences as well; due to the increased number of quadrature nodes required as the angles tend to [math] in the universal quadrature rules. In order for the universal quadrature rules to remain efficient, they are generated for the range $(\theta_{1},\theta_{2})\in(\frac{\pi}{12},2\pi-\frac{\pi}{12})$. Regions with narrower angles should be handled on case to case basis and region with careful discretization of the boundary coupled with special purpose quadrature rules which account for the specific singular behavior of the solutions in the vicinity of triple junctions. In fig. 7, the top right missing corner corresponds to $\theta_{3}\in(0,\pi/12)$.

Referring to fig. 7, we observe that the condition number of the discrete linear systems varies mildly as we vary the angles $\theta_{1},\theta_{2}$, with a maximum condition number of $2.8$. The discontinuity in the plot is explained by the different choice of geometries for regions I, IV.

6.4 Application - Polarization computation

In this section, we demonstrate the efficiency of our approach for computing polarization tensors for a perturbed hexagonal lattice with cavities. The polarization computation corresponds to the following particular setup of the triple junction problem, $\mu_{i}=1$, $f_{i}=0$, $\nu_{i}=\varepsilon_{i}$, where $\varepsilon_{i}$ denotes the permittivity of the medium, and $g_{1}(\boldsymbol{x})=(\varepsilon_{\ell(i)}-\varepsilon_{r(i)})n_{1}(\boldsymbol{x})$ or $g_{2}(\boldsymbol{x})=(\varepsilon_{\ell(i)}-\varepsilon_{r(i)})n_{2}(\boldsymbol{x})$, where $\boldsymbol{x}=(x_{1},x_{2})\in\Gamma_{i}$, $\boldsymbol{n}(\boldsymbol{x})=(n_{1}(\boldsymbol{x}),n_{2}(\boldsymbol{x}))$, and $\varepsilon_{\ell(i)},\varepsilon_{r(i)}$ are the conductivities of the regions on either side of the edge $\Gamma_{i}$. If $u_{1}$ is the solution corresponding to $g_{1}$ and $u_{2}$ is the solution corresponding to $g_{2}$, then the polarization tensor $P$ is the $2\times 2$ matrix given by

[TABLE]

Note that in this particular setup, we only need to solve the problem corresponding to the operator $\mathcal{K}_{\text{neu}}$, as the solution $\sigma$ for $\mathcal{K}_{\text{dir}}[\sigma]=0$ is $\sigma=0$. Let $\rho_{1},\rho_{2}$ denote the solutions of eq. 15 corresponding to boundary data $g_{1}$ and $g_{2}$ respectively. Using properties of the single layer potential, the integrals of the polarization tensor can be expressed in terms of $\rho$ as

[TABLE]

We compare the efficiency of our approach to RCIP which to the best of our knowledge is the state of the art method for such problems. The geometry is generated using a regular hexagonal lattice inside the unit square whose vertices are perturbed in a random direction by a tenth of the side length, and the permittivity $\varepsilon$ is region $i$ is given by $10^{c_{i}}$ where $c_{i}$ is a uniform random number between $[-1,1]$. The choice of parameters for the problem setup is identical to the setup in section 11 in [34].

We discretize the geometry with $3$ panels on each edge of roughly equal size, and the reference solution is computed using $5$ panels on each edge. The geometry contains $10688$ vertices, $15855$ edges, and $5189$ regions. There are $1395240$ degrees of freedom for the coarse discretization (approximately $88$ degrees of freedom per edge) and $1902600$ degrees of freedom for the reference solution. These discretizations required $131$ iterations for GMRES to converge to a relative residual of $10^{-16}$, and the absolute error in the polarization tensor when compared to the reference solution is $5.1\times 10^{-12}$. In comparison, RCIP using approximately $71$ degrees of freedom per edge in a hexagonal lattice with $5293$ inclusions obtained an accuracy of $2\times 10^{-14}$ in computing the $2,2$ entry of the polarization matrix and required $105$ GMRES iterations to converge.

The polarization tensor for this configuration, correct to 13 significant digits, is given by

[TABLE]

Remark 6.2**.**

We note that the performance of our approach is close to current state of the art methods such as RCIP [22]. In our examples, further improvements in speed can be achieved using additional compression techniques to reduce the degrees of freedom in the resulting linear system [35, 36].

7 Concluding remarks and future work

In this paper we analyze the systems of boundary integral equations which arise when solving the Laplace transmission problem in composite media consisting of regions with polygonal boundaries. Our discussion is focused on the particular case of composite media with triple junctions (points at which three distinct media meet) though our analysis extends to higher-order junctions in a natural way.

We show that under some restrictions the solutions to the boundary integral equations corresponding to a triple junction is well-approximated by a linear combination of powers $t^{\beta_{j}}$ where $t$ denotes the distance from the corner along the edge, and the $\beta_{j},$ $j=1,2,\dots$ is a countable collection of real numbers obtained by solving a certain equation depending only on the material properties of the media and the angles at which the interfaces meet.

In addition to the theoretical interest of the result, our analysis also enables an easy construction of near-optimal discretizations for triple junctions. In particular, RCIP which is the leading method for solving electrostatic problems on multiple junction interfaces, requires approximately $71$ discretization nodes per edge to compute solutions to near machine precision accuracy, where as our proposed discretization achieves an accuracy of $5\times 10^{-12}$ using roughly $88$ discretization nodes per edge. Finally, we illustrate the properties of this discretization with a number of numerical examples.

The results of this paper admit a number of natural extensions and generalizations. Firstly, the analysis outlined in this paper extends almost immediately to junctions involving greater numbers of media. However, the construction of an efficient Nyström discretization of higher-order junctions requires special care since the solutions to corresponding integral equations are not $\mathbb{L}^{2}$ functions on the boundary. (in fact the solutions are known to be $\mathbb{L}^{1}$ functions on the boundary [21]). Secondly, with a small modification a similar analysis should be possible for boundary integral equations arising from triple junction problems for other partial differential equations such as the Helmholtz equation, Maxwell’s equations, and the biharmonic equation. This line of inquiry is being vigorously pursued and will be reported at a later date.

Finally, a similar approach will also work for generating discretizations of triple junctions in three dimensions. This is particularly valuable since geometric singularities in three-dimensions can often result in prohibitively large linear systems. Accurate discretization with few degrees of freedom would greatly improve the size and complexity of systems which could be simulated.

8 Acknowledgments

The authors would like to thank Alex Barnett, Charles Epstein, Leslie Greengard, Vladimir Rokhlin, and Kirill Serkh for many useful discussions.

Appendix A Analysis of $\mathcal{K}_{\text{dir}}$

First we present the proof of theorem 3.1. In order to do so, we require the following technical lemma which describes the double layer potential defined on a straight line segment with density $s^{\beta}$ at an arbitrary point near the boundary. Here $s$ is the distance along the segment.

Lemma A.1**.**

Suppose that $\Gamma$ is an edge of unit length oriented along an angle $\theta$, parameterized by $s(\cos{(\theta)},\sin{(\theta)})$, $0<s<1$. Suppose that $\boldsymbol{x}=t(\cos{(\theta+\theta_{0})},\sin{(\theta+\theta_{0})})$ (see fig. 9) where $0<t<1$, and $\boldsymbol{x}\not\in\Gamma$. Suppose that $\sigma(s)=s^{\beta}$ for $0<s<1$, where $\beta\geq 0$. If $\beta$ is not an integer, then

[TABLE]

If $\beta=m$ is an integer, then

[TABLE]

In the following lemma, we compute the potential $\mathcal{K}_{\text{dir}}[{\boldsymbol{v}}t^{\beta}]$, in the vicinity of triple junction with angles $\theta_{1},\theta_{2},\theta_{3}$, material parameters ${\boldsymbol{d}}=(d_{(1,2)},d_{(2,3)},d_{(3,1)})$, where ${\boldsymbol{v}}\in\mathbb{R}^{3}$ and $\beta$ are constants (see fig. 2).

Lemma A.2**.**

Consider the geometry setup of the single vertex problem presented in section 3. For a constant vector ${\boldsymbol{v}}\in\mathbb{R}^{3}$, suppose that the density on the edges is of the form

[TABLE]

If $\beta$ is not an integer, then

[TABLE]

where $\mathcal{A}_{\text{dir}}$ is defined in eq. 21 and

[TABLE]

If $\beta=m$ is an integer, then

[TABLE]

where

[TABLE]

Proof.

The result follows from repeated application of lemma A.1 for computing $\mathcal{D}_{(l,m):(i,j)}\sigma_{(i,j)}$. ∎

The proof of theorem 3.1 then follows immediately from lemma A.2.

We now turn our attention to the proof of lemma 3.1, which provides a construction of $\beta,{\boldsymbol{v}}$ satisfying the conditions of theorem 3.1. In order to do that, we first observe that if one of $a,b$, or $c$ is [math], then the expression of $\det{\mathcal{A}_{\text{dir}}}$ simplifies significantly, and there exists an explicit construction of $\beta$ satisfying $\det{\mathcal{A}_{\text{dir}}}(a,b,\beta)=0$. Recall that we use interchangeably use the following variables for the material properties $(a,b,c)=(d_{3,1},d_{1,2},d_{2,3})$. Having established the existence of analytic $\beta,{\boldsymbol{v}}$ on a 1-D manifold which is a subset $(a,b)\in(-1,1)^{2}$, we now analytically continue these values of $\beta,{\boldsymbol{v}}$ to carve out the open region $S$ on which $\beta,{\boldsymbol{v}}$ can be analytically extended. This proof is discussed in section A.1.

A.1 Existence of $\beta,{\boldsymbol{v}}$ satisfying theorem 3.1

The determinant of the matrix $\mathcal{A}_{\text{dir}}(a,b,\beta)$ is given by

[TABLE]

where $c=-(a+b)/(1+ab)$, and

[TABLE]

Given the formula above, for all $(a,b)\in(-1,1)^{2}$ when $\beta=m\geq 0$ is an integer, $\det{\mathcal{A}_{\text{dir}}(a,b,\beta)}=0$. When $m\neq 0$, the matrix $\mathcal{A}_{\text{dir}}$ has rank-2, since the matrix is similar to an anti-symmetric matrix and is not identically zero. The null vector ${\boldsymbol{v}}$ of $\mathcal{A}_{\text{dir}}(a,b,m)$ is given by ${\boldsymbol{v}}_{m}=[\sin{(m\theta_{3})},\sin{(m\theta_{1})},\sin{(m\theta_{2})}]^{T}$, i.e., the pair $(m,{\boldsymbol{v}}_{m})$ always satisfies eq. 21. When $\beta=0$, $\mathcal{A}_{\text{dir}}(a,b,\beta)=0$ and hence for any ${\boldsymbol{v}}\in\mathbb{R}^{3}$, the pair $\beta,{\boldsymbol{v}}$ satisfies eq. 21. Based on this observation we set

[TABLE]

We now turn our attention to constructing the remaining $\beta_{i,j}$, the corresponding vectors ${\boldsymbol{v}}_{i,j}$, and their regions of analyticity $S_{i,j}$, $i=1,2,\ldots$, $j=1,2$. From eq. 51, the remaining values of $\beta_{i,j}$ as a function of the material parameters $(a,b)$ are defined implicitly via the roots of the equation $\alpha(a,b,c(a,b);\beta_{i,j}(a,b))=0$, where $c=-(a+b)/(1+ab)$ and $\alpha$ is defined in eq. 52.

It turns out that the implicit solutions $\beta(a,b)$ of $\alpha(a,b,c(a,b);\beta(a,b))=0$, are known when $a=0$, $b=0$, or $c=0$. This gives us an initial value for defining $\beta_{i,j}$ in order to apply the implicit function theorem, and extend it to a region containing the segments $a=0$, $b=0$, or $c=0$. Given this strategy, let $R_{1},\dots,R_{6}\subset(-1,1)\times(-1,1)$ be defined as follows (see fig. 10)

[TABLE]

In the following, we will consider only the segment $R_{1}$; and construct an open region $S_{i,j}^{1}\subset(-1,1)^{2}$ which contains $R_{1}$ on which we define a family of functions $\beta_{i,j}(a,b):S^{1}_{i,j}\to\mathbb{R}$, $j=1,2$, which satisfy the conditions of lemma 3.1. Analogous results hold for the open sets containing the remaining segments $R_{2},R_{3},\ldots R_{6}$ with almost identical proofs. The region of analyticity for $\beta_{i,j}$ is then given by $S_{i,j}=\cup_{k=1}^{6}S_{i,j}^{k}$.

Definition A.1**.**

For $(a,0)\in R_{1}$ and $i=1,2\dots$ let $\beta_{i,1}(a,0)$ be the solution to the equation

[TABLE]

such that

[TABLE]

Similarly, for $i=1,2,\dots$ let $\beta_{i,2}(a,0)$ be the solution to the equation

[TABLE]

such that

[TABLE]

The existence of $\beta_{i,j}$ for $i=1,2,\dots$ and $j=1,2$ satisfying these conditions is guaranteed by the following  lemma A.3, proved in [28].

Lemma A.3**.**

Suppose that $\delta\in\mathbb{R}$, $0<|\delta|<1$ and $\theta\in(0,2\pi)$ and $\theta/\pi$ is irrational. Consider the equations

[TABLE]

Then there exist an countable collection of functions $z_{i}^{\pm}(\delta),$ $i=1,2,\dots$ such that

$\sin^{2}(\pi z_{i}^{\pm}(\delta))=\delta^{2}\sin^{2}(z_{i}^{\pm}(\delta)(\pi-\theta))$* for all $\delta\in[0,1],$ and $i=1,2,\dots$* 2. 2.

the functions $z_{i}^{\pm}$ are analytic in $(0,1)$, 3. 3.

$\lim_{\delta\to 0}z_{i}^{\pm}(\delta)=i,$** 4. 4.

$z_{i}^{+}(\delta)>i$* and $z_{i}^{-}(\delta)<i$ for all $\delta\in(0,1)\,.$*

The following lemma extends the domain of definition of the functions $\beta_{i,j}$, $j=1,2$, to some open subset $S_{i,j}^{1}$ containing $R_{1}$.

Lemma A.4**.**

Suppose $\theta_{1},$ $\theta_{2}$ and $\theta_{3}$ are positive numbers summing to $2\pi$, and $\theta_{1}/\pi$, $\theta_{2}/\pi$, and $\theta_{3}/\pi$ are irrational numbers. Suppose that $\beta_{i,j}$ are defined as above for $i=1,2,\dots$ and $j=1,2$. For $a\in(0,1)$, the function $\beta_{i,j}$ satisfies

[TABLE]

Moreover, there exists a unique extension of $\beta_{i,j}$ to an analytic function of $(a,b)$ on an open neighborhood $R_{1}\subset S_{i,j}^{1}\subset(-1,1)^{2}$ which satisfies

[TABLE]

Proof.

We begin by observing that for $j=1,2,$ $\beta_{i,j}$ satisfies

[TABLE]

[TABLE]

Upon multiplication by

[TABLE]

and using eq. 66 we get

[TABLE]

which does not vanish for all $a>0$. Thus by the implicit function theorem, there exists an analytic extension of $\beta_{i,j}$ to a neighborhood $(a,b)\in R_{1}\subset S_{i,j}^{1}\subset(-1,1)^{2}$ which satisfies $\alpha(a,b,c(a,b);\beta_{i,j})=0$. ∎

The following theorem establishes the analyticity of the null vectors of $\mathcal{A}_{\text{dir}}(a,b;\beta)$ in a neighborhood of $R_{1}$ when $\beta=\beta_{i,j}.$

Theorem A.1**.**

For each $j=1,2$, and $i=1,2,\dots$, the matrix $\mathcal{A}_{\text{dir}}(a,b,\beta_{i,j})$ defined in eq. 21 has a null-vector ${\boldsymbol{v}}_{i,j}$ whose entries are analytic functions of $(a,b)$ on $S_{i,j}^{1}$.

Proof.

Since $\beta_{i,j}$ is such that the matrix $\mathcal{A}_{\text{dir}}(a,b,\beta_{i,j})$ is singular, it has a null vector ${\boldsymbol{v}}_{i,j}$. Moreover, as long as $(a,b)\neq(0,0)$ and $\beta_{i,j}$ is not an integer, the matrix $\mathcal{A}_{\text{dir}}$ has rank at least $2$. Thus [math] is an eigenvalue of $\mathcal{A}_{\text{dir}}(a,b,\beta_{i,j}(a,b))$ with multiplicity $1$ for all $(a,b)\in S_{i,j}^{1}$. Since the entries of the matrix $\mathcal{A}_{\text{dir}}$ are analytic functions of $(a,b)$, we conclude that the entries of ${\boldsymbol{v}}_{i,j}$ are analytic on $S_{i}^{1}$. ∎

Finally, each $S_{i,j}^{k}$ is an open subset containing the segments $R^{k}$, $k=1,2,\ldots 6$. Then $S_{i,j}=\cup_{k=1}^{6}S_{i,j}^{k}$ is an open subset of $(-1,1)^{2}$ containing $\cup_{k=1}^{6}R_{k}$. Thus, for any finite $N$, $|\cap_{i=0}^{N}\cap_{j=0}^{2}S_{i,j}|>0$.

A.2 Completeness of density basis

Recall that for any $\beta,{\boldsymbol{v}}$ which satisfy the conditions of theorem 3.1, and $\sigma={\boldsymbol{v}}t^{\beta}$, the potential $\mathcal{K}_{\text{dir}}[\sigma]$ corresponding to any of these densities is an analytic function. In order to show that, the potential corresponding to a particular collection of $\beta,{\boldsymbol{v}}$ span all polynomials of a fixed degree on all the three edges meeting at the triple junction, we explicitly write down the linear map from the coefficients of the density in the ${\boldsymbol{v}}t^{\beta}$ basis to the coefficients of Taylor series of the potentials on each of the edges using lemma A.2. We then observe that this mapping is invertible along the line segments corresponding to $a=0$, $b=0$, or $c=0$, and since the mapping is an analytic function of the parameters $(a,b)$, it must also be invertible in an open region containing the segments $a=0$, $b=0$ or $c=0$. This part of the proof is discussed in section A.2.

For any integer $N>0$, let $\hat{S}_{N}$, denote the common region of analyticity of $\beta_{i,j},{\boldsymbol{v}}_{i,j}$, $j=0,1,2$, $i=0,1,2\ldots N$, i.e. $S_{N}=\cup_{k=1}^{6}S^{k}_{N}$, where $S_{N}^{k}=\cap_{i=0}^{N}\cap_{j=0}^{2}S_{i,j}^{k}$. By construction, $R_{j}\subset S^{j}_{N}$ for all $N$. We now prove the result theorem 3.2 in one of the components of $S_{N}$, say $S_{N}^{1}$. The proof for the other components follows in a similar manner.

Let ${\boldsymbol{p}}_{i}=[p_{i,0},p_{i,1},p_{i,2}]^{T}$, and suppose that

[TABLE]

Then, using lemma A.2, since $\beta_{i,j},{\boldsymbol{v}}_{i,j}$ are such that $\mathcal{A}_{\text{dir}}(a,b,\beta_{i,j})\cdot{\boldsymbol{v}}_{i,j}=0$, the potential corresponding to this density on the boundary $(\Gamma_{(1,2)},\Gamma_{(2,3)},\Gamma_{(3,1)})$ is given by

[TABLE]

where $B_{i,j}$ are the $3\times 3$ matrices given by

[TABLE]

Let ${\bf B}$ denote the $3(N+1)\times 3(N+1)$ matrix whose $3\times 3$ blocks are given by $B_{i,j}$, $i,j=0,1,2,\ldots N$.

Recall that on $R_{1}\subset S_{N}^{1}$, $b=0$, $\beta_{i,0}=i$, $\beta_{i,1}$, satisfies $\sin{(\pi\beta_{i,1})}=-a\sin{(\beta_{i,1}(\pi-\theta_{3}))}$, $\beta_{i,2}$ satisfies $\sin{(\pi\beta_{i,2})}=a\sin{(\beta_{i,2}(\pi-\theta_{3}))}$ , $i=0,1,2\ldots$, and the corresponding vectors ${\boldsymbol{v}}_{i,j}$, $i=0,1,2\ldots$, $j=0,1,2$, are given by

[TABLE]

where

[TABLE]

Furthermore, the matrices ${\bf C}$ and $\boldsymbol{C}_{\text{diag}}$ defined in eq. 48, and eq. 50 respectively, also simplify to

[TABLE]

and

[TABLE]

Let $u_{(1,2),i},u_{(2,3),i},u_{(3,1),i}$ denote the coefficient of $|t|^{i}$ in the Taylor expansions of $u_{(1,2)},u_{(2,3)},u_{(3,1)}$ respectively. Let $P$ denote the permutation matrix whose action is given by

[TABLE]

Then along $R_{1}$, the matrix $PBP^{T}$ is demonstrated in fig. 11.

The matrices $D_{1},D_{2}$ are diagonal and are given by

[TABLE]

The matrices $C_{1},C_{2}$ are Cauchy matrices whose entries are given by

[TABLE]

Since we have assumed $\theta_{1}/\pi,\theta_{2}/\pi,\theta_{3}/\pi$, to be irrational, we note that $\eta_{i}>0$ and that $\sin{(m\theta_{3})}\neq 0$ for all $m\neq 0$. Thus, the diagonal matrices $D_{1},D_{2}$ are invertible. Furthermore on $(a,0)$, neither of $\beta_{i,1}$ or $\beta_{i,2}$, take on integer values lemma A.3. Thus, the Cauchy matrices $C_{1},C_{2}$ are invertible.

Let $T$ denote the bottom-right $2(N+1)\times 2(N+1)$ block. Then from the structure of $P{\bf B}P^{T}$ and the fact that the diagonal matrix $D_{1}D_{2}$ is invertible, it is clear that $B$ is invertible if and only if $T$ is invertible.

Remark A.1**.**

The matrix $T$ is the mapping from the coefficients of the singular basis of solutions for the transmission problem with angle $\pi\theta_{3}$ and material parameter $a$ to the corresponding coefficients of the Taylor expansion of the potential on the edges $(2,3),(3,1)$. The invertibility of $T$ follows from the analysis in [28]. We present the proof here in terms of the notation used in this paper.

Upon applying an appropriate permutation matrix $P_{2}$ to $T$ from the right and the left, we note that

[TABLE]

The matrix $P_{2}TP_{2}^{T}$ is invertible if and only if it’s bottom right $2N\times 2N$ corner is invertible. Let $I_{N}$ denote the $N\times N$ identity matrix, then the bottom right corner of $P_{2}TP_{2}^{T}$ factorizes as as

[TABLE]

which is clearly invertible since the matrices $D_{1},C_{1},C_{2}$ are invertible.

Finally, using all of these results, it follows that the matrix ${\bf B}$ is invertible for all $(a,0)=R_{1}$. Since all of the quantities involved are analytic, on every compact subset of $S^{1}_{N}$, we conclude that the matrix ${\bf B}$ is invertible in an open neighborhood $R_{1}\subset\tilde{S}^{1}_{N}\subset S^{1}_{N}$. By construction $|\tilde{S}^{1}_{N}|>0$.

Appendix B Analysis of $\mathcal{K}_{\text{neu}}$

All the proofs for the analysis of $\mathcal{K}_{\text{neu}}$ are similar to the corresponding proofs of $\mathcal{K}_{\text{dir}}$. We only present the analogs of lemmas A.1 and A.2.

In the following lemma we present the directional derivative of a single layer potential defined on straight line segment with density $s^{\beta}$ at an arbitrary point near the boundary. Here $s$ is the distance along the segment, at an arbitrary point near the boundary.

Lemma B.1**.**

Suppose that $\Gamma$ is an edge of unit length oriented along an angle $\pi\theta$, parameterized by $s(\cos{(\theta)},\sin{(\theta)})$, $0<s<1$. Suppose that $\boldsymbol{x}=t(\cos{(\theta+\theta_{0})},\sin{(\theta+\theta_{0})}$ and $\boldsymbol{n}=(-\sin{(\theta+\theta_{0})},\cos{(\theta+\theta_{0})})$ (see fig. 9) where $0<t<1$, and $\boldsymbol{x}\not\in\Gamma$. Suppose that $\sigma(s)=s^{\beta-1}$ for $0<s<1$, where $\beta\geq 1/2$. If $\beta$ is not an integer, then

[TABLE]

If $\beta=m$ is an integer, then

[TABLE]

In the following lemma, we compute the potential at a triple junction with angles $\pi\theta_{1},\pi\theta_{2},\pi\theta_{3}$, and material parameters ${\boldsymbol{d}}=(d_{(1,2)},d_{(2,3)},d_{(3,1)})$ (see fig. 2).

Lemma B.2**.**

Consider the geometry setup of the single vertex problem presented in section 3. For a constant vector ${\boldsymbol{v}}\in\mathbb{R}^{3}$, suppose that the density on the edges is of the form

[TABLE]

If $\beta$ is not an integer, then

[TABLE]

where $\mathcal{A}_{\text{neu}}$ is defined in eq. 25, and ${\bf C}({\boldsymbol{d}},k)$ is defined in eq. 48. If $\beta=m$ is an integer, then

[TABLE]

where $\boldsymbol{C}_{\text{diag}}$ is defined in eq. 50.

Proof.

The result follows from repeated application of lemma B.1 for computing $\mathcal{D}^{\ast}_{(l,m):(i,j)}\sigma_{(i,j)}$. ∎

The proof of theorem 3.3 then follows immediately from lemma B.2.

In the following lemma, we prove that $\beta_{i,j},S_{i,j}$, $i=1,2,\ldots$, $j=0,1,2$, defined in section A.1 satisfy $\beta_{i,j}(a,b)>1/2$ for all $(a,b)$ in an open subset $T_{i,j}\subset S_{i,j}$.

Lemma B.3**.**

Suppose that $\beta_{i,j},S_{i,j}$, $i=1,2,\ldots$, $j=0,1,2$, are as defined in section A.1. Then there exists an open subset $T_{i,j}\subset S_{i,j}$, such that $\beta_{i,j}(a,b)>1/2$ for all $(a,b)\in T_{i,j}$. Moreover for any $N>0$, $\cap_{i=1}^{N+1}\cap_{j=0}^{2}|T_{i,j}|>0$.

Proof.

Since $\beta_{i,0}=i$, the statement is trivially true with $T_{i,j}=(-1,1)^{2}$. Since $\beta_{i,j}=z_{i}^{\pm}(\delta,\theta)$ on $a=0,b=0$, or $c=0$, for appropriate parameters $\delta,\theta$, we conclude that $\beta_{i,j}>1/2$, on $a=0,b=0$, or $c=0$, for $i=1,2,\ldots$, $j=1,2$. Since $\beta_{i,j}$ are analytic on $S_{i,j}$, there exists an open subset containing the segments $a=0,b=0$, or $c=0$, which we denote by $T_{i,j}$, such that $\beta_{i,j}(a,b)>1/2$ for all $(a,b)\in T_{i,j}$. Since each $T_{i,j}$ is an open subset of $(-1,1)^{2}$, containing $\cup_{k=1}^{6}R_{k}$, we conclude that $|\cap_{i=1}^{N+1}\cap_{j=0}^{2}T_{i,j}|>0$.

∎