Computations of eigenvalues and resonances on perturbed hyperbolic surfaces with cusps

TL;DR

This paper introduces a fast method for computing the scattering matrix on hyperbolic surfaces with cusps using the Neumann-to-Dirichlet map, enabling detailed analysis of resonances and arithmetic properties.

Contribution

The paper presents a novel, efficient approach to compute scattering matrices from FEM-derived Neumann-to-Dirichlet maps on hyperbolic surfaces with cusps.

Findings

Accurate computation of scattering matrices from FEM data.

Ability to analyze resonance behavior under conformal perturbations.

Identification of arithmetic surfaces via numerical experiments.

Abstract

In this paper we describe a simple method that allows for a fast direct computation of the scattering matrix for a surface with hyperbolic cusps from the Neumann-to-Dirichlet map on the compact manifold with boundary obtained by removing the cusps. We illustrate that even if the Neumann-to-Dirichlet map is obtained by a Finite Element Method (FEM) one can achieve good accuracy for the scattering matrix. We give various interesting examples of how this can be used to investigate the behaviour of resonances under conformal perturbations or when moving in Teichm\"{u}ller space. For example, based on numerical experiments we rediscover the four arithmetic surfaces of genus one with one cusp. This demonstrates that it is possible to identify arithmetic objects using FEM.

Figures

Tables

Table 1. Table 1 : Eigenvalues λ = 1 4 + t 2 𝜆 1 4 superscript 𝑡 2 \lambda=\frac{1}{4}+t^{2} for the space of odd functions on A 0 subscript 𝐴 0 A_{0} . Data from [ 4 ] for comparison

Computed $t$	Data from [4]
9.53369	9.5336…
12.1730	12.1730…
14.3585	14.3585…
16.1381	16.1380…
16.6443	16.6442…
18.1809	18.1809…
19.4847	19.4847…

Table 2. Table 2 : Resonances for B r subscript 𝐵 𝑟 B_{r}

Computed resonances for				Poles of (13)–(15)
$B_1=A_0$	$B_{1/\sqrt{2}}$	$B_{1/\sqrt{3}}$	$B_{0.5001}$
0.2500 + 7.0674 $i$	0.2499 + 7.0676 $i$	0.2501 + 7.0681 $i$	0.2499 + 7.0707 $i$	${\zeta }_1/2\approx$ 0.2500 + 7.0674 $i$
0.2500 + 10.5110 $i$	0.2499 + 10.5116 $i$	0.2502 + 10.5129 $i$	0.2484 + 10.5182 $i$	${\zeta }_2/2\approx$ 0.2500 + 10.5110 $i$
0.2500 + 12.5054 $i$	0.2504 + 12.5063 $i$	0.2497 + 12.5093 $i$	0.2516 + 12.5168 $i$	${\zeta }_3/2\approx$ 0.2500 + 12.5054 $i$
0.2500 + 15.2125 $i$	0.2495 + 15.2145 $i$	0.2489 + 15.2165 $i$	0.2470 + 15.2297 $i$	${\zeta }_4/2\approx$ 0.2500 + 15.2124 $i$
0.2501 + 16.4676 $i$	0.2495 + 16.4700 $i$	0.2491 + 16.4747 $i$	0.2528 + 16.4889 $i$	${\zeta }_5/2\approx$ 0.2500 + 16.4675 $i$
0.2500 + 18.7931 $i$	0.2496 + 18.7960 $i$	0.2497 + 18.8023 $i$	0.2419 + 18.8246 $i$	${\zeta }_6/2\approx$ 0.2500 + 18.7931 $i$
0.2499 + 20.4594 $i$	0.2507 + 20.4635 $i$	0.2424 + 20.4683 $i$	0.2570 + 20.4880 $i$	${\zeta }_7/2\approx$ 0.2500 + 20.4594 $i$
		-0.0001 + 2.8597 $i$		$\pi i/\log \left(3\right)\approx$ 2.8596 $i$
	-0.0000 + 4.5325 $i$			$\pi i/\log \left(2\right)\approx$ 4.5324 $i$
		-0.0006 + 8.5796 $i$		$3\pi i/\log \left(3\right)\approx$ 8.5788 $i$
	0.0002 + 13.5984 $i$			$3\pi i/\log \left(2\right)\approx$13.5971 $i$
		0.0007 + 14.3001 $i$		$5\pi i/\log \left(3\right)\approx$14.2980 $i$

Table 3. Table 3 : Embedded eigenvalues λ = 1 4 + t 2 𝜆 1 4 superscript 𝑡 2 \lambda=\frac{1}{4}+t^{2} for the space of even functions on B r subscript 𝐵 𝑟 B_{r} . All eigenvalues have multiplicity one. ∗ denotes eigenvalues for the so called old-forms missed in [ 37 ]

$B_1=A_0$		$B_{1/\sqrt{2}}$		$B_{1/\sqrt{3}}$
Computed	Data from	Computed	Data from	Computed	Data from
$t$	[4]	$t$	[17, 37]	$t$	[17, 37]
13.7798	13.7797…	8.92297	8.92288	5.09885	5.09874
17.7387	17.7386…	10.9206	10.9204	8.03918	8.03886
19.4237	19.4847…	13.7802	13.7798∗	9.74450	9.74375
		14.6855	14.6852	11.3470	11.3464
		16.4044	16.4041	11.8906	11.8900
		17.7394	17.7386∗	13.1362	13.1351
		17.8788	17.8780	13.7810	13.7798∗
		19.1261	19.1254	14.6278	14.6262
		19.4245	19.4235∗	15.8012	15.7995
				16.2727	16.2710
				16.7384	16.7362
				17.5021	17.5006
				17.7413	17.7385∗
				18.6501	18.6474
				18.9662	18.9626
				19.4268	19.4235∗
				19.8997	19.8961

Table 4. Table 4 : Resonances for C ℓ , 0 subscript 𝐶 ℓ 0 C_{\ell,0} when ℓ ℓ \ell is a special length. The actual computed values for C arccosh ⁡ ( 5 ) , 0 subscript 𝐶 arccosh 5 0 C_{\operatorname{arccosh}\left(5\right),0} may differ by one in the last digit from those shown in the first column

Computed resonances for			Poles of (18)–(22)
$C_{\mathrm{arccosh}\left(2\right),0}$	$C_{\mathrm{arccosh}\left(3\right),0}$	$C_{\mathrm{arccosh}\left(9\right),0}$
and
$C_{\mathrm{arccosh}\left(5\right),0}$
0.2500 + 7.0678 $i$	0.2498 + 7.0680 $i$	0.2500 + 7.0677 $i$	${\zeta }_1/2\approx$ 0.25+ 7.0674 $i$
0.2498 + 10.5130 $i$	0.2501 + 10.5127 $i$	0.2496 + 10.5117 $i$	${\zeta }_2/2\approx$ 0.25 + 10.5110 $i$
0.2507 + 12.5071 $i$	0.2495 + 12.5079 $i$	0.2495 + 12.5082 $i$	${\zeta }_3/2\approx$ 0.25 + 12.5054 $i$
0.2508 + 15.2170 $i$	0.2496 + 15.2183 $i$	0.2488 + 15.2183 $i$	${\zeta }_4/2\approx$ 0.25 + 15.2124 $i$
0.2505 + 16.4737 $i$	0.2497 + 16.4742 $i$	0.2499 + 16.4745 $i$	${\zeta }_5/2\approx$ 0.25 + 16.4675 $i$
		-0.0000 + 1.9520 $i$	$\pi i/\log \left(5\right)\approx$ 1.9520 $i$
0.0000 + 2.8596 $i$			$\pi i/\log \left(3\right)\approx$ 2.8596 $i$
-0.0001 + 4.5324 $i$	-0.0001 + 4.5326 $i$		$\pi i/\log \left(2\right)\approx$ 4.5324 $i$
		-0.0000 + 5.8562 $i$	$3\pi i/\log \left(5\right)\approx$ 5.8559 $i$
-0.0007 + 8.5797 $i$			$3\pi i/\log \left(3\right)\approx$ 8.5788 $i$
		-0.0005 + 9.7608 $i$	$5\pi i/\log \left(5\right)\approx$ 9.7599 $i$
0.0041 + 13.6006 $i$	-0.0034 + 13.6003 $i$		$3\pi i/\log \left(2\right)\approx$ 13.5971 $i$
		-0.0022 + 13.6671 $i$	$7\pi i/\log \left(5\right)\approx$ 13.6639 $i$
0.0009 + 14.3046 $i$			$5\pi i/\log \left(3\right)\approx$ 14.2980 $i$

Table 5. Table 5 : Resonances for C ℓ , 1 / 2 subscript 𝐶 ℓ 1 2 C_{\ell,1/2} when ℓ ℓ \ell is a special length

Computed resonances for			Poles of (23)
$C_{2\mathrm{arccosh}\left(\frac{3}{2}\right),1/2}$	$C_{2\mathrm{arccosh}\left(2\right),1/2}$	$C_{2\mathrm{arccosh}\left(3\right),1/2}$
0.2499+7.0681 $i$	0.2499+7.0678 $i$	0.2499+7.0681 $i$	${\zeta }_1/2\approx$ 0.25+7.0674 $i$
0.2503+10.5123 $i$	0.2501+10.5128 $i$	0.2500+10.5121 $i$	${\zeta }_2/2\approx$ 0.25+10.5110 $i$
0.2499+12.5078 $i$	0.2512+12.5067 $i$	0.2492+12.5079 $i$	${\zeta }_3/2\approx$ 0.25+12.5054 $i$
	0.0+4.5324 $i$		$\pi i/\log \left(2\right)\approx$ 4.5324 $i$
	0.0003+13.6004 $i$		$3\pi i/\log \left(2\right)\approx$ 13.5971 $i$

Table 6. Table 6 : Embedded eigenvalues λ = 1 4 + t 2 𝜆 1 4 superscript 𝑡 2 \lambda=\frac{1}{4}+t^{2} and their multiplicities μ ​ ( t ) 𝜇 𝑡 \mu(t) for C ℓ , 0 subscript 𝐶 ℓ 0 C_{\ell,0} and C ℓ , 1 / 2 subscript 𝐶 ℓ 1 2 C_{\ell,1/2} when ℓ ℓ \ell is a special length. The actual computed values for C arccosh ⁡ ( 5 ) , 0 subscript 𝐶 arccosh 5 0 C_{\operatorname{arccosh}\left(5\right),0} and C 2 ​ arccosh ⁡ ( 3 ) , 1 / 2 subscript 𝐶 2 arccosh 3 1 2 C_{2\operatorname{arccosh}\left(3\right),1/2} may differ in the last digit from those shown in the first and fourth columns, resp. The last two eigenvalues for C arccosh ⁡ ( 3 ) , 0 subscript 𝐶 arccosh 3 0 C_{\operatorname{arccosh}\left(3\right),0} are shown for comparison with those for B 1 / 2 subscript 𝐵 1 2 B_{1/\sqrt{2}} in Table 3 . A subset of these eigenvalues are eigenvalues for the groups Γ ~ 0 ​ ( 5 ) subscript ~ Γ 0 5 \widetilde{\Gamma}_{0}(5) , Γ ~ 0 ​ ( 6 ) subscript ~ Γ 0 6 \widetilde{\Gamma}_{0}(6) , and Γ 3 superscript Γ 3 \Gamma^{3} . These are in good agreement with those computed in [ 22 ] and [ 33 ]

