Computation of free boundary minimal surfaces via extremal Steklov eigenvalue problems

Recently Fraser and Schoen showed that the solution of a certain extremal Steklov eigenvalue problem on a compact surface with boundary can be used to generate a free boundary minimal surface, i.e., a surface contained in the ball that has (i) zero mean curvature and (ii) meets the boundary of the ball orthogonally (doi:10.1007/s00222-015-0604-x). In this paper, we develop numerical methods that use this connection to realize free boundary minimal surfaces. Namely, on a compact surface, $\Sigma$, with genus $\gamma$ and $b$ boundary components, we maximize $\sigma_j(\Sigma,g) \ L(\partial \Sigma, g)$ over a class of smooth metrics, $g$, where $\sigma_j(\Sigma,g)$ is the $j$-th nonzero Steklov eigenvalue and $L(\partial \Sigma, g)$ is the length of $\partial \Sigma$. Our numerical method involves (i) using conformal uniformization of multiply connected domains to avoid explicit parameterization for the class of metrics, (ii) accurately solving a boundary-weighted Steklov eigenvalue problem in multi-connected domains, and (iii) developing gradient-based optimization methods for this non-smooth eigenvalue optimization problem. For genus $\gamma =0$ and $b=2,\dots, 9, 12, 15, 20$ boundary components, we numerically solve the extremal Steklov problem for the first eigenvalue. The corresponding eigenfunctions generate a free boundary minimal surface, which we display in striking images. For higher eigenvalues, numerical evidence suggests that the maximizers are degenerate, but we compute local maximizers for the second and third eigenvalues with $b=2$ boundary components and for the third and fifth eigenvalues with $b=3$ boundary components.


