The peak heat flux conjecture for the first Dirichlet eigenmode of convex planar domains

THE PEAK HEA T FLUX CONJECTURE F OR THE FIRST DIRICHLET EIGENMODE OF CONVEX PLANAR DOMAINS ZIJIAN W ANG, JEREMY G. HOSKINS, MANAS RACHH, AND ALEX H. BARNETT Abstract. In this paper, we study the scale-in v arian t quan tit y G (Ω) = ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) λ 1 , where u 1 is the first L 2 -normalized Dirichlet Laplace eigenfunction of a Euclidean domain Ω and λ 1 is its eigenv alue. This is related to the p eak b oundary heat flux in the long time limit. F or con v ex domains we pro v e that ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) is upp er-bounded by a (domain-independent) constant m ultiple of λ 1 . Using lay er potentials, we derive shap e-deriv ativ e formulae for effi- cient gradient computations. When combined with high-order Nystr¨ om discretization, a fast boundary integral equation solver, and eigen v alue ro otfinding, this allows us to n umerically op- timize G o v er a class of rounded polygonal discretized domains. Based on extensive numerical experiments, we then conjecture that, ov er the set of conv ex domains, G is maximized by the semidisk, with the p eak flux at the cen ter of the diameter. T o lend analytical supp ort to this conjecture, we prov e that the semidisk is a critical point of G under infinitesimal perturbations of its circular arc. 1. Introduction W e study the behavior of the Dirichlet eigenfunctions of a con v ex planar domain Ω, the non- trivial solutions to (1)  − ∆ u = λu, in Ω , u = 0 , on ∂ Ω , where ∆ is the t wo-dimensional Laplace operator. It is well known that if Ω is op en, connected, and bounded, then there exists an infinite sequence of Dirichlet eigenv alues 0 < λ 1 < λ 2 ≤ λ 3 ≤ . . . , with associated L 2 (Ω)-normalized eigenfunctions u 1 , u 2 , u 3 , · · · ∈ L 2 (Ω). Hassell and T ao show ed that the L 2 norms of their b oundary normal deriv ativ es ∂ n u 1 , ∂ n u 2 , . . . gro w lik e √ λ k as k → ∞ [23]. Suc h eigenfunctions (the “mo des of a drum”) arise in a wide range of ph ysical contexts, including heat flow, quan tum mechanics, Bro wnian motion, acoustics, optics, and structural engineering. How ev er, in many applications, the first eigenfunction, i.e., the one with the smallest eigen v alue, determines the long-time b eha vior of the system. This motiv ates the follo wing question: How big c an the normal derivative of u 1 b e on the b oundary of Ω , and what shap e, if any, maximizes it? Under rescaling of Ω b y a linear factor, sa y α , then ∂ n u 1 scales lik e α − 2 , and this is also true for λ 1 . Th us p erhaps the simplest quantit y that depends only on the shap e (not the size) of Ω is ∂ n u 1 /λ 1 . One could then ask for the solution, or even the existence of solutions, to the follo wing optimization problem: argmax Ω ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) λ 1 , (ZW) Department of Applied and Computational Mathematics, Y ale Universit y , New Hav en, CT 06511, USA (JH) Department of Statistics and CCAM, University of Chicago, IL 06037, USA (MR) Department of Mathematics, Indian Institute of T echnology Bom bay , Pow ai, Mumbai 400076, India (AB) Center for Computational Mathematics, Flatiron Institute, New Y ork, NY 10010, USA E-mail addresses: zijian.wang@yale.edu , jeremyhoskins@uchicago.edu , mrachh@iitb.ac.in , abarnett@ flatironinstitute.org Date: Marc h 18, 2026. 1 2 ZW, JH, MR, AND AB − 1 .0 − 0.5 0.0 0.5 1.0 0 . 0 0 . 5 1 . 0 (a) u 1 for Ω the s emidisk (0,0) 0.0 0.5 1.0 1.5 s 0 2 4 𝜕 n u 1 ( s ) / 𝜆 1 0. 0 0. 1 0. 2 0. 3 (b) scaled normal derivative of u 1 on 𝜕Ω Figure 1. (a) The L 2 (Ω)-normalized first Diric hlet eigenfunction u 1 of the unit semidisk, shown using the color scale to the righ t. (b) Its b oundary deriv ative plotted vs s , the counterclockwise arc-length from the b ottom-left corner ( − 1 , 0) of Ω in panel (a). The b oundary function is divided by λ 1 to make it scale- in v ariant. Its maximum ov er s , denoted b y G (Ω), is sho wn as a dot in (b), and o ccurs at s = 1 (corresp onding to the origin in panel (a)). Conjecture 3 is that there is no other conv ex domain Ω that exceeds this maxim um of C ∗ . where the optimization is ov er all b ounded conv ex domains Ω ⊂ R 2 . Remark 1 (A physical interpretation) . If U ( x, t ) satisfies the homo gen e ous he at e quation ∂ t U = ∆ U in Ω , with a zer o-temp er atur e c ondition on ∂ Ω , then in the long time limit U ( x, t ) ∼ c 1 e − λ 1 t u 1 ( x ) , so that the distribution of he at flux thr ough the b oundary is pr op ortional to ∂ n u 1 . The question is then: what shap e with a given de c ay r ate pumps the most flux thr ough some b oundary p oint? 1 This has a similar flavor to the r e c ently solve d [16] hot sp ots conjecture c onc erning the long-time b ehavior for c onvex domains inste ad with Neumann (insulating) b oundary c onditions. Despite its apparen t simplicit y , relatively little seems to be known about this problem. Here w e pro v e a partial analytic result and provide numerical evidence motiv ating a conjecture on the extremizing shap e. In particular, w e establish the following upp er b ound. Theorem 2. L et Ω b e a b ounde d c onvex domain in the plane. If λ 1 is its first Dirichlet eigenvalue and u 1 the c orr esp onding L 2 -normalize d eigenfunction then ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) ≤ C λ 1 , (2) wher e C is a c onstant indep endent of the domain. The following conjecture is based on extensive numerical exp erimen ts. Conjecture 3. If we define G (Ω) := ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) λ 1 , then sup Ω ∈D G (Ω) = max Ω ∈D G (Ω) = C ∗ , (3) wher e D is the sp ac e of al l b ounde d c onvex domains in the plane, and the maximum is c onje ctur e d to b e attaine d when Ω is a semidisk. The optimal c onstant C ∗ is given by C ∗ := 1 √ π j 1 , 1 | J 0 ( j 1 , 1 ) | ≈ 0 . 3655840228073865 , wher e J n is the Bessel function of the first kind of or der n , and j n,i is the i -th p ositive r o ot of J n . 1 T o make a precise physical interpretation, one would need a set-up in which c 1 is constant. Constant initial conditions U ( · , 0) ≡ 1, for example, give c 1 = ∥ ∂ n u 1 ∥ L 1 ( ∂ Ω) /λ 1 and hence a different optimization problem. PEAK HEA T FLUX CONJECTURE 3 Fig. 1 illustrates this. The conjecture is surprising for tw o reasons: 1) the extremal shap e has corners, in contrast to the extremizers for other kno wn lo w-lying eigen-quantities [2, 43, 44], and 2) the semidisk is also a shap e that ac hieves maximal gro wth in ∥ ∂ n u k ∥ L ∞ ( ∂ Ω) in the high-fr e quency limit k → ∞ , and at the same point, since the p oin t is self-fo c al for ra ys under billiard reflections in Ω. While this semiclassical ( k → ∞ ) concept of self-fo cusing is irrelev ant at k = 1, the coincidence is intriguing. Since ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) is a maxim um ov er b oundary points and is in v ariant under rigid motions of the domain, we may translate and rotate Ω so that a b oundary p oint at which ∂ n u 1 is maximized lies at the origin (0 , 0) and the inw ard unit normal points v ertically upw ard (as in Fig. 1). W e denote the set of all domains satisfying these prop erties as D 0 and define the follo wing function F : D 0 → R + b y F (Ω) := ∂ n u 1 (0 , 0) λ 1 , (4) where ∂ n u 1 ( · ) denotes the inw ard normal deriv ative of u 1 . W e can restate Conjecture 3 as sup Ω ∈D 0 F (Ω) = max Ω ∈D 0 F (Ω) = | √ π j 1 , 1 J 0 ( j 1 , 1 ) | − 1 . (5) As a first step to wards this conjecture, the follo wing theorem establishes that the semidisk is a critical point of the functional F if w e restrict the v ariations to the circular arc of the boundary . In particular, w e consider infinitesimal boundary perturbations whic h fix the bottom straigh t piece of the b oundary , and obtain the following result. Theorem 4. L et V : [0 , π ] → R b e a Lipschitz c ontinuous function with V (0) = V ( π ) = 0 , and define the p ar ametric domain Ω t := { ( r cos θ , r sin θ ) ∈ R 2 | 0 ≤ r < 1 + tV ( θ ) , 0 < θ < π } , wher e t is a r e al p ar ameter. Then, d dt     t =0 F (Ω t ) = 0 . (6) Remark 5 (Conv exity) . The r e ader may wonder why our c onje ctur e r estricts to c onvex domains. The answer is that the functional G (Ω) , e quivalently F (Ω) , is unb ounde d when nonc onvex domains ar e al lowe d. A n example is the se ctor of angle β > π , namely in p olars Ω β := { ( r , θ ) ∈ R 2 | r < 1 , 0 < θ < β } , which has a reentran t corner at the origin. We have u 1 ( r , θ ) = cJ ν ( j ν, 1 r ) sin ν θ , wher e ν = π /β < 1 , so using J ν ( z ) ∼ ( z / 2) ν / Γ( ν + 1) as z → 0 we get that the normal derivative along θ = 0 is ∂ n u 1 = r − 1 ∂ θ u 1 ∼ C r ν − 1 , as r → 0 , which is unb ounde d. As β → π + , Ω β b e c omes arbitr arily (Hausdorff ) close to the semidisk, but G (Ω β ) r emains infinite. The remainder of the pap er is organized as follows: In Section 2, we review relev ant b ounds on eigenfunctions and b oundary normal deriv atives, as w ell as connections to classical extremal- domain problems and the shap e-deriv ativ e viewp oin t. In Section 3, w e pro ve Theorem 2, then, after gathering some general results on shap e deriv ativ es of eigenfunctions, Theorem 4. In Section 4, w e present analytic results and qualitative examples for a simple domain class, motiv ating the conjectured extremizer. In Section 5, w e describ e the boundary in tegral formulation and the shap e-deriv ative calculations used to assemble gradients of the discretized ob jective. In Section 6, w e report n umerical experiments based on these gradients and discuss the observed extremizing shap es. W e conclude in Section 7 with a summary and directions for future work. 2. Back ground and rela ted w ork 2.1. Bounds on eigenfunctions and b oundary normal deriv atives. The study of Laplacian eigenfunctions and eigenv alues has a long history , reflecting their ubiquity in applications. Indeed, one could trace this history back at least to the 18th cen tury and the analysis of wa ves on a string. Since then there has b een a v ast b ody of w ork on asymptotic and non-asymptotic b eha vior for a wide v ariet y of b oundary conditions, both in Euclidean space and on manifolds. Of particular 4 ZW, JH, MR, AND AB in terest has b een the relationship to geometric features of the domain. Classic examples in v olv- ing the eigenv alues include W eyl’s law relating the growth of eigen v alues to the volume of the domain [64], and Kac’s celebrated pap er “Can y ou hear the shap e of a drum?” [31]. A wonderful review of the geometric effects on eigenfunctions is provided by Greb enk ov and Nguy en [20]. One vein of research of particular interest in applications is quan tifying v arious norms of eigen- functions, and the lo cations of their maxima or high-amplitude regions. Suc h work has a v ery differen t flav or for lo w-frequency eigenfunctions (such as u 1 ) than it do es in the high-frequency limit, where semiclassical and microlo cal analysis are often inv olved. Most b ounds on u 1 relate v arious domain norms. Pa yne–Rayner show ed that in R 2 , ∥ u 1 ∥ L 2 (Ω) ≤ √ λ 1 √ 4 π ∥ u 1 ∥ L 1 (Ω) with equality for a disk [45]. This was subsequently extended to higher dimensions by the same authors [46] and simplified b y Kohler-Jobin [34] with the in tro duction of an auxiliary problem. W ang et al. further generalize the R n ≥ 3 case to compact minimal surface with weakly connected Lipsc hitz b oundary [62]. W ork has been done to establish analogous b ounds for other op erators, e.g. p -Laplacians [1] and pseudo-Laplacians [41]. F or L ∞ norms, Pa yne and Stakgold show ed that for conv ex planar domains [48] π 2 | Ω | ∥ u 1 ∥ L 1 (Ω) ≤ ∥ u 1 ∥ L ∞ (Ω) . (7) Moreo ver, they prov ed a point wise b ound sho wing that the maxim um of u 1 dep ends on its location in Ω: u 1 ( x ) ≤ d ( x, ∂ Ω) √ λ 1 | Ω | ∥ u 1 ∥ L 1 (Ω) . (8) V an den Berg obtained an explicit inradius-based L ∞ b ound for L 2 (Ω)-normalized Diric hlet eigen- functions on b ounded, op en, connected domains in R d [61]. A complementary line of research studies the lo cation of maximizing p oin t for the solutions to Schr¨ odinger equations (which special- izes to Laplacian eigenfunctions when the p oten tial V is constant.) Rac hh and Steinerb erger [54] established that for simply connected planar domains, the distance from the maximizing p oin t to the b oundary is low er bounded b y ( | V | ∞ ) − 1 / 2 up to a universal constant factor, whic h translates to λ − 1 / 2 in the context of Laplacian eigenfunctions. This result is extended to fractional Schr¨ odinger equations with appropriate exp onen ts [6]. W e refer the reader to [20] for further extensions of these inequalities and related results. One conv enient tool for deriving such inequalities is Dirichlet Gr e en ’s function for the domain, whic h satisfies the PDE ( − ∆ x G ( x, y ) = δ ( x − y ) , x, y ∈ Ω , G ( x, y ) = 0 , x ∈ ∂ Ω . (9) Using the identit y u k ( x ) = λ k Z Ω G ( x, y ) u k ( y ) d y , (10) Moler and Pa yne [40] show ed, among other things, that for Diric hlet eigenfunctions | u k ( x ) | ≤ λ m ∥ u k ∥ L ∞ (Ω) w ( x ) , k = 1 , 2 , . . . (11) where w ( x ) = R Ω G ( x, y ) d y is the torsion function , a m uch studied quan tity in its o wn right dating bac k to the w ork of St V enan t. More recently , v an den Berg show ed that if Ω is an op en, b ounded and connected se t in d dimensions and u 1 is the L 2 (Ω)-normalized first Dirichlet eigenfunction, then ∥ u 1 ∥ L ∞ (Ω) ≤ 2 2 − d 2 π d/ 4 Γ( d/ 2)  j d − 2 2 , 1  d − 2 2    J d 2  j d − 2 2 , 1     ρ − d/ 2 , (12) PEAK HEA T FLUX CONJECTURE 5 where ρ is the inradius [61]. In the plane ( d = 2), this reduces to ∥ u 1 ∥ L ∞ (Ω) ≤ 1 √ π J 1 ( j 0 , 1 ) ρ − 1 . (13) W e now turn to norm b ounds that hold for all eigenfunctions. The main question is ho w suc h norms gro w with resp ect to λ K , as k → ∞ , and thus results tend to in v olve semiclassical or microlo cal analysis. The maximal b ound on the sup norm was pro v ed in 1952 by Levitan [37] (see also Av akumovic [5] and H¨ ormander [29]), stating that ∥ u k ∥ L ∞ (Ω) ≤ c (Ω) λ 1 / 4 k , k = 1 , 2 , . . . (14) for some constant c dep ending on Ω but indep endent of k . This pow er of 1 4 is sharp, being ac hieved for instance by the disk: the radially-symmetric subsequence u 0 ,m ( r , θ ) = J 0 ( j 0 ,m r ) / √ π J 1 ( j 0 ,m ), m = 1 , 2 , . . . , has a peak at r = 0 growing as p j 0 ,m from large-argumen t asymptotics of J 1 , and one recalls λ 0 ,m = j 2 0 ,m . Note that all ra ys launched from the disk cen ter return, and at equal times: it is a self-fo cus. Indeed, in 2002, Sogge–Zelditc h [57] show ed that this maximal gro wth can only happ en with a self-fo cus. F or b oundary deriv atives of eigenfunctions, for ∂ Ω smo oth, Hassell–T ao [23] show ed the L 2 -norm b eha vior c (Ω) λ 1 / 2 k ≤ ∥ ∂ n u k ∥ L 2 ( ∂ Ω) ≤ C (Ω) λ 1 / 2 k , k = 1 , 2 , . . . . (15) One w ould exp ect dimensionally that the sup norms of the b oundary deriv atives hav e one p o wer of √ λ k more than the sup norms of u k themselv es in (14), and indeed the growth b ound in 2D, ∥ ∂ n u k ∥ L ∞ ( ∂ Ω) ≤ C (Ω) λ 3 / 4 k , k = 1 , 2 , . . . , (16) is kno wn to b e sharp for C ∞ domains, a result summarized in [65, Cor. 2.1] (although the smo oth- ness is not made explicit there), and prov en for C ∞ manifolds with ∂ Ω conca v e (this excludes planar domains) more recen tly by Sogge–Zelditc h [58, Cor. 1.2]. Zelditch [65, § 5.4] men tions that for planar domains (16) it is ac hieved for the semidisk or any (presumably smaller) sector of the disk. F urthermore, paralleling results for the interior, this sharp growth is known only to b e ac hieved if there is a self-fo cal b oundary p oin t [58, 65]. In this paper, we fo cus on the first eigenfunction and on the b oundary quantit y ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) , aiming to understand its largest p ossible size (after the scale-inv arian t normalization by λ 1 ). Note that this problem is closely related in spirit to w ork on the torsion function and gradient b ounds at “failure points.” In [30], numerical metho ds were used to prop ose candidates for extremizing do- mains in that con text. Interestingly , for that problem the shap e of the extremizer, if an extremizer exists at all, is still an op en question. 2.2. Extremal domains and eigen v alue problems. Extremal domains related to sp ectral prop erties under fixed constraints hav e been w ell-explored due to their impact in physical science and engineering. One of the most famous example is Lord Rayleigh’s conjecture that the disk minimizes the first Diric hlet Laplacian eigenv alue within all planar domains of a fixed area [55]. These isop erimetric problems are known for b eing easy to state yet hard to prov e [24]. Almost 30 y ears later, this was prov ed by F ab er and Krahn [35] using a Sch warz rearrangement argument. Recen tly , it has b een strengthened to sharp quantitativ e b ounds via F raenk el asymmetry [11, 12]. Numerous w orks ha ve inv estigated extremal domains in a v ariety of settings. One direction is to v ary the underlying PDE problem. F or Neumann boundary conditions, the first Diric hlet eigen v alue is zero with the constant eigenfunction. Szeg¨ o in t w o dimensions and W einberger in n dimensions prov ed that the ball is the minimizing shap e for the second Neumann eigenv alue [60, 63]. F or Robin b oundary conditions with p ositiv e b oundary parameter, Bossel sho wed that the disk minimizes the first Diric hlet eigen v alue in R 2 and Daners extended it in to higher dimensions [9, 15]. F or negativ e boundary parameters, F reitas et al. prov ed that the ball is a maximizer of the first eigen v alue when the boundary parameter is not to o negativ e in R 2 [19]. Moreov er, they sho wed that for sufficiently negative b oundary parameter and general R n , the ball is not the maximizer. This serves as the one of v ery first examples where the extremal domain is not a ball, in contrast to the ma jority of earlier F ab er–Krahn type results. 6 ZW, JH, MR, AND AB Another direction of w ork is where the constrain ts for the admissible shap es are v aried. If we only consider p olygons, P´ oly a and Szeg¨ o pro v ed that regular n -gon minimizes the first Dirichlet Laplacian eigen v alue when n = 3 , 4 using a Steiner’s symmetrization, which fails at n = 5 [51]. The case when n = 5 remains an op en problem. Going beyond simply connected domains, Pa yne et al. show ed that the ann ulus is the extremal domain within doubly connected domains, with Diric hlet/Neumann conditions on the outer/inner b oundary resp ectiv ely [47]. Finally , one could v ary the ob jective. P´ oly a show ed, among other things, that the union of tw o disks minimizes the second Dirichlet eigen v alue [50] under fixed area. When a con vexit y constrain t is added, it is known that the optimal shap e con tains t wo straight lines, but is not a stadium, although it is numerically close to one (see Oudet [44], who also minimized λ k at fixed area, for k = 3 , 4 , . . . , 10). Kennedy extended P´ olya’s result to Robin eigenv alues with positive b oundary parameters [33]. Rather than optimizing a sp ecific eigen v alue, Ashbaugh et al. prov ed that the ball maximizes the r atio of the second to first Diric hlet eigen v alue [3]. Numerical optimization studies (using smo oth radial F ourier series represen tations, and where the eigen v alue problem was solv ed by particular [43] or fundamental solution [2] metho ds) hav e found the planar domains that extremize a v ariety of such small-index sp ectral ratios; the domains seem to be smooth. Inspired b y the recently established gradien t bound of the torsion function, as well as the nontrivial (conjectured) extremal domain for that problem [30], here w e study in the maximum magnitude of the boundary normal deriv ative of the first Diric hlet eigenfunction. Rather than enforcing an area constraint, we normalize by the first eigenv alue to mak e our ob jective scale-inv ariant. 2.3. Shap e deriv ativ es via la yer p oten tials. Boundary integral equations provide an efficient framew ork for computing Laplacian eigenfunctions and th us ob jective functions like the one in (4), see [66], for example. La yer potentials are fundamen tal building blo c ks for computing eigenfunc- tions using boundary integral equations [14]. How ever, optimizing methods for (4) tend to require gradien t information for faster conv ergence, w hic h in turn requires F r ´ echet deriv atives of lay er p oten tials. In [52], Potthast established F r ´ echet differentiabilit y of b oundary integral operators in the con text of in v erse acoustic scattering, which pro vides theoretical supp ort as w ell as to ols for conv ergence analysis of iterated metho ds [27, 28]. In particular, P otthast show ed the differen- tiabilit y of single and double lay er p oten tials as they appear in the Brakhage-W erner form ulation of the Helmholtz equation [10, 14]. Subsequently , the he also extended the results to Neumann b oundary conditions [53] the results of whic h are relev ant for ev aluating shap e deriv atives of (4). 3. Proofs of Theorem 2 and Theorem 4 The pro of of Theorem 2 is a straightforw ard application of the Green’s function metho d, com- bined with an estimate on the normal deriv ative of the torsion function in con v ex domains. Pr o of of The or em 2. W e begin by observing that if G is the Dirichlet Green’s function for Ω defined in (9), then u 1 ( x ) = λ 1 Z Ω G ( x, x ′ ) u 1 ( x ′ ) dx ′ . (17) Up on taking the gradient of b oth sides, we obtain ∇ u 1 ( x ) = λ 1 Z Ω ∇ x G ( x, x ′ ) u 1 ( x ′ ) dx ′ . (18) Without loss of generality , let u 1 ≥ 0 . F urthermore, for x ∈ ∂ Ω − n ( x ) · ∇ u 1 ( x ) ≥ 0 , − n ( x ) · ∇ G ( x, x ′ ) ≥ 0 , (19) where n ( x ) is the outw ard p oin ting normal. It follows that ∥ n · ∇ u 1 ∥ L ∞ ( ∂ Ω) ≤ λ 1 ∥ u 1 ∥ L ∞ (Ω) max x ∈ ∂ Ω Z Ω n ( x ) · ∇ x G ( x, x ′ ) dx ′ . (20) PEAK HEA T FLUX CONJECTURE 7 The in tegral on the righ t-hand side of the previous expression is the normal deriv ative of the torsion function w for Ω , which itself satisfies the following PDE ( ∆ w = − 1 , in Ω , w = 0 , on ∂ Ω . (21) An inequalit y of Sp erb (see Keady and McNabb as w ell) [32, 59] gives, for con vex domains, the follo wing b ound for the maxim um of the gradient of the torsion function, ∥∇ w ∥ L ∞ (Ω) ≤ ρ, (22) where ρ once again denotes the inradius of Ω. Up on substitution of this b ound into the previous inequalit y , and using v an den Berg’s inequalit y (13), w e obtain ∥ n · ∇ u 1 ∥ L ∞ ( ∂ Ω) ≤ 1 √ π J 1 ( j 0 , 1 ) λ 1 , (23) whic h completes the proof. □ W e note that C from the ab o ve proof is C ≈ 1 . 08676163613127, and thus about three times larger than the optimal C ∗ discussed in the introduction. 3.1. Rate of c hange of eigenfunction normal deriv ative data with resp ect to shap e deformations. T o address Theorem 4, we first need t wo lemmas for the rate of c hange of the normal deriv ative data of the i th Dirichlet eigenmo de u i of a domain undergoing a general b ound- ary deformation, which may be of independent in terest. In Section 3.2, these will then b e applied to the mo de i = 1 of the semidisk to pro v e the theorem. Let Ω t ⊂ R 2 b e a domain with piecewise smo oth b oundary , whose shap e changes with a parameter t ∈ R . In particular, let V be the out ward normal “v elo cit y” on ∂ Ω with resp ect to c hanging t , at t = 0, as in Theorem 4. Let u i b e the i th normalized Dirichlet eigenmo de of Ω t , as in (1). The classical Hadamard form ula for the eigen v alue rate of c hange is (24) ˙ λ i = − Z ∂ Ω V ( ∂ n u i ) 2 ds, where we use the standard notation that a dot ab o ve a quantit y denotes ∂ t , i.e., its partial with resp ect to the deformation parameter. Here ds is the arc-length elemen t on ∂ Ω. F or this, see Grinfeld [22]. He derives formulae (his Sec. 1.4) for eigenfunction c hanges only in cases where the resulting formula is lo cal on the b oundary; ho w ever, the normal-deriv ativ e data case that we need b elo w is not cov ered. Recall the calculus of moving surfaces (e.g., see Grinfeld [22, § 3.2]). F or a function f defined in the closure of Ω there are tw o t yp es of time-deriv ative at x ∈ ∂ Ω: the total deriv ativ e D t f ( x ) := ( d/dt ) f ( x ( t )) in which x ( t ) mov es with the b oundary (this is often written D/D t , or δ /δ t in Grinfeld), and the usual partial ∂ t f ( x ) = ˙ f ( x ) in whic h the spatial lo cation x is held fixed as t c hanges. They are related by the c hain rule [22, (42)]: (25) D t f = ˙ f + V ∂ n f on ∂ Ω . As an application, if u i is the ( t -dep enden t) Dirichlet eigenfunction of the deforming domain, then u i | ∂ Ω = 0 for all t , so that (25) gives (see [22, (70)]) the partial deriv ativ e (26) ˙ u i = − V ∂ n u i on ∂ Ω . The Thomas rule (e.g. see Lemma 10) gives the rate of c hange of the unit normal vector n ( x ) at a p oin t x ∈ ∂ Ω moving with the normal deformation V , due to the lo cal rotation rate of the surface, (27) D t n = −∇ ∂ Ω V , where ∇ ∂ Ω is the surface gradient op erator, equal to τ ∂ s in our case of d = 2 dimensions, τ b eing the unit tangen t v ector in the increasing s sense. W e will also need the lo cal mean curv ature of ∂ Ω at t = 0, namely κ = 1 /R , where R is the signed lo cal radius of curv ature (the sign is p ositiv e for the circle; note that [22, p.664] uses the opp osite sign). 8 ZW, JH, MR, AND AB Lemma 6. L et u i b e a Dirichlet eigenfunction of a domain Ω whose t -dep endent shap e change is sp e cifie d by the b oundary normal velo city function V . Then the total and p artial derivatives of its b oundary data ar e r elate d by (28) D t ( ∂ n u i ) = ∂ n ˙ u i − V κ∂ n u i on ∂ Ω . Pr o of. W e follow [22, Sec. 4.2] but for the Dirichlet case. The pro duct rule for the total deriv ativ e giv es D t ( ∂ n u i ) = D t ( n · ∇ u i ) | ∂ Ω = n · D t ∇ u i + ( D t n ) · ∇ u i . Applying (25) to each comp onen t in the first D t , and (27) to the second D t , gives (29) D t ( ∂ n u i ) = n · ∇ ˙ u i + V n · ∇∇ u · n − ∇ ∂ Ω V · ∇ u i , where ∇∇ denotes the Hessian matrix. F or a general Helmholtz solution (∆ + λ ) u = 0 in R d w e ha ve [22, (74)] n · ∇∇ u · n = − λu − ∆ ∂ Ω u − κ∂ n u, on ∂ Ω , where ∆ ∂ Ω is the surface Laplace-Beltrami op erator (noting the sign change from Grinfeld). This can b e prov ed easily in the d = 2 case via lo cal p olar co ordinates. Since u i additionally has zero Diric hlet b oundary condition this simplifies to n · ∇∇ u i · n = − κ∂ n u i . Applying this to (29), and noting that ∇ u i is purely normal, gives (28). □ W e next show that ˙ u i solv es a BVP in the static (fixed) domain Ω. Lemma 7. L et u i b e a simple normalize d Dirichlet eigenfunction of the b ounde d domain Ω under- going shap e change with b oundary velo city function V , as ab ove. Then its r ate of change ψ := ˙ u i is the unique solution to the c onstr aine d inhomo gene ous Dirichlet Helmholtz interior BVP (∆ + λ i ) ψ = − ˙ λ i u i =: f in Ω (30) ψ = − V ∂ n u i =: g on ∂ Ω (31) ( ψ , u i ) L 2 (Ω) = 0 . (32) Pr o of. The PDE comes from the taking ∂ t of (∆ + λ i ) u i = 0. The b oundary condition is simply (26). By themselves this pair (30)–(31) is an on-resonance (secular) inhomogeneous Diric hlet Helmholtz BVP , whose null-space is th us the span of u i . Thus a solution ψ exists if and only if the data f and g satisfy a standard compatibility condition (see, e.g., [39, Thm. 4.10]). This condition comes from Green’s iden tity applied to the data and the n ull-space vector u i (i.e., solution to the homogeneous adjoint BVP), Z Ω [ ψ (∆ + λ i ) u i − u i (∆ + λ i ) ψ ] dx = Z ∂ Ω ( ψ ∂ n u i − u i ∂ n ψ ) ds. Applying (30), using (∆ + λ i ) u i = 0, and using that u i = 0 on ∂ Ω, the condition is ( u i , f ) L 2 (Ω) + Z ∂ Ω g ∂ n u i ds = 0 . F or f and g defined in (30)–(31), since u i is normalized, the condition holds by the Hadamard form ula (24). This assures existence of a (nonunique) solution to the BVP (30)–(31). The or- thogonalit y condition (32), whic h mak es ψ unique, comes from applying the product rule to the t -deriv ative of the normalization R Ω u 2 i dx = 1 to give ( u i , ˙ u i ) Ω = 0. (See [22, (68)]; a subtlet y is that the t - dependent fundamen tal theorem of calculus for volume in tegrals is needed, which in the Diric hlet case do es not in tro duce an y new terms.) □ PEAK HEA T FLUX CONJECTURE 9 Finally , it is conv enient to split the solution of (30)–(32) into a particular solution ˜ ψ plus a homogeneous solution w . Namely , if ˜ ψ solves (33) (∆ + λ i ) ˜ ψ = − ˙ λ i u i in Ω with arbitrary b oundary condition, and then w is the unique solution to (∆ + λ i ) w = 0 in Ω (34) w = − V ∂ n u i − ˜ ψ on ∂ Ω (35) ( w , u i ) L 2 (Ω) = − ( ˜ ψ , u i ) L 2 (Ω) , (36) then ψ = ˜ ψ + w solves (30)–(32), and is th us the desired eigenfunction partial ˙ u i . 3.2. Application to u 1 for the semidisk: pro of of Theorem 4. The first Dirichlet eigen- function of the semidisk 0 < r < 1, 0 < θ < π , is, in p olar co ordinates, (37) u 1 ( r , θ ) = − 2 √ π J ′ 1 ( k 1 ) J 1 ( k 1 r ) sin θ , where the normalization constant (negativ e so that u 1 is nonnegative) follo ws from, eg, [42, (10.22.37)], and its eigenv alue is λ 1 = k 2 1 = j 2 1 , 1 . Prop osition 8. L et Ω b e the semidisk, u 1 and λ 1 b e as ab ove, then for any ˙ λ 1 ∈ R , (38) ˜ ψ ( r , θ ) = − ˙ λ 1 √ π k 1 J ′ 1 ( k 1 ) r J 0 ( k 1 r ) sin θ satisfies the i = 1 c ase of (33) , i.e., is a p articular solution. Pr o of. The Laplacian in p olar co ordinates is ∆ = 1 r ∂ r ( r ∂ r · ) + 1 r 2 ∂ θθ . By separation of v ariables, ˜ ψ ( r , θ ) = U ( r ) sin θ where U solv es the inhomogeneous radial ODE (39) L 1 [ U ]( r ) := 1 r ( r U ′ ) ′ − 1 r 2 U ( r ) + λ 1 U ( r ) = 2 ˙ λ 1 √ π J ′ 1 ( k 1 ) J 1 ( k 1 r ) , and v anishes at r = 0 so that ˜ ψ is regular at the origin. By direct differentiation, L 1 [ r J 0 ( k 1 r )] = 2 ∂ r J 0 ( k 1 r ) = − 2 k 1 J 1 ( k 1 r ), so that U ( r ) = − ˙ λ 1 / ( √ π k 1 J ′ 1 ( k 1 )) · r J 0 ( k 1 r ) satisfies (39). □ W e now specialize to a b oundary normal deformation function V ( θ ) on the semidisk, and which v anishes on the straigh t part of ∂ Ω, as in Theorem 4. In this case, since ∂ n u 1 ( θ ) = − (2 k 1 / √ π ) sin θ on the semidisk r = 1, and ds = dθ , Hadamard (24) b ecomes (40) ˙ λ 1 = − 4 k 2 1 π Z π 0 V ( θ ) sin 2 θ dθ , whic h ma y b e inserted into (38) to define ˜ ψ . It then only remains to solve the homogeneous on-resonance BVP (34)–(36); for this we use what is (p erhaps confusingly) known as the metho d of p articular solutions . Prop osition 9. With V an arbitr ary smo oth function on 0 < θ < π , r e c al ling (40) and ˜ ψ given by (38) , the unique i = 1 homo gene ous solution to (34) – (36) takes the form of the F ourier-Bessel sine series (41) w ( r , θ ) = ∞ X ℓ =1 c ℓ J ℓ ( k 1 r ) sin ℓθ, for some c o efficients c ℓ ∈ R , and wher e, in p articular, c 1 = 0 . 10 ZW, JH, MR, AND AB Pr o of. Each term in (41) solv es the Helmholtz equation (34) in the semidisk with zero b oundary condition on the straight b ottom b oundary . Thus each c ℓ , ℓ = 2 , 3 , . . . is uniquely determined b y matc hing F ourier sine series coefficients of the b oundary data on the righ t-hand side of (35). Uniqueness here holds because J ℓ ( k 1 )  = 0, ℓ = 2 , 3 , . . . ; see [42, § 10.21(i)]. The formula for these c ℓ are not hard to write, but are not needed. Y et, b ecause J 1 ( k 1 ) = 0, c 1 is undetermined by (35), and is instead fixed b y (36). Inserting (38) and (37), dropping prefactors for simplicity , one sees that the right-hand side of (36) is a m ultiple of Z π 0 Z 1 0  r J 0 ( k 1 r ) sin θ · J 1 ( k 1 r ) sin θ  r drdθ = π 2 Z 1 0 r 2 J 0 ( k 1 r ) J 1 ( k 1 r ) dr = 0 , using a Bessel function identit y [42, (10.22.7)] with µ = 0 and ν = 1. Th us ( w , u 1 ) L 2 (Ω) = 0, so c 1 = 0. □ W e finally apply all of the ab o ve results. Let 0 denote the origin of the semidisk; it is easy to c heck using (37) that this p oint is the global maximum of ∂ n u 1 on ∂ Ω. Since there is no surface deformation at the origin, D t ∂ n u i ( 0 ) = ∂ n ˙ u i ( 0 ), by Lemma 6. Then by the quotient rule, to show d F /dt = 0 as claimed in Theorem 4, it is sufficient to show the v anishing of (42) λ 1 ∂ n ˙ u 1 ( 0 ) − ˙ λ 1 ∂ n u 1 ( 0 ) , for an arbitrary deformation V on the semidisk. How ever, by Lemma 7, ∂ n ˙ u 1 ( 0 ) = ∂ n ψ ( 0 ) = ∂ n ˜ ψ ( 0 ) + ∂ n w ( 0 ) = − ˙ λ 1 √ π k 1 J ′ 1 ( k 1 ) + 0 where we used Prop. 8 for the first term, and see that the second term v anishes b y Prop. 9, since J ′ ℓ (0) = 0, ℓ = 2 , 3 , . . . while c 1 = 0. Y et using (37) and J 0 (0) = 1 shows that the second term of (42) is − ˙ λ 1 · − 2 k 1 √ π J ′ 1 ( k 1 ) . Recalling λ 1 = k 2 1 sho ws that the t wo terms in (42) cancel, completing the pro of of the theorem. 4. Anal ytic resul ts for rect angular domains Before in tro ducing the n umerical supp ort for the conjectured extremizer, we first examine F in (3) (and the resulting optimal constant) when D is restricted to rectangles. The purp ose of this brief analysis is tw o-fold: a) firstly to illustrate the inv ariance of F under scaling, and b) to demonstrate that the most symmetric shape in the class do es not optimize the function, thus hin ting that the circle need not necessarily optimize F o ver all conv ex domains. Let R a,b b e a rectangle with side lengths a ≥ b > 0. Recall that its first eigenfunction and eigen v alue are given b y u ( x, y ) = sin π x a sin π y b , λ = π 2 ( a − 2 + b − 2 ) . Direct computation shows that the Euclidean norm of its gradien t ac hiev es maxim um at b oundary: (43) |∇ u | 2 = π 2  a − 2 cos 2 π x a sin 2 π y b + b − 2 sin 2 π x a cos 2 π y b  = π 2 a 2 b 2  a 2 X − ( a 2 + b 2 ) X Y + b 2 Y  , where X := sin 2 π x a and Y := sin 2 π y b tak e v alues in [0 , 1]. This, together with the assumption that a ≥ b , implies that ∥∇ u ∥ L ∞ ( ∂ Ω) = π b . Another straightforw ard computation sho ws that ∥ u ∥ 2 L 2 (Ω) = Z a 0 Z b 0 sin 2 ( π x/a ) sin 2 ( π y /b ) dydx = ab 4 . PEAK HEA T FLUX CONJECTURE 11 Com bining the ab o ve computations, w e hav e ϕ ( R a,b ) = π b π 2 ( a − 2 + b − 2 ) q ab 4 = 2 π  ( a − 2 + b − 2 ) a 1 / 2 b 3 / 2  − 1 . By comparing the exponents, w e see that ϕ ( R ra,r b ) = ϕ ( R a,b ) for an y scaling factor r > 0. Let α := a b ≥ 1 b e the side ratio of the rectangle. Using scale-inv ariance, we obtain ϕ ( R a,b ) = ϕ ( R α, 1 ) = 2 π 1 ( α − 2 + 1) α 1 / 2 = 2 π α 3 / 2 α 2 + 1 . Finally , it follows that α = √ 3 is the optimal rectangle, for which F ≈ 0 . 362794816. 5. Shape deriv a tives and gradient comput a tion This section describ es an integral equation based approach for ev aluating F and its shap e deriv atives. Our starting p oin t is a b oundary integral characterization of Diric hlet eigenpairs (a related approach is describ ed in [66]), follow ed by the computation of the shap e deriv atives of F extending the results from [27, 52]. 5.1. Boundary integral equations and lay er potentials. Let Ω ⊂ R 2 b e a bounded C 2 domain with b oundary ∂ Ω. Let k > 0 and supp ose u  = 0 solv es the Helmholtz equation (∆ + k 2 ) u = 0 in Ω , u = 0 on ∂ Ω , (44) so that k 2 is a Diric hlet eigen v alue. Let x, y ∈ ∂ Ω, let dS denote arc-length measure, and let n ( x ) denote the outw ard unit normal. The free-space Green’s function for (∆ + k 2 ) in R 2 is given by G k ( x, y ) := i 4 H (1) 0 ( k | x − y | ) , (45) where H 1 0 ( z ) is the Hank el function of the first kind of order 0. Consider the following lay er p oten tial op erators defined on the b oundary ∂ Ω (46) S k [ σ ]( x ) := Z ∂ Ω G k ( x, y ) σ ( y ) ds y , S ′ k [ σ ]( x ) := p . v . Z ∂ Ω ∂ n x G k ( x, y ) σ ( y ) ds y , S ′′ k [ σ ]( x ) := f . p . Z ∂ Ω ∂ 2 n x G k ( x, y ) σ ( y ) ds y , D k [ σ ]( x ) := p . v . Z ∂ Ω ∂ n y G k ( x, y ) σ ( y ) ds y , D ′ k [ σ ]( x ) := f . p . Z ∂ Ω ∂ n x ∂ n y G k ( x, y ) σ ( y ) ds y , with p . v . denoting Cauch y principal v alue and f . p . denoting Hadamard finite part [39]. In a sligh t abuse of notation, we use S k [ σ ]( x ) to also denote the operator m apping σ on the b oundary ∂ Ω to a solution of the Helmholtz equation in Ω. Using Green’s iden tities, the Dirichlet eigenfunction u can be represented via the single-lay er p oten tial S k with density σ := − ∂ n u | ∂ Ω , i.e. u = S k [ ∂ n u | ∂ Ω ]. Using standard jump conditions [14] and taking the limit (47) σ ( x ) = ∂ n u ( x ) = lim h → 0 + n ( x ) · ∇ u ( x − hn ( x )) , w e get that σ satisfies the equation 1 2 σ ( x ) − S ′ k [ σ ]( x ) = 0 , x ∈ ∂ Ω . (48) F ollowing the pro cedure in [66], the Diric hlet eigenv alue can b e obtained by solving the nonlinear eigen v alue problem of determining k suc h that I / 2 − S ′ k has a nontrivial n ull space. 12 ZW, JH, MR, AND AB 5.2. Curv e parametrization and normal p erturbations. Let γ : [0 , L ] → R 2 parameterize ∂ Ω and let τ denote the unit tangent, n the outw ard unit normal, and κ the signed curv ature. In this section w e denote the smo oth scalar deformation field b y ν on ∂ Ω, consider the normal p erturbation γ ε ( s ) := γ ( s ) + εn ( s ) ν ( s ) , (49) th us ε denotes what was called t in Sec. 3.1. W e denote by δ ( · ) the deriv ative at ε = 0, i.e. δ f := ∂ ε | ε =0 f . (50) Let s b e an arbitrary smo oth parameter on ∂ Ω, so that dS = | ˙ γ ( s ) | ds, τ := ˙ γ | ˙ γ | . (51) F or simplicity , for functions and op erators defined on ∂ Ω (such as la yer p oten tials, σ, n, τ , κ, ν ), we write f ( s ) to denote their v alues at p oin t γ ( s ). F or a scalar function ν on ∂ Ω, w e write dν dτ := ∂ τ ν for its tangential deriv ative. W e will use the following standard geometric identities (see, e.g., [17] for a more comprehensive treatment). Lemma 10. Under the normal p erturb ation γ ε = γ + εnν , we have δ ( | ˙ γ | ) = κν | ˙ γ | , and δ n = − τ dν dτ (Thomas rule) . (52) Pr o of. Let R denote clo c kwise rotation b y π / 2 and define the unit normal b y n := | ˙ γ | − 1 R ˙ γ (so Rτ = n and Rn = − τ ). Differen tiating γ ε = γ + εnν with respect to s giv es ∂ s γ ε = ˙ γ + ε ( ˙ n ν + n ˙ ν ), so δ | ˙ γ | = ˙ γ | ˙ γ | · ( ˙ n ν + n ˙ ν ) = τ · ˙ n ν = κν | ˙ γ | . (53) F or the normal v ariation, differen tiating n = | ˙ γ | − 1 R ˙ γ yields δ n = − δ | ˙ γ | | ˙ γ | 2 R ˙ γ + 1 | ˙ γ | R ( ˙ n ν + n ˙ ν ) = − τ dν dτ , (54) using R ˙ γ = | ˙ γ | n , R ˙ n = κ | ˙ γ | n , and dν dτ = ˙ ν / | ˙ γ | . □ 5.3. Shap e deriv ative and v ariations of la y er p otentials. Let x ∗ = (0 , 0) ∈ ∂ Ω denote the b oundary p oin t at which we ev aluate the normal deriv ative used for the ob jective function F , as defined in (4). The b oundary in tegral equation (48) determines σ (and hence u ) only up to scale. Th us, we express the functional in terms of the unnormalized eigenfunction as F (Ω) = ∂ n u ( x ∗ ) k 2 p N (Ω) = σ ( x ∗ ) k 2 p N (Ω) , (55) where N (Ω) := ∥ u ∥ 2 L 2 (Ω) . Since w e only consider p erturbations fixing x ∗ and the lo cal b oundary orien tation, we hav e δ ( σ ( x ∗ )) = ( δ σ )( x ∗ ). Differentiating (55) gives δ F = ( δ σ )( x ∗ ) k 2 √ N − 2( δ k ) σ ( x ∗ ) k 3 √ N − σ ( x ∗ ) δ N 2 k 2 N 3 / 2 . (56) Th us we require δk , δσ , and δ N . The v ariations δ k and δ σ follow from differentiating the boundary in tegral equation (48). Since S ′ k dep ends on b oth the b oundary and the wa ven umber k , differenti- ating pro duces t wo operator v ariations: the shap e deriv ative δ ( S ′ k ) and the deriv ative with respect to the wa ven umber k denoted by ( δ k ) ∂ k S ′ k . 0 =  1 2 I − S ′ k  [ δ σ ] − δ ( S ′ k )[ σ ] − ( δ k ) ( ∂ k S ′ k )[ σ ] . (57) F ormulas for δ ( S ′ k ) and ∂ k S ′ k are shown in the lemma b elo w and can b e obtained via standard calculations [38]. PEAK HEA T FLUX CONJECTURE 13 Lemma 11. L et ν b e smo oth on ∂ Ω and c onsider the p erturb ation γ ε = γ + εnν . Defining S τ ,k [ σ ]( x ) := p . v . Z ∂ Ω ∂ τ x G k ( x, y ) σ ( y ) ds y , (58) we have δ ( S ′ k ) [ σ ] = ν S ′′ k [ σ ] + D ′ k [ ν σ ] − dν dτ S τ ,k [ σ ] + S ′ k [ κν σ ] (59) ∂ k S ′ k [ σ ] = ik 4 Z ∂ Ω H 1 0 ( k | x − y | )( x − y ) · n ( x ) σ ( y ) ds y . (60) Pr o of. W e first compute ∂ k S ′ k . Straightforw ard differentiation of the kernel giv es us ∂ n x G k ( x, y ) = ik 4 H 1 1 ( k | x − y | ) ( x − y ) · n x | x − y | . (61) Using identit y [42, Eq. 10.6.6], w e hav e ∂ k  k H 1 1 ( k | x − y | )  = | x − y | H 1 0 ( k | x − y | ) . (62) Substituting this into (61) gives the desired result. Next, we compute δ ( S ′ k ) [ σ ]. Fix a smooth parametrization γ : [0 , L ] → ∂ Ω. F or simplicity , w e write S ′ k [ σ ]( t ) for the op erator ev aluated at the target p oin t γ ( t ) (and similarly for D ′ k , S ′′ k , S τ ,k ). The p erturbed k ernel is S ′ k,ε [ σ ]( t ) := p . v . Z L 0 n ε ( t ) · ∇ x G k  γ ε ( t ) , γ ε ( s )  σ ( s ) | ˙ γ ε ( s ) | ds. (63) Straigh tforward computation gives us (64) δ ( S ′ k ) [ σ ]( t ) = Z L 0  ν ( t ) ∂ 2 n x G k ( γ ( t ) , γ ( s )) + ν ( s ) ∂ n x ∂ n y G k ( γ ( t ) , γ ( s ))  σ ( s ) | ˙ γ ( s ) | ds + Z L 0  δ n ( t ) · ∇ x G k ( γ ( t ) , γ ( s ))  σ ( s ) | ˙ γ ( s ) | ds + Z L 0 ∂ n x G k ( γ ( t ) , γ ( s )) σ ( s ) δ | ˙ γ ( s ) | ds, Applying Lemma 10, the last t w o integrals b ecome (65) Z L 0  δ n ( t ) · ∇ x G k ( γ ( t ) , γ ( s ))  σ ( s ) | ˙ γ ( s ) | ds = − dν dτ ( t ) S τ ,k [ σ ]( t ) Z L 0 ∂ n x G k ( γ ( t ) , γ ( s )) σ ( s ) δ | ˙ γ ( s ) | ds = S ′ k [ κν σ ] Rewriting the first integral in (64) in terms of S ′′ k , D ′ k yields (59). □ Although S ′′ k and D ′ k are h yp ersingular, the com bination in (59) has a kernel that is bounded on the diagonal. T o see this, we first separate the leading order singularity of the Helmholtz Green’s function by writing: G k ( x, y ) = − 1 2 π log ρ + k 2 8 π Q k ( ρ ) + R k ( ρ ) , Q k ( ρ ) := ρ 2 log ρ, ρ := | x − y | (66) where R k is C 2 . Fix the target parameter t and set r := t − s, γ r := γ ( t ) − γ ( s ) , ν r := ν ( t ) − ν ( s ) . Then as r → 0, γ r = r ˙ γ ( t ) + O ( r 2 ) , | γ r | 2 = r 2 | ˙ γ ( t ) | 2 + O ( r 3 ) , ν r = r ˙ ν ( t ) + O ( r 2 ) , (67) and ˙ ν ( t ) = | ˙ γ ( t ) | dν dτ ( t ) b y definition of the tangential deriv ative. Moreo ver, since γ r is tangen t to first order, n ( t ) · γ r = O ( r 2 ), and hence ∂ n x log | γ r | = O (1) . 14 ZW, JH, MR, AND AB The leading order singularit y in (59) arises from differen tiating the singular part − 1 2 π log | x − y | t wice in normal directions, whic h produces a | γ r | − 2 singularit y . A t the level of leading terms, one has ∂ 2 n x log | γ r | = 1 | γ r | 2 + O (1) , ∂ n x ∂ n y log | γ r | = − 1 | γ r | 2 + O (1) , (68) so the first tw o terms in (59) contribute, up to b ounded remainders, (69) − 1 2 π  ν ( t ) ∂ 2 n x log | γ r | + ν ( s ) ∂ n x ∂ n y log | γ r |  = − 1 2 π ν ( t ) − ν ( s ) | γ r | 2 + O (1) = − 1 2 π ν r | γ r | 2 + O (1) = − 1 2 π 1 r | ˙ γ ( t ) | dν dτ ( t ) + O (1) . On the other hand, ∂ τ x log | γ r | = τ ( t ) · γ r | γ r | 2 = 1 r | ˙ γ ( t ) | + O (1) , (70) so the normal-rotation term − dν dτ ( t ) S τ ,k [ σ ] cancels the 1 r singularit y coming from ν ( t ) S ′′ k [ σ ] + D ′ k [ ν σ ]. The remaining pieces (including S ′ k [ κν σ ]) are bounded, since they in volv e only first deriv atives of log | γ r | and the C 2 remainder R k . Note that Q ′′ k ( ρ ) is O (log ρ ), which is b ounded when combined with the ν ( t ) − ν ( s ) factor, following essentially the same calculations as in (69). Th us the kernel of the combination in (59) is b ounded on the diagonal. 5.4. Computing δ k and δ σ . Let k = k ε and σ = σ ε satisfy (48) on the perturb ed b oundary ∂ Ω ε . Differentiating (48) yields  1 2 I − S ′ k  [ δ σ ] = δ ( S ′ k )[ σ ] + ( δ k ) ( ∂ k S ′ k )[ σ ] . (71) Let µ b e a left nullv ector of 1 2 I − S ′ k and normalize it by ⟨ µ, σ ⟩ L 2 ( ∂ Ω) = 1. T aking the L 2 ( ∂ Ω) inner pro duct of (71) with µ and using ( 1 2 I − S ′ k ) ∗ µ = 0 yields 0 =  µ, δ ( S ′ k )[ σ ]  L 2 ( ∂ Ω) + ( δ k )  µ, ( ∂ k S ′ k )[ σ ]  L 2 ( ∂ Ω) . (72) Assuming ⟨ µ, ( ∂ k S ′ k )[ σ ] ⟩ L 2 ( ∂ Ω)  = 0 (e.g. if the first eigen v alue is simple), w e obtain δ k = − ⟨ µ, δ ( S ′ k )[ σ ] ⟩ L 2 ( ∂ Ω) ⟨ µ, ( ∂ k S ′ k )[ σ ] ⟩ L 2 ( ∂ Ω) . (73) Once δ k is known, equation (71) determines δ σ up to adding a multiple of σ . T o fix this, we add the rank one integral op erator µ ⟨ σ, ·⟩ that imp oses the orthogonality constraint ⟨ σ, δ σ ⟩ L 2 ( ∂ Ω) = 0 , (74) whic h yields a unique solution of (71). 5.5. Normalization and deriv ative of the ob jectiv e. T o ev aluate δ N in (56), we use a stan- dard identit y due to Rellich [56] that for an y Dirichlet eigenfunction (44) with eigenv alue k 2 , 2 k 2 Z Ω u 2 dx = Z ∂ Ω ( ∂ n u ) 2 ( x · n ) ds. (75) Then, recalling σ = − ∂ n u , we hav e N (Ω) = ∥ u ∥ 2 L 2 (Ω) = 1 2 k 2 Z ∂ Ω σ 2 ( x · n ) ds. (76) Differen tiating (76) and using Lemma 10 together with δ ( x · n ) = δ x · n + x · δ n = ν − ( x · τ ) dν dτ , (77) PEAK HEA T FLUX CONJECTURE 15 w e obtain δ N = − 2 δ k k N + 1 k 2 Z ∂ Ω σ δ σ ( x · n ) ds + 1 2 k 2 Z ∂ Ω σ 2  ν − ( x · τ ) dν dτ + κν ( x · n )  ds. (78) Com bining (56), (73), (71), and (78) giv es the deriv ative of F under normal p erturbation ν . 5.6. Non-normal p erturbations. In this subsection, we consider a non-normal deformation field V ( s ). The perturb ed curve b ecomes γ ε ( s ) := γ ( s ) + εV ( s ) . (79) Note that V ( s ) = n ( s ) ν ( s ) corresp onds to the normal p erturbation scenario that we studied earlier. F ollowing the notation used in Lemma 10, we obtain the deriv ativ es of | ˙ γ | and n using similar calculations. Lemma 12. Under the non-normal p erturb ation γ ε ( s ) := γ ( s ) + εV ( s ) , we have δ ( | ˙ γ | ) = ˙ γ · ˙ V | ˙ γ | , δ n = R ˙ V | ˙ γ | − n ˙ γ · ˙ V | ˙ γ | 2 . (80) Next, we compute the deriv ative of S ′ k . Lemma 13. The derivative of S ′ k is given by (81) δ ( S ′ k )[ σ ]( t ) = Z L 0 ( U 1 ( t, s ) + U 2 ( t, s ) + U 3 ( t, s )) σ ( s ) | ˙ γ ( s ) | ds U 1 ( t, s ) : = R dV dτ ( t ) · ∇ x G k ( γ ( t ) , γ ( s )) U 2 ( t, s ) : = n ( t ) ·  ∇ 2 x G k ( γ ( t ) , γ ( s ))( V ( t ) − V ( s ))  U 3 ( t, s ) : = ∂ n x G k ( γ ( t ) , γ ( s ))  τ ( s ) · dV dτ ( s ) − τ ( t ) · dV dτ ( t )  . Pr o of. Recall that S ′ k [ σ ]( t ) = p.v. Z L 0 ( n ( t ) · A ( t, s )) σ ( s ) | ˙ γ ( s ) | ds, A ( t, s ) : = ∇ x G k ( γ ( t ) , γ ( s )) (82) Using the pro duct rule, we hav e δ ( S ′ k )[ σ ]( t ) = Z L 0  ( δ n ( t )) · A ( t, s ) + n ( t ) · δ ( A ( t, s )) + ( n ( t ) · A ( t, s )) δ | ˙ γ ( s ) | | ˙ γ ( s ) |  σ ( s ) | ˙ γ ( s ) | ds. (83) By straightforw ard calculations and an application of the c hain rule, w e hav e ∂ ε ∇ x G k ( γ ε ( t ) , γ ε ( s )) =  ∇ 2 x G k ( γ ε ( t ) , γ ε ( s ))  ∂ ε γ ε ( t ) + ( ∇ y ∇ x G k ( γ ε ( t ) , γ ε ( s ))) ∂ ε γ ε ( s ) . (84) Ev aluating at ε = 0 gives us δ A ( t, s ) = ∇ 2 x G k ( γ ( t ) , γ ( s ))( V ( t ) − V ( s )) . (85) This, together with Lemma 12, giv es the desired result. □ W e need to establish b oundedness of the k ernel δ ( S ′ k ) on the diagonal. The argumen t is essen- tially the same as the one used in the normal case. In the non-normal case, the 1 /r cancellation is b et ween U 1 and U 2 . W e define r and γ r as b efore, and let V r = V ( t ) − V ( s ). W e consider the logarithmic part of the free-space Green’s function. Straightforw ard computations give us ∇ x log | γ r | = γ r | γ r | 2 , ∇ 2 x log | γ r | = I | γ r | 2 − 2 γ r γ T r | γ r | 4 , (86) where T denotes vector transp ose and I is the identit y matrix. F or U 1 , we hav e R ˙ V | ˙ γ | · γ r | γ r | 2 = R ˙ V · ( r ˙ γ + O ( r 2 )) r 2 | ˙ γ | 3 = − ˙ V · n r | ˙ γ | 2 + O (1) , (87) 16 ZW, JH, MR, AND AB where the last equality uses the definition of the rotation R . F or U 2 , we hav e n T  I | γ r | 2 − 2 γ r γ T r | γ r | 4  V r = n · V r | γ r | 2 − 2( n · γ r )( γ r · V r ) | γ r | 4 . (88) Since n · ˙ γ = 0, b oth the numerator and the denominator of the second term are O ( r 4 ). Moving to the first term, we hav e n · V r | γ r | 2 = r n · ˙ V + O ( r 2 ) r 2 | ˙ γ | 2 + O ( r 3 ) = n · ˙ V r | ˙ γ | + O (1) . (89) This cancels with the 1 /r singularity in (87). Same as the normal case, it is straightforw ard to v erify b oundedness of the remaining terms. T o rearrange δ S ′ k in la y er p oten tials, it is conv enient to decomp ose V into normal and tangen tial contributions: V = αn + β τ . (90) Straigh tforward computation gives us δ S ′ k [ σ ] =  − dα dτ S τ ,k [ σ ] + S ′ k [ κασ ] + αS ′′ k [ σ ] + D ′ [ ασ ]  +  β d dτ S ′ k [ σ ] − S ′ k  dβ dτ σ  , (91) where the first part agrees with the normal case. Finally , we differentiate the normalization term. Lemma 14. The derivative of the normalization term N is given by δ N = − 2 δ k k N + 1 k 2 Z ∂ Ω σ δ σ ( x · n ) ds + 1 2 k 2 Z ∂ Ω σ 2  V · n + x ·  R dV dτ  ds (92) Pr o of. W e start with the Rellic h identit y stated in (75). Recall that N = 1 2 k 2 Z L 0 σ ( s ) 2 B ( s ) ds, B ( s ) := ( γ ( s ) · n ( s )) | ˙ γ ( s ) | . (93) An application of the pro duct rule giv es us δ N = − 2 δ k k N + 1 2 k 2 Z L 0  2 σ ( s ) δ σ ( s ) B ( s ) + σ ( s ) 2 δ B ( s )  ds. (94) In the normal p erturbation case, w e formulate δ k and δ σ using δ ( S ′ k ). Therefore, we do not need to recompute these t wo terms. It remains to compute δ B , which follo ws from the iden tities in Lemma 12. Observe that (95) δ B = ( δ γ · n ) | ˙ γ | + ( γ · δ n ) | ˙ γ | + ( γ · n ) δ | ˙ γ | = | ˙ γ | ( V · n ) + γ · R ˙ V | ˙ γ | − n ˙ γ · ˙ V | ˙ γ | 2 ! + ( γ · n ) ˙ γ · ˙ V | ˙ γ | = ( V · n ) | ˙ γ | + γ · R ˙ V . Substituting this back in to (94) completes the pro of. □ Remark 15. In our numeric al exp eriments, we appr oximate a domain Ω using a r ounde d p olygon with N vertic es ( r i cos θ i , r i sin θ i ) i =1 , 2 ,...N . T o p erform optimization on F , we r e quir e the gr adient with r esp e ct to the r adial p ar ameters r i . We start by evaluating δ F for 2 N non-normal p erturb a- tions, c orr esp onding to the tr anslation of e ach vertex in the x and y dir e ctions. F or e ach vertex, we pr oje ct its x and y derivatives onto the (cos θ i , sin θ i ) dir e ction to obtain the derivative with r esp e ct to r i . PEAK HEA T FLUX CONJECTURE 17 5.7. F ast ev aluation of shap e deriv atives. The ev aluation of the shap e deriv ativ e of F requires the ev aluation of the v ariations δ k , δ σ , and δ N . In order to ev aluate δ k we need to compute δ S ′ k [ σ ], and ∂ k S ′ k [ σ ]. The standard Helmholtz fast m ultip ole metho d (FMM) in tw o dimensions can be used for the fast ev aluation of δ S ′ k [ σ ]. Ho wev er, the ev aluation of ∂ k S ′ k [ σ ] requires small mo difications as the k ernel does not satisfy the Helmholtz equation. A simple calculation sho ws that (96) ∂ k S ′ k [ σ ] = ik 4 Z ∂ Ω ( x − y ) · n ( x ) H 1 0 ( k | x − y | ) σ ( y ) ds . This kernel can be ev aluated using three Helmholtz FMM computations. Let (97) ϕ 1 ( x ) = ik 4 Z ∂ Ω y 1 H 1 0 ( k | x − y | ) σ ( y ) ds , ϕ 2 ( x ) = ik 4 Z ∂ Ω y 2 H 1 0 ( k | x − y | ) σ ( y ) ds , ϕ 3 ( x ) = ik 4 ( x · n ( x )) Z ∂ Ω H 1 0 ( k | x − y | ) σ ( y ) ds , eac h of which is computable using the standard Helmholtz FMM. Then (98) ∂ k S ′ k [ σ ]( x ) = ϕ 3 ( x ) − n 1 ( x ) ϕ 1 ( x ) − n 2 ϕ 2 ( x ) . T urning to the computation of δ σ , w e need to solve an integral equation with the op erator 1 / 2 I − S ′ k . As observ ed in Remark 15, the v ariation δ F is computed one perturbation direction at a time. Moreov er, the computation of the ob jective function F also requires solving the same in tegral equation with different data. Thus when computing δ F for domains described by N radii, O ( N ) solves with the op erator 1 / 2 I − S ′ k are required, and fast direct solvers tend to b e the metho d of c hoice in this setting. In particular w e use the recursive sk eletonization approach [26] for obtaining a compressed approximation of the inv erse. Finally , the ev aluation of δ N is straightforw ard as it requires an in tegral of a smo oth function on the b oundary ∂ Ω, once δ σ , and δ k are av ailable. 6. Numerical Experiments W e appro ximate domains in D 0 b y polygons parameterized in p olar coordinates. F or N ≥ 2, w e define equally spaced angles θ i := i − 1 N − 1 π for i = 1 , . . . , N . Given radii r = ( r 1 , . . . , r N ) ∈ R N > 0 , w e set p i ( r ) := ( x i , y i ) := ( r i cos θ i , r i sin θ i ) , i = 1 , . . . , N , (99) and let ˜ Ω N ( r ) b e the N -gon obtained by connecting p 1 ( r ) , . . . , p N ( r ) in order and closing the b oundary with the segment [ p N ( r ) , p 1 ( r )] ⊂ { y = 0 } . The domain Ω N ( r ) is then obtained by rounding ˜ Ω N ( r ) follo wing the pro cedure in [18] with a minor mo dification discussed below. Remark 16. The p olygonal domain is r ounde d to simplify the c omputation of the Dirichlet eigen- value, which is obtaine d by finding the r o ots of a discr etization of the F r e dholm determinant of I − 2 S ′ k [66] . The F r e dholm determinant is wel l-define d for se c ond kind F r e dholm op er ators of tr ac e class to which I − 2 S ′ k b elongs when the b oundary is at le ast C 2 . The numeric al pr o c e dur e discusse d her e c an b e extende d to p olygonal domains by inste ad c omputing the eigenvalue via the zer os of f ( k ) = 1 / ( v · ( I − 2 S ′ k ) − 1 u ) for some smo oth functions u and v define d on ∂ Ω , se e [13, 36] , for example. The rounding pro cedure in [18] replaces the p olygon in the vicinit y of every vertex by an affine transformation of a smooth approximation of the graph of f ( x ) = | x | which agrees with | x | outside the interv al [ − h, h ]. The rounding parameter h at vertex p i is chosen to b e α min( | p i − p i +1 | , | p i − p i − 1 | ) for some α < 1 / 2. How ever, the resulting b oundary deformations V corresp onding to the x and y p erturbations of the vertex p i will b e non-smo oth whenever | p i − p i +1 | = | p i − p i − 1 | . W e instead use a smo oth version of the min function giv en by smin( a, b ) = 1 log (exp (1 /a ) + exp (1 /b ) − 1) , 18 ZW, JH, MR, AND AB whic h satisfies min( a, b ) / 2 ≤ smin( a, b ) ≤ min( a, b ). In all the exp erimen ts in this section α = 0 . 1. W e define the discretized ob jective f N : R N > 0 → R , f N ( r ) := F (Ω N ( r )) , (100) where F is the functional in (4). W e discuss details of n umerical ev aluations of the ob jective in Section 6.1. F or a fixed N , we maximize f N b y gradient ascen t starting from some initialization to find the optimal rounded N -gon. The gradient is assem bled from shape deriv atives of the individual terms in F and then conv erted to deriv atives with resp ect to the parameters r i , as describ ed in Section 5. 6.1. Ob jective computation. F or eac h parameter v ector r , we compute the first Diric hlet eigen- v alue λ 1 of Ω N ( r ) together with the v alue of the functional F (Ω N ( r )). W riting λ 1 = k 2 1 , our computation is based on the b oundary integral form ulation described in Section 5. In particu- lar, if u 1 is the first eigenfunction and σ := − ∂ n u 1 denotes its inw ard normal deriv ative (with n the outw ard unit normal), then σ lies in the n ullspace of the second-kind op erator 1 2 I − S ′ k 1 . Th us k 1 is c haracterized b y the existence of a nontrivial solution of  1 2 I − S ′ k  [ σ ] = 0. Since S ′ k is compact for smo oth Γ = ∂ Ω, the corresp onding F redholm determinan t v anishes precisely at Diric hlet eigenv alues. Using this fact, w e locate k 1 b y searching for the smallest positive ro ot of this determinant [8, 66]. T o do so n umerically , we use the MA TLAB softw are pack age c hunkIE [4]. It rounds the corners of a p olygon and splits the rounded p olygon’s boundary into panels. Each panel is represen ted b y 16th-order Gauss–Legendre no des. W e discretize 1 2 I − S ′ k b y a high-order Nystr¨ om metho d, ev aluate the resulting determinant at Chebyshev p oin ts in a brack et [ a, b ] for k 1 , and extract the smallest real ro ot using a Cheb yshev interpolant [49]. T o accelerate the rep eated determinant ev aluations, w e use the recursive sk eletonization routines from FLAM [25] and fmm2d [21] via c hunkIE. Once k 1 is known, we compute a corresp onding discrete nontriv al null v ector σ , then normalize it using the Rellich identit y (76), N (Ω) := ∥ u 1 ∥ 2 L 2 (Ω) = 1 2 k 2 1 Z Γ ( x · n ) σ ( x ) 2 d s x , (101) whic h w e appro ximate by b oundary quadrature. Finally , we ev aluate σ at the distinguished b oundary p oin t x ∗ = (0 , 0) b y high-order interpolation on the panel containing x ∗ and compute the ob jectiv e via (55): f N ( r ) = F (Ω N ( r )) = σ ( x ∗ ) k 2 1 p N (Ω N ( r )) . (102) 6.2. Results. W e first fix N = 16 and maximize f 16 from sev eral initial p olygons, as sho wn in Figure 2. T o test robustness under discretization refinement, we increase the n um b er of vertices during the optimization. Starting from N = 8, we alternate b et w een optimizing at fixed N until ∥∇ f N ( r ) ∥ 2 ≤ η and refining the polygon b y inserting additional v ertices. W e terminate refinement once N reac hes a prescrib ed target N target . See Algorithm 1 for details of our v ertex refinement strategy . Algorithm 2 details the gradient ascent pro cedure, including the choice of step size with tw o lo cal ev aluations as well as fallback strategy . See Figure 3 for details on a large-scale exp erimen t using Ver texRefinement , which result with a semidisk-lik e shap e with N = 113 v ertices. W e use a large enough K so that GradientAscent is not limited by the maxim um n umber of iterations. PEAK HEA T FLUX CONJECTURE 19 Figure 2. Optimization results with N = 16: Each ro w represen ts a different initialization. Going from top to b ottom, w e hav e circle, square and triangle initializations. Each column is a different time step during optimization. Algorithm 1 Ver texRefinement Require: target num b er of vertices N target , gradient tolerance η > 0, max iterations K 1: Initialize N ← 8 and radii r ∈ R N > 0 . 2: while true do 3: r ← GradientAscent ( r, η , K ). 4: if N ≥ N target then 5: break 6: end if 7: Let θ i := i − 1 N − 1 π and p i := ( r i cos θ i , r i sin θ i ) for i = 1 , . . . , N . 8: Set ˜ r 2 i − 1 ← r i for i = 1 , . . . , N . 9: for i = 1 , . . . , N − 1 do 10: ˜ p 2 i ← [ p i , p i +1 ] ∩ { ( ρ cos θ i + θ i +1 2 , ρ sin θ i + θ i +1 2 ) : ρ > 0 } . 11: ˜ r 2 i ← ∥ ˜ p 2 i ∥ 2 . 12: end for 13: Set r ← ˜ r and N ← 2 N − 1. 14: end while Figure 3. A large scale exp erimen t that starts with N = 8 vertices, and gradu- ally increases the n umber of v ertices exceeds N target = 100 with gradien t threshold set to η = 5 × 10 − 6 . Under vertex addition strategy in Ver texRefinement , the resulting p olygon has N = 113 vertices. 20 ZW, JH, MR, AND AB Algorithm 2 GradientAscent Require: initialization r (0) ∈ R N > 0 , gradient tolerance η > 0, max iterations K 1: P arameter: step size h = 10 − 2 2: for k = 0 , 1 , . . . , K − 1 do 3: g ← ∇ f N ( r ( k ) ) 4: if ∥ g ∥ 2 ≤ η then 5: break 6: end if 7: Set direction d ← g / ∥ g ∥ 2 8: f 0 ← f N ( r ( k ) ) 9: f + ← f N ( r ( k ) + hd ) and f − ← f N ( r ( k ) − hd ) 10: a ← f + + f − − 2 f 0 2 h 2 and b ← f + − f − 2 h 11: α ∗ ← − b 2 a 12: if α ∗ < 0 then 13: r ( k +1) ← r ( k ) + g 14: else 15: r ( k +1) ← r ( k ) + α ∗ d 16: end if 17: end for 18: return updated radii r PEAK HEA T FLUX CONJECTURE 21 7. Conclusion and future work W e studied the size of the normal deriv ative of the first Diric hlet eigenfunction on the boundary of a conv ex planar domain. W e prov ed a universal upp er b ound of the form ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) ≤ C λ 1 (Theorem 2). Motiv ated by numerical exp erimen ts, w e form ulated Conjecture 3, predicting that the semidisk maximizes the scale-inv arian t quantit y G (Ω) = ∥ ∂ n u 1 ∥ L ∞ ( ∂ Ω) /λ 1 among b ounded con vex planar domains. As a partial result, w e prov ed that the semidisk is a critical p oint of the p oin twise functional F under a natural class of perturbations of the semicircular arc (Theorem 4). W e also described a b oundary in tegral form ulation and derived the shape deriv ativ es needed to assemble gradients of the discretized ob jectiv e efficiently . The resulting gradient ascent exp er- imen ts, across a range of initializations and discretization refinements, consistently conv erged to semidisk-lik e shap es, providing further evidence for Conjecture 3. There are several natural directions for future work. First, one would like to establish existence of an extremizer for G in a suitable top ology . Since G is scale-inv arian t, one ma y restrict to a class of conv ex domains that is compact in the Hausdorff metric; by the Blasc hk e selection theorem [7] suc h classes are sequentially compact. Thus, a key remaining step is to prov e an appropriate con tinuit y (or upp er semi-contin uity) statement for G under Hausdorff conv ergence. Second, it w ould b e interesting to obtain a more rigorous c haracterization of the maximizer by proving lo cal optimalit y of the semidisk (e.g. via second-v ariation analysis and conv exity constraints). Finally , the same methodology can b e applied to related optimization problems, such as other norms of ∂ n u 1 on ∂ Ω, different eigenmo des, or alternativ e b oundary conditions. A cknowledgment W e are grateful for discussions with Andrew Hassell, Charlie Epstein, and Stefan Steinerb erger. J.G.H. was supp orted in part by a Sloan Research F ellowship. The Flatiron Institute is a division of the Simons F oundation. 22 ZW, JH, MR, AND AB References [1] Angelo Alvino, Vincenzo F erone, and Guido T rombetti, On the pr op erties of some nonline ar eigenvalues , SIAM Journal on Mathematical Analysis 29 (1998), no. 2, 437–451. [2] Pedro R. S. Antunes and Pedro F reitas, Numerical optimization of low eigenvalues of the Dirichlet and Neu- mann Laplacians , J. Optim. Theory Appl. 154 (2012), 235–257. [3] Mark S Ashbaugh and Rafael D Benguria, A sharp b ound for the ratio of the first two eigenvalues of Dirichlet Laplacians and extensions , Annals of Mathematics 135 (1992), no. 3, 601–628. [4] T ravis Askham, Manas Rac hh, Michael O’Neil, Jerem y Hoskins, Daniel F ortunato, Shidong Jiang, F redrik F ryklund, T ristan Go odwill, Hai Y ang W ang, and Hai Zh u, chunkIE: a MA TLAB inte gr al equation toolb ox , 2024. V ersion 1.0.0, av ailable at https://chunkie.readthedocs.io/ . [5] V o jisla v G Av akumo vi´ c, ¨ ub er die Eigenfunktionen auf geschlossenen Riemannschen Mannigfaltigkeiten , Math- ematische Zeitschrift 65 (1956), no. 1, 327–344. [6] Anup Biswas, L o cation of maximizers of eigenfunctions of fr actional Schr¨ odinger’s equations , Mathematical Physics, Analysis and Geometry 20 (2017), no. 4, 25. [7] Wilhelm Blaschk e, Kr eis und kugel , V eit Leipzig, 1916. [8] F olkmar Bornemann, On the numeric al evaluation of Fre dholm determinants , Mathematics of Computation 79 (2010), no. 270, 871–915. [9] M-H Bossel, Membr anes ´ elastiquement li´ ees: extension du th ´ eor` eme de R ayleigh–F ab er–Kr ahn et de l’in´ egalit´ e de Chee ger , Comptes rendus de l’Acad´ emie des sciences. S´ erie 1, Math ´ ematique 302 (1986), no. 1, 47–50. [10] Helmut Brakhage and P eter W erner, ¨ Ub er das Dirichletsche A ußenr aumprob lem f¨ ur die Helmholtz sche Schwingungsgleichung , Archiv der Mathematik 16 (1965), no. 1, 325–329. [11] Lorenzo Brasco, Guido De Philippis, and Bozhidar V elic hko v, F ab er–Kr ahn inequalities in sharp quantitative form (2015). [12] Lorenzo Brasco and Aldo Pratelli, Sharp stability of some sp e ctral ine qualities , Geometric and F unctional Analysis 22 (2012), no. 1, 107–135. [13] H Cheng, W Crutchfield, M Do ery, and Leslie Greengard, F ast, ac cur ate inte gr al e quation methods for the analysis of photonic crystal fib ers i: The ory , Optics Express 12 (2004), no. 16, 3791–3805. [14] David Colton and Rainer Kress, Integr al e quation metho ds in sc attering the ory , SIAM, 2013. [15] Daniel Daners, A F ab er–Kr ahn ine quality for Robin pr oblems in any sp ac e dimension , Mathematisc he Annalen 335 (2006), no. 4, 767–785. [16] Jaume de Dios Pon t, Convex sets c an have interior hot sp ots , 2024. [17] Michel C Delfour and J-P Zol´ esio, Shapes and geometries: metrics, analysis, differ ential c alculus, and opti- mization , SIAM, 2011. [18] Charles L Epstein and Michael O’Neil, Smo othe d c orners and sc atter e d waves , SIAM Journal on Scientific Computing 38 (2016), no. 5, A2665–A2698. [19] Pedro F reitas and David Krejˇ ci ˇ r ´ ık, The first Robin eigenvalue with ne gative b oundary p ar ameter , Adv ances in Mathematics 280 (2015), 322–339. [20] Denis S Greb enk ov and B-T Nguy en, Ge ometric al structure of Laplacian eigenfunctions , SIAM Rev. 55 (2013), no. 4, 601–667. [21] Leslie Greengard and Vladimir Rokhlin, A fast algorithm for particle simulations , Journal of computational physics 73 (1987), no. 2, 325–348. [22] P Grinfeld, Hadamar d’s formula inside and out , Journal of optimization theory and applications 146 (2010), 654–690. [23] Andrew Hassell and T erence T ao, U pp er and lower b ounds for normal derivatives of Dirichlet eigenfunctions , Math. Res. Lett. 9 (2002), no. 2–3, 289–305. MR1909646 (2003k:58047) [24] Antoine Henrot, Extremum pr oblems for eigenvalues of el liptic op erators , Springer, 2006. [25] Kenneth L Ho, FLAM: F ast linear algebr a in MA TLAB—algorithms for hier ar chical matric es , Journal of Open Source Software 5 (2020), no. 51, 1906. [26] Kenneth L Ho and Leslie Greengard, A fast dir e ct solver for structure d line ar systems by r e cursive skeletoniza- tion , SIAM Journal on Scientific Computing 34 (2012), no. 5, A2507–A2532. [27] T Hohage and C Schormann, A Newton-type metho d for a transmission pr oblem in inverse scattering , Inv erse Problems 14 (1998), 1207–1227. [28] Thorsten Hohage, Lo garithmic c onver gence r ates of the iter atively re gularize d Gauss–Newton metho d for an inverse potential and an inverse sc attering problem , Inv erse problems 13 (1997), no. 5, 1279. [29] Lars V alter H¨ ormander, On the Riesz me ans of sp e ctr al functions and eigenfunction exp ansions for el liptic differ ential op erators , Matematika 12 (1968), no. 5, 91–130. [30] Jeremy G Hoskins and Stefan Steinerb erger, T owar ds optimal gr adient bounds for the torsion function in the plane , The Journal of Geometric Analysis 31 (2021), 7812–7841. [31] Mark Kac, Can one he ar the shap e of a drum? , The American Mathematical Monthly 73 (1966), no. 4P2, 1–23. [32] Grant Keady and Alex McNabb, The elastic torsion pr oblem: solutions in c onvex domains , NZ Journal of Mathematics 22 (1993), no. 43-64, 30. PEAK HEA T FLUX CONJECTURE 23 [33] James Kennedy, An isop erimetric ine quality for the se c ond eigenvalue of the Laplacian with Robin b oundary c onditions , Pro ceedings of the American Mathematical So ciet y 137 (2009), no. 2, 627–633. [34] Marie-Th´ er` ese Kohler-Jobin, Sur la pr emi` er e fonction pr opr e d’une membr ane: une extension ` a N dimensions de l’in´ egalit´ e isop´ erim´ etrique de Payne-R ayner , Zeitschrift f¨ ur angewandte Mathematik und Physik ZAMP 28 (1977), 1137–1140. [35] Edgar Krahn, ¨ Ub er eine von R ayleigh formulierte Minimaleigenschaft des Kr eises , Mathematisc he Annalen 94 (1925), no. 1, 97–100. [36] Jun Lai and Shidong Jiang, Se c ond kind inte gr al e quation formulation for the mode c alculation of optic al wave guides , Applied and Computational Harmonic Analysis 44 (2018), no. 3, 645–664. [37] Boris Moiseevich Levitan, On the asymptotic b ehavior of the sp e ctral function of a self-adjoint differ ential e quation of the se c ond or der , Izvestiy a Rossiisk oi Ak ademii Nauk. Seriya Matematichesk ay a 16 (1952), no. 4, 325–352. [38] Lihan Liu, Jingqiu Cai, and Y ongzhi Stev e Xu, R e gularize d Newton iter ation metho d for a p enetrable c avity with internal me asur ements in inverse sc attering pr oblem , Mathematical Metho ds in the Applied Sciences 43 (2020), no. 5, 2665–2678. [39] William Charles Hector McLean, Str ongly el liptic systems and b oundary inte gr al e quations , Cambridge uni- versit y press, 2000. [40] CB Moler and LE Pa yne, Bounds for eigenvalues and eigenve ctors of symmetric op er ators , SIAM Journal on Numerical Analysis 5 (1968), no. 1, 64–70. [41] J Mossino, A gener alization of the Payne–Rayner isop erimetric ine quality , Bolletino della Unione Mathematica Italiana 2 (1983), no. 3, 335–342. [42] F rank W. J. Olver, Daniel W. Lozier, Ronald F. Boisvert, and Charles W. Clark (eds.), NIST handb o ok of mathematic al functions , Cambridge Universit y Press, 2010. http://dlmf.nist.gov . [43] Braxton Osting, Optimization of spe ctr al functions of Dirichlet–Laplacian eigenvalues , J. Comput. Phys. 229 (2010), no. 22, 8578–90. [44] ´ Edouard Oudet, Numeric al minimization of eigenmo des of a membr ane with r esp e ct to the domain , ESAIM: Control, Optimisation and Calculus of V aritions 10 (2004), 315–335. [45] Lawrence E Payne and Margaret E Rayner, An isop erimetric ine quality for the first eigenfunction in the fixe d membr ane pr oblem , Zeitschrift f¨ ur angewandte Mathematik und Physik ZAMP 23 (1972), 13–15. [46] , Some isop erimetric norm b ounds for solutions of the Helmholtz e quation , Zeitschrift f ¨ ur angewandte Mathematik und Physik ZAMP 24 (1973), 105–110. [47] Lawrence Edward Pa yne and Hans F W ein b erger, Some isop erimetric ine qualities for membr ane fre quencies and torsional rigidity , Journal of Mathematical Analysis and Applications 2 (1961), no. 2, 210–216. [48] LE Payne and I Stakgold, On the me an value of the fundamental mode in the fixed membr ane pr oblem , Applicable Analysis 3 (1973), no. 3, 295–306. [49] Rodrigo B Platte and Lloyd N T refethen, Chebfun: a new kind of numerical c omputing , Progress in industrial mathematics at ecmi 2008, 2010, pp. 69–87. [50] G P´ olya, On the char acteristic fr e quencies of a symmetric membrane , 1955. [51] George P´ olya and G´ ab or Szeg¨ o, Isop erimetric ine qualities in mathematic al physics (1951). [52] Roland P otthast, Fr´ echet differ entiability of b oundary inte gr al op er ators in inverse acoustic sc attering , Inv erse Problems 10 (1994), 431–447. [53] , F r´ echet differ entiability of the solution to the ac oustic Neumann sc attering pr oblem with r esp e ct to the domain (1996). [54] Manas Rachh and Stefan Steinerb erger, On the lo c ation of maxima of solutions of Schr¨ odinger’s e quation , Communications on Pure and Applied Mathematics 71 (2018), no. 6, 1109–1122. [55] John William Strutt Baron Rayleigh, The the ory of sound , V ol. 2, Macmillan, 1896. [56] F ranz Rellich, Darstel lung der Eigenwerte von ∆ u + λu = 0 dur ch ein Randinte gral , Math. Z. 46 (1940), 635–636. [57] Christopher D. Sogge and Stev e Zelditch, Riemannian manifolds with maximal eigenfunction gr owth , Duk e Mathematical Journal 114 (2002), no. 3, 387 –437. [58] , Sup norms of Cauchy data of eigenfunctions on manifolds with c onc ave b oundary , Comm unications in Partial Differential Equations 42 (2017), no. 8, 1249–1289. [59] Ren´ e Peter Sp erb, Maximum principles and their applic ations (1981). [60] G´ ab or Szeg¨ o, Ine qualities for certain eigenvalues of a membr ane of given ar e a , Journal of Rational Mechanics and Analysis 3 (1954), 343–356. [61] M V an Den Berg, On the L-infinity norm of the first eigenfunction of the Dirichlet Laplacian , Poten tial Analysis 13 (2000), no. 4, 361–366. [62] Qiaoling W ang and Changyu Xia, Isoperimetric b ounds for the first eigenvalue of the laplacian , Zeitschrift f ¨ ur angewandte Mathematik und Physik 61 (2010), no. 1, 171–175. [63] Hans F W einberger, An isop erimetric ine quality for the n-dimensional fr ee membr ane pr oblem , Journal of Rational Mechanics and Analysis 5 (1956), no. 4, 633–636. [64] H. W eyl, Ueb er die asymptotische V erteilung der Eigenwerte , Nachric hten v on der Gesellschaft der Wis- senschaften zu G¨ ottingen, Mathematisch-Ph ysik alische Klasse 1911 (1911), 110–117. [65] Steve Zelditch, Bil liar ds and b oundary tr ac es of eigenfunctions , 2003, pp. 1–22. 24 ZW, JH, MR, AND AB [66] Lin Zhao and Alex Barnett, R obust and efficient solution of the drum pr oblem via Nystr¨ om appr oximation of the Fre dholm determinant , SIAM Journal on Numerical Analysis 53 (2015), no. 4, 1984–2007.