Stress fields around two pores in an elastic body: exact quadrature domain solutions

Analytical solutions are given for the stress fields, in both compression and far-field shear, in a two-dimensional elastic body containing two interacting non-circular pores. The two complex potentials governing the solutions are found by using a conformal mapping from a pre-image annulus with those potentials expressed in terms of the Schottky–Klein prime function for the annulus. Solutions for a three-parameter family of elastic bodies with two equal symmetric pores are presented and the compressibility of a special family of pore pairs is studied in detail. The methodology extends to two unequal pores. The importance for boundary value problems of plane elasticity of a special class of planar domains known as quadrature domains is also elucidated. This observation provides the route to generalization of the mathematical approach here to finding analytical solutions for the stress fields in bodies containing any finite number of pores.


Introduction
There is a rich tradition in the field of linear plane elasticity of making use of analytic function theory, coupled with conformal mapping techniques, to find the stress distributions in two-dimensional elastic a role [20]. The present work appears to be the first to point out, and showcase, the role of quadrature domains in problems of plane elasticity. Importantly, our methodology also extends to finding the stress distributions around any finite number of pores in an unbounded medium. This matter is discussed in more detail in §7.

Function theory in an annulus
To solve the linear elasticity problems of interest here, we need to introduce some convenient functions in the annulus ρ < |ζ | < 1. All those needed can be constructed from the basic function (2.1) Standard techniques [21] can be used to confirm that this infinite product is convergent for all finite ζ = 0 and 0 < ρ < 1. It is easy to verify directly from the definition (2.1)-see appendix A for details-that P(ρ 2 ζ ) = −ζ −1 P(ζ ), P(ζ −1 ) = −ζ −1 P(ζ ). (2.2) An important feature of the function P(ζ ) is that it vanishes when ζ = 1, and also at all points ζ = ρ 2n , where n is any non-zero integer. We also introduce which is the logarithmic derivative of P(ζ ) multiplied by ζ , and the prime notation denotes differentiation with respect to ζ . It has a simple pole at ζ = 1 with residue +1, that is, near ζ = 1, + locally analytic function. L(ζ ) has a second-order pole at ζ = 1 with strength −1, that is, near ζ = 1, L(ζ ) = − 1 (ζ − 1) 2 + locally analytic function.
(2.8) L(ζ ) is an example of a loxodromic function [22]: a function H(ζ ) is defined to be a loxodromic function if it is meromorphic everywhere inside (and on the boundary of) some fundamental annulus ρ ≤ |ζ | < ρ −1 and satisfies the functional relation H(ρ 2 ζ ) = H(ζ ). (2.9) The annulus ρ ≤ |ζ | < ρ −1 is called fundamental because, given the singularity structure of H(ζ ) in this annulus, its singularities in all other annuli filling out the complex ζ -plane follow from (2.9). For L(ζ ), because its only singularity in the fundamental annulus (or, strictly speaking, on its boundary) is a second-order pole at ζ = 1, and because it satisfies (2.9), then L(ζ ) qualifies as a loxodromic function.  Locations of the three poles (crosses) and three zeros (small circles) for the conformal mapping (3.1). In the annulus ρ < |ζ | < 1, which is the pre-image of the elastic body, the mapping has a pole at √ ρ and a zero at − √ ρ. It also has two simple poles, at √ ρ −1 e ±iθ , in the unphysical annulus 1 < |ζ | < ρ −1 as well as two simple zeros at ± √ ρ −1 .
The function K(ζ ) in (2.3), although not itself loxodromic, will be our main tool later in constructing loxodromic functions having simple pole singularities.

