{'fraction_non_alphanumeric': 0.03602162162162162, 'fraction_numerical': 0.008648648648648649, 'mean_word_length': 4.426610348468849, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wave focusing by submerged islands and gravitational analogues

We study water waves propagating over a smooth obstacle in a fluid of varying depth, motivated by the observation that submerged islands in the ocean act as effective lenses that increase the amplitude and destructive power of tsunami waves near focal points. We show that islands of substantial height (compared to the water depth) lead to strong focusing in their immediate vicinity, and generate caustics of either cusp or butterfly type. We highlight similarities and differences with focusing of (high-frequency) gravitational waves by a neutron star. In the linear regime, the comparison is made precise through an effective-spacetime description of the island-fluid system. This description is then put to practical use: we identify caustics by solving the Raychaudhuri equation (a transport equation) along rays of the effective metric. Next, the island-fluid scattering processes are examined in detail (i.e.~deflection angle, phase shifts, scattering amplitudes) using numerical simulations and analytical techniques, including the eikonal approximation and its generalisation in the form of the Gaussian beam approximation. We show that the techniques capture the key features of the simulations. Finally, we extend the eikonal approximation to the dispersive regime, demonstrating that the essential features are robust in dispersive settings. This paves the way for future exploration in a controlled laboratory set-up.


I. INTRODUCTION
Wave propagation in inhomogeneous media is a vast and fascinating subject, in which even everyday systems give rise to a variety of beautiful phenomena such as coronae, rainbows and glories. These effects are associated with critical points and caustics at which the ray-optics treatments of wave propagation break down [1,2]. In this paper, we consider surface water waves propagating over a fluid with a varying depth, a problem common in oceanography and coastal engineering, and its associated critical phenomena.
In Ref. [3], Berry considered the propagation of water waves over smooth obstacles, showing that a submerged island acts as a lens for surface waves. The amplification from lensing multiplies the devastating power of tsunami waves, and this amplification is greatest near the focal point at the cusp of the caustic. Shallow islands create focal points that are located relatively far from the submerged island. In this region, analytical solutions for the wave profile near the focus point may be found by means of the paraxial approximation [3], where all rays contributing to the focusing are assumed to propagate parallel to one another.
In this work, we move beyond the case of shallow submerged islands to consider more substantial obstacles with a typical relative height h 0 /h ∞ ≈ 1/2, and focal points that are close to, or even on top of, the submerged island. The paraxial approximation is invalid in this case, and we turn instead to methods including geometric optics, numerical simulations, and the Gaussian beam approximation.
The key features of the process are illustrated in Fig. 1. A wave approaches the island; as the water depth decreases, its speed of propagation decreases also, causing focussing and an associated increase in the wave amplitude (Fig. 1A). In passing over the island, the wavefront slows down until, at the focal point, it develops a pair of cusps. These cusps in the wavefront move outwards along a caustic wedge (red dashed line in Fig. 1B). Inside the caustic wedge, the wave now has leading and secondary fronts. At the focal point, there is substantial amplification; just inside the caustic wedge we anticipate the amplification and diffraction effects commonly associated with rainbows [2].
Remarkably, the focussing of water waves by a submerged island is qualitatively similar, in many respects, to the focussing of gravitational waves by an astrophysical compact object, such as a star or neutron star [4,5]. The gravitational potential well of a compact object causes gravitational time dilation, which slows down the wavefront, and consequently generates a focal point and cusp caustic. For our Sun, this focal point is circa 550 a.u. (astronomical units) distant from its centre. On the other hand, for a neutron star the focal point is expected to be either inside the star, or close to its crust [4,5]. In this picture, one may think of the Sun as analogous to the shallow submerged island considered by Berry [3], and a neutron star as analogous to the more substantial obstacles considered here.
This analogy extends further. Under certain physical assumptions, water waves propagating in a Neutron star A C B FIG. 1. Qualitative description of wave scattering with a focusing object. Panel A depicts a schematic of the underwater island along the axis y = 0. The blue curves represent the profile of the waves as passing over the submerged obstacle. The wavelength shortens as the water depth decreases; this is followed by an increase in amplitude. Panel B, reproduced from [4], shows rays and waves in a neutron star spacetime. The blue and purple lines depict the geometrical wavefronts which are retarded by the strong gravitational field of the neutron star. The red dashed curve is the caustic traced out by the cusps in the wavefront. Panel C shows the behavior of rays (in black) and geometrical wavefronts (in red) propagating over the submerged island. Note the striking similarities between panel B and C. background fluid flow are described by the equation for a massless scalar field on an effective curved space-time. This observation is the basis of the field of analogue gravity [6,7], which offers laboratory platforms to experimentally investigate fundamental phenomena occurring in curved space-times, such as Hawking radiation [8][9][10], superradiance [11] or cosmological particle production [12,13]. Here, we show that an underwater island will result in an effective space-time similar to the one of a dense astrophysical object, such as a neutron star. Moreover, the analogy offers a set of tools that can be applied to both the hydrodynamical system as well the gravitational one. In addition, the analogy offers the possibility to experimentally investigate processes occurring around compact bodies and to test mathematical methods employed to described them.
The paper is structured as follow. In Section II, we derive the governing equation for linear shallow water waves propagating over an submerged obstacle and show that it can be seen as the wave equation for massless scalar fields on an effective curved spacetime. We then study the geodesics of this analogue space-time and show that our system exhibits rainbow scattering and focusing at the origin of the cusp caustic. In Section III, we investigate the propagation of waves in our system, both numerically and through semi-analytical methods. In particular, we explore the high-frequency regime using the eikonal (or geometrical) approximation as well as the Gaussian beam approximation which allows us to evolve our geometrical wave through the caustic. In order to pave the way for experimental investigations of the effects described, we consider the impact of dispersive effects by means of the geometrical ray approximation. This is done in Section IV. Finally, Section V concludes our work and discusses possible experimental realisation of our systems. Two appendices contains further technical details on the analytical solution in the case of a parabolic island and on the born approximation to estimate the scattering cross-section.

A. Governing equations
We consider a stationary, irrotational and inviscid fluid, described by a scalar velocity potential φ such that the velocity of the fluid is given by v = ∇φ. The surface of the fluid is a 2-dimensional surface with coordinates x = (x, y). An obstacle is placed under the surface such that the fluid depth varies in space. The obstacle is assumed to have a shallow gradient such that derivatives of the fluid depth can be neglected. Small deformations, δh, of the fluid interface are related to perturbations of the velocity potential φ( x) via ∂ t cosh(ih( ∇ 2D ))φ = −gδh. Here ∇ 2D is the two-dimensional gradient ∇ 2D = e x ∂ x + e y ∂ y (we will omit the subscript in the following). Under these assumptions and neglecting non-linearities, the velocity potential obeys the following wave equation [14]: When considering shallow surface waves, i.e., wavelengths significantly longer than the fluid depth, this reduces to the non-dispersive wave equation where ∆ ≡ ∇ · ∇, and c(x) is the propagation speed of the wave given by In the following, we will keep the discussion as general as possible by not specifying a particular profile for the underwater island. Later, to obtain and visualize quantitative results, we will specify the water depth as with the default values h ∞ = 0.04 m, h 0 = 0.02 m, r 0 = 0.3 m and n = 4. These parameters represent an underwater island which can be set up in existing water tanks [11]. This particular class of profile also allows for application of specific analytical techniques, such as the Born approximation. In Appendix C we present a full analytical calculation for a parabolic island,

B. Rays and the eikonal approximation
For a general submerged island profile, the wave equation (2) is in general not solvable in closed form. Nevertheless, approximation methods can be applied to gain insight and to make predictions. A simple but effective method is the eikonal approximation, in which the wave is described in terms of a coherent collection of rays [15]. As well as being applicable in fluid-mechanical systems, the eikonal approximation has also been successfully applied in analogue-gravity settings in order to describe in wave-vortex interaction scenario [16,17] and to describe light-ring mode emission in such systems [18,19].
The eikonal approximation rests on the assumption that the phase of the wave varies rapidly in comparison with its amplitude. One seeks solutions to the wave equation Eq. (2) of the form where S and A are the local phase and amplitude of the wave, and is an order-counting parameter. Expanding the wave equation in a hierarchical fashion in powers of −1 yields a system of equations for the local phase and amplitude as asymptotic series in . The leading order term in this expansion is the eikonal equation, that determines the leading term in the expansion of the phase, S 0 . The lines of constant phase S 0 are the wavefronts.
The eikonal equation has a natural interpretation as a Hamilton-Jacobi equation, corresponding to the Hamiltonian H( x, t; k, ω) that determines the rays of the system (see [15,20] for detailed discussions). The Hamiltonian is obtained by making the substitution (∂ t S 0 , ∇S 0 ) → (−ω, k) in Eq. (7), to obtain The condition H = 0 is nothing other than the dispersion relation.
The rays are parametrized curves ( x(τ ), t(τ )) that are the solutions of Hamilton's equations. In the scattering scenario, we consider a congruence of rays that impinge from infinity, x → ∞, which are asymptotically parallel and in phase, such that the incoming wavefronts are straight lines orthogonal to the rays. Practically, we solve the Hamilton's equation in Cartesian coordinates to find the trajectories of the rays x(τ ) as well as the variation of the momentum k( x(τ )) along the rays. In a time-independent system the associated frequency ω is constant. Explicitly, we solve the following system of equationsẋ where the dot represents a derivative with respect to the ray parameter τ (i.e.ẋ = dx/dτ , etc.). Once the ray and its associated momentum have been computed numerically, we can reconstruct the eikonal phase along the trajectory. This is done using the definition of the momentum ∇S 0 = k. The eikonal wavefronts are then found as constant phase lines across a congruence of rays.
In the case of an axisymmetric submerged island, it is convenient to work with cylindrical coordinates (r, θ), in which the Hamiltonian becomes where (ω, m) = ( dt dτ , r 2 c 2 (r) dθ dτ ) represent the frequency and azimuthal number which are constants of motion, and k r is conjugate to r. Figure 2 depicts a congruence of rays incident on the submerged island from right infinity.

The amplitude and the caustic
To find the leading-order term in the amplitude of the wave, A 0 , one extends the hierarchical expansions to sub-leading order in . For a stationary system, of interest here, this equation is FIG. 2. Rays passing over an island form a caustic. A congruence of rays (grey curves) approaches from the right. The solid black line corresponds to a ray passing inside the rainbow wedge, and its associated impact parameter b and deflection angle Θ(b). The dashed black line shows the rainbow ray, with impact parameter b r , which is maximally deflected (see also Fig. 3). The red curve highlights the caustic where neighbouring rays cross; this is located using the Raychauduri equation (24). The inset panel reveals the structure of the butterfly caustic.
This equation implies that the flux of wave action A 2 0 ∇S 0 is conserved along a tube of rays. Using the Hamilton-Jacobi equation and the stationarity condition, the amplitude along a ray is where d( x(τ )) is the cross-section of the tube of rays. The cross section of the tube vanishes at points where neighbouring rays meet; here the leading-order amplitude A 0 diverges, and the asymptotic expansion breaks down. The caustic is the set of all such points.

Deflection, geometrical scattering and the rainbow angle
The ray trajectories provide insight into the scattering of waves with the submerged island. Of particular interest is the deflection angle, Θ(b), for an incident ray with impact parameter b = mc ∞ /ω (see Fig. 2). The deflection angle is obtained by integrating dθ/dr using Hamilton's equation to get where r 0 is a turning point satisfying The deflection angle Θ(b) is shown in Fig. 3. The blue curve corresponds to the non-dispersive case, and is given by Eq. (14). The purple, yellow, and red curves correspond to the deflection angle in the dispersive case at 2, 3 and 4 Hz respectively (see IV for a discussion of dispersive effects). In every case the deflection angle has an extremum at b = b r , defining the rainbow angle θ r = |Θ(b r )|. Note that dispersion leads to smaller deflection.
The geometrical scattering length is defined as the density of rays passing into the unit angle dθ. Rays with impact parameter b are deflected by an angle Θ(b) while the rays with impact parameter b + db are deflected by an angle Θ(b) + (dΘ/db)db. Consequently, the geometrical scattering length is As shown in Fig. 3 and Fig. 2, there is a ray with impact parameter b r that is maximally deflected, such that Θ (b r ) = 0 (and Θ (b r ) > 0). This is known as a rainbow ray [2]. This ray determines the rainbow angle of the caustic wedge, θ r = |Θ(b r )|. Formally, the geometrical scattering length diverges at this angle, indicating the breakdown of the method, as expected at any caustic feature.
C. The effective spacetime and the Raychauduri equation

The effective spacetime
In accordance with the analogue-gravity paradigm, Eq. (2) can be rewritten as a Klein-Gordon equation for a scalar field propagating on a curved spacetime as [6], Here x µ are spacetime coordinates (with x 0 = t and x i = ( x) i ), and g is the determinant of the covariant metric tensor g µν whose inverse is given in contravariant form by In Eq. (17), we have used the Einstein summation convention which implies that repeated indices are summed over. Indices are lowered with g µν and raised with g µν . It is important to note that Eq. (17) is nothing more than a formal rewriting of the wave equation (2). The effective spacetime description is a tool that offers a new point of view on the problem of wave propagation in inhomogeneous media.

Geometric optics and the null geodesics
In a similar fashion to Sec. II B, one may now seek an asymptotic solution to Eq. (17) in the short-wavelength regime. Inserting the ansatz into (17) yields at leading order the eikonal equation where k µ ≡ ∇ µ S 0 , and ∇ µ denotes the covariant derivative on the effective spacetime. By taking a derivative of the eikonal, one obtains the geodesic equation The integral curves x µ (λ) that satisfy dx µ /dλ = k µ are geodesics of the effective spacetime; and these geodesics are null by virtue of (20). In summary, the rays of the eikonal of the physical system of Sec. II B correspond with the null geodesics in the effective spacetime. The geodesics may be found from the Hamiltonian H(x µ , k µ ) = 1 2 g µν k µ k ν , which derives from interpreting the eikonal equation (20) as the corresponding Hamilton-Jacobi equation. By extending the expansion to sub-leading order, one finds that the amplitude A 0 is governed by a transport equation

Focusing and the Raychaudhuri's equation
We shall now show that the effective spacetime description yields a practical method for computing the caustic, that is, the set of points where neighbouring rays meet. This is done by solving a transport equation for the expansion scalar ϑ associated with a bundle (or congruence) of geodesics, that quantifies the way neighbouring geodesics converge (or diverge). This transport equation is known as Raychaudhuri's equation. Following [21], the Raychaudhuri equation for this system is where k µ is tangent to a null geodesic, and R µν is the Ricci tensor, a tensor which describes the local curvature of the effective spacetime (see Appendix A for further details). For an arbitrary submerged island, the Ricci tensor is given (see Appendix A) in Cartesian coordinates by In the case of an axisymmetric island, one can express the coefficient of the Ricci tensor in polar coordinates as To find the caustic from the Raychauduri equation in practice, we seek points along geodesics at which ϑ → −∞. After the change of variables ϑ = u u , where u = u(λ) is some function and u ≡ du dλ , finding the caustic reduces to finding points at which u is zero (providing that its derivative is well behaved). The Raychaudhuri equation becomes

D. Characterising the caustic
The caustic separates the (x, y)-plane in two regions with distinct wave profiles (see Fig. 2). In the region outside of the caustic, each point is connected to the initial wavefront by a single ray, and the phase function is single-valued. Conversely, in the region inside the caustic, each point is connected to the initial wavefront by more than one ray, which results in a multi-valued phase function, and constructive/destructive interference effects.
In the case of interest here, the caustic belongs to the butterfly caustic class. Inside the primary region of the butterfly caustic, that is between the rainbow rays, but outside of the star-shaped region, each point is reach by 3 different rays; while inside the star-shape region itself (Fig. 2, inset), points may be reached by at most 5 different rays [22]. The multivaluedness of the phase means that the geometrical wavefronts are not smooth inside the caustic; rather, they develop cusps and self-intersection as they move through the caustic. Figure 4 shows the caustic structure for several different island profiles, in the family of Eq. 4. Changing the parameter n changes the shape of the caustic. In particular, the caustic changes from a single cusp to a butterfly as n increases. We now seek to understand this phenomenon at a deeper level.
The butterfly caustic observed in our system can be understood as a projection of an hypersurface in a four-dimensional space of control parameters onto the two-dimensional (x, y) plane of the water surface. The particular symmetry of the underwater island fixes one of the four parameters (C 3 = 0, see below). The caustic is then found by taking a two-dimensional slice through a three-dimensional space. This can be explored by varying the island profile. In particular, by varying the parameter n in Eq. (4), governing the fall off of the island, we shift the height of the slice through the surface, revealing the transition from a cusp caustic to the characteristic star shape of the butterfly caustic. Figure 5 illustrates this point by showing the caustic shape as a function of n, which generates a 2D surface in a 3D space.
In catastrophe theory, caustics are found from a generating function, φ(s; C), where s = (s i ) are state variables and C = (C i ) are control parameters. Essentially, the variables (s i ) parametrize the rays going from an initial surface to a point characterised by the variables (C i ). Note that the control parameters will include the coordinates of the end point but may also contain other parameters governing the media in which rays propagate. The caustic is then found by looking for singularities of the gradient map, from s-space to C-space, defined by the condition vanishes, that is One of the main result of catastrophe theory is that caustics are divided into equivalence classes, and all elements of a class can be deformed into one another smoothly. Each equivalence class is described by a generating function, in the form of a standard polynomial. The butterfly caustic has a single state parameter, s, and its standard polynomial is [22] From Eqs. (25) and (26), we can express C 1 and C 2 as functions of (s, C 3 , C 4 ). For a fixed value of C 3 , and C 4 , the caustic is found as a curve parametrised by s in the (C 1 , C 2 ) plane. The transition from the cusp to the butterfly caustic can be seen by fixing C 3 = 0 and varying C 4 for some positive to some negative value. For C 4 > 0, the caustic in the (C 1 , C 2 ) plane is a cusp caustic, while for C 4 < 0 it is a butterfly caustic. The chrysalis point, from which the butterfly emerges, is at C 4 = 0.
From the previous description we can relate (C 1 , C 2 ) to the Cartesian coordinates (x, y) and C 4 to the parameter n governing the island profile. In particular, C 4 will be a monotonically decreasing function of n, such that C 4 (n < n chrysalis ) > 0 and C 4 (n > n chrysalis ) < 0.

E. The Gaussian Beam Approximation
As we have seen, approximating the amplitude of the wave in the vicinity of caustics is impossible by means of the geometrical approximation, which predicts an infinite amplitude at the focus point. One can deal with such infinities in some cases by modifying the ray method and introducing special functions. This is the case in the Airy treatment of the rainbow scattering for example, which can be applied to a variety of caustic problems in one spatial dimension. This method is based on obtaining an approximate wave equation in the vicinity of the caustic. Another method of curing the singularities of the geometrical description is to reintroduce some "wave flesh onto the classical bones" [23]. Instead of considering the wave as a congruence of rays, one can consider it as being a collection of beams. Each beam is centred on the underlying rays, and the amplitude is not located exactly on the ray but on a Gaussian profile transverse to the ray. This Gaussian beam can spread and focus depending on the inhomogeneities of the media. For this reason, the method is called the Gaussian beam approximation [24,25]. More precisely, the wave equation (2) is reduced to a local wave equation around each ray. The essential ingredients of the Gaussian beam approximation are outlined below.
First, define an orthonormal coordinate basis ( e s , e n ) adapted to each ray, where e s is tangent to the ray and e n is transverse, such that e s . e n = 0. A point in the vicinity of the ray is located with coordinate (s, n) where s is the arc-length along the ray and n the displacement along e n . In this coordinate system, we look for solution to the wave amplitude of the Gaussian form where Γ(s) = p(s) q(s) and (p(s), q(s)) are unknown functions which obey the following system of differential equations: Here we are interested in the amplitude of the ray passing through the caustic, which we take to be located on the y = 0 axis. For such central ray, the arc-length, s, is precisely equal to its displacement from the initial point hence s = x 0 − x. Now, along the central ray, we take the normal to be n = (0, 1), and therefore the normal coordinates can be substituted directly as y, because dy dn = 1. In Cartesian coordinates (p, q) now obeys the following differential equations In the traditional ray method, q, p and Γ are real and q vanishes at the caustic points which causes the amplitude to diverge. In order to avoid the singularity, p and q must be complex valued, and we can think of the Gaussian beam as a collection of complex rays [26]. Since the coefficients in Eq. (30) are real, it follows that the initial condition must be complex. It turns out that the evolution of the Gaussian is characterised by a single complex parameter, which can be represented via two real numbers. Those two real numbers can be interpreted as the initial effective half-width of the Gaussian beam and the distance between the initial position of the beam and the point where the half-width is minimum (i.e., the distance between the initial point and the caustic) [24]. Since we are only interested in the amplification at the caustic, the initial half-width can be normalised and we are left with a single parameter controlling the evolution of the Gaussian beam. To finalise setting-up the initial condition, we choose the minimum half-width of the Gaussian beam to be located at the focus point determined using the eikonal approximation. As shown in Sec. II C, surface waves propagating over the submerged island obey a massless scalar field equation in a curved spacetime, Eq. (17). We consider here a monochromatic solution, φ = Re e iωt φ ω ( x), which can be constructed from a sum over partial waves, as follows: Here, the radial profiles φ ω,m (r) satisfy the following radial equation Far from the obstacle, the depth of the water tends to a constant, h(r) → h ∞ and the wave equation (32) takes the form with c 2 ∞ = ω 2 /gh ∞ , which admits the solutions with c ∞ k ∞ = ω. A out/in ω,m represent the (complex) amplitude of radially outgoing/ingoing modes.
We are interested in the scattering of an incident monochromatic plane wave. Hence, we seek solutions which far from the island, are a superposition of a plane wave propagating towards x → −∞ and a radially-outgoing scattered component, The function f ω (θ) is the scattering amplitude, which can be expressed as a partial-wave sum, The rotational invariance of the island profile implies that the phase shifts are symmetric, δ m = δ −m , and we can rewrite the expansion of the scattering amplitude as a m cos(mθ), with a m = 1 2iπk By decomposing the incoming plane wave onto the azimuthal basis as, e ikx = m i m J m (kr)e imθ , and using the asymptotic form at infinity of the Bessel functions of the first kind, J m (kr), the phase shifts are found in terms of the ingoing/outgoing amplitudes in Eq. (34), In the case of a submerged obstacle, the phase shifts are real (δ m ∈ R) since there is no absorption or amplification and hence |A out /A in | = 1.

Geodesic phase shifts
As we have just seen, in the high frequency/large m limit, waves can be described as a collection of rays. Hence, we can express in this limit the phase shifts in terms of the properties of the rays. In the large-m regime, there is the following correspondence between the deflection angle (14) and the phase shifts, We show in Appendix B that, in the large m limit, and for the general shape of underwater island we consider in Eq. (4), the deflection function behaves as Θ = O(m −n ), where n governs the fall off of the obstacle at infinity. This implies that, for n = 1, the phase shifts behaves as δ m = O(m n−1 ) and will tend to zero as m → ∞. Conversely, in the case of gravitational body of mass M , the deflection angle falls off as ∼ 4M/b, which leads to a logarithmic behavior of for the phase shifts, and a Rutherford-like (θ −4 ) divergence in the scattering cross section in the forward direction.

A. Numerical Method
We solved the 1D wave equation (32) with a numerical method to obtain the radial profiles of the partial waves φ ω,m (r). We start with appropriate initial conditions at r = 1, by selecting initial values φ ω,m ( ) and φ ω,m ( ) which determined from the solution to the approximate wave equation at small radius, This equation admits the regular solution φ ω,m ∼ r 1/2+|m| . We then integrate numerically Eq. (32) using the NDSolve solver in Mathematica into the far-field. The ingoing and outgoing amplitudes A in/out ω,m in Eq. (34) are obtained by matching the numerical solution onto the generalized series solutions, The coefficients b j are obtained from the recurrence relation that is obtained by substituting the generalized series solution (41) into Eq. (32) and expanding order-by-order in powers of 1/r. Finally, we calculate the phase shifts δ m by inserting the ingoing and outgoing amplitudes in Eq. (38). From the set of radial profiles, we reconstruct a monochromatic plane wave incident on the submerged island by matching the relative phases of each mode and requiring that it satisfies Eq. (35) in the far field. Explicitly, the incident plane wave is reconstructed by evaluating the following sum When reconstructing the incident plane wave, as well as the scattering cross-section from the partial wave expansion, we numerically compute the radial profile of a finite number of m modes, up to m ≤ M , and then use Bessel functions, which are the solutions to the wave equation for large m, to complete the sum up to second cut-off, M 2, which is greater than M . The value of M is chosen such that the phase shift obtained numerically satisfy |1 − e 2iδm | < 10 −3 , and M 2 such that the Bessel function are negligible in the region of interest.
B. Results Figure 6 shows the phase shifts obtained numerically, and the comparison with values from Born approximation in Appendix B, for different frequencies: f = 2Hz, f = 4Hz, and f = 6Hz. We can see the qualitative agreement between the approximate estimate and the numerical value of the phase-shifts. Figure 7 shows an incident plane wave of frequency f = 4Hz encountering an underwater island. The focusing of the incident plane wave and the resulting increase in amplitude is clearly visible in  Fig. 7, and for this frequency the amplification factor is ∼ 3. Interference effects are visible within the rainbow wedge. The eikonal wavefronts (red lines) closely track the undulations in the full numerical solution. The eikonal wavefronts were found by starting from an initial wavefront on the right hand side of Fig. 7 and then ray-tracing using Hamilton's equations. The island parameters used in the partial wave expansion and the eikonal methods are identical and are given in Eq. (4). Figure 8 shows the scattering cross section of an incident plane wave on the submerged obstacle for the frequencies f = 2Hz,f = 4Hz, f = 6Hz, and f = 8Hz. The rainbow angle θ r is indicated by the vertical dashed line. Inside the rainbow angle (θ < θ r ), the scattering cross section has an oscillatory behaviour, whereas outside the angle it has an exponential fall-off, indicatiing a shadow region. This is the typical behaviour in rainbow scattering [1,4,5]. The scattering cross section shown in Fig. 8 was computed using the partial wave expansion given by Eq. (37) and summed over 0 < m < 300. The phase shift are obtained numerically for m ≤ 80, and via the Born approximation for m > 80. The noise visible in the scattering cross section of the high frequency for large angle arises from the discrepancy between the numerical phase-shift and the Born approximation for m ∼ 80. The partial sum is computed up to N = 300 with the first 80 modes solved numerically and the contribution for the higher m is estimated using the Born approximation. The vertical dashed curve represent the location of the rainbow angle θ r . For angles θ < θ r , the scattering cross section has an oscillatory behavior; for θ > θ r it decays exponentially.
C. Amplitude at the caustic Figure 9 shows the wave profile along the y = 0 axis for a wave with angular frequency ω = 8π. The blue curve shows the numerical profile obtained by the partial wave expansion method, and the red curve shows the wave profile obtained analytically after approximating the underwater island with a finite size parabola (see Appendix.C). The purple/yellow envelope depicts the amplitude of the wave in the Gaussian beam/eikonal approximation, respectively. The eikonal approximation of the amplitude diverges at the focus point, whereas the Gaussian beam remains finite and provide an overall good description of the profile of the wave and its amplification. to the original wave amplitude. Second, there is a smaller decrease in amplitude as the waves approach the obstacle. This decrease in amplitude is also present in the eikonal approximation, and can be understood from Eq. (13); it is due to the fact that the propagation speed decreases faster than distance between neighbouring rays. Finally, we observe a dip in the amplitude downstream from the caustic cusp. The minimum (near x = −1 in Fig. 9) is seen in the numerical simulation, as well as in the parabola approximation. Qualitatively, this dip is due to destructive interference between multiple wavefronts inside the rainbow wedge. It is not present in the Gaussian beam profile, nor in the eikonal profile, because the contribution from secondary rays has not been included.

A. From geodesics to rays
Thus far, we have focused on the case of water waves with a linear dispersion relation, which allowed us to establish a precise analogy between our system and that of gravitational waves propagating through compact bodies. Of course, this linear dispersion relation is an approximation (one amongst many) to the physical system, and it is well known that surface water waves are subject to dispersive effects. Considering the impact of dispersion in a general setting is beyond the scope of this paper; however, if we restrict attention to the behaviour of high-frequency waves, then by means of the eikonal approximation, it is possible to extend the analysis of Sec. II B to the dispersive regime.
It was shown in Ref. [16] that the trajectories of the "particles" that make up the eikonal waves can be obtained from the dispersive Hamiltonian given by with k = k 2 x + k 2 y . The eikonal trajectories, which previously corresponded to geodesics of the effective space-time, are now called rays, since one cannot define a single effective metric for which all the rays would be the geodesics. We can find the rays, and thus the eikonal wavefronts by solving Hamilton's equations, as before. Figure 10 shows the eikonal wavefronts for the dispersive and non-dispersive systems, for an incoming wave of frequency f = 3Hz. We may draw two general observations from this figure. First, that the qualitative behaviour between the two regimes is similar, that is, we still observe the presence of a rainbow ray (see also Fig. 3), and a time delay leading to a focusing of the wave. Second, that the dispersion, or more precisely the subluminal dispersion of water waves, leads to a smaller deflection of the rays by the underwater island. This behaviour can be understood in terms of geometrical optics and Snell's law n D θ D = n L θ L (for small angles), where θ D,L is the refracted angle and n D,L the refractive index for the dispersive and linear regime respectively. Since the refractive index is inversely proportional to the wave speed, we have that n D > n L which implies that θ D < θ L . A consequence is that, in the dispersive case, the rainbow angle is narrower, and the focal cusp is further from the scattering centre.

V. CONCLUSION
In the preceding sections, we have analysed the scattering of surface water waves passing over submerged obstacles. As Berry found [3], underwater islands act as lenses which focus and amplify Dispersive effects lead to a narrowing of the rainbow wedge, and a cusp focal point that is located further away from the obstacle. incident waves. In the short-wavelength limit, this leads to the formation of caustics near which the power of the wave (or of the tsunami) is enhanced. We have shown here that a submerged island of significant height relative to the water depth will generate strong focusing immediately downstream (see Fig. 2). Moreover, we find that for substantial islands, the caustics formed can be of cusp or butterfly type, with the latter emerging for steep-sided islands (see Fig. 4).
We simulated the wave scattering process in the linear regime using a partial-wave decomposition into m-modes (Sec. II F). The key features of the scattering patterns are described by several semi-analytical techniques. The eikonal approximation yields wavefronts that are a good match to those seen in the simulation (see Fig. 7). However, the eikonal amplitude diverges at the caustic, and the eikonal short-wavelength assumption breaks down where neighbouring rays intersect. To circumvent this issue, we adopted the Gaussian beam approximation, which yields a valid prediction for the wave profile across the focal point (the cusp). In Fig. 9, we observed robust qualitative agreement between the Gaussian beam approximation and the numerical simulation. To augment these approximation schemes, we also derived an exact analytical expression for waves propagating over an island of parabolic shape (see Appendix C).
Remarkably, this hydrodynamical system has features in common with strong-gravity systems in astrophysics. In the absence of dispersion and dissipation, the equations governing water waves propagating over a fluid with varying height are mathematically equivalent to those governing a massless scalar field in an effective 2D spacetime (see Sec. II C). The effective spacetime is qualitatively similar to (a 2D slice through) the spacetime of a massive, dense body, such as a neutron star. Consequently, the wave scattering patterns possess similar features. For example, in both cases we expect rainbow scattering, which is the manifestation of constructive/destructive interference effects on the inside of the rainbow wedge associated with a maximally-deflected ray. We have characterised this effect in our system and observed its presence in our numerical simulation (see Fig. 8). By comparison with the results of Refs. [4,5] in the gravitational context, we find that water waves passing over a submerged island closely resemble gravitational waves focussed by the spacetime curvature of a neutron star; with the most important differences arising from the dimensionality of the systems (2D vs 3D).
The analogue-gravity description is more than a mathematical curiousity, however. We showed that by solving the Raychaudhuri equation in the effective spacetime, we can locate and characterise the caustic (see Sec. II D). The Raychaudhuri equation is most familiar in relativistic context in the context of (Penrose-Hawking) singularity theorems. It is a transport equation that describes the rate of change of the cross-sectional area of a congruence of rays. The caustic is the set of points where that cross-sectional area passes through zero.
Finally, we considered the effect of dispersion on the geometrical picture of wave propagation (Sec. IV). We established that the key features of the scattering process are modified but not eliminated by dispersive effects. In Fig. 10, we observe a similar convergence of rays leading to focusing downstream of the obstacle, and the rainbow scattering is still present; the main difference here is the narrowing of the rainbow angle. Consideration of dispersive effects is an important step towards an experimental realisation of a neutron-star analogue in a wavetank, and direct measurements of the predicted phenomena.
In the results presented in this work, we have selected physical parameters (h ∞ , h 0 r 0 , n and f 0 ) that we anticipate are relevant to (future) wavetank experiments. Any implementation in physical media would offer the possibility to investigate regimes which are not fully modelled here. For example, obtaining an accurate description of the focusing of dispersive waves in the vicinity of the caustic remains an open problem. One could certainly imagine that the Gaussian beam approximation used here could be extended into the dispersive regime, but such an extension lies beyond the scope of this paper. Another interesting regime to explore, experimentally and theoretically, is the one arising from Berry's original proposal, namely the formation of non-linear waves as a result of the focusing process. As the wave focuses, its amplitude will locally increase. This amplification may result in a breakdown of the linear description studied in this paper and one may expect non-linear processes to come into play, such as the generation of solitary or rogue waves [27][28][29].
From the (inverse) metric in Eq. (18), we can explicitly compute the Ricci tensors, as follows. Since the metric is time-independent it implies that all time components of the Christoffel symbols vanish. Defining f = 1 2h(r)r dh dr , then the connection form A µ ν = Γ µ ρν dx ρ in (t, x, y) coordinates (suppressing the time components) is: The curvature form is F = dA + A ∧ A, and A ∧ A vanishes in our case, leaving F = dA (d is the usual exterior derivative, and the wedge product of two matrices is matrix multiplication with wedge products on each component). Consequently, A calculation, best illustrated with index notation, confirms that Performing a contraction of the µ and ρ indices yields the Ricci tensor in the form R µν dx µ dx ν = (∂ x (xf ) + ∂ y (yf ))dx 2 + (∂ x (xf ) + ∂ y (yf ))dy 2 .
By comparing this form of the solution with the requirement that the wave be a superposition of a plane wave and a purely outgoing waves, we can express the phase shift in the Born approximation as: Since I is real, we can see that |e 2iδ B m | = 1. Using the underwater island profile (4), we can expand the propagation speed (3) as a power series in 1/x: where A = c 2 0 /c 2 ∞ and expand the potential U as: Using this expansion with n = 4 we can get the various contribution to the integral I, by using the properties of the spherical Bessel function [32]: with . (B16)

Born approximation to the deflection angle
Using the Born approximation, we can estimate the deflection angle in the large m limit m. Using (B12), we deduce that the Born phase shift is given by: (B17) In the large m limit, we have that I ≈ I 1 and I 1. Using the Taylor expansion of arctan near 0, we have that: Finally using the relation between the deflection angle and the phase shift Θ = 2dδ m /dm, we can estimate the deflection angle in the Born approximation as:

Application to the scattering cross-section
The Born approximation gives us a way to estimate the phase-shift by evaluating the integral (B11). The full potential U can be used in the integral or one can substitute the leading order expansion of U for large x in order to analytically get the large m limit to the phase-shift. These approximations can be used in combination with the numerical method detailed in III A to accurately compute the scattering cross-section.
The scattering cross-section can be computed from the partial wave expansion by evaluating the infinite sum given in Eq.(37). This sum can be split into three component which will be evaluated using different method: f ω (θ) = a m cos(mθ) = S 1 + S 2 + S 3 .
(B20) S 1 , corresponding the low-m contribution, can be evaluated by using the phase-shift coefficients obtained by solving numerically the wave equation.
S 2 , corresponding the mid-m contribution, can be evaluated the phase shifts via the Born approximation by numerically integrating (B11) and using the full potential U . S 3 , corresponding the large-m contribution, can be estimated by using the large-m limit analytic expression for the phase-shifts in the Born approximation with Eq.(B16). Indeed S 3 is given as a combination of the Hurwitz-Lerch transcendent Φ(z, s, a) which satisfy: Φ(z, s, a) = ∞ k=0 z k (k + a) s .
By including high order terms in the large-m limit of the Born phase shifts, one can increase the accuracy of the estimate for S 3 . Figure 11 shows the contribution of the different sums in computing the scattering cross section. We can see that the biggest gain in accuracy comes from including S 2 which consists in approximating the phase shifts using the Born approximation with the full potential U . FIG. 11. Contribution to the scattering cross section for f = ω 2π = 8Hz of the various sums S 1 , S 2 and S 3 . Here the cut off between the sums are N 1 = 80 and N 2 = 280.
which can be rewritten as: where k 0 = ω 0 /c 0 , c 2 0 = gh 0 and b = B/h 0 . This equation is a hypergeometric equation and its solutions are given by hypergeometric function 2 F 1 : with s = b(k 2 0 + bm 2 ). From the asymptotic expression of the hypergeometric functions when r → 0, and the boundary condition at the origin, we deduce that A 2 = 0. Therefore the physical solution is: 2. Region II : r ≥ R In the region r ≥ R, the speed of the waves is a constant and the wave equation has the form The solutions to his equation are given in terms of Bessel functions: From the asymptotic expansion of the Bessel functions, we can expression the in and out coefficient in terms of B 1 and B 2 as This implies that the phase shift can be simply written as

Matching condition
The solutions in region I and region II should agree at r = R. This implies that: where F = 2 F 1 b|m|−s 2b , b|m|+s 2b , 1 + |m|, bR 2 , J m = J m (k ∞ R) and Y m = Y m (k ∞ R). Similarly, their derivatives must also be continuous. This implies that The two matching conditions gives us a relation between B 1 and B 2 : where the coefficient α is fully determined by the parameters of the system. Explicitly with Inserting this into the expression for the phase shift, we get that e 2iδm = α − i α + i . (C20)