$C_{\mathrm{arccosh}\left(2\right),0}$		$C_{\mathrm{arccosh}\left(3\right),0}$		$C_{\mathrm{arccosh}\left(9\right),0}$		$C_{2\mathrm{arccosh}\left(\frac{3}{2}\right),1/2}$		$C_{2\mathrm{arccosh}\left(2\right),1/2}$
and						and
$C_{\mathrm{arccosh}\left(5\right),0}$						$C_{2\mathrm{arccosh}\left(3\right),1/2}$
$t$	$\mu \left(t\right)$	$t$	$\mu \left(t\right)$	$t$	$\mu \left(t\right)$	$t$	$\mu \left(t\right)$	$t$	$\mu \left(t\right)$
2.42507	1	2.89100	2	2.00968	1	2.95648	2	2.89101	2
3.24141	1	3.25000	1	3.44480	1	3.53606	2	3.25001	1
3.97879	1	4.13811	2	3.70334	1	3.70339	1	4.13811	2
4.54850	1	4.36806	1	4.13245	1	4.51375	2	4.36809	1
4.64665	1	4.95729	2	4.65694	1	5.50420	2	4.95731	2
4.94791	1	5.62822	1	4.89729	1	5.81512	2	5.62824	1
5.09888	1	6.02334	2	5.05935	1	5.87951	1	6.02335	2
5.19203	1	6.22332	1	5.34525	1	6.62069	1	6.22332	1
5.35557	1	6.69430	2	5.87949	1	6.64683	2	6.69441	2
6.12073	1	7.22111	1	6.05422	1	6.78381	2	7.22571	1
		⋮
		8.92338	1
		⋮
		10.9213	1

Table 7. Table 7 : Resonances for D 𝐷 D

0.2500 + 7.0675 $i$	\rdelim}3*[${\zeta }_1/2\approx$ 0.25+7.0674 $i$]
0.2500 + 7.0676 $i$
0.2499 + 7.0680 $i$
0.2497 + 10.5117 $i$	\rdelim}3*[${\zeta }_2/2\approx$ 0.25+10.5110 $i$]
0.2496 + 10.5121 $i$
0.2500 + 10.5125 $i$
0.2504 + 12.5060 $i$	\rdelim}3*[${\zeta }_3/2\approx$ 0.25+12.5054 $i$]
0.2506 + 12.5062 $i$
0.2496 + 12.5078 $i$
0.2510 + 15.2152 $i$	\rdelim}3*[${\zeta }_4/2\approx$ 0.25+15.2124 $i$]
0.2489 + 15.2153 $i$
0.2490 + 15.2154 $i$
0.2503 + 16.4695 $i$	\rdelim}3*[${\zeta }_5/2\approx$ 0.25+16.4675 $i$]
0.2504 + 16.4699 $i$
0.2494 + 16.4724 $i$
0.2463 + 18.7973 $i$	\rdelim}3*[${\zeta }_6/2\approx$ 0.25+18.7931 $i$]
0.2454 + 18.7983 $i$
0.2503 + 18.7984 $i$
-0.0000 + 4.5324 $i$	\rdelim}2*[$\pi i/\log \left(2\right)\approx$ 4.5324 $i$]
-0.0000 + 4.5325 $i$
0.0002 + 9.0654 $i$	\rdelim}2*[$2\pi i/\log \left(2\right)\approx$ 9.0647 $i$]
0.0002 + 9.0657 $i$
-0.0005 + 13.5988 $i$	\rdelim}2*[$3\pi i/\log \left(2\right)\approx$ 13.5971 $i$]
-0.0006 + 13.5997 $i$
-0.0000 + 18.1367 $i$	\rdelim}2*[$4\pi i/\log \left(2\right)\approx$ 18.1294 $i$]
0.0018 + 18.1369 $i$

Equations

$$y^{-2}(\mathrm{d}{x}^2+\mathrm{d}{y}^2).$$

$$⟨f,g⟩=\int f(z)\overline{g(z)}{y}^{-2}\mathrm{d}x\mathrm{d}y.$$

$$\Delta ={\Delta }_{\mathbb{H}}=-y^2\left({\partial }_x^2+{\partial }_y^2\right)$$

$${\int }_a^{\infty }{\int }_{-1/2}^{+1/2}|f(z){|}^2{y}^{-2}\mathrm{d}x\mathrm{d}y$$

$$X=M{\cup }_{\partial Z}Z,Z=Z_1\bigsqcup \cdots \bigsqcup Z_p,Z_k=(\mathbb{R}/\mathbb{Z})\times [a_k,\infty )$$

$$P=\Delta +V(x),$$

$$\lambda =s(1-s),$$

$$s=\frac{1}{2}+\mathrm{i}t,$$

$$\lambda =\frac{1}{4}+t^2.$$

$$(P-\lambda {)}^{-1}=(P-s(1-s){)}^{-1},$$

$${E_j(z,s)|}_{Z_k}={\delta }_{j,k}{y}_k^s+C_{j,k}(s)y_k^{1-s}+T_{j,k}(z_k,s),$$

$$\overline{C(s)}=C^{\ast }(s)=C(\overline{s}),$$

$$C(s)C(1-s)=\text{mathbbm}1,$$

$$(\Delta -s(1-s))f(z)=0$$

$$\left(-y^2\frac{d^2}{d{y}^2}+4{\pi }^2{m}^2{y}^2-s(1-s)\right)f_m(y)=0.$$

$$T_{j,k}(z,s)=\sqrt{y}\sum _{m\in \mathbb{Z}\backslash \{0\}}{a}_{m,j,k}(s)K_{\mathrm{i}t}(2\pi |m|y)e^{2\pi \mathrm{i}mx},$$

$$\left(P-s_r(1-s_r)\right)f_r=0,$$

$${f_r|}_{Z_k}=a_{r,k}{y}_k^{1-s_r}+R_{r,k}(z_k),$$

$$⟨P\psi ,\varphi ⟩-⟨\psi ,P\varphi ⟩={\int }_{\partial M}\left(\psi (z)\frac{\partial \overline{\varphi }}{\partial \mathbf{n}}(z)-\frac{\partial \psi }{\partial \mathbf{n}}(z)\overline{\varphi }(z)\right)\mathrm{d}z$$

$${\int }_{\partial M}\left(\psi (z)\frac{\partial \overline{\varphi }}{\partial \mathbf{n}}(z)-\frac{\partial \psi }{\partial \mathbf{n}}(z)\overline{\varphi }(z)\right)\mathrm{d}z=\sum _{k=1}^p{\int }_{\partial M_k}\left(\psi (z)\frac{\partial \overline{\varphi }}{\partial \mathbf{n}}(z)-\frac{\partial \psi }{\partial \mathbf{n}}(z)\overline{\varphi }(z)\right)\mathrm{d}z.$$

$${\int }_{\partial M_k}\left(\psi (z)\frac{\partial \overline{\varphi }}{\partial \mathbf{n}}(z)-\frac{\partial \psi }{\partial \mathbf{n}}(z)\overline{\varphi }(z)\right)\mathrm{d}z={\int }_{-\frac{1}{2}}^{\frac{1}{2}}{\left(\psi (z)a_k\frac{\partial \overline{\varphi }}{\partial y}(z)-a_k\frac{\partial \psi }{\partial y}(z)\overline{\varphi }(z)\right)|}_{y=a_k}\frac{1}{a_k}\mathrm{d}x.$$

$$P_{\text{Neu}}{\Phi }_j={\lambda }_j{\Phi }_j,$$

$${\frac{\partial {\Phi }_j}{\partial \mathbf{n}}|}_{\partial M}=0,$$

$$\begin{array}{rl}(P-\lambda )\psi & =0,\text{ in }M,\\ {\frac{\partial \psi }{\partial \mathbf{n}}|}_{\partial M}& =f.\end{array}$$

$${\mathcal{N}}^M(s)f:=\psi {|}_{\partial M},$$

$${\mathcal{N}}^M(s)f=\sum _j\frac{1}{{\lambda }_j-s(1-s)}⟨f,{\varphi }_j{⟩}_{L^2(\partial M)}{\varphi }_j,$$

$$\left({\mathcal{N}}^M(s)-{\mathcal{N}}^M(s_0)\right)f=\sum _j\frac{s_0(1-s_0)-s(1-s)}{({\lambda }_j-s(1-s))({\lambda }_j-s_0(1-s_0))}⟨f,{\varphi }_j{⟩}_{L^2(\partial M)}{\varphi }_j.$$

$$f(z)=\sum _{m\in \mathbb{Z}}{f}_m(y)e_m(x),$$

$$L_{\mathrm{cusp}}^2(Z^a)=\{f\in L^2(Z^a)\mid f_0(y)=0a.e.\}.$$

$$(\Delta -\lambda )\psi =0,$$

$${-a\frac{\partial \psi }{\partial y}|}_{y=a}=f.$$

$${\mathcal{N}}^{Z^a}(s)(f-\mathbf{av}(f)):=\psi {|}_{y=a},$$

$$\mathbf{av}(f):={\int }_{-1/2}^{1/2}f(x)\mathrm{d}x.$$

Full text

1. 1 Introduction and setup

2. 2 Plan of the paper and discussion of the results

3. 3 The Neumann-to-Dirichlet operator on $\partial M$

4. 4 The Neumann-to-Dirichlet operator on cusps

5. 5 The relation between the Neumann-to-Dirichlet operator and the scattering matrix

6. 6 Numerical computation of the scattering matrix, resonances and embedded eigenvalues

7. 6.1 Scattering matrix and resonances

8. 6.2 Error estimates for the scattering matrix

9. 6.3 Embedded eigenvalues

10. 7 Examples and numerical studies

11. 7.1 General set-up

12. 7.2 Benchmarking

13. 7.3 $A_{\phi}$. The modular domain and conformal perturbations

14. 7.3.1 Description of the surface

15. 7.3.2 Known properties of the spectrum

16. 7.3.3 Numerical results

17. 7.4 $B_{r}$. Artin’s billiard

18. 7.4.1 Description of the surface

19. 7.4.2 Known properties of the spectrum

20. 7.4.3 Numerical results

21. 7.5 $C_{\ell,\tau}$. Hyperbolic surfaces of genus one with one cusp

22. 7.5.1 Description of the surface

23. 7.5.2 Known properties of the spectrum

24. 7.5.3 Numerical results, case 1: $\tau=0$, varying $\ell$

25. 7.5.4 Numerical results, case 2: $\tau=\frac{1}{2}$, varying $\ell$

26. 7.5.5 Numerical results, case 3: $\ell=2\operatorname{arccosh}\left(\frac{3}{2}\right)\approx 1.924847$, varying $\tau$

27. 7.5.6 Numerical results, case 4: equal length geodesics, varying $\tau$ and $\ell=\ell^{*}(\tau)$

28. 7.6 $D$. The hyperbolic surface of genus zero with three cusps

29. 7.6.1 Description of the surface

30. 7.6.2 Numerical results

31. Acknowledgements

32. A Tables of resonances and eigenvalues