Introduction
Recently, A. Fraser and R. Schoen discovered a rather surprising connection between an extremal Steklov eigenvalue problem and the problem of generating free boundary minimal surfaces in the Euclidean ball [FS11,FS13,FS15]. These findings have been further developed [FTY14,FS19,GL20] and were recently reviewed in [Li19]. In this paper, we develop numerical methods to further investigate this connection. We first briefly review some of these previous results before stating the contributions of the present work.
Here, for fixed surface Σ with genus γ and b boundary components, we consider the dependence of the j-th Steklov eigenvalues on the metric, i.e., the mapping g → σ j (Σ, g). It is known that for any smooth Riemannian metric g, we have the following upper bound on the j-th Steklov eigenvalue in terms of the topological invariants γ and b, (2) σ j (Σ, g) L(∂Σ, g) ≤ 2π(γ + b + j − 1) ∀j ∈ N.
Free boundary minimal surfaces. Denote the closed n-dimensional Euclidean unit ball by B n := {x ∈ R n : |x| ≤ 1} and the (n − 1)-dimensional unit sphere by S n−1 = ∂B n . Let M ⊂ B n be a d-dimensional submanifold with boundary ∂M = M ∩ S n−1 . We say that M is a free boundary minimal submanifold in the unit ball if (i) M has zero mean curvature and (ii) M meets S n−1 orthogonally along ∂M.
When d = 2, we call M a free boundary minimal surface in the unit ball or, more simply, a free boundary minimal surface. For a good visual aid to understanding the definition of free boundary minimal surfaces (and a peak at the results of this paper), we recommend the reader take a look at the free boundary minimal surfaces displayed in Figures 13 and 14.
Fraser and Schoen's connection. Fraser and Schoen observed that a d-dimensional submanifold M ⊂ B n with boundary ∂M = M ∩ S n−1 is a free boundary minimal surface if and only if the coordinate functions x i , i = 1, . . . , n restricted to M are Steklov eigenfunctions with eigenvalue σ = 1. Furthermore, they showed the following theorem.
Theorem 1.1 ( [FS13]). Let Σ be a compact surface with boundary. Suppose that g 0 is a smooth metric on Σ attaining the supremum in (3) for some j ∈ N. Let U be the n-dimensional eigenspace corresponding to σ j (Σ, g 0 ). Then, there exist independent Steklov eigenfunctions u 1 , . . . , u n ∈ U which give a (possibly branched) conformal immersion u = (u 1 , · · · , u n ) : Σ → B n such that u(Σ) is a free boundary minimal surface in B n and, up to rescaling of the metric, u is an isometry on ∂Σ.
Theorem 1.1 gives a method for using the solution of (3) to compute free boundary minimal surfaces. The simplest such example is the equatorial disk, obtained as the intersection of B 3 with any two-dimensional subspace of R 3 . This can be constructed from Weinstock's result that inequality in (2) with j = 1, γ = 0, and b = 1 is attained only by the round disk, D [Wei54]. In this case, for the eigenvalueσ 1 (0, 1) = 2π, we have the two-dimensional eigenspace given by span{x, y}.
The equatorial disk is given as the map u : D → R 2 , defined by u(x, y) = x y .
For genus γ = 0 and b = 2 boundary components, the extremal metric is rotationally invariant and the corresponding free boundary minimal surface is the critical catenoid. We will discuss this example further in Section 3. For genus γ = 0 and b ≥ 3 boundary components, the extremal metric is not known explicitly, but it is known that the corresponding free boundary minimal surface is embedded in B 3 and star-shaped with respect to the origin [FS13]. In [GL20], the authors used homogenization methods to construct surfaces that have large first Steklov eigenvaluẽ σ 1 . In particular, free boundary minimal surfaces of genus γ = 0 with particular symmetries (e.g., symmetries of platonic solids) were constructed numerically. The authors proved that the first nonzero Steklov eigenvalue, σ 1 , of these surfaces is 1 and emphasized that it is not known whether these surfaces have extremal first eigenvalues among all surfaces with the same genus and number of boundary components. We will compare our results to these surfaces in Section 5.
In [FTY14], Fan, Tam, and Yu extended the study of (3) to higher values of j on the cylinder (γ = 0, b = 2) among rotationally symmetric conformal metrics. They obtained different results for even and odd eigenvalues. They showed that the maximum of theσ 2j−1 , j ∈ N among all rotationally symmetric conformal metrics on the cylinder is achieved by the j-fold covering of the critical catenoid immersed in R 3 . The maximum ofσ 2 is not attained. The maximum of theσ 2j for j ≥ 2 among all rotationally symmetric conformal metrics on the cylinder is achieved by the j-fold covering of the critical Möbius band. These results will be further discussed in Section 3 and further compared to our computed surfaces in Section 5.
Results and outline. In this paper, we develop computational methods for solving the extremal Steklov eigenvalue problem (3) and thus generating free boundary minimal surfaces via Theorem 1.1. This approach is used to realize free boundary minimal surfaces beyond the known examples of equatorial disks, the critical catenoid, the critical Möbius band, and their higher coverings discussed above.
In Section 2, we explain how the conformal uniformization of multiply connected domains can be used to significantly reduce the complexity of the general Steklov eigenproblem (1) and extremal Steklov eigenproblem (3). The argument relies on two ingredients: (1) The uniformization result that for a smooth, compact, connected, genus-zero Riemannian surface with b boundary components, (Σ, g), there exists a conformal mapping f : (Σ, g) → (Ω, ρI), where Ω is a disk with b − 1 holes and ρI is a conformally flat metric.
(2) The composition v • f of a function v with a conformal map f is harmonic if and only if v is harmonic.
Let D = {x ∈ R 2 : |x| ≤ 1} be the unit disk and be a punctured unit disk with b − 1 holes, This argument implies that it is sufficient to consider the family of (flat!) Steklov eigenproblems, where ∆ is the Laplacian on Ω, ∂ n is the outward normal derivative, and ρ > 0 is a density function. The extremal Steklov eigenvalue problem (3) for genus γ = 0 is transformed tõ Here,σ j = σ j L, σ j is the j-th nontrivial eigenvalue satisfying (4), and L = ∂Ωc,r ρ(x) dx is the total length of ∂Ω c,r . The first two constraints simply state that the holes are contained in the domain and are pairwise disjoint.
In Section 3, we explicitly solve the Steklov eigenvalue problem on a rotationally symmetric annulus (i.e., γ = 0, b = 2, c 1 = 0, and ρ constant on each boundary component) and describe the critical catenoid and its higher coverings in detail. These Steklov eigenvalues and corresponding free boundary minimal surfaces will be used to verify our computational methods.
In Section 4, we develop numerical methods for computing Steklov eigenvalues satisfying (4) on multiply connected domains, computing the solution to the optimization problem (5), and the computation of free boundary minimal surfaces from the Steklov eigenfunctions. In brief, we use the method of particular solutions to compute Steklov eigenvalues, gradient-based interior point methods for the optimization problem, and compute the mapping to a surface by minimizing a particular energy. These methods build on previous computational methods for extremal eigenvalue problems on Euclidean domains, including minimizing Laplace-Dirichlet eigenvalues over Euclidean domains of fixed volume or perimeter [Oud04, Ost10, AF12, OK13, OK14, AO17, BBG17], maximizing Steklov eigenvalues over two-dimensional Euclidean domains of fixed volume [AKO17,BBG17]. These methods have recently been extended to more general geometric settings. In particular, [KLO17] maximized Laplace-Beltrami eigenvalues over conformal classes of metrics with fixed volume and compact Riemannian surfaces of fixed genus (γ = 0, 1) and volume.
In Section 5, we present the results of numerous computations. For genus γ = 0 and b = 2, . . . , 9, 12, 15, 20 boundary components, we numerically solve the extremal Steklov problem (5) for the first eigenvalue. We include figures displaying the optimal punctured disks and three linearlyindependent eigenfunctions associated to the first eigenvalue, as well as tabulate the values of the obtained Steklov eigenvalues. We also plot the associated free boundary minimal surfaces, which are visually striking. Finally, in Section 5, we also present results for maximizing higher eigenvalues. Here, numerical evidence suggests that the maximizers are degenerate, but we compute local maximizers for the second and third eigenvalues with b = 2 boundary components and for the third and fifth eigenvalues with b = 3 boundary components. For brevity, we were only able to report the results for selected values of b and j; the results of additional computations can be found onÉ. Oudet's website [Oud20], along with gifs.
We conclude in Section 6 with a discussion.