Conformal mapping from an annulus
Consider the special class of unbounded elastic bodies exterior to two finite area pores described by conformal maps from the annulus ρ < |ζ | < 1 having the form where R, θ and ρ are real constants. (There is a minor abuse of notation here in that z is taken to denote both the complex coordinate in the physical plane and the conformal mapping function, but this should not cause confusion.) From its form as a ratio of products of P functions, the poles and zeros of the function (3.1) are clear: because P(ζ ) vanishes simply when ζ = 1 then, restricting to the poles and zeros in the annulus ρ < |ζ | < 1/ρ, we see that there are three simple zeros at − √ ρ, ± √ ρ −1 and three simple poles at √ ρ, √ ρ −1 e ±iθ . Figure 2 shows a schematic illustrating the pole and zero locations of (3.1) in the annulus ρ < |ζ | < ρ −1 . Under the mapping (3.1), the circle |ζ | = ρ (called C 1 ) and the circle |ζ | = 1 (called C 0 ) are transplanted to the pore boundaries, the point ζ = √ ρ maps to z = ∞ and ζ = − √ ρ maps to the origin z = 0. Not all choices of parameters are physically admissible though: we must restrict to parameter choices for which the image of the annulus ρ < |ζ | < 1 under the mapping is one-toone. A necessary (but not sufficient) condition is that z (ζ ) = 0 for ρ < |ζ | < 1. For all choices of admissible parameters, it can be verified that the images of the circles |ζ | = ρ, 1 are rotations of each other by angle π about the origin. To see this, note that properties (2.1) of P(ζ ) can be used to show that, for all ζ = 0, (a) ( b) Figure 3. (a,b) Typical two-pore shapes described by the mapping (3.1) for θ = π/3 (a) and θ = 2π/3 (b) and ρ = 0.005, 0.05, 0.1 and 0.2 (superposed) with centroids fixed. The elastic body is the unbounded doubly connected region exterior to the pores. For small area and θ = 2π/3, the pores are near ellipses with semi-axis ratio 1 : 3.
It is easy to check that, owing to the reflectional symmetry of its pole and zero locations about the real ζ -axis, the conformal mapping (3.1) satisfies which is consistent with the reflectional symmetry of the image domain about the real z-axis. Of crucial importance for our development is the observation that, on use of (2.2), it can be checked that (3.1) satisfies the functional relation for all choices of the parameters R, θ and ρ. Because (3.1) can also be seen to be meromorphic for all ζ = 0 then z(ζ ) is a loxodromic function. It is this feature of the conformal mappings that causes the domains whose images they represent to be called quadrature domains (see §7 for more discussion). Given that z(ζ ) satisfies (3.5), it should be clear that, if z(ζ ) is meromorphic, it is enough to study its poles and zeros in the annulus ρ < |ζ | < ρ −1 because its poles and zeros in all other annuli will then follow from (3.5). In this way, the poles and zeros depicted in figure 2 are the only ones with which we need be concerned (even though, in fact, the function z(ζ ) has a countable infinity of zeros and poles). The analysis to follow pertains for any choices of the three real parameters ρ, θ and R. However, in the example calculations given later, we focus on the particular choice This restriction reduces the pore shapes studied to a two-parameter family. Even more, we will normalize the domains to be such that the centroid of the two pores are at ±1 on the real axis, thereby further reducing the class of shapes to a one-parameter family. This condition on the centroids fixes R as a function of ρ and θ. Mathematically, ρ then serves as the natural governing parameter; altering it corresponds physically to varying the area of the two symmetric pores. As ρ → 0, the pore areas tend to zero; as ρ increases, the area of each grows and the separation between closest points on the pores decreases.
To give an idea of the broad class of two-pore shapes described by formula (3.1), figure 3 shows a superposition of possible pore shapes for the values ρ = 0.005, 0.05, 0.1 and 0.2 and for the two choices θ = π/3 and θ = 2π/3. For the latter choice of θ the transition, as ρ increases from zero, from small elliptical pores (having semi-axes in the ratio 1 : 3) to near-touching pores exhibiting several high curvature boundary portions is clear. Figure 4 shows the pore configuration for ρ = 0.5 and θ = 2π/3. An interesting feature is that the two 'outer' boundaries of the two pores, when considered together, appear to make up an almost closed circumscribing curve that is very   . Pore shapes with ρ = 0.5 and θ = 2π/3. The pores are close to touching at three near cusps. The two 'outer' pore boundaries together appear to form a single circumscribing circular boundary (but with two small gaps). close to circular. As ρ increases further, the three near-cusps on each pore boundary draw closer together until they almost touch.
It should be emphasized that, while we focus here on the special choice (3.6), generalizing the following analysis to any other value of θ is completely straightforward.
The special family of shapes we have just chosen affords us an opportunity to study stress interaction effects between two highly non-circular pores. Remarkably, those stress fields can be found in explicit analytical form for the two physically interesting scenarios examined in § §5 and 6.