1. 1 Fundamental domain of a cusp in $\mathbb{H}$
2. 2 Part of a cusp isometrically embedded into $\mathbb{R}^{3}$
3. 3 A surface of genus one with two hyperbolic cusps
4. 4 Fundamental domain for the modular surface
5. 5 $c(x,y)$** as a function on the modular domain**
6. 6 Trajectories of eight selected resonances of $A_{\phi_{q}}$ as $q$ changes
7. 7 **The reduced modular domain for Artin’s billiard **$B_{1/\sqrt{2}}$
8. 8 **Trajectories of four selected resonances for **$B_{r}$
9. 9 **Resonances for $B_{r}$ with special values of **$r$
10. 10 Fundamental domain for a hyperbolic surface of genus one with one cusp
11. 11 **Trajectories of four selected resonances for **$C_{\ell,0}$
12. 12 **Sequence of symmetry reductions of $C_{\ell,0}$, **$\ell=\operatorname{arccosh}(3)$
13. 13 **Resonances and embedded eigenvalues for $C_{\ell,0}$ with special lengths **$\ell$
14. 14 **Trajectories of four selected resonances for **$C_{\ell,1/2}$
15. 15 **Resonances for $C_{\ell,1/2}$ for special lengths **$\ell$
16. 16 Fundamental domain for a hyperbolic surface of genus zero with three cusp
17. 17 **Resonances for **$D$

1. 1 **Eigenvalues $\lambda=\frac{1}{4}+t^{2}$ for the space of odd functions on **$A_{0}$
2. 2 **Resonances for **$B_{r}$
3. 3 **Embedded eigenvalues $\lambda=\frac{1}{4}+t^{2}$ for the space of even functions on **$B_{r}$
4. 4 Resonances for $C_{\ell,0}$ when $\ell$ is a special length
5. 5 Resonances for $C_{\ell,1/2}$ when $\ell$ is a special length
6. 6 Embedded eigenvalues $\lambda=\frac{1}{4}+t^{2}$ and their multiplicities $\mu(t)$ for $C_{\ell,0}$ and $C_{\ell,1/2}$ when $\ell$ is a special length
7. 7 **Resonances for **$D$

\listofvideos

1 Introduction and setup

Suppose that $\mathbb{H}=\{z=x+\mathrm{i}y\mid y>0\}$ is the upper half-plane with the hyperbolic (of constant curvature $-1$) metric

[TABLE]

The Riemannian measure is then $y^{-2}\mathrm{d}x\mathrm{d}y$ and the $L^{2}$-inner product is given by

[TABLE]

The metric Laplace operator

[TABLE]

is essentially self-adjoint with domain $C^{\infty}_{0}(\mathbb{H})$, and later on we do not distinguish notationally operators and their closures, if there is no danger of confusion.

The map $(x,y)\mapsto(x+1,y)$ is an isometry of the upper half-space , and the quotient of the set $\mathbb{H}_{a}=\{z=x+\mathrm{i}y\mid y>a\}$ by this isometry results in a so-called hyperbolic cusp with height $a>0$. Thus, such a cusp $Z^{a}$ is topologically equivalent to $\mathbb{S}^{1}\times[a,\infty)$ and it is equipped with a metric of constant negative curvature.

Figure 1 shows a fundamental domain that becomes $Z^{a}$ when the parallel sides are identified. Of course the space of smooth functions $C^{\infty}(Z^{a})$ on $Z^{a}$ can be identified with smooth functions on $\mathbb{H}$ periodic in $x$ (with period one) and similarly, $L^{2}(Z^{a})$ can be identified with the set of measurable functions $f(z)$ on $\mathbb{H}$, periodic in $x$ (with period one), such that the $L^{2}$-norm

[TABLE]

is finite. We will in the following use these identifications without further mention. The Neumann Laplace operator $\Delta_{Z^{a}}$ on the cusp $Z^{a}$ is obtained by imposing Neumann boundary condition on the operator $\Delta=-y^{2}\left(\partial_{x}^{2}+\partial_{y}^{2}\right)$ on the boundary $({\mathbb{R}}/{\mathbb{Z}})\times\{a\}$. This operator is self-adjoint and has spectrum consisting of an absolutely continuous part $[\frac{1}{4},\infty)$ and of a discrete set of non-negative eigenvalues with finite dimensional eigenspaces (see for example [28] or [21] and references there).

Suppose that $X$ is a complete two-dimensional Riemannian manifold (or orbifold with finitely many isolated orbifold singularities) that is either a hyperbolic cusp or a disjoint union of hyperbolic cusps outside of a compact region. Thus, we are assuming that $X$ has the form

[TABLE]

such that the Riemannian metric $g$ on $X$ restricted to a neighbourhood of the cusp $Z_{k}$ is the hyperbolic metric defined above (see Figure 3).

Now assume that $P$ is a formally self-adjoint differential operator of Laplace type on $X$ acting on functions (which means that $P=-g^{ij}\partial_{i}\partial_{j}+\text{lower order terms}$), and let $\Delta$ be the Laplace operator acting on functions on $X$. Thus, $P-\Delta$ is a first order operator and we will assume that $P-\Delta$ is compactly supported away from each $Z_{k}$. The simplest example would be

[TABLE]

where $V\in C^{\infty}(X)$ is a potential that is supported in the interior of $M$. However, we do not want to exclude more general cases here. Since $X$ is complete the operator $P$ is essentially self-adjoint on $C^{\infty}_{0}(X)$.

Remark 1.1.

All our formulae and conclusions hold true with the obvious modifications if $M$ has additional boundary components and/or conical singularities away from cusps $Z$, and appropriate elliptic boundary conditions are imposed there to make $P$ a self-adjoint operator.

Manifolds with such cusps were considered and analysed in [28] and [29] and the gluing constructions for the heat kernel carry over to our setting. The structure of the spectrum and the generalised eigenfunctions can also be inferred from the meromorphic continuation of the resolvent. This approach can be found for example in [11]. In the following we summarise the known results.

As the Neumann Laplace operator on $Z^{a}$, the operator $P$ has spectrum consisting of the absolutely continuous part $[1/4,\infty)$ of multiplicity $p$ and, maybe, eigenvalues of finite multiplicity. As usual the resolvent $(P-\lambda)^{-1}$ is often more conveniently described using other parameters $s$ and $t$ which are related to the spectral parameter $\lambda$ in the following way,

[TABLE]

The set on which the resolvent is naturally defined as a meromorphic function with values in the space of bounded operators is $\mathbb{C}\backslash[\frac{1}{4},\infty)$ in terms of $\lambda$, the half-plane $\operatorname{Re}(s)>1/2$ in terms of $s$, and the lower half-space in terms of the parameter $t$. The resolvent,

[TABLE]

viewed as an operator from $L^{2}_{\text{comp}}(X)$ to $H^{2}_{\text{loc}}(X)$, admits a meromorphic continuation as a function of $s$ to the entire complex plane with poles of finite rank (that is, all the negative Laurent coefficients are finite rank operators). These poles correspond to eigenvalues and so-called scattering resonances.

The generalised eigenfunction $E_{j}(z,s)$ of the operator $P$, attached to the cusp $Z_{j}$, can be constructed from the resolvent and therefore admits a meromorphic continuation to $\mathbb{C}$ as a function of $s$. When restricted to $Z_{k}$ it is of the form

[TABLE]

where $T_{j,k}(z,s)$ is in $L^{2}(Z_{k})$. Here $z_{k}=x_{k}+\mathrm{i}\,y_{k}$ denotes the coordinates on the cusp $Z_{k}$. Both $C(s)$ and $T(z,s)$ are meromorphic matrix-valued functions of $s$ in the entire complex plane. The matrix-valued meromorphic function $C(s)$ is defined by (1) and is normally referred to as the scattering matrix. It satisfies the relations

[TABLE]

implying that it is unitary on the absolutely continuous spectrum.

Since $Z^{a}$ has a natural $\mathbb{S}^{1}$-action we can decompose, in the case $P=\Delta$, the solutions of

[TABLE]

into the Fourier modes $f(z)=\sum\limits_{m\in\mathbb{Z}}f_{m}(y)e^{2\pi\mathrm{i}mx}$ that satisfy

[TABLE]

For $m=0$ and $s\neq\frac{1}{2}$ this implies that $f_{0}(y)$ is a linear combination of $y^{s}$ and $y^{1-s}$. For $m\not=0$ the general solution of this ODE can be expressed in terms of Bessel functions. Then we obtain

[TABLE]

where convergence is in $C^{\infty}(Z^{a})$. Here $K_{\nu}$ is the modified Bessel $K$-function of order $\nu$.

Poles of the scattering matrix $C(s)$ are called resonances. Resonances correspond to poles of the generalised eigenfunctions $E_{j}(z,s)$ and the coefficient of the lowest term in the Laurent expansion at a resonance $s_{r}$ is proportional to a function $f_{r}\in C^{\infty}(X)$ such that

[TABLE]

and

[TABLE]

where $R_{r,k}$ is exponentially decaying as $y_{k}\to\infty$. The function $f_{r}$ is sometimes referred to as the resonant state at the resonance $s_{r}$.

2 Plan of the paper and discussion of the results

The main aim of this paper is to demonstrate that the domain decomposition using the Neumann-to-Dirichlet map leads to a simple and fast numerical scheme allowing to compute the scattering matrix on spaces with hyperbolic cusps.

The paper is structured as follows.

In Sections 3 and 4 we construct the Neumann-to-Dirichlet maps on the compact part of a hyperbolic surface and on the cusps, respectively. In Section 5 we show that the scattering matrix can be extracted from the Neumann-to-Dirichlet operator of a compact part of a hyperbolic surface with cusps by means of simple linear algebra methods. In particular, if a numerical approximation of the Neumann-to-Dirichlet map at a spectral point is provided by any method, fast and standard linear algebra routines can be used to extract the scattering matrix. In Section 6 we show that in fact standard finite element methods are already sufficient to calculate the scattering matrix, and hence the scattering resonances, with good accuracy if the spectral parameter is not too large. Various examples of constant but also non-constant curvature are treated and discussed in detail in Section 7. We compare them to known values for arithmetic surfaces as computed for example by Winkler ([37]) and Hejhal ([17]). Since our method is extremely fast and flexible we were able to produce moving pictures that show how scattering resonances move with conformal perturbations or in Teichmüller space. Figures illustrating this are included in Section 7. In particular, in genus one case we identified several surfaces for which the scattering matrix is expressible in terms of the Riemann zeta function. These surfaces correspond to the four arithmetic surfaces known to exist in genus one with one cusp. It seems that these arithmetic surfaces are the only ones (up to isomorphism) for which the resonances are lined up along critical lines.

For surfaces of constant negative curvature there are direct fast converging methods that allow the computation of embedded eigenvalues and scattering resonances. For example Hejhal’s algorithm can be used to compute embedded eigenvalues with extreme accuracy (see for example [5], see also [4]), and is used to compute large numbers of high lying eigenvalues (for example [22] for arithmetic examples). Variations have also been used to track resonances (for example [9, 1, 3]). Our approach is different in that it treats the compact part as a black-box and also allows for perturbations away from constant curvature. The correspondence between the scattering matrix and the Neumann-to-Dirichlet map can be used to relate number theoretic questions to transmission problems. This approach was taken independently in [8] in the context of quotients of hyperbolic space by Fuchsian groups and leads to a reformulation of the Riemann hypothesis in terms of transmission eigenvalues.

We would like to point out that numerical instabilities leading to spurious eigenvalues or eigenvalues being missed seem to be absent in our approach. We give several tables comparing our results to known computations in arithmetic constant curvature situations.

3 The Neumann-to-Dirichlet operator on $\partial M$

The operator $P$ is a formally self-adjoint elliptic differential operator on $X$ that coincides with the Laplace operator near the boundary of $M$. Therefore, we have Green’s formula

[TABLE]

for all $\phi,\psi\in C^{\infty}(M)$. In our case the boundary $\partial M$ is a disjoint union of components $\partial M_{k}=\partial Z_{k}$ each of which is isometric to the circle. We therefore have