The Euclidean Steklov eigenproblem
In Section 2.1, we explain how the conformal uniformization of multiply connected domains can be used to significantly reduce the complexity of the general Steklov eigenproblem (1) and extremal Steklov eigenvalue problem (3) to obtain the Euclidean Steklov eigenproblem and (4) and extremal Steklov eigenvalue problem (5), respectively. In Section 2.2, we also compute the eigenvalue derivatives with respect to the density and shape parameters and discuss optimality conditions for the extremal Steklov eigenvalue problem (5).
2.1. Conformal uniformization of multiply-connected surfaces and the Steklov eigenproblem. The uniformization theorem for compact, genus-zero Riemann surface without boundary states that such surfaces can be conformally mapped to the Riemann sphere. Here, we use a generalization of this result for multiply-connected surfaces; see [Hen86,Theorem 17.1b], [GL99], [ZYZ + 09], and [JGHW18,p.123].
Theorem 2.1 ( [GL99]). Suppose (Σ, g) is a smooth, compact, connected, genus-zero Riemann surface with b boundary components. Then Σ can be conformally mapped to a unit disk with b − 1 circular holes. That is, there exists a punctured unit disk with b − 1 holes, Ω c,r = D \ ∪ b−1 i=1 D i , and a conformal map f : (Σ, g) → (Ω c,r , ρI), where ρI is a conformally flat metric. Furthermore, two such mappings differ by a Möbius transformation.
Remark 2.2. The uniqueness of the conformal map up to a Möbius transformation means that it is possible to center one of the holes at the origin and center another hole on the positive x-axis. Thus, fixing these three parameters, the dimension of the parameter space of hole centers and radii , is 1 for b = 2 and 3b − 6 for b ≥ 3, which is the dimension of the conformal module.
We now sketch a brief derivation of (4) from (1). Let f : (Σ, g) → (Ω c,r , ρI) be a conformal mapping. It is well-known that v = u•f : Σ → R is harmonic if and only if u : Ω c,r → R is harmonic [Olv17]. This justifies (4a). We show (4b) on a flat domain for simplicity. Write Remark 2.2 shows that our parameterization of Ω c,r is over-complete, as the following example further demonstrates.
Example 2.3. Denote Ω 1 as an eccentric annulus with boundaries |z − c 1 | < r 1 and |z| < 1 and Ω 2 as an concentric annulus r 2 < |x| < 1 where c 1 , r 1 , r 2 are real numbers and x, z ∈ C. A conformal mapping h : Ω 1 → Ω 2 is given by az where a and r 2 are determined by mapping c 1 + r 1 , c 1 − r 1 to r 2 , −r 2 and satisfy .
In this example, z = h(x) = x+a 1+ax , and In Figure 1, the mapping is shown for c 1 = r 1 = 1 4 and the resulting a = r 2 = 2 − √ 3. Thus, the eccentric annulus Ω 1 with boundary density ρ = 1 has the same Steklov spectrum as the concentric annulus Ω 2 with boundary density ρ(x) given above. In particular, this example shows that the decomposition of perturbations of a metric into conformal and non-conformal directions is not equivalent to either changing (c 1 , r 1 ) or ρ, respectively. While changing ρ is a conformal perturbation, a change in (c 1 , r 1 ) gives a perturbation to the metric that has components in both the conformal and non-conformal directions.
The following two examples illustrate what happens to the boundary density ρ when Σ becomes "pinched".
Example 2.4. We consider the conformal mapping h : D → Ω α from the unit disk |x| ≤ 1 to the Hippopede domain, Ω α , AK19]. When α = 1, this is the identity mapping on the unit disk and as α → 0 + , it maps a unit disk to two "kissing" disks. In the left and center panels of Figure 2, the mapping is shown for α = 1 10 . Here, we compute, Let x = e iθ . In the right panel of Figure 2, we plot ρ α (θ) for α = 1 50 , 1 10 , and 1 5 . We observe that ρ α (θ) becomes singular as α → 0 + at θ = π 2 and 3π 2 . In Example 2.4, the density is singular at two points. The following example illustrates how the density function can become singular at a single point.
Example 2.5. We consider a radius r 1 := 0.833 disk, Ω 1 = {|x| < r 1 }, (see Figure 3(c)) and a domain Ω 2 consisting of the union of two disks with radii r 1 and 1 and a 'neck' of width 2α (see Figure 3(a)). Define the conformal mapping h : Ω 1 → Ω 2 as the composition of the two functions Figure 2. A conformal mapping from a Hippopede shape to a unit disk. See Example 2.4.
The constants a, b, c are chosen as Figure 3(a) and (b). The constant β is chosen so that Ω 1 maps to a unit disk and the zero in Ω 1 maps to −i(1 − β). See Figure 3(b) and (c).When β is small, this function maps points which are uniformly distributed on ∂Ω 1 to points that accumulate near −i on the unit disk. The boundary density, ρ, can be obtained via the product rule, As shown in Figures 3(d), the density reaches a large value at θ = π 2 . Figure 3(e) shows the detail profile of the density function about one.
2.2. Eigenvalue derivatives with respect to the density and shape parameters. In this section, we consider σ andσ = σL as a function of ρ and the shape Ω c,r . We first compute the derivatives with respect to ρ.
We describe below optimality conditions when the multiplicity of the optimized eigenvalue is greater than one. Our formulation is highly inspired by previous articles [ESI + 07, FS15,BO16]. We first need the following regularity result; see [LP15, Theorem 3.2].
We can now evaluate directional derivatives based on previous parametrizations: Lemma 2.8. Let σ be an eigenvalue of multiplicity p > 1 of the weighted Steklov system (4) for some nonnegative boundary density ρ. Denote by E σ the corresponding eigenspace. Let ρ ε = ρ+εδρ be a perturbation of ρ for some δρ ∈ L 2 (∂Ω) . Let (σ i (ε)) 1≤i≤p and (u i (ε)) 1≤i≤p be some smooth parametrizations as the ones given by Lemma 2.7. Then σ i = d dε σ i (ε)| ε=0 are the eigenvalues of the quadratic form q δρ defined on E σ ⊂ L 2 (∂Ω, ρ) by Differentiationg this equality with respect to ε and evaluating at ε = 0 gives Thus, with v = u j (0) and using (7) replacing i per j and v by u i (0), we obtain Moreover, the σ i are eigenvalues of this quadratic form.
We can now establish optimality conditions with respect to the boundary density in case of multiple eigenvalues.
Proposition 2.9. Let j ≥ 1 and Ω a smooth domain of R 2 . Assume a nonnegative ρ ∈ L 2 (∂Ω) maximizes the product σ j (ρ)L(ρ) among all nonnegative functions of L 2 (∂Ω) where L(ρ) = ∂Ω ρ dx and σ j (ρ) is the j-th eigenvalues of system (4). If σ j (ρ) is of multiplicity p > 1 and E σ j its eigenspace, there exists a basis of p functions u 1 , . . . , u p of E σ j which satisfy Proof. The proposition is an almost direct consequence of Lemma 2.8 and of Hahn-Banach separation theorem. Consider the convex hull K = Co u 2 , u ∈ E σ j . We want to prove that the function identically equal to one belongs to K. If it is not the case, by Hahn-Banach theorem applied to the finite dimensional normed vector subspace of C 1 (∂Ω) spanned by K and 1, there exists a function δρ ∈ C 1 (∂Ω) such that ∂Ω δρ dx > 0 and which satisfies, for all u ∈ E σ j , ∂Ω u 2 δρ dx ≤ 0.
This last inequality asserts that the quadratic form q δρ on E σ j has nonnegative eigenvalues. Thus, both the p eigenvalues and the weighted length increase in the direction of δρ. As a consequence, for ε small enough, the product of σ j (ρ + εδρ)L(ρ + εδρ) is strictly greater than σ j (ρ)L(ρ) due to the strict inequality of the separation result which contradicts the optimality.
To compute the derivatives of σ andσ with respect to the centers c and radii r, we first compute the shape derivative with respect to perturbations of the boundary of Ω c,r . This result extends a result in [DKL14, AKO17, BBG17] to ρ = 1.
Proposition 2.10. Consider the perturbation x → x+τ v. Then a simple (unit-normalized) Steklov eigenpair (σ, u) satisfies the perturbation formula wheren is the outward unit normal vector,t denotes the tangential direction, and where κ is the signed curvature of the boundary. We also have L = ∂Ω κρ(v ·n) − ρ t (v ·t) dx.
Differentiating the normalization equation, ∂Ω ρu 2 dx = 1, we have that where κ is the curvature of the boundary and ρ = −∇ρ · v. Extending ρ constantly in the normal direction, we have ρ + (v ·n)ρ n = −ρ t (v · t) where t denotes the tangential direction. We then have that Combining this with (9), we obtain the desired result.
Using Proposition 2.10, we can now compute the derivatives of σ andσ for the domain Ω c,r = D \ ∪ b−1 i=1 D i with respect to a center c i and radius r i of D i as follows. To compute the derivative with respect to r i , we choose a perturbation v so that v ·n = −1 and v ·t = 0 on ∂D i .
Then, noting that κ = −1/r i , we obtain To compute the derivative with respect to c i , we take two perturbations v of the form v ·n = cos θ and v ·t = sin θ on ∂D i and v ·n = sin θ and v ·t = − cos θ on ∂D i , to obtain (11) ∇ c i σ = ∂Ω |∇u| 2 − 2ρ 2 σ 2 u 2 + σ r i ρu 2 cos θ sin θ + σρ t u 2 sin θ − cos θ dx.
Remark 2.11. In [FS15], a detailed study of perturbations to the metric yield two conditions for a maximal Steklov eigenvalue. The first comes from the study of perturbations in "conformal directions" and, as in Proposition 2.9, result in the existence of eigenfunctions {u j } n j=1 such that the map U = [u 1 | · · · |u n ] : Ω → B n satisfies U (∂Ω) ⊂ S n−1 . The second condition comes from the study of non-conformal perturbations of the metric and give that the map U : Ω → B n has isothermal coordinates, i.e., satisfies Since a change in the parameters (c, r) gives a perturbation to the metric that has components in both the conformal and non-conformal directions (see Remark 2.2 and Example 2.3), this second condition is nontrivial to obtain from (10) and (11).