Complex potentials in plane elasticity
The stresses and displacements under plane strain or plane stress conditions of a two-dimensional material can be expressed in terms of two complex potentials φ(z) and ψ(z). The displacements are given by where the displacement components are denoted by (u, v), G is the shear modulus of the host material, and the so-called Kolosov constant [1] κ = 3 − 4ν (for plane strain) and κ = (3 − ν)/(1 + ν) (for plane stress), with ν being the Poisson ratio. The potentials φ(z) and ψ(z) can be found by solving boundary value problems of different type involving the surface tractions. We now consider two such problems separately.

Two uniformly pressurized pores
Consider an infinite elastic body containing two pores, with no stresses acting at infinity, and with a uniform hydrostatic pressure of magnitude p = 1 acting along the pore boundaries. We will find the stress distribution around the pores subject to this uniform hydrostatic loading. Our results are a generalization, to two non-circular pores, of the work of Ekneligoda & Zimmerman [10] who studied single isolated pores. In [10], it is shown that for problems in which the traction is specified along the contour then For us, there is an important difference to the treatment in [10] arising from the fact that the elastic medium we consider is doubly connected; the analysis of multiply connected geometries is discussed in detail in the standard texts [1,2]. Owing to some additive freedoms in the specification of the potentials φ(z) and ψ(z), a constant of integration that would ordinarily appear in (5.2) has there been set to zero with impunity. However, for problems with multiple boundaries, only one such constant can be taken to vanish leaving others to be determined.
For the problem of two symmetric pores considered here, it turns out that the relevant conditions on the pore boundaries can be taken to be where γ is some real constant to be found and, clearly, neither of the integration constants have been set equal to zero here (unless it turns out that γ = 0-the determination of γ is discussed below). This is because the symmetry of the geometrical configuration with respect to rotation by π around the origin can be used to argue that φ(z) and ψ(z) can be chosen to be odd functions: which fixes all additive degrees of freedom in φ(z) and ψ(z) and necessitates the appearance of ±γ in (5.3). The fact that the two constants are negatives of each other follows from (5.4) and the symmetry of the pores. We now introduce the functions Because it is a one-to-one conformal mapping, z(ζ ) is analytic everywhere in the annulus ρ < |ζ | < 1 except for the required simple pole at ζ = √ ρ which maps to z = ∞. But because φ(z) and ψ(z) are analytic in the elastic medium, and decay as z → ∞, the composed functions Φ(ζ ) and Ψ (ζ ) turn out to be analytic everywhere in the annulus ρ < |ζ | < 1.
In terms of the functions (5.6), the complex conjugate of the first boundary condition in (5.3) takes the formΦ where we have used the fact thatζ = 1/ζ on C 0 . The complex conjugate of the second boundary condition in (5.3) becomes where we have used the fact thatζ = ρ 2 /ζ on C 1 . On use of property (3.5) of the conformal mapping function, (5.8) becomes Subtraction of (5.7) and (5.9) now leads tō This is a key observation that underlies the success of our approach. It should be clear from the sequence of steps above that (5.10) is a consequence of the condition (3.5) pertaining to the class of quadrature domains to which we have restricted attention.
As a point of interest, meromorphic functions satisfying the second of the functional relations in (5.10)-which reduces to (2.9) when γ = 0-are sometimes called quasi-loxodromic functions.
In view of (5.10), the analyticity properties of Φ(ζ ) throughout the complex plane will be determined by its analyticity properties in the annulus ρ < |ζ | < 1/ρ. From the complex conjugate of (5.7), we find But, Φ(ζ ) and Ψ (ζ ) are both known to be analytic in ρ < |ζ | < 1, so that are analytic in the annulus 1 < |ζ | < 1/ρ (recall also that, for a one-to-one conformal mapping, z (ζ ) = 0 for ρ < |ζ | < 1). Hence, in the annulus 1 < |ζ | < 1/ρ, (5.11) implies that the only possible singularities of Φ(ζ ) are those inherited from z(ζ ). However, (3.1) reveals that z(ζ ) has just two simple pole singularities there at the points e ±iθ / √ ρ. It is then natural to propose that A and B can be found by equating the residues of the simple poles at e ±iθ / √ ρ on each side of equation (5.11). Indeed, it is found that where a is such that near ζ = e −iθ / √ ρ, where we have used the fact that The constant C follows from ensuring that φ(0) = 0 as required in (5.5) is satisfied. Because ζ = − √ ρ is the pre-image of z = 0, then In summary, we have found that the complex potentials for two pores under hydrostatic loading are where A is given in (5.15).
With the stress fields determined in this way, it is possible to compute the pore compressibility C pc [10] defined by and is the initial area of pore 1 in figure 1 (taken to be the image of |ζ | = 1). These quantities were found in [10] for the class of isolated pores studied there. For plane strain, it is preferable to consider the non-dimensionalized quantity GC pc 1 − ν (5.25) because this turns out to be purely a function of the pore geometry. Zimmerman [9] showed this in the single pore case using algebraic manipulations based on a conformal mapping representation, but it pertains to the multiple pore case too as we show in appendix C using more general arguments. Appendix C gives details of how, on substituting for Ψ (ζ ) from (5.20) and use of the fact that κ + 1 = 4(1 − ν) for plane strain, to establish that implying that the quantity of interest (5.25) is The right-hand side of (5.27) is purely a function of the pore geometry and is independent of the Poisson ratio.
To compute the values of the right-hand side of (5.27), whose integrand is known explicitly in terms of the special functions P(ζ ), K(ζ ) and L(ζ ), we can make use of convenient infinite sum representations of K(ζ ) and L(ζ ) are given in (A 9) of appendix A. Truncated versions of these can be used for numerical computation. P(ζ ) can similarly be evaluated easily by truncating the infinite product (2.1) that defines it. The explicit integral in (5.27) can be evaluated, to exponential accuracy, by the trapezoidal rule.   .27), against pore area A for θ = 2π/3. Here, ρ is used as a parameter: for each ρ, the value of R is chosen that fixes the centroids of the pores to be at ±1. The centroid position C of pore 1 is readily computed from the formula where we have used (3.4). The parameter R can be chosen, as an explicit function of ρ and θ, so that C = 1 for pore 1. The conformal map (3.1) is then fully determined and the pore area A can be found from (5.24). As ρ → 0, so that the pore area tends to zero, the shapes of the pores to which we have restricted attention tend to ellipses elongated in the y-direction with major axis three times larger than the minor axis; this can be seen by simply inspecting plots of the images under the mapping (3.1), or by analysis of this formula as ρ → 0. This feature can be seen from figure 3 where the small-area pores appear to be elliptical in shape. In this limit, the graph tends to the value 10/3 in agreement with the results of Zimmerman [9] for an isolated ellipse under a hydrostatic loading. He showed that for the conformal mapping to the exterior of an ellipse of the form and, for an ellipse with semi-axes in the ratio 1 to 3, then a 1 = − 1 2 . This feature, which serves as a check on our more general mathematical approach, is to be expected, because in this small area limit, the typical lengthscale of each pore is small compared with their separation meaning that they become increasingly unaware of each other's presence. Then, each pore acts like an isolated pore under uniform hydrostatic pressure.
As the pore area increases from zero the interaction effects and shape deformations lead to a net decrease in compressibility. However, eventually, as the pore areas increase further the interaction effects lead to an increase of the compressibility above the starting value of 10/3 relevant for small pore areas. The limit ρ → 1 − is a singular mathematical limit and must be treated with care; the very steep vertical slope of the graph in figure 5  noted that for the validity of our two-pore analysis, we require ρ to be strictly less than unity on topological grounds; moreover, the rate of convergence of the infinite product (2.2) deteriorates as ρ → 1 − . In all calculations used to plot figure 5, we allowed ρ to increase as high as 0.9 and it is interesting that, as ρ becomes larger, the two pores almost merge in such a way that the circumscribed boundary tends to that of a single circular pore. Figure 5 shows that the associated compressibility appears to tend to the value 4 and that behaviour might also be expected: while the compressibility of a purely circular pore is well known to be 2 [9,10] the apparent factorof-2 discrepancy arises only, because we have normalized the pore compressibility with respect to the area A of a single pore. However, it is clear from the geometry of figure 4 that this is precisely half the area of the full circumscribed semicircle marked out by the 'exterior part' of each pore boundary.
In both limits ρ → 0 and ρ → 1 − , the analytical solutions for two pores therefore agree with the relevant single pore results lending us confidence in the analysis. The compressibility of other classes of pore shapes can be studied using a similar analysis. This is discussed further in §7.