[TABLE]

Given a particular boundary component $\partial M_{k}$ we can choose coordinates $(x,y)$ such that the cusp $Z_{k}$ corresponds to $S^{1}\times[a_{k},\infty)$. In this case $\frac{1}{a_{k}}dx$ is the natural Riemannian measure induced by the metric on the boundary and $a_{k}\frac{\partial}{\partial y}$ is the unit normal vector field. We therefore have

[TABLE]

We can hence construct another self-adjoint operator $P_{\text{Neu}}$ on $L^{2}(M)$ by restricting $P$ to $M$ and imposing Neumann boundary conditions at the boundary $\partial M$. Since $P_{\text{Neu}}$ is self-adjoint and elliptic there exists an orthonormal basis in $L^{2}(M)$ consisting of smooth eigenfunctions $(\Phi_{j})_{j\in\mathbb{N}}$ such that

[TABLE]

where $\lambda_{1}\leq\lambda_{2}\leq\ldots\to\infty$ are the corresponding eigenvalues.

If $\lambda\in\mathbb{C}$ is not a Neumann eigenvalue then for each $f\in C^{\infty}(\partial M)$ there exists a unique function $\psi\in C^{\infty}(M)$ such that

[TABLE]

The so-called Neumann-to-Dirichlet operator $\mathcal{N}^{M}(s):C^{\infty}(\partial M)\to C^{\infty}(\partial M)$ is defined as

[TABLE]

where $\psi\in C^{\infty}(M)$ is the solution of (4).

Separating between the different boundary components the Neumann-to-Dirichlet map can also be thought of as a matrix of operators $\mathcal{N}^{M}_{kj}(s):C^{\infty}(\partial M_{j})\to C^{\infty}(\partial M_{k})$. It is well known that $\mathcal{N}^{M}(s)$ is a pseudodifferential operator of order $-1$ whose full symbol depends only on the germ of the metric near the boundary (see [25] in case $s=0$, but the proof given there works in general). In particular the off-diagonal terms of the matrix $\mathcal{N}^{M}_{kj}(s)$ are smoothing operators and the diagonal ones are pseudodifferential operators of order $-1$ acting on $C^{\infty}(\partial M_{j})$. Using Green’s formula one easily obtains

[TABLE]

where $\phi_{j}=\Phi_{j}|_{\partial M}$ are the restrictions of the Neumann eigenfunctions $\Phi_{j}$ to the boundary $\partial M$ of $M$ and the sum converges in $H^{1/2}(\partial M)$ (see [24]). Taking differences one obtains

[TABLE]

This converges in $H^{3/2}(\partial M)$ uniformly with respect to the $H^{2}(\partial M)$-norm of $f$. In particular, $\mathcal{N}^{M}(s)$ is a meromorphic family of pseudodifferential operators of order $-1$ with first order poles at $s_{j}$ that are related to the Neumann eigenvalues $\lambda_{j}$ of $P_{\text{Neu}}$ by $\lambda_{j}=s_{j}(1-s_{j})$. The family of operators $\mathcal{N}^{M}(s)$ is hence completely determined by the data $(\phi_{j},\lambda_{j})_{j\in\mathbb{N}}$.

4 The Neumann-to-Dirichlet operator on cusps

Since the $Z^{a}$ admits an $S^{1}$-action the space $L^{2}(Z^{a})$ every function $f\in L^{2}(Z^{a})$ may be decomposed into Fourier modes

[TABLE]

where $e_{m}(x)=e^{2\pi\mathrm{i}mx}$. The functions with vanishing zero Fourier coefficients form a sub-space in $L^{2}(Z^{a})$, the so called cuspidal functions

[TABLE]

The orthogonal complement $L^{2}_{0}(Z^{a})$ of $L^{2}_{\mathrm{cusp}}(Z^{a})$ is then the space of functions that do not depend on $x$. This space is canonically isomorphic to $L^{2}((a,\infty),y^{-2}dz)$. The Neumann Laplace operator leaves both spaces invariant. Its restriction to $L^{2}_{0}(Z^{a})$ has absolutely continuous spectrum $[\frac{1}{4},\infty)$ and the restriction to $L^{2}_{\mathrm{cusp}}(Z^{a})$ has purely discrete spectrum consisting of eigenvalues of finite multiplicity accumulating at $\infty$. If $\lambda=s(1-s)$ is not a eigenvalue of the Neumann Laplace operator on $L^{2}_{\mathrm{cusp}}(Z^{a})$ then for each $f\in L^{2}(S^{1})$ with $\int\limits_{-1/2}^{1/2}f(x)\mathrm{d}x=0$ there exists a unique function $\psi\in L^{2}(Z^{a})$ such that

[TABLE]

We will define the cuspidal Neumann-to-Dirichlet operator $\mathcal{N}^{Z^{a}}(s):C^{\infty}(S^{1})\to C^{\infty}(S^{1})$ as

[TABLE]

where

[TABLE]

This operator has an explicit description in terms of Bessel functions. Namely, it follows directly from the expansion into Fourier modes that for any $m\not=0$ we have

[TABLE]

and $\mathcal{N}^{Z_{k}}(s)e_{0}=0$. Since the boundary of $M$ consists of a disjoint union of components $\partial M_{k}$ we can assemble the Neumann-to-Dirichlet operator to an operator $\mathcal{N}^{c}(s)$ acting on $L^{2}(\partial M)=\bigoplus_{k=1}^{p}L^{2}(\partial M_{k})$ by defining

[TABLE]

In the same way the averaging operator can be assembled to a map $\mathord{\bf av}\!:L^{2}(\partial M)\to L^{2}(\partial M)$.

5 The relation between the Neumann-to-Dirichlet operator and the scattering matrix

The generalised eigenfunctions $E_{j}(z,s)$ form a meromorphic family of functions satisfying $(P-\lambda)E_{j}(z,s)=0$ on all of $M$. Hence,

[TABLE]

On the other hand the restriction of $E_{j}(z,s)$ to each cusp has an expansion of the form (1) with a decaying tail term. We therefore have

[TABLE]

This means in particular that $\frac{\partial}{\partial\mathbf{n}}E_{j}(z,s)|_{\partial M}$ is in the kernel of the map

[TABLE]

Note that the averaging map ${\mathord{\bf av}}\!:L^{2}(\partial M)\to L^{2}(\partial M)$ is the orthogonal projection onto the space of locally constant functions $L^{2}_{0}(\partial M)$ on $\partial M$. This space is naturally identified with $\mathbb{C}^{p}$, the $k$-th component being identified with the function value on the boundary component $\partial M_{k}$.

Theorem 5.1.

Suppose that $s\not=\frac{1}{2}$ is a complex number that is not a pole of $\mathcal{N}^{M}(s)$ or $\mathcal{N}^{c}(s)$, and not a pole of the scattering matrix $C(s)$. Suppose furthermore that $s(1-s)$ is not an $L^{2}$-eigenvalue of $P$. Then the kernel of the map $(\mathbbm{1}-\mathord{\bf av})\mathcal{N}^{M}(s)+\mathcal{N}^{c}(s)$ is $p$-dimensional and spanned by $\{\frac{\partial}{\partial\mathbf{n}}E_{k}(z,s)|_{\partial M}\mid 1\leq k\leq p\}$.

Proof.

The assumptions imply that the generalised eigenfunctions $E_{j}(z,s)$ exist at $s$. We have already shown that $\frac{\partial}{\partial\mathbf{n}}E_{j}(z,s)|_{\partial M}$ is in the kernel of $(\mathbbm{1}-\mathord{\bf av})\mathcal{N}^{M}(s)+\mathcal{N}^{c}(s)$. Moreover, any non-zero linear combination of $E=\sum_{j}c_{j}E_{j}$ such that $\frac{\partial}{\partial\mathbf{n}}E_{j}(z,s)|_{\partial M}=0$ will give rise to an $L^{2}$-Neumann eigenfunction on $Z$ by taking the non-zero part of its Fourier expansion. Since we excluded Neumann eigenvalues on the cusp by the requirement that $s$ is not a pole of $\mathcal{N}^{c}$, the functions $\frac{\partial}{\partial\mathbf{n}}E_{j}(z,s)|_{\partial M}$ are linearly independent. Now suppose that $g\in L^{2}(\partial M)$ is in the kernel of $(\mathbbm{1}-\mathord{\bf av})\mathcal{N}^{M}+\mathcal{N}^{c}$. Both $(\mathbbm{1}-\mathord{\bf av})\mathcal{N}^{M}$ and $\mathcal{N}^{c}$ are elliptic pseudodifferential operators of order $-1$ and their principal symbols coincide. Hence, their sum is elliptic too and, by elliptic regularity, $g\in C^{\infty}(\partial M)$. This means that there is a function $F_{M}\in C^{\infty}(M)$ and a function $F_{Z}\in L^{2}(Z)\cap C^{\infty}(M)$ such that

[TABLE]

These equations imply that the functions $F_{M}$ and $F_{Z}$, when expanded into Fourier modes, have the same non-zero Fourier coefficients on each cusp. Hence, for each cusp $Z_{k}$ there exist coefficients $a_{k}$ and $b_{k}$ such that the function

[TABLE]

has the same Fourier expansion as $F_{M}$ on $Z_{k}$. Therefore, we can construct a globally defined function $F\in C^{\infty}(X)$ which agrees with $F_{M}$ on $M$, such that

[TABLE]

Now use Green’s identity on a cut-off domain $M_{R}$ obtained by cutting off the cusps $Z_{k}$ at $y=R$ and use the fact that the tail term is exponentially decaying:

[TABLE]

Define $E(z):=\sum\limits_{k=1}^{p}a_{k}E_{k}(z,s)$. Then the above implies that $F-E\in L^{2}$. Since $s(1-s)$ was assumed not to be an $L^{2}$-eigenvalue we conclude that $F=E$. ∎

Theorem 5.2.

Suppose that $s\not=\frac{1}{2}$ is a complex number that is not a pole of $\mathcal{N}^{M}(s)$ or $\mathcal{N}^{c}(s)$, and not a pole of the scattering matrix $C(s)$. Suppose furthermore that $s(1-s)$ is not an $L^{2}$-eigenvalue of $P$. Let $V$ be the kernel of the map $(\mathbbm{1}-\mathord{\bf av})\mathcal{N}^{M}(s)+\mathcal{N}^{c}(s)$ and define the maps

[TABLE]

Then the map $(s-1)Q_{2}+Q_{1}$ is invertible and

[TABLE]

where $A$ is the diagonal matrix $A=\operatorname{diag}(a_{1},a_{2},\ldots,a_{p})$.

Proof.

By the previous theorem we can use the basis $\left\{\phi_{k}:=\frac{\partial}{\partial\mathbf{n}}E_{k}(z,s)|_{\partial M}\mid 1\leq k\leq p\right\}$ in $V$ to check the invertibility of the map $(s-1)Q_{2}+Q_{1}$ and the formula for the scattering matrix. By the expansion (1) we have

[TABLE]

Therefore, $\left(((s-1)Q_{2}+Q_{1})\phi_{j}\right)_{k}=(2s-1)\delta_{j,k}a_{k}^{s}$ and the right hand side is a non-singular matrix. Moreover, $\left((sQ_{2}-Q_{1})\phi_{j}\right)_{k}=(2s-1)C_{j,k}(s)a_{k}^{1-s}$. The formula (7) is immediately implied by this. ∎

In the case of one cusp the above reduces to a generalised eigenvalue problem.

Corollary 5.3.