Steklov eigenvalues of rotationally symmetric annuli and the critical catenoid
Here, we discuss the Steklov eigenvalues of rotationally symmetric annuli, the critical catenoid, and coverings of the critical catenoid. These results are also discussed in [FS11,FTY14] using cylindrical coordinates, but it useful to review these computations and have them written in annular coordinates for comparison and discussion; see also [Mar14,Dit04].
For each k = 1, 2, . . ., there are also eigenfunctions that are oscillatory in θ of the form u(r, θ) = (Ar k + Br −k ){cos kθ, sin kθ}, where A, B are constants. Here, the brackets indicate that we can choose either cos or sin; the corresponding eigenvalue has multiplicity two. Using the boundary conditions we obtain the 2 × 2 generalized eigenproblem, This is equivalent to the eigenproblem k ρ 1 sρ s sinh(−k log s) from which one obtains the real positive eigenvalues σ k,± = k 2ρ 1 sρ s coth(−k log s) ρ 1 + sρ s ± (ρ 1 + sρ s ) 2 − 4ρ 1 sρ s tanh 2 (−k log s) .
In Figure 4(left), for ρ s /ρ 1 = 11.01609, we display the length-normalized Steklov eigenvalues for various values of s. The eigenvalue corresponding to the radially symmetric eigenfunction is plotted in red. The thin vertical line indicates the value s = 0.090776. For this value of s, the first Steklov eigenvalue has multiplicity three and length-normalized eigenvalueσ = 10.47478. In Figure 4(right), we plot contours of two of the eigenfunctions; the third can be obtained by rotating the image of the lower eigenfunction by π 2 . 3.2. Extremal eigenvalues for rotationally symmetric annuli. We consider the extremal eigenvalue problem for rotationally symmetric annuli, max s,ρs,ρ 1σ j ,σ j := σ j L.
3.2.1. The first eigenvalue. We first consider j = 1. By the symmetry of ρ 1 and sρ s , we obtain the optimality condition sρ s = ρ 1 =: ρ.
In this case, we have the two length-normalized eigenvalues and associated L 2 (∂Ω, ρ)-normalized eigenfunctions The two values of σL are equal when s is the unique solution of the transcendental equation The solution is approximately given by s = 0.090776. We now consider the map U : A s → B 3 , defined by Note that this map has coordinates that are linear combinations of the above eigenfunctions. One can check that these are isothermal coordinates, i.e., and satisfy U (∂A s ) ⊂ S 2 ⊂ R 3 , i.e., Furthermore, it is not difficult to check that U (A s ) is the critical catenoid. That is, is a catenoid and the critical catenoid is the catenoid with α = α * = β 2 + cosh 2 β − 1 2 where β = − log √ s ≈ 1.19968 is the unique solution of β = coth β. It is known that the critical catenoid is a free boundary minimal surface [FS15].