Two pores subject to far-field shear
As a second example of the versatility of our construction, consider the same class of elastic bodies with two pores but now with an imposed state of pure shear stress at infinity. The following results generalize, to two non-circular pores, the work of Ekneligoda & Zimmerman [11].
Following Ekneligoda & Zimmerman [11], we consider two problems, called problems 1 and 2. In problem 1, the far-field shear is taken to be aligned with the real and imaginary axes, so that τ xy = τ yx = τ and the relevant complex potentials are simplỹ The stress state of problem 1 has the required form at infinity, but produces an unwanted traction on the boundaries of the two pores. Problem 2, which is now defined in the elastic medium containing the two pores, is taken to have zero stress at infinity and has boundary tractions that, when added to those of problem 1, cancel them out. By linearity, the solution for a body with two traction-free pores and the required shear stress at infinity is given by the sum of the solutions to problems 1 and 2. The potentialsφ(z) andψ(z) in (6.1) solve problem 1, but it remains to solve problem 2, whose potentials will be denoted by φ(z) and ψ(z). In this situation, the boundary conditions on the two pores boundaries are now Again, we can argue on the basis of symmetry that φ(z) and ψ(z) are odd, and the constants in (6.2) are chosen to be consistent with this symmetry. On introduction of the functions defined in (5.6), the complex conjugate of the first boundary condition in (6.2) isΦ where we have used the fact thatζ = 1/ζ on C 0 . The complex conjugate of the second boundary condition in (6.2) isΦ On use of the key property (3.5) of the conformal mapping function, (6.4) becomes Hence, (6.3) and (6.5) again lead to the result The complex conjugate of (6.3) implies that By inspection, we see that Φ(ζ ) now inherits the singularity structure of both z(ζ ) andz(1/ζ ) in the annulus 1 < |ζ | < 1/ρ: it, therefore, has three simple poles located at 1/ √ ρ, e ±iθ / √ ρ. It is natural to pose that for some complex constants A, B, C and D that satisfy The appearance of the three K-functions means that Φ(ζ ) has the required simple pole singularities; the constraint (6.9) derives from the requirement that Φ(ζ ) satisfies (6.6). The constants A, B and C are determined by equating residues of the left-and right-hand sides of (6.7) at the three poles at 1/ √ ρ, e ±iθ / √ ρ. After some straightforward algebra to compute residues, it is found that where b is such that near ζ = 1/ √ ρ we havē √ ρ + analytic function, (6.11) (a formula for b is given in appendix B) and X is the same quantity defined in (5.17). From (6.9), we find which is purely imaginary. With A, B and C determined in (6.10), D is determined from (5.5): Finally, with both Φ(ζ ) and z(ζ ) now known explicitly, Ψ (ζ ) is given by (6.3): (6.14) In summary, we have found the stress fields for Problem 2 associated with two pores with an imposed state of pure shear at infinity to involve the two complex potentials and