Assume $X$ has one cusp, i.e. $p=1$, and suppose that $s\not=\frac{1}{2}$ is a complex number that is not a pole of $\mathcal{N}^{M}(s)$ or $\mathcal{N}^{c}(s)$, and not a pole of the scattering matrix $C(s)$. Then either the pair $(\mathcal{N}^{M}(s)+\mathcal{N}^{c}(s),\mathord{\bf av})$ has precisely one generalised eigenvalue $G(s)$, or $C(s)=\frac{s}{s-1}a^{2s-1}$. In the former case the scattering matrix can be computed from this eigenvalue as

[TABLE]

6 Numerical computation of the scattering matrix, resonances and embedded eigenvalues

6.1 Scattering matrix and resonances

Theorems 5.1 and 5.2 yield an extremely simple and fast algorithm to compute the scattering matrix, resonances or eigenvalues for the situation described above. In this section we will assume that $s\not=\frac{1}{2}$ is a complex number that is not a pole of $\mathcal{N}^{M}(s)$ or $\mathcal{N}^{c}(s)$, and not a pole of the scattering matrix $C(s)$. In the following we take $(e_{m})_{m\in\mathbb{Z}}$ to be the orthonormal basis of $L^{2}({\mathbb{R}}/{\mathbb{Z}},dx)$ consisting of Fourier modes $e_{m}(x)=e^{2\pi\mathrm{i}mx}$. Since each boundary component $\partial M_{k}$ can be identified with a circle, this gives an orthonormal basis $(e_{m,k})_{m\in\mathbb{Z},\;k=1,\ldots,p}$ in $L^{2}(\partial M)$. We will write $(e_{\alpha})_{\alpha\in I}$ where the index set for $\alpha=(\alpha_{1},\alpha_{2})$ is $I:=\mathbb{Z}\times\{1,\ldots,p\}$.

The boundary data $(\phi_{j},\lambda_{j})_{j\in\mathbb{N}}$ of Neumann eigenvalues can be used to compute the matrix elements of the Neumann-to-Dirichlet operator $\mathcal{N}^{M}(s)$ using (5) and the Fourier expansion in the basis $(e_{\alpha})$ giving

[TABLE]

Using (6), convergence in (8) is accelerated if we compute the matrix elements $\mathcal{N}^{M}_{\alpha\beta}(s_{0})$ directly at a single particular value $s_{0}$, cf. [24]. Then

[TABLE]

and the series in (9) converges more rapidly than the one in (8). The acceleration trick may be repeated if one computes directly $\mathcal{N}^{M}_{\alpha,\beta}(s_{j})$ for several particular values of $s_{j}$.

The matrix elements of $\mathcal{N}^{Z^{a}}(s)$ are simply

[TABLE]

Moreover, $\mathord{\bf av}_{\alpha,\beta}=\delta_{\alpha_{1},0}\delta_{\beta_{1},0}$.

We would then like to find the $p\times p$ matrix $G(s)$ such that

[TABLE]

The idea of the numerical approximation is of course to truncate this Fourier basis and approximate the above matrices by finite matrices by considering only $0\leq|\alpha_{1}|,|\beta_{1}|\leq J$ for some large integer $J$. We denote by $\widetilde{N}^{M}$, $\widetilde{N}^{c}$, and $\widetilde{\mathord{\bf av}}$ the finite matrices obtained from truncating the Fourier expansion at $J$. Then these matrices are $(2J+1)p\times(2J+1)p$ matrices with complex entries.

We use the finite element method to compute the Neumann boundary data $(\phi_{j},\lambda_{j})_{j\in\mathbb{N}}$ in terms of the numbers $\lambda_{j}$ and the Fourier modes $\langle\phi_{j},e_{\alpha}\rangle$. Once these data are obtained a finite element approximation to $\widetilde{N}^{M}(s)$ can be computed very quickly for arbitrary $s\in\mathbb{C}$ in a given compact subset of the complex plane. The matrices $\widetilde{N}^{c}$ can be computed very fast using a well known continued fraction expansion for the Bessel $K$-function [6, Section 17],

[TABLE]

where we use Gauss’ notation

[TABLE]

In order to compute the scattering matrix numerically in the above approximation we proceed as follows. By Theorem 5.1 the operator $T(s):=(\mathbbm{1}-\mathord{\bf av})\mathcal{N}^{M}(s)+\mathcal{N}^{c}(s)$ has a $p$-dimensional kernel spanned by $\frac{\partial}{\partial\mathbf{n}}E_{k}(z,s)|_{\partial M}$. We compute the cut off approximation

[TABLE]

If $J$ is large enough this matrix will have precisely $p$ small singular values. We can therefore perform a singular value decomposition to construct an orthonormal system of singular vectors $(v_{1},\ldots,v_{p})$ with small singular values. The system of vectors $(\widetilde{\mathord{\bf av}}(v_{1}),\ldots,\widetilde{\mathord{\bf av}}(v_{p}))$ determines a $p\times p$ matrix $\widetilde{Q}_{1}$. Similarly the vectors $(\widetilde{\mathord{\bf av}}\widetilde{N}^{M}(v_{1}),\ldots,\widetilde{\mathord{\bf av}}\widetilde{N}^{M}(v_{p}))$ determine a $p\times p$ matrix $\widetilde{Q}_{2}$. Since the set of invertible maps is open the matrix $(s-1)\widetilde{Q}_{2}+\widetilde{Q}_{1}$ is invertible if the approximation is good enough. By Theorem 5.2 we then get a numerical approximation of the scattering matrix by

[TABLE]

As before $A$ is the diagonal matrix $A=\mathrm{diag}(a_{1},a_{2},\ldots,a_{p})$.

Since resonances are poles of $C(s)$ and we have the functional equation $C(s)C(1-s)=\mathbbm{1}$, the resonances are precisely the zeros of the determinant of $C(1-s)$.

6.2 Error estimates for the scattering matrix

In this section we will show that in principle the error in the computation can be made rigorous if the exterior and interior Neumann-to-Dirichlet maps are obtained by a method with rigorous errors. Let us start assuming that we have a mechanism at our disposal to estimate the first Sobolev norm of $\mathcal{N}^{M}(s)\Phi-\Psi$ for given smooth functions $\Phi$ and $\Psi$. This depends on a chosen method of computation of the Neumann-to-Dirichlet map.

On $L^{2}(\partial M)$ we have the orthonormal basis $(e_{m,k})_{m\in\mathbb{Z},\;k=1,\ldots,p}$. We define the Fourier multiplier $q:H^{s}(\partial M)\to H^{s-1}(\partial M)$ by $qe_{m,k}=\left(|m|+1\right)e_{m,k}$. The operator $q$ is a first order pseudodifferential operator and can also be expressed in terms of the Laplace operator on the boundary. For concreteness we fix the $H^{s}$-norm on $\partial M$ as $\|\Phi\|_{H^{s}(\partial M)}=\|q^{s}\Phi\|_{L^{2}(\partial M)}$.

In the following we will assume that $s\in\mathbb{C}$ is fixed such that the assumptions of Theorem 5.2 hold: $s\not=\frac{1}{2}$ is not a pole of $\mathcal{N}^{M}(s)$ or $\mathcal{N}^{c}(s)$, not a pole of the scattering matrix $C(s)$ and $s(1-s)$ is not an eigenvalue of $P$. Then

[TABLE]

is a zero order elliptic pseudodifferential operator. In particular [math] is not in the essential spectrum of $T^{\prime*}T^{\prime}$. This implies that the self-adjoint operator $|T^{\prime}|$ has [math] as a multiplicity $p$ eigenvalue and a spectral gap in the sense that the spectrum is contained in $\{0\}\cup[K_{1},\infty)$ for some $K_{1}>0$.

Our numerical approximation takes place in the finite dimensional subspace $\mathcal{W}_{J}$ of functions $f$ that have a finite Fourier expansion of the form

[TABLE]

As before $J$ is a sufficiently large integer. The method will then usually find an orthonormal set vectors $(v_{1},\ldots,v_{p})$ in $\mathcal{W}_{J}$ such that

[TABLE]

If $P_{0}$ is the orthogonal projection onto the $p$-dimensional kernel of $|T^{\prime}|$ it follows that

[TABLE]

Applying the numerical approximation of $\mathcal{N}^{M}(s)$ we obtain another set of vectors $(w_{1},\ldots,w_{p})$ in the subspace. Given an error estimate on the Dirichlet-to-Neumann map as assumed we will get a bound of the form

[TABLE]

The approximations $\widetilde{Q}_{1}$ and $\widetilde{Q}_{2}$ of the maps $Q_{1}$ and $Q_{2}$ can be though of as finite rank operators with range in the subspace $\mathcal{W}_{J}$ that vanish on the orthogonal complement of $\mathcal{W}_{J}$. Recall that $\widetilde{Q}_{1}=(\widetilde{\mathord{\bf av}}(v_{1}),\ldots,\widetilde{\mathord{\bf av}}(v_{p}))$ and $\widetilde{Q}_{2}=(\widetilde{\mathord{\bf av}}(w_{1}),\ldots,\widetilde{\mathord{\bf av}}(w_{p}))$. If we choose $(P_{0}v_{1},\ldots,P_{0}v_{p})$ as a basis in the kernel of $T$ to describe $Q_{1}$ and $Q_{2}$ we obtain

[TABLE]

where the norms are the operator norms of the respective matrices. Let

[TABLE]

and

[TABLE]

Then, assuming $\epsilon_{1}\,K_{2}<1$, we obtain

[TABLE]

Collecting everything we can now estimate the error of the approximated scattering matrix

[TABLE]

as

[TABLE]

If the error for the Neumann-to-Dirichlet map is known, and the norm $\|\mathcal{N}^{M}(s)\|$ and the spectral gap $K_{1}$ can be estimated or computed, the error for the scattering matrix can be explicitly bounded. In principle this makes it possible to use interval arithmetics to rigorously prove interval bounds for the scattering matrix. The scattering matrix is holomorphic in the resolvent set and by a classical theorem of Hurwitz the zeros of uniform approximations of $C(s)$ converge to zeros of $C(s)$. A quantitative version of this is given in [31]. This allows to estimate the error of the approximation of the computed resonances. Note that using a finite truncation of (6) approximates the Neumann-to-Dirichlet operator in the correct norm.

6.3 Embedded eigenvalues

Since we assumed $P$ was self-adjoint any eigenvalues will have to be on the real line. There are two classes of eigenvalues: those below the continuous spectrum and those embedded into the continuous spectrum. We will refer to the eigenvalues $\lambda<1/4$ as small eigenvalues and the eigenvalues $\lambda\geq 1/4$ as the embedded eigenvalues. Embedded eigenvalues correspond to real values of $t$ and therefore the real part of $s$ for these eigenvalues will always be $1/2$. As a consequence the zero modes of the Fourier expansion of these eigenfunctions in the cusp has to vanish. We therefore make the following observation.

Theorem 6.1.

The embedded eigenvalues away from the poles of $\mathcal{N}^{M}$ and $\mathcal{N}^{c}$ are exactly those values of $\lambda=s(1-s)\in\left[\frac{1}{4},+\infty\right)$ for which there exists a non-zero vector $f\in C^{\infty}(S^{1})$ such that $\mathord{\bf av}(f)=0$ and

[TABLE]