Higher eigenvalues.
For larger values of j, we numerically solve (13). In Figure 5, we plot the value ofσ j as a function of s and ρ s /ρ 1 for j = 1, . . . , 6. The maximum value of σ j L is indicated and data for the maximum values is also tabulated. Observe that for j = 1, 3, . . . , 6, we have that sρ s = ρ 1 .
For odd j = 2m − 1, m ∈ N, from the results of Fan, Tam, and Yu [FTY14], we have that the extremum is attained at the crossings of the two length-normalized eigenvalues with associated L 2 (∂Ω, ρ)-normalized eigenfunctions The two values of σL are equal when s is the unique solution of the transcendental equation We obtainσ 2m−1 = mσ 1 , for m ≥ 1. The extremal metric is achieved by the m-fold cover of the critical catenoid. For even j, Fan, Tam, and Yu [FTY14] show the following. For j = 2, the extremal value is not attained among rotationally symmetric annuli and for even j ≥ 4, the extremal value is attained. For m ≥ 2, we have σ 2m L = 4mπ tanh(

Computational Methods
In Section 2, we described how conformal maps could be used to reduce the general Steklov eigenproblem (1) to the Euclidean Steklov eigenproblem (4). In this section, we describe the computational methods used to solve the Euclidean Steklov eigenproblem (4), optimization methods used to solve the extremal eigenvalue problem (3), and methods for computing the minimal surface from the Steklov eigenfunctions.
4.1. Solving the Euclidean Steklov eigenproblem (4). We use the method of particular solutions to solve the Steklov eigenproblem (4). This method for multiply-connected Laplace problems was recently discussed in [Tre18]. The methods rely on the following Theorem.
Theorem 4.1 (Logarithmic Conjugation Theorem [Tre18]). Suppose Ω is a finitely connected region, with K 1 , . . . , K N denoting the bounded components of the complement of Ω. For each j, let a j be a point in K j . If u is a real valued harmonic function on Ω, then there exist an analytic function f on Ω and real numbers c 1 , . . . , c N such that Let M ∈ N * and consider some fixed punctured disk Ω c,r . Based on Theorem 4.1, we define the finite basis B to approximate solutions of eigenvalue problem (4) as the union of the harmonic rescaled real and imaginary parts of the functions For instance, we rescaled the basis polynomial Re 1 (z−c 2 ) 3 by a factor r 3 2 so that this basis function takes values of order 1 on the second circle. Consider now (p l ) 1≤l≤L a uniform sampling with respect to arc length of ∂Ω c,r . Using B, we approximate solutions of eigenvalue problem (4b) by the solution of the non symmetric square generalized eigenvalue problem 1≤l≤L, φ∈B and B = (φ(p l )) 1≤l≤L, φ∈B .
Example 4.2. To illustrate the complexity of the approach to obtain a fine approximation of eigenvalues, we considered a circular domain with four holes and L = 5000 points; see Figure 6(left). We evaluated the first six nontrivial eigenvalues with a high number of B elements for M = 50. In Figure 6(right), you can observe the evolution of the error with respect to M for M taking values from 2 to 10. Taking the converged values as an approximation of the exact ones, in this specific example, it can be observed that with M = 10 the error is already smaller than 10 −8 . Here, the first nontrivial eigenvalue has multiplicity two, so the curves are almost indistinguishable.  Table 1. The first ten nontrivial Steklov Eigenvalues, σ j , of the Hippopede domain, Ω α , for α = 0.1, 0.06, 0.04. The last column are the values, known analytically, that appear in the limit as α → 0.
Example 4.3. We now consider a geometric convergence study related to Example 2.4; see also Figure 2. Using the mapping from the unit disk to the Hippopede domain, Ω α , we study the limit as α → 0. Our computations are performed on the unit disk with non-constant density, ρ, as given in Example 2.4. In the limit, the density becomes singular, and the purpose of this example is to illustrate that a weakness of our numerical method is that we cannot accurately compute eigenvalues of pinched domains (α → 0) or, equivalently, if the density is singular. The results are displayed in Table 4.1. The values for the disjoint union of two radius 0.5 disks, obtained in the limit α → 0, are given in the rightmost column of Table 4.1. We note a very slow convergence of the eigenvalues as α → 0.

4.2.
Optimization methods for extremal Steklov eigenvalues (5). We used gradient-based optimization methods to solve the extremal Steklov eigenvalue problem (5). We first describe our parameterization of the boundary 4.2.1. Parameterizing the geometry. Let ρ ∈ L ∞ (∂Ω c,r ) be the boundary density and denote the restriction of ρ to the i-th disk boundary by Finally, denote D k := D and ρ k the restriction of ρ to ∂D k . Thus, if Ω c,r has b boundary components, the geometry is described by the parameters Since ∂D(c i , r i ) ∼ = S 1 , we expand each ρ i in the truncated Fourier series From Remark 2.2, it would be possible to center one of the holes at the origin and another on the positive x-axis. However, we found that the representation of the boundary density ρ for finite basis size (finite N ) was better without fixing these centers. 4.2.2. Gradient based optimization methods. As in [AKO17], to handle multiple eigenvalues, we trivially transform (5) into the following problem max t (16a) We approximated the positivity constraint ρ ≥ 0 by imposing the positivity on all L sample points, This approximation leads to linear inequalities with respect to the coefficients (A i,l , B i,l ) only. We also augment the previous optimization problem with the geometrical constraints in (5) by imposing the (few) quadratic constraints on the variables (c i , r i ) 1≤i≤k−1 : Using the derivatives computed in (6), (10), and (11), together with the interior point method implemented in [BNW06], we solved (16). All results of section 5, have been obtained with the following parameters: M = 30 (maximal order of basis elements), L = 10 4 (number of sampling points) and at most 5, 000 iterations to reach a first order optimality condition criteria to a relative precision of 10 −6 . Observe that in all cases, we were able to recover the multiplicity three of the optimal eigenvalue up to 6 digits.
In our implementation, the computational cost is proportional to the number of connected components of the boundary. For instance, one hour of computation on a standard laptop was required to obtain the desired precision for three boundary components.

Computing the free boundary minimal surface from the Steklov eigenfunctions.
At this point we assume that we have successfully solved the extremal Steklov problem (5) and want to use Theorem 1.1 to compute the associated free boundary minimal surface using the Steklov eigenfunctions.
Let σ denote the optimal eigenvalue and assume that it has multiplicity n. Define the mapping v = [v 1 , . . . , v n ] : Ω → R n , where {v i } n i=1 is some choice of basis for the n-dimensional eigenspace. For A ∈ R n , we consider the map u A : Ω → R n , defined by x ∈ Ω.
To identify the matrix A, so that u A : Ω → R n satisfies (17), we construct the objective function (18) where W (u) = 1 4 (|u| 2 − 1) 2 . We then minimize J(A) over A ∈ R n×n . In all experiments in section 5, using this selection process, we were able to obtain three eigenfunctions which take values in the sphere on ∂Ω to an absolute pointwise error bounded by 10 −3 . Moreover, since we have a parameterization of the surface, using the well-known analytic formula, we were able to compute the mean curvature of the surfaces, which in all cases was bounded by 10 −2 . The mean curvature and the Gaussian curvature are plotted on the free boundary minimal surface at [Oud20]. Additionally, the angle that the boundary makes with the normal vector to the sphere is less than one degree.         Table 2. For different number of boundary components b, we report the value of the first nontrivial normalized Steklov eigenvalueσ 1 = σ 1 L, the value obtained by [GL20], and the configuration of the centers of the boundary components. For b = 8 and b = 20, our configuration of boundary components differs from [GL20], so the values should not be directly compared (indicated with an asterisk). Figure 15. Convex polytopes associated to the center of mass of boundary connected components of minimal surfaces in the ball. First row. Four (first plot) and six boundary connected components (the two remaining plots). Second row. Two views of a square antiprism associated to a minimal surface with a boundary made of height connected components and an icosahedron associated to a minimal surface with twelve connected components in its boundary (last plot).

Numerical solutions of the extremal Steklov eigenvalue problem and the corresponding free boundary minimal surfaces
In this section, we describe the solutions for the extremal Steklov eigenvalue problem (5), for various number of boundary components (BC), b, and eigenvalue number, j, and the corresponding free boundary minimal surfaces (FBMS). 5.1. First nontrivial eigenvalue (k = 1). We first consider the first nontrivial eigenvalue (k = 1) for varying numbers of BC, b = 2, . . . , 9, 12, 15, 20. In each case, the multiplicity of the extremal eigenvalue is three, as expected [FS15]. In Figure 7, we plot the optimal punctured disks, Ω c,r , for b = 2, . . . , 9 and b = 12 BC. In Figures 8 and 9, we plot three linearly independent eigenfunctions associated to the first eigenvalue on their respective punctured disk for b = 2, 3, 4, 5 BC. For these values of b, the corresponding optimal densities are plotted in Figures 10, 11, and 12. In Figures 13  and 14, we plot the corresponding (approximate) FBMS in the ball for b = 3, 4, 5, 12, 15 BC. In all cases, the BC of the FBMS are positioned at very symmetric locations, as further illustrated in Figure 15. Values ofσ 1 and additional information about these configurations are recorded in Table 4.3. Additional figures, including gifs, can be found at [Oud20] and were not included here for brevity. Figure 16. Six linearly independent eigenfunctions associated to the third eigenvalue for three boundary components.
We now make a few more detailed remarks for the problem with the various number of BC, b, considered, especially for values of b that are related to the platonic solids. For some values of b, we also compare to the FBMS discussed in [GL20].
For b = 2, we recover the critical catenoid, the known FBMS [FS15] that we also discussed in Section 3. Note that in Figure 7 the hole is centered within the disk and in Figure 10, the density is constant on each BC. The eigenfunctions plotted in Figure 8 exhibit symmetries and are explicitly given in Section 3; see Figure 4(right).
For b = 3, the FBMS has BC positioned with centers on an equilateral triangle inscribed on a great circle of the sphere; see Figure 13. Interestingly, the holes in the domain, Ω c,r , are slightly asymmetrically configured; see Figure 7. The densities plotted in Figure 10 do not exhibit symmetry. The eigenfunctions plotted in Figure 8 do not exhibit symmetries, but this could be a result of our (arbitrary) choice within the three dimensional eigenspace.
For b = 4, the FBMS has BC positioned with centers at the vertices of a regular tetrahedron; see Figure 13. This is further illustrated in Figure 15, where the BC are overlaid on a regular tetrahedron. A similar minimal surface was computed in [GL20] and the value ofσ 1 is within 10 −4 ; see Table 4.3. In Figure 7, the holes in the domain, Ω c,r , are slightly asymmetrically configured. In Figure 11, the density on the outer boundary is nearly constant and the densities on the inner boundaries are similar to each other. There is no clear structure to the eigenfunctions potted in Figure 9.
For b = 5, the FBMS has BC positioned with centers at the vertices of a triangular bipyramid; see Figure 14. In Figure 7, the holes in the domain, Ω c,r , are not only asymmetrically configured, but the radii of the holes vary. In Figure 11, the density on the outer boundary is nearly constant and the densities on the inner boundaries are similar to each other. Again, the eigenfunctions plotted in Figure 9 do not appear to be structured.
For b = 6, the FBMS has BC positioned with centers at the vertices of a regular octahedron; see Figure 14. This is further illustrated in Figure 15, where the BC are overlaid on a regular Figure 17. Nine first linearly independent eigenfunctions associated to the fifth eigenvalue for three boundary components. octahedron. Again, a similar minimal surface was computed in [GL20] and the value ofσ 1 is within 5 × 10 −3 ; see Table 4.3. In Figure 7, the holes in the domain, Ω c,r , are slightly asymmetrically configured; there is a small hole near the origin and four holes of equal radii roughly centered at the vertices of a square. In Figure 11, the density on the outer boundary is nearly constant and the densities on the inner boundaries are similar to each other.
For b = 7, the FBMS has BC positioned at the vertices of a pentagonal bipyramid. Figures of the FBMS can be found at [Oud20]. In Figure 7, the domain, Ω c,r , has a small (uncentered) hole surrounded by five holes.
For b = 8, the FBMS has BC positioned at the vertices of a square antiprism; see [Oud20] and Figure 15. Interestingly, we obtainσ 1 ≈ 16.4954 for this surface, which is larger than the value obtained for the FBMS with BC at the vertices of a cube, as discussed in [GL20], with valuẽ σ 1 ≈ 16.4249; see Table 4.3. In Figure 7, the domain, Ω c,r , has three smaller holes surrounded by four larger holes. For b = 9, the FBMS has BC positioned at the vertices of a triaugmented triangular prism. Figures of the FBMS can be found at [Oud20]. In Figure 7, the domain, Ω c,r , has three smaller holes surrounded by five larger holes.
For b = 12, the FBMS has BC positioned at the vertices of a regular icosahedron; see Figure 14. This is further illustrated in Figure 15, where the BC are overlaid with a regular icosahedron. A similar minimal surface was computed in [GL20] and the value ofσ 1 is within 10 −4 ; see Table 4.3. In Figure 7, the domain, Ω c,r , have one small uncentered hole, surrounded by five medium-sized holes, surrounded by five larger holes.
For b = 15 the FBMS is plotted in Figure 14. The FBMS has BC that are positioned with centers with triangular symmetry.
For b = 20 the FBMS has irregularly located BC; a figure can be found at [Oud20]. Interestingly, we obtainσ 1 ≈ 19.7076 for this surface, which is larger than the value obtained for the FBMS with BC at the vertices of a regular dodecahedron, as discussed in [GL20], with valueσ 1 ≈ 19.57189; see Table 4.3.

5.2.
Higher eigenvalues (j ≥ 2). Here, we consider the extremal Steklov eigenvalue problem (5), for higher eigenvalues,σ j , j ≥ 2. Less in known in this case and, in particular, the multiplicity of the optimal eigenvalue, and hence the dimension in which the FBMS exists, is unknown.
We recall from [FTY14] (see also Section 3) that by maximizing σ j for odd j among rotationally symmetric annuli yields an j+1 2 covering of the critical catenoid, a FBMS with b = 2 boundary components and j-th normalized Steklov eigenvalue, σ j = j + 1 2σ 1 , j odd.
We first consider b = 2 BC and eigenvalue j = 2. In this case, the density ρ on the outer boundary of the punctured disk becomes degenerate and resembles the ρ discussed in Example 2.5 and displayed in Figure 3. We believe that this ρ corresponds to the critical catenoid glued to a disc, but this is difficult to resolve using our numerical method; see Example 4.3. For other higher eigenvalues, we see similar phenomena for some initializations of ρ. However, there are a few values of eigenvalue number j and BC b, that give interesting local maximizers and are very robust with respect to the initialization.
For b = 3 BC the FBMS obtained by maximizing the j = 3 and j = 5 eigenvalues are displayed in Figure 18. If Figures 16 and 17, the first few eigenfunctions are plotted in the optimal domains, Ω c,r . The eigenvalues obtained areσ 3 = 23.6659 andσ 5 = 34.5317. Note that, again, these are local maximizers since larger eigenvalues can be obtained by gluing two or four balls to the surface attained by maximizing the first eigenvalue with b = 3 BC, to obtain eigenvaluesσ 3 = 12.0120 + 2 · 2 · π ≈ 24.5784 andσ 5 = 12.0120 + 2 · 4 · π ≈ 37.1447.

Discussion
In this paper, we developed computational methods to maximize the length-normalized j-th Steklov eigenvalue,σ j (Σ, g) := σ j (Σ, g)L(∂Σ, g) over the class of smooth Riemannian metrics, g on a compact surface, Σ, with genus γ and b boundary components. Our numerical method involves (i) using conformal uniformization of multiply connected domains to avoid explicit parameterization for the class of metrics, (ii) accurately solving a boundary-weighted Steklov eigenvalue problem in multi-connected domains, and (iii) developing gradient-based optimization methods for this nonsmooth eigenvalue optimization problem. Using the connection due to Fraser and Schoen [FS15], the solutions to this extremal Steklov eigenvalue problem for various values of b boundary components are used to generate free boundary minimal surfaces.
In hindsight, it may have been better to perform these computations on a punctured sphere rather than a punctured disk, as a punctured disk distinguishes one boundary (the 'outer' one). In particular, by considering a punctured sphere, it may be that the holes appear more symmetrically than for a punctured disk; see Figure 7.
Beyond further exploring higher eigenvalues j and higher numbers of boundary components b, there are a number of interesting extensions of this work. In particular, we would be very interested to compute extremal Steklov eigenvalues on the Möbius band, torus, and other higher genus surfaces and use the associated eigenfunctions to generate free boundary minimal surfaces. We're also interested in related extremal eigenvalue problems, involving convex combinations of Steklov eigenvalues or Steklov eigenvalues for the p-Laplacian.