Quadrature domains
Where did formula (3.1) come from? And how can other formulae, describing other pore shapes, be derived for which the analysis of this paper also works?
To answer these questions, we remark that the foregoing analysis has implicitly demonstrated the role played, in problems of linear plane elasticity, by a special class of doubly connected planar domains produced as the images of conformal mappings of the loxodromic form (3.1). This family of elastic bodies is just one example of a much broader class of planar domains known as (multiply connected) quadrature domains [18,20]. Other quadrature domains will similarly give rise, in principle, to solutions for the elastic stress problems considered here that are expressible in closed form. The question of how to construct multiply connected quadrature domains has been a topic of much recent research interest [18,20]. That constructive theory led to formula (3.1).
To explain this in more detail, note that we have focused here on the example of two pores as this case affords the clearest insights into the associated function theory. But, all the same mathematical ideas pertain to the case of any finite distribution of pores. If M + 1 pores are present, with M > 1, then to generalize our analysis the surrounding elastic medium should be taken to be an unbounded (M + 1)-connected quadrature domain. Then, for the conformal mapping, a suitable choice of pre-image ζ domain is the unit ζ -disc now with M smaller circular discs excised [23,24]; the annulus ρ < |ζ | < 1 used here is the unit disc with just a single, concentric, circular disc excised. It is possible to associate, to any such circular pre-image domain, a function called the Schottky-Klein prime function [23], often denoted as a function of just two variables, ω(ζ , α), even though it also depends on the choice of the circular preimage domain (i.e. the geometry of the pre-image circles). To within a constant of proportionality, the function P(ζ /α), which has clearly played a crucial role in our analysis, is exactly the Schottky-Klein prime function [23] associated with the concentric annulus ρ < |ζ | < 1. Crowdy & Marshall [24] have given a detailed account of how to construct conformal mappings to multiply connected quadrature domains using the Schottky-Klein prime function as a building block. Those conformal mappings have the property of being automorphic functions; in the two-pore case, those automorphic functions are precisely the loxodromic functions considered here.
Indeed, Crowdy & Marshall [25] have shown that the very same class of mappings (3.1) explored here for these plane elasticity problems also provide analytical solutions of the twodimensional Euler equations of fluid dynamics describing two rotating vortex patches. They also showed how that work can be extended to the case of any number of co-rotating vortex patches in equilibrium [26]; to do so, use is made of multiply connected quadrature domains and the function theory based on the Schottky-Klein prime function just described [23,24]. In the same way, the mathematics of Crowdy & Marshall [26] can be exploited to solve plane elasticity problems for elastic bodies with multiple pores by extension of the methodology of the present paper. This is left as the topic for future work.
It is important to point out a significant mathematical fact: Gustafsson [27] has shown that the class of multiply connected quadrature domains is dense (in an appropriate sense [27]) in the class of all multiply connected planar domains (with sufficiently smooth boundary components). To get a flavour of this in the doubly connected case, we note that generalizations of (3.1) that are also loxodromic functions (and, therefore, correspond to doubly connected quadrature domains) are given by (see [24] for more details): where α maps to z = 0, whereas β maps to z = ∞ (with ρ < |α|, |β| < 1) and, in order that the function is loxodromic, the other poles {r k |k = 1, . . . , N} and zeros {s k |k = 1, . . . , N} of the mapping (7.1) must be in the annulus 1 < |ζ | < 1/ρ and satisfy the single condition By taking N sufficiently large, and by appropriate choices of the parameters R, ρ, α, β, {r k , s k |k = 1, . . . , N}, approximate conformal mappings of this type to any unbounded doubly connected region can be found, in principle. Then, the analysis of this paper can be repeated for the new mappings (7.1) to find, also in terms of the Schottky-Klein prime function, analytical expressions for the complex potentials for the stress fields.
In this way, it is possible, in principle at least, to find approximating analytical solutions to the stress problems of plane elasticity in any multiply connected region by identifying a satisfactory quadrature domain approximant and generalizing the solution scheme laid out here. Of course, as quadrature domains with greater numbers of poles are needed to approximate the domains this becomes less feasible in practice. Nevertheless, a broad class of new analytical solutions to the basic problems of plane elasticity is essentially available. Certainly, elastic bodies corresponding to low order quadrature domains-that is, those with associated conformal mappings which are meromorphic functions with a small number of poles-can be studied explicitly, as we have done here.