Computations of embedded eigenvalues face the problem that it is not possible to numerically distinguish between an embedded eigenvalue and a resonance that is close to the spectrum. Rigorous error estimates that guarantee the existence of an embedded eigenvalue therefore always need some additional information about the geometry or, in the constant curvature case, arithmetic nature of the surface (see for example [5]). The mathematically rigorous numerical part of this work is mostly about the computation of the scattering matrix and of resonances. We therefore only briefly sketch how one detects embedded eigenvalues or resonances close to the spectrum. We are looking for vectors $v$ that satisfy $\widetilde{\mathord{\bf av}}(v)=0$ and for which $\left(\widetilde{N}^{M}+\widetilde{N}^{c}\right)v$ is small. For numerical stability the $\mathrm{QR}$-decomposition $\widetilde{B}(s)=\widetilde{\mathbf{Q}}(s)\widetilde{\mathbf{R}}(s)$ of the matrix

[TABLE]

is performed. This matrix maps $\mathbb{C}^{(2J+1)p}$ to $\mathbb{C}^{(2J+1)p}\oplus\mathbb{C}^{(2J+1)p}\oplus\mathbb{C}^{(2J+1)p}\oplus\mathbb{C}^{(2J+1)p}$. Let ${\mathbf{P}}$ be the projection onto the first two summands. We are looking to find values of $s$ for which there exists a vector $v$ for which ${\mathbf{P}}\widetilde{B}(s)v$ is very small whereas $(1-p)\widetilde{B}(s)v$ is not. Since $\widetilde{\mathbf{R}}$ is invertible these are exactly the small singular values of the matrix ${\mathbf{P}}\widetilde{\mathbf{Q}}(s)$. Thus, our method of finding embedded eigenvalues is to plot the smallest singular value of ${\mathbf{P}}\widetilde{\mathbf{Q}}(s)$ as a function of $s=\frac{1}{2}+\mathrm{i}\sqrt{\lambda-\frac{1}{4}}$. If the smallest singular value is close to zero for some $s=\frac{1}{2}+\mathrm{i}\sqrt{\lambda-\frac{1}{4}}$ this amounts to a small spectral gap $K_{1}$, i.e. a small $(p+1)$-st singular value of $T^{\prime}$. Hence the error estimate for the scattering matrix near such a point becomes much worse, reflecting the fact that we may also have a resonance close to the spectrum. The existence of resonances near the spectrum or embedded eigenvalues and completeness of a list of computed values can heuristically be verified using Turing’s method and variants of the Weyl law that have been proved in this context [28]. Versions of the Weyl law with error estimates are available in the constant curvature case, see e.g. [34].

7 Examples and numerical studies

7.1 General set-up

In this section we will consider several examples of manifolds with cusps. We will be focusing on the Laplace operator acting on functions, i.e. in all the examples we will have $P=\Delta$. These examples are divided into groups as follows.

• $A_{\phi}$

the modular domain with its constant curvature metric changed by a conformal factor $\mathrm{e}^{\phi}$. This family is parametrised by smooth functions $\phi$ on the modular surface.

• $B_{r}$

a triangular domain that is sometimes referred to as Artin’s billiard and that interpolates between the Hecke triangular surfaces. This family is parametrised by a real number $r>\frac{1}{2}$.

• $C_{\ell,\tau}$

the surfaces of genus $(1,1)$ and constant curvature, i.e. the punctured torus. The Teichmüller space of genus $(1,1)$ has dimension $2$ and therefore this family is parametrised by a length parameter $\ell>0$ and a twist parameter $\tau\in[0,1)$.

• $D$

the unique hyperbolic surface of genus zero with three cusps.

In all these examples we decomposed the surface into compact part $M$ and a cusp-part. The method allows the freedom of choosing a cut-off parameter $a$. In the examples below $a$ was usually chosen in the interval $[0.3,2]$, depending on the geometry. Note that choosing significantly higher values of $a$ decreases the accuracy of a Neumann-to-Dirichlet map approximation, and choosing a small $a$ creates meshing problems due to a “narrow” compact part $M$. Experiments indicate that the dependence of computed eigenvalues and resonances upon a choice of $a$ in a suitable subset of the above interval is negligible.

To compute a numerical approximation of the Neumann-to-Dirichlet map we use the accelerated expansion (9). We used the finite element framework FreeFEM++ ([14] and [13]) to compute the Neumann-to-Dirichlet map at some point $s_{0}$ and to compute the boundary data of the first 1000 Neumann eigenvalues. On the boundary Fourier modes up to $|m|=40$ were used. In the FEM implementation we used discretisation with up to 200 points on the boundary of the compact part. The Neumann-to-Dirichlet map on the cusps is computed using (10) and (11).

The corresponding data were expressed in terms of Fourier modes on the boundary and imported into a Mathematica script that directly computed the scattering matrix by the method described before. Since then the scattering matrix was available as a numerical function, we used Newton’s root finding algorithm to locate zeros of its determinant. The functional equation (3) was then used to determine the scattering resonances. The poles and the zeros of the scattering matrix are located in the half-planes $\operatorname{Re}{s}>\frac{1}{2}$ and $\operatorname{Re}{s}<\frac{1}{2}$, respectively. Unless resonances are very close to the spectrum there are therefore no issues due to poles and zeros being close together. One can thus use the argument principle and contour integration to count the number of resonances in a region bounded away from the spectrum. We have found that in practice Newton’s root finding algorithm finds all resonances away from the spectrum in a fast reliable manner. This is due to the well-behaved analytic properties of the scattering matrix. To locate and track resonances that are very close to the spectrum we start from a perturbation of the surface and then use predictive algorithms based on polynomial extrapolation to follow the path of the resonance. This way even the resonances that seem to have high order touching points with the spectrum could be tracked.

The following conventions are assumed in all videos and graphs:

• •

Resonances and eigenvalues are traced in a part of $(\operatorname{Re}s,\operatorname{Im}s)\in\left(-\infty,\frac{1}{2}\right]\times[0,+\infty)$ quadrant of the $s$-plane.

• •

The resonances very close to the continuous spectrum $\frac{1}{2}+\mathrm{i}\left[0,+\infty\right)$, as well as embedded eigenvalues, are shown in blue. Not all embedded eigenvalues are shown.

• •

The resonances on the critical line $\frac{1}{4}+\mathrm{i}[0,+\infty)$ or very close to it are shown in red.

• •

The resonances on the line $\operatorname{Re}s=0$ or very close to it are shown in green.

• •

The eigenvalue at $s=0$ is never shown.

• •

In graphs showing the trajectories of resonances, the starting points of the trajectories are marked by a disk, and the end points by a square.

• •

In geometry figures, the arcs shown in the same colour and line type are identified; solid black lines indicate Neumann conditions imposed on the arcs.

7.2 Benchmarking

In the case of the modular surface $A_{0}$ the scattering matrix can be expressed in terms of the Riemann zeta function, see (12), and we could compare and compute the relative error of our approximation. The scattering matrix $C(s)$ computed for $s=\frac{1}{2}+\mathrm{i}t$, $t\in[0.,30.]$ (this amounts to the interval $[0.25,900.]$ in the spectrum) had a maximal relative error of about $0.25\%$. On the interval $[0,10]$ for $t$ we even obtained a maximal relative error not exceeding $0.004\%$. Similar errors hold on the critical line. We note that these approximations are surprisingly good considering that a finite element approximation was used. The finite element method and subsequent computations were carried out with double precision.

As shown in Subsection 6.2 the $(p+1)$-th singular value of $q\widetilde{T}$ is a measure of the spectral gap. Away from resonances close to the real line or embedded eigenvalues the size of the first $p$ smallest singular values of $q\widetilde{T}$ compared to the $(p+1)$-th small singular value was extremely small (typically of an order of a double precision rounding error) in our computations. Hence, using the terminology of Subsection 6.2 the numerical estimate for $K_{1}^{-1}\delta_{1}$ was very small and the theoretical error was dominated by $\delta_{2}$ which stems from the approximation of the Neumann-to-Dirichlet map. In our case most of the errors are due to the FEM approximation and decrease with mesh refinement.

The computational cost of building the scattering matrices using FEM realisation of Neumann-to-Dirichlet maps is relatively low if $\operatorname{Im}s\leq 30$ (this of course depends on the FEM implementation and the number of eigenvalues used). The real runtime costs actually occur when we look for complex roots and poles of the scattering matrix and trace individual resonance dependence on the parameters.

7.3 $A_{\phi}$. The modular domain and conformal perturbations

7.3.1 Description of the surface

This surface can be obtained from the domain

[TABLE]

by gluing along the boundary as follows. The sides $x=-\frac{1}{2}$ and $x=\frac{1}{2}$ are identified by means of the parallel translation $x\mapsto x+1$. The circular arc $\{(x,y)\in\mathbb{H}\mid x^{2}+y^{2}=1,0\leq x\leq\frac{1}{2}\}$ is identified with $\{(x,y)\in\mathbb{H}\mid x^{2}+y^{2}=1,-\frac{1}{2}\leq x\leq 0\}$ using the map $x\mapsto-x$. This results in a hyperbolic surface with one cusp and two orbifold singularities at the points $(0,1)$ and $(1/2,\sqrt{3}/2)$, the latter identified with $(-1/2,\sqrt{3}/2)$.

This surface can be decomposed into a compact part and a cusp as indicated in Figure 4. It can also be obtained directly as a quotient $X=\mathrm{PSL}(2,\mathbb{Z})\backslash\mathbb{H}$ as the above described domain is a fundamental domain of the $\mathrm{PSL}(2,\mathbb{Z})$ action, and the boundary components are identified using the maps $z\mapsto z+1$ and $z\mapsto-\frac{1}{z}$ (see, for example, [21] for an introduction).

While the metric $y^{-2}(dx^{2}+dy^{2})$ has constant curvature $-1$ we can consider a function $\phi$ which is compactly supported in the interior of the shaded region and change the metric by a conformal factor $\mathrm{e}^{\phi}$ to $\mathrm{e}^{\phi(x,y)}y^{-2}(dx^{2}+dy^{2})$. If

[TABLE]

this conformal transformation leaves the volume of $X$ unchanged. The surface equipped with this modified metric will in general have non-constant curvature.

7.3.2 Known properties of the spectrum

Case $\phi=0$:** In the case of constant curvature $-1$ ($\phi=0$) this surface is arithmetic. It has infinitely many embedded eigenvalues (the so-called Maass-eigenvalues) satisfying a Weyl law as shown by Selberg [32] using his trace formula. The scattering matrix $C(s)$ can be computed explicitly and equals**

[TABLE]

where $\Lambda(s)=\pi^{-\frac{s}{2}}\Gamma(\frac{s}{2})\zeta(s)$ and $\zeta(s)$ is the Riemann zeta function [15], [20] . Moreover, the Maass eigenvalues have been computed with great accuracy and verified by a rigorous algorithm [5], see also [4]. This surface therefore provides an excellent test for our method.

General case:** In case $\phi$ is non-constant (i.e. if curvature is non-constant) one expects at least some of the embedded eigenvalues to dissolve and become resonances ([30]). This has become known as the Sarnak-Phillips conjecture. Similarly the resonances will move away from the critical line.**

7.3.3 Numerical results

Case $\phi=0$:** Since the surface is symmetric with respect to the transformation $x\mapsto-x$ one can use symmetry reduction and consider the space of even and odd functions. These are functions on**

[TABLE]

satisfying either Dirichlet (odd functions) or Neumann (even functions) boundary conditions at the boundary. The spectrum on the space of odd functions is purely discrete, and there are no resonances. Several first eigenvalues on the space of odd functions, and their comparison with the results of ****[4]**** are presented in Table 1.

The results for the space of odd functions are below in Section 7.4.3.

The curve in moduli space:** We chose the the family of the conformal factors**

[TABLE]

with parameter $q$ in the interval $q\in[-2.,2.]$. One can sum over the group $\mathrm{PSL}(2,\mathbb{Z})$ to make this conformal factor a function on the surface. For numerical purposes the additional terms introduced in that way are however irrelevant as they are below working (double) precision. Note that the resulting family of metrics has constant curvature precisely at $q=0$. Moreover, the volume is constant along this curve in the moduli space of metrics. We computed the Neumann-to-Dirichlet data at $200$ points in the parameter interval on the cutoff surface with boundary at $a=2.2$, with $100$ discrete points on the boundary, as well as $600$ eigenvalues and their boundary data, and Fourier modes with $m$ between $-15$ and $15$. One can then trace the resonances as they move along the curve, see Video 1.

The same computation was performed using the family of conformal factors

[TABLE]

We omit the results which are very similar.

7.4 $B_{r}$. Artin’s billiard

7.4.1 Description of the surface

Given a positive parameter $r>\frac{1}{2}$ this surface can be obtained from the domain

[TABLE]

by gluing along the boundary as follows. The sides $x=-\frac{1}{2}$ and $x=\frac{1}{2}$ are identified by means of the parallel translation $x\mapsto x+1$. The circular arc $\{(x,y)\in\mathbb{H}\mid x^{2}+y^{2}=r^{2},0\leq x\leq\frac{1}{2}\}$ is identified with $\{(x,y)\in\mathbb{H}\mid x^{2}+y^{2}=r^{2},-\frac{1}{2}\leq x\leq 0\}$ using the map $x\mapsto-x$. This results in hyperbolic surface with one cusp and two conical singularities. Since the surface is symmetric with respect to the transformation $x\mapsto-x$ one can use symmetry reduction and consider the space of even and odd functions. These are functions on

[TABLE]

see Figure 7, satisfying Dirichlet or Neumann boundary conditions at the boundary. Since the spectrum on the space of odd functions is pure discrete we consider here only the spectrum on the subspace of even functions.

7.4.2 Known properties of the spectrum

There are various cases when this surface can be obtained as a quotient of $\mathbb{H}$ by a Hecke triangle group $G_{q}$, $q\geq 3$, namely when $r^{-1}=2\cos\frac{\pi}{q}$. In particular, for $r=\frac{1}{\sqrt{3}}$, $r=\frac{1}{\sqrt{2}}$, and $r=1$, the resulting surfaces are arithmetic and correspond to the surfaces obtained from the Hecke triangle groups $G_{q}$ in the cases $q=6,4,3$, respectively (see ****[26]****). These are the only arithmetic cases.

They have infinitely many embedded eigenvalues satisfying Weyl’s law. In each of the above three cases the scattering matrix, and hence the resonances, can be expressed explicitly in terms of the zeros of the Riemann $\zeta$-function. This has been done explicitly in ****[19]**** but the formulae can also be deduced using the known expressions for the congruence subgroups (****[16]**** and ****[20]****). The results are

[TABLE]

Because of the different choice of cusp-width our formulae differ by a factor $3^{1/2-s}$ and $2^{1/2-s}$ respectively from ****[19]**** in the first two cases. These surfaces were recently investigated in the context of the Sarnak-Phillips conjecture by Hillairet and Judge, who proved that for generic $r$ there are no eigenvalues (****[18]****) in the subspace of even functions.

7.4.3 Numerical results

We have computed resonances for $1000$ equidistant points in the parameter range $r\in[0.54,1.20]$ and tracked them, see Video 2, and also Figure 8 for selected resonances.

The resonances and embedded eigenvalues for the special arithmetic cases are shown in Figure 9.

We have also investigated the case $r=0.5001$, which is close to the limiting case $r=1/2$. Since in this case another layer of continuous spectrum appears one expects resonances to accumulate near the spectrum as $r\to\frac{1}{2}$. Apart from these resonances clustering around the spectrum we find stable ones that seem to converge to half the Riemann zeros, see Figure 9.

The numerical values of resonances and embedded eigenvalues for four special cases, and comparison with theoretical predictions and known results are collected in Tables 2 and 3.

7.5 $C_{\ell,\tau}$. Hyperbolic surfaces of genus one with one cusp

7.5.1 Description of the surface

The Teichmüller space for genus one surfaces of constant negative curvature $-1$ and one cusp is two dimensional and can be parameterised by the two Fenchel-Nielsen coordinates $\ell>0$ and $\tau\in[0,1)$. The parameter $\ell$ is the length of a primitive closed geodesic and the angle $\tau$ is the twist parameter along this geodesic. Given the above two parameters we have an explicit description of the corresponding surface of genus one with one cusp as follows.

For a fixed $\ell>0$, let $\alpha>0$ be the angle $\alpha=2\arctan(\tanh{\frac{\ell}{4}})=\arcsin(\tanh{\frac{\ell}{2}})$. Then the fundamental domain of the surface with Fenchel-Nielsen coordinates $(\ell,\tau)$ is the non-compact domain $D$ bounded by the following oriented geodesic arcs, $\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4},\gamma_{5},\gamma_{6},\gamma_{7}$,

[TABLE]

for brevity, we use the complex coordinate $x+\mathrm{i}y$ on $\mathbb{H}$.

Note that $\gamma_{2},\gamma_{4}$ and $\gamma_{6}$ are the images of $\gamma_{3},\gamma_{5}$ and $\gamma_{7}$ respectively under the reflection about the $y$-axis $x+\mathrm{i}y\mapsto-x+\mathrm{i}y$. Figure 10 depicts the fundamental domain decomposed into a compact part (darker shading) and a cusp (lighter shading). The surface $C_{\ell,\tau}$ is formed as follows. The infinite geodesic $\gamma_{6}$ is identified with $\gamma_{7}$ using the hyperbolic translation $z\mapsto z-1$. The geodesic arc $\gamma_{2}$ is identified with $\gamma_{3}$ using the hyperbolic motion along $\gamma_{1}$. Once these identifications are completed both $\gamma_{1}$ and $\gamma_{4}\cup\gamma_{5}$ become closed boundary geodesics of length $\ell$. These boundary components can be glued together as follows. First shift all points on $\gamma_{1}$ by $\tau\ell$. Then use hyperbolic translation along $\gamma_{2}$ and $\gamma_{3}$ to map the geodesic onto $\gamma_{4}\cup\gamma_{5}$. The resulting surface $C_{\ell,\tau}$ is a surface of genus one with one cusp such as the one depicted in Figure 2. The light-shaded region in Figure 10 gives a hyperbolic surface $M$ of genus one with horocyclic boundary.

In addition to the closed geodesics of length $\ell$, $C_{\ell,\tau}$ has another closed geodesics of length

[TABLE]

The two lengths are equal whenever

[TABLE]

or equivalently when $\ell=\ell^{*}(\tau)$ is the positive solution of

[TABLE]

7.5.2 Known properties of the spectrum

In general the Laplace operator on the surface will have simple continuous spectrum and may have embedded eigenvalues. The expectation is however, that these embedded are generically absent. There are several special cases for which the surface $C_{\ell,\tau}$ is symmetric, and therefore a symmetry reduction can be performed. We will single out and discuss several particular families.

7.5.3 Numerical results, case 1: $\tau=0$, varying $\ell$

With the twist parameter $\tau$ fixed the only remaining parameter is the length parameter $\ell$. Since in this case the twist is zero, the curve $\gamma_{2}$ becomes a closed simple geodesic of length

[TABLE]

on the resulting surface. As the function $\ell^{\prime}(\ell,0)$ is monotone decreasing in $\ell$, $\ell^{\prime}(\ell^{\prime}(\ell,0),0)\equiv\ell$, and also $\ell=\ell^{\prime}(\ell,0)$ when $\ell=\ell^{*}(0)=\operatorname{arccosh}(3)\approx 1.762747$, our parametrisation of $C_{\ell,0}$ is not unique: namely, the surfaces $C_{\ell,0}$ and $C_{\ell^{\prime}(\ell,0),0}$ are always isometric. Therefore it only makes sense to track resonances for $\ell\leq\ell^{*}(0)$. We have nevertheless analysed some special values of $\ell>\ell^{*}(0)$ to verify that our numerical results do not depend on the choice of parametrisation.

We have tracked the resonances in the interval $\ell\in[1.2,\operatorname{arccosh}(3)]$, see Video 3, and also Figure 11 for the trajectories traced by four selected resonances.

In this interval there are several special lengths. The numerically found resonances and the first ten embedded eigenvalues for these special lengths are in Tables 4 and 6, resp.

$\ell=\operatorname{arccosh}(2)\approx 1.316957$

In this special case the scattering matrix takes the form

[TABLE]

Indeed, one can conjugate the generators of the corresponding Fuchsian group into the following three matrices

[TABLE]

These clearly generate a subgroup $\Gamma$ of the arithmetic group $\widetilde{\Gamma}_{0}(6)$, that is the group generated by $\Gamma_{0}(6)$ together with its Atkin-Lehner (Fricke) involutions. Since the Atkin-Lehner involutions act transitively on the four cusps of $\Gamma_{0}(6)\backslash\mathbb{H}$ the domain $\widetilde{\Gamma}_{0}(6)\backslash\mathbb{H}$ has only one cusp. The scattering matrix for $\Gamma_{0}(N)$ has been computed in ****[16]**** and ****[20]****. If $N$ is square-free then according to Hejhal ****[16, Vol 2, p 536]**** the full scattering matrix equals

[TABLE]

where

[TABLE]

The vector $\left(\begin{matrix}1\\ 1\end{matrix}\right)$ is an eigenvector of $M_{q}(s)$ with eigenvalue $\frac{1+q^{1-s}}{1+q^{s}}$. The Atkin-Lehner involutions act transitively on the cusps for square-free $N$. Therefore, in this case, $\widetilde{\Gamma}_{0}(N)$, the group generated by the Atkin-Lehner involutions and $\Gamma_{0}(N)$, gives rise to a quotient with one cusp. The scattering matrix for $\widetilde{\Gamma}_{0}(N)$ must then be the restriction of the scattering matrix to functions invariant under the Atkin-Lehner involutions. The invariant functions correspond to the the span of the vector $\bigotimes\limits_{q\mid N}\left(\begin{matrix}1\\ 1\end{matrix}\right)$ in this representation of the scattering matrix. This vector is an eigenvector of $\bigotimes\limits_{q|N}M_{q}(s)$ with eigenvalue $\prod\limits_{q|N}\frac{1+q^{1-s}}{1+q^{s}}$. Summarising, in the case of square-free $N$, the scattering matrix for $\widetilde{\Gamma}_{0}(N)$ is given by

[TABLE]

This has also been obtained in ****[22, Lemma 5]****. Since $\widetilde{\Gamma}_{0}(6)$ acts on the fundamental domain for our group but leaves the cusp invariant, the scattering matrix for $\Gamma$ must be the same, apart from the extra factor $2^{1-2s}$ appearing because of the cusp width $2$. We refer to ****[36]**** for details of this argument.

Figure 13 shows the computed resonances for $\ell=\operatorname{arccosh}(2)$.

$\ell=\operatorname{arccosh}(3)\approx 1.762747$

In this case $\alpha=\dfrac{\pi}{4}$, so all the boundary arcs of $C_{\ell,0}$ have the same radius $\dfrac{\sqrt{2}}{8}$. Also, the length $\ell^{\prime}(\ell,0)$ of the second distinguished closed geodesic $\gamma_{2}$ coincides with $\ell$. We can carry out the following sequence of symmetry reductions, see Figure 12.