Discussion
The stress fields around a family of two-pore configurations in an elastic body have been found in analytical form by invoking some novel mathematical ideas from quadrature domain theory.
One parameter in the description of the shapes studied here is the so-called conformal modulus ρ [23], which governs the degree of 'interaction' of the two pores. Our class of solutions, and its various generalizations, affords an opportunity to study, in an explicit way, how pore interactions affect quantities of physical interest such as compressibility, shear compliance and other measures of excess strain in a body owing to the presence of a pair of cavities. For illustration, we have investigated in detail the compressibility of a particular family of two pore configurations.
An alternative approach to the two-pore problem considered here is via purely numerical schemes, such as boundary integral methods. However, such methods would have difficulty in accurately resolving the stress fields in the vicinity of the high curvature, close-to-touching regions exhibited by the family of pore shapes studied here. Our analytical solutions should serve as valuable benchmarks for any such numerical formulations.
Finally, it is not without import to remark that the mathematical approach expounded here also works, of course, for simply connected quadrature domains. It is a basic fact [18,28] that the boundaries of simply connected quadrature domains can be parametrized by rational function conformal mappings. The function theory in that case becomes significantly easier with the solutions for Φ(ζ ) and Ψ (ζ ) also being rational functions. But then the approach we have demonstrated here reduces to ideas that are already well known in solid mechanics; see, for example, the review on the bending of thin plates using the so-called 'rational function method' by Hasebe & Wang [29]. What we have shown here is that if one wants to extend the rational function method to multiply connected regions, then the natural mathematical generalization is to consider the class of multiply connected quadrature domains; the boundaries of those domains can be parametrized by so-called automorphic functions [23,24,26] (i.e. the loxodromic functions of this paper) which are the natural generalizations of rational functions to higher-genus Riemann surfaces. The ideas of this paper therefore form the basis of what one might call the 'loxodromic function method' for doubly connected regions or the 'automorphic function method' for higher connected regions approximated by multiply connected quadrature domains.