First of all the domain $X=C_{\ell,0}$ has a reflection symmetry $x\mapsto-x$. It follows that we have a natural decomposition of $L^{2}(X)$ into invariant subspaces for $\Delta$ consisting of even and odd functions. On the subspace of odd functions the spectrum of $\Delta$ is discrete since it is part of the space of cups forms, and the Eisenstein series are all even. The subspace of even functions corresponds to the space of functions on half of the domain $\{z=x+\mathrm{i}y\in X\mid x\geq 0\}$ satisfying Neumann boundary conditions along $\gamma_{6}$, $\gamma_{2}$ and along the $y$-axis part of the boundary of the resulting domain and periodic boundary conditions that identify $\gamma_{4}$ with the right half of $\gamma_{1}$. Now we have another symmetry $x\mapsto\dfrac{1}{4}-x$. Again the space of odd functions is contained in the space of cusp forms and the even part corresponds to considering the domain $\left\{z=x+\mathrm{i}y\in X\mid 0\leq x\leq\dfrac{1}{4}\right\}$ with Neumann boundary conditions everywhere along its boundary. The resulting domain has yet another symmetry $x\mapsto\dfrac{1}{4}-x$. The Laplace operator on the space of even functions on this domain corresponds the Laplace operator on the domain $\left\{z=x+\mathrm{i}y\in D\mid 0\leq x\leq\dfrac{1}{8}\right\}$ with Neumann boundary conditions everywhere along its boundary. Since the symmetry reduction of Artin’s billiard for $r=\dfrac{1}{\sqrt{2}}$ leads, after scaling by a factor $\dfrac{1}{4}$, to an isometric domain this shows that the continuous spectral subspace of $C_{\operatorname{arccosh}(3),0}$ and that of $B_{\frac{1}{\sqrt{2}}}$ are unitarily equivalent and the scattering matrices as well as the resonances coincide up to a scaling factor. One therefore has

[TABLE]

The above discussion also shows that the discrete spectrum consists of several parts, each belonging to mixed Dirichlet-Neumann problems on certain domains.

Figure 13 shows the computed resonances for $\ell=\operatorname{arccosh}(3)$.

In a similar way as before one can conjugate the generators of the corresponding Fuchsian group into

[TABLE]

which is a subgroup of the arithmetic group $\widetilde{\Gamma}_{0}(2)$. Equation (21) can therefore also be derived from (20) in the same way as before.

$\ell=\operatorname{arccosh}(5)\approx 2.29243$

This case is isometric to the case $\ell=\operatorname{arccosh}(2)$ since these two lengths are related by (17), see also the discussion following that formula.

$\ell=\operatorname{arccosh}(9)\approx 2.88727$

This case is isometric to the case $\ell=\left(\operatorname{arccosh}\left(\frac{3}{2}\right)\right)\approx 0.962424$ which lies outside our computed range. In this case the scattering matrix is given by

[TABLE]

The generators of the corresponding Fuchsian group can be conjugated to

[TABLE]

These therefore generate a subgroup $\Gamma$ of $\widetilde{\Gamma}_{0}(5)$. The surface $\widetilde{\Gamma}_{0}(5)\backslash\mathbb{H}$ has one cusp. The group $\widetilde{\Gamma}_{0}(5)$ acts on our surface $G\backslash\mathbb{H}$ and the action fixes the cusp. This implies that the scattering matrices of $G$ and of $\widetilde{\Gamma}_{0}(5)$ coincide modulo a possible factor coming from the normalisation of the cusp-width. In the same way as before, equation (20) (see also ****[2, equation (5)]****) gives the formula (22).

Figure 13 shows the computed resonances for $\ell=\operatorname{arccosh}(9)$.

7.5.4 Numerical results, case 2: $\tau=\frac{1}{2}$, varying $\ell$

We have tracked the resonances in the interval $\ell\in[1.12485,2.72485]$, see Video 4, and also Figure 14 for the trajectories traced by four selected resonances.

For twist $\tau=\frac{1}{2}$ we find the following special lengths. The numerically found resonances and the first ten embedded eigenvalues for these special lengths are in Tables 5 and 6, resp.

$\ell=2\operatorname{arccosh}\left(\frac{3}{2}\right)\approx 1.924847$

One can check by direct computation that for this particular $\ell$ the twist parameter $\tau=\frac{1}{2}$ is the unique twist for which the length of the second simple closed geodesic generating the fundamental group coincides with $\ell$. Any hyperbolic surface of genus one with one cusp that possesses two simple closed curves of that length that intersect in one point only will therefore be isometric to this surface. In particular, it is isometric to the arithmetic one punctured torus described by Cohn in ****[10]**** and by Gutzwiller in ****[12]****. The scattering matrix is known to be equal to

[TABLE]

where the extra factor $6^{1-2s}$ relative to ****[12]**** is because the cusp width in ****[12]**** was chosen to be $6$ rather than one. Its scattering resonances coincide with the one for the modular domain and are there directly related to the non-trivial zeros of the Riemann zeta function. Figure 15 shows the computed resonances for $\ell=2\operatorname{arccosh}\left(\frac{3}{2}\right)$. The form of the scattering matrix (23) can also be derived as follows. The generators of the Fuchsian group can be conjugated into

[TABLE]

which is a subgroup of $\mathrm{PSL}(2,{\mathbb{R}})$. Therefore, $\mathrm{PSL}(2,{\mathbb{R}})$ acts on our surface and fixes the cusp. Hence, the scattering matrix coincides with that of the modular domain up to a factor $6^{1-2s}$, since the generator $\left(\begin{matrix}1&6\\ 0&1\end{matrix}\right)$ yields a cusp of width $6$.

The special values close to the critical line and to the imaginary line are compared with theoretical prediction of (23) in Table 5.

We list some embedded eigenvalues for the twist parameter $\tau=\frac{1}{2}$ in Table 6. Note that some of the double eigenvalues coincide with those for the group $\Gamma^{3}$ from ****[33]****. Additionally, some embedded eigenvalues for the twist parameter $\tau=\frac{1}{2}$ have been computed in ****[23]****, however the authors have missed quite a few embedded eigenvalues in their list. They correctly identify two multiplicity two eigenvalues at 2.95648 and 4.51375, but do for example miss the multiplicity two eigenvalue at about 3.53606 and the simple eigenvalue at about 3.70339, cf. Table 6. We have performed a heuristic check using Weyl’s law and Turing’s method and our list appears to be complete.

$\ell=2\operatorname{arccosh}(2)\approx 2.6339157$

This case can be shown to be isometric to the surface $C_{\operatorname{arccosh}{3},0}$ by computing the generators, and our independent numerical results are in full agreement.

Figure 15 shows the computed resonances for $\ell=2\operatorname{arccosh}(2)$.

$\ell=2\operatorname{arccosh}(3)\approx 3.525494$

This surface is isometric to the one with $\ell=2\operatorname{arccosh}\left(\frac{3}{2}\right)$, and our independent numerical results are in full agreement.

7.5.5 Numerical results, case 3: $\ell=2\operatorname{arccosh}\left(\frac{3}{2}\right)\approx 1.924847$, varying $\tau$

The dynamics of resonances is shown in Video 5.

7.5.6 Numerical results, case 4: equal length geodesics, varying $\tau$ and $\ell=\ell^{*}(\tau)$

The dynamics of resonances is shown in Video 6.

Remark 7.1.

There are precisely four isomorphism classes of smooth arithmetic surfaces of genus one with one cusp (see [27], and also [35]). One can use the generators for the four surfaces $C_{\operatorname{arccosh}(2),0}$, $C_{\operatorname{arccosh}(3),0}$, $C_{\operatorname{arccosh}(9),0}$, $C_{2\operatorname{arccosh}\left(\frac{3}{2}\right),\frac{1}{2}}$ and identify them, using [35, Theorem 4.1], with the four known arithmetic cases. We discovered these special parameters by looking for values of the Fenchel-Nielsen parameters for which the resonances are all along critical lines. The numerical data and the location of the scattering poles then allowed us to conjecture formulae for the scattering matrix. We are very grateful to Andreas Strömbergsson, who saw the relation to $\widetilde{\Gamma}_{0}(N)$ from the formulae and was willing to share his expertise on Atkin-Lehner theory. This made it possible to provide proofs for the corresponding formulae (18) and (22).

7.6 $D$. The hyperbolic surface of genus zero with three cusps

7.6.1 Description of the surface

This surface is unique up to isometry and can be constructed as follows. Take the domain in the upper half space with boundary given by the the four curves $\gamma_{1},\gamma_{2},\gamma_{3}$ and $\gamma_{4}$ (see Figure 16). Here $\gamma_{1}$ and $\gamma_{2}$ are the two half-circles of radius $\frac{1}{2}$ centered at $z=\frac{1}{4}$ and $z=-\frac{1}{4}$ respectively. The curves $\gamma_{3}$ and $\gamma_{4}$ are the half lines perpendicular to the real axis originating from $z=-\frac{1}{2}$ and $z=\frac{1}{2}$ respectively. The surface is obtained by identifying $\gamma_{1}$ and $\gamma_{2}$, as well as $\gamma_{3}$ and $\gamma_{4}$. The three cusps are then located at $z=0$, $z=\frac{1}{2}$, and at infinity. The surface can also be obtained as a quotient of the upper half space by the subgroup $\Gamma_{0}(4)$ in $\mathrm{PSL}(2,{\mathbb{R}})$ which is generated by the matrices $\left(\begin{matrix}1&1\\ 0&1\end{matrix}\right)$ and $\left(\begin{matrix}1&0\\ 4&1\end{matrix}\right)$. The cusps at $z=0$ and $z=\frac{1}{2}$ can be removed from the surface by cutting along a horocycle (see Figure 16) and one then obtains two cusps. Each cusp is isometric to a standard cusp of some height. Removing the three cusps in this way one remains with a compact surface with three boundary components. This corresponds to the darker shaded region in Figure 16. Note that the points $z=\frac{1}{2}$ and $z=-\frac{1}{2}$ (belonging to the compactification of the hyperbolic plane) are identified.

7.6.2 Numerical results

Here the scattering matrix is a $3\times 3$ matrix, and our algorithm computes this reliably. In order to find the resonances we locate the zeros of the determinant of the scattering matrix and make use of the functional equation. Numerically we find that the resonances in this case are of multiplicity three at half the non-trivial roots of the Riemann zeta function, with additional resonances of multiplicity two at the points $\frac{\mathrm{i}k\pi}{\log{2}},\;k\in\mathbb{Z}\backslash\{0\}$, see Figure 17 and Table 7. Our root finding algorithm finds roots very close to one another in the case of multiplicities. It factors out an already detected root from the function and is therefore able to detect other roots close to the already found one. We can not distinguish numerically between true multiplicities and resonances that are very close to one another. What we find numerically is in excellent agreement with the known value of the scattering matrix for $\Gamma_{0}(4)$ ****[7]****:

[TABLE]

and the resonances $\frac{\mathrm{i}k\pi}{\log{2}},\;k\in\mathbb{Z}\backslash\{0\}$ are again due to the rational factors in the scattering matrix.

Acknowledgements

We are very grateful to Andreas Strömbergsson for important comments on a draft version of this paper (see Remark 7.1). We also thank an anonymous referee for many valuable comments and suggestions on an earlier version of the paper.

Appendix A Tables of resonances and eigenvalues