Statistical mechanics of specular reflections from fluctuating membranes and interfaces
Amir Azadi David R. Nelson

























 








Department of Physics, Harvard University, Cambridge, MA 02138, USA
Wolfram Research, Somerville, MA 02144, USA



Department of Physics, Harvard University, Cambridge, MA 02138, USA
School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA



We study the density of specular reflection points in the geometrical optics limit when light scatters off fluctuating interfaces and membranes in thermodynamic equilibrium. We focus on the statistical mechanics of both capillary-gravity interfaces (characterized by a surface tension) and fluid membranes (controlled by a bending rigidity) in thermodynamic equilibrium in two dimensions. Building on work by Berry, Nye, Longuet-Higgins and others, we show that the statistics of specular points is fully characterized by three fundamental length scales, namely, a correlation length ξ, a microscopic length scale ℓ and the overall size L of the interface or membrane. By combining a scaling analysis with numerical simulations, we confirm the existence of a scaling law for the density of specular reflection points, n_spec, in two dimensions, given by n_spec∝ℓ^-1 in the limit of thin fluctuating interfaces with the interfacial thickness ℓ≪ξ_I. The density of specular reflections thus diverges for fluctuating interfaces in the limit of vanishing thickness and shows no dependance on the interfacial capillary-gravity correlation length ξ_I. Although fluid membranes under tension also exhibit a divergence in n_spec∝(ξ_Mℓ)^-1/2, the number of specular reflections in this case can grow by decreasing the membrane correlation length ξ_M.
 

Statistical mechanics of specular reflections from fluctuating membranes and interfaces
    David R. Nelson
    December 30, 2023
=======================================================================================





§ INTRODUCTION

The intricate dancing pattern of bright points and curves which can be seen as sea surface reflections or as caustics at the bottom of a swimming pool, in the presence of sunlight, echoes the geometry of rippling water surfaces. These specular points and caustic reflections were known historically to the ancient Greeks, and were sketched and appreciated by Leonardo da Vinci <cit.>.
Provided light wavelengths are short compared to ripple sizes, they obey the principles of classical geometric optics as embodied in Fermat's principle <cit.>. However, systematic studies of the statistics of the specular points and optical caustics which encode the geometrical features of fluctuating surfaces are relatively recent. Among the theoretical and experimental studies <cit.>, of particular interest is the pioneering work by Munk and Cox <cit.> who developed a method to measure the roughness of the sea surface agitated with wind driven waves, thus determining the effect of wind on the mean square surface slope from the spatial distribution of specular reflections. 
  
The distribution of specular reflections is also connected to problems in condensed matter and stochastic processes which include the statistics of zero crossings of random functions, and the noise currents embodied in the Shot noise <cit.>.
The statistics of maxima, minima and saddle points, for two dimensional membranes and interfaces is related to the density of specular reflection points of a stochastic scalar field, in one and two
dimensions. The generalization to a variety of physical systems is the subject of extensive works by Berry, Upstill  <cit.>, Longuet-Higgins <cit.> and others inspired in part by the catastrophe theory of caustic formation by light <cit.>. In Refs. <cit.>, Halperin et al., have studied related statistical properties of the zeros of the n-dimensional vector fields in the d-dimensional space. These points for n=d (or curves in n=d-1), characterize the topological singularities in the orientation of a vector field <cit.>. 
Here we apply related ideas to specular reflections of various types membranes and interfaces in thermodynamic equilibrium. By “interface" we mean a boundary between a liquid and gas phase, where the restoring forces are gravity and surface tension. We use “membrane" to denote the system such as lipid bilayers with aqueous phase above and below, where the dominant forces are a bending rigidity and an effective surface tension. 


An analysis of specular reflections off an undulating water surface illuminated by a distant light source such as the sun <cit.> presents fascinating problems linked to the rich statistical dynamics of capillary and gravity waves. With wind-driven wave excitations, generated by the non-equilibrium dynamics of the atmosphere, one must deal with complex nonlinear mode-couplings and energy transfers across length scales associated with driven wave turbulence <cit.>. In this paper, we investigate two simpler models of reflecting surfaces, associated with interfaces and membranes in thermodynamic equilibrium in two dimensions. Exact results are readily generated, which can be checked by straightforward computer simulations. In this sense, our investigation is analogous to the “absolute equilibrium" models of homogeneous, isotropic turbulence, which are interesting in their own right and can sometimes provide insights into more complex problems, such as the direction of turbulent energy cascades <cit.>.

Our first system, designed to model specular reflections off, say, air-water interfaces, assumes that the configurations f(x⃗), of a single-valued interface height profile in d-dimensions (see Fig. (<ref>(a)) are governed by an equilibrium probability distribution
P_int[f(x⃗)]∝exp[-F_int[f(x⃗)]/k_BT], with an equilibrium interfacial free energy given by 

    F_int[f(x⃗)]=1/2∫ d^d-1x [σ|∇⃗ f(x⃗)|^2+ρ_0 g f^2(x⃗)].
		
Here, the first term arises from the gradient expansion of the contribution from a surface tension σ, 
F_s=σ∫ d^d-1 x √(1+ |∇⃗ f(x⃗)|^2) and the second from the gravitational potential energy, ρ_0g∫_0^h(x⃗)f'df'=1/2ρ_0 g h^2(x⃗), integrated over the surface of an incompressible liquid with height h(x⃗) and mass density ρ_0, in equilibrium with a vapor phase of negligible density. Here g is the gravitational constant and f(x⃗) is the deviation of the fluid height from its equilibrium value h_0. The interfacial length scale ξ_I=√(σ/ρ_0 g) marks the boundary between capillary and gravity wave excitations <cit.>. A similar long-wavelength description describes solid-vapor interfaces above the roughening transition <cit.>.

Our second model describes specular reflections off a membrane, e.g., a thermally fluctuating lipid bilayer, suspended in water across a hole of fixed cross-sectional area with zero osmotic pressure difference (see Fig. (<ref>(b)). We use the Monge representation to describe a nearly flat membrane embeded in d dimensions located at r⃗(x_1,...,x_d-1)=[x_1,...,x_d-1, f(x⃗)]. In the absence of this constraint, fluctuating membranes are controlled by a bending free energy, F_b=1/2κ∫ d^d-1x |∇^2 f(x⃗)|^2, characterized by a bending rigidity κ. We assume an incompressible lipid bilayer, and implement the constraint of fixed total membrane area, A_0 via a Lagrange multiplier α, and thus describe the statistical mechanics by a total free energy, F_tot=F_b-α∫ d^d-1 x√(1+|∇⃗ f(x⃗)|^2). Upon expanding the square root in this integral and neglecting a constant term, the probability of a particular membrane configuration is given by 
P_mem[f(x⃗)]∝exp[-F_mem[f(x⃗)]/k_BT], where 

    F_mem[f(x⃗)]=1/2∫ d^d-1x [κ|∇^2 f(x⃗)|^2+α|∇⃗ f(x⃗)|^2].
	
The constant α is fixed by the condition that 

    ∫ d^dx⟨√(1+|∇⃗ f(x⃗)|^2)⟩=A_0>A,
	
where the brackets indicate a α-dependent thermal average over P_mem[f(x⃗)] in an ensemble specified by Eq. <ref>. Here A_0 is the area an incompressible membrane would have at zero temperature and A is the projected area, shown in Fig. 1b. 
In the small gradient approximation, assumed throughout this paper, α is thus determined by,
∫ d^d-1x⟨ |∇⃗ f(x⃗)|^2⟩=A_0-A
(If A<A_0, α is negative and the membrane can be subject to a buckling instability; in this paper, we assume α>0). Note that the effect of the constraint of fixed area leads to a tension-like contribution to the effective free energy. Eq. (<ref>) now defines a membrane length scale ξ_M=√(κ/α). The dominant restoring force for fluctuations with wavelengths λ less than ξ_M is the bending rigidity, while the tension term dominates for λ>ξ_M. Similar free energies arise for one-dimensional polymers in the wormlike-chain approximation with a pulling force applied to the ends <cit.>.





In this article, by combining analytics and numerical simulations, we study the interplay between the geometry and statistical mechanics of specular reflection with material properties that describe membranes and interfaces in thermal equilibrium.
We begin by calculating the number of zeros of a random Gaussian scalar field, f(x) (or more generally the number of times f(x)=y, some fixed height), controlled by probability distributions such as Eqs. (<ref>) and (<ref>) using path integral methods. A generalization of this method allows us to calculate the distribution of specular points, generated by a light source far from the plane of the membrane or interface, as seen by a distant observer, due to reflections from a random surface. 

In this paper, we focus for simplicity on d=2, i.e. one-dimensional membranes and interfaces. We hope to publish a paper on specular reflections off two-dimensional surfaces in d=3, including an experimental test using reflecting undulating surfaces generated by a three-dimensional printer, in the future.  

In the following sections, we focus first on the density of specular points for the special case where the light source and the observer both reside on a single line perpendicular to the average plane of the reflecting surface. We then adapt this calculation, using the paraxial approximation <cit.>, to the case of grazing angles, directions far from the average surface normal where the observer is nevertheless remains far from the reflecting interfaces or membranes.
These calculations allow us to explore how the density of specular points is related to the correlation lengths, ξ_I and ξ_M, introduced above to characterize the height-height correlation functions of membranes and interfaces. 

After treating membranes and interfaces in two dimensions, we explore certain geometrical features of the fluctuating surfaces and membranes such as average mean curvature, which is closely related to the density of specular reflection points. The interface and membrane correlation lengths ξ_I and ξ_M determine the density of specular reflection points, when combined with a microscopic length scale ℓ and an overall interface or membrane size L. Having tabulated the results of specular reflection theory, we then test the validity of the scaling laws for the density of specular points with Monte-Carlo simulations of interfaces and membranes with varying correlation lengths. 





§ EQUILIBRIUM THERMODYNAMICS OF FLUCTUATING INTERFACES AND MEMBRANES



 §.§ Free energy and height correlations in two dimensions

The configuration of a one-dimensional fluctuating interface or membrane (curve) in two dimensions can be described by the height function f(x), which we assume is single valued, that characterizes the deviation from the flat state. A long wavelength free energy that can be used to describe both interfaces and membranes is given by, 

    F_1d   =   1/2∫ dx [c(∇^2 f)^2+b(∇⃗ f)^2+a f^2].

If we set c=0, and take a to be proportional to the gravitational constant and b to the surface tension, we obtain the capillary-gravity interface model of Eq. (<ref>).
Note that the first term in eq. (<ref>) ∼ c(∇^2 f)^2, could then be introduced in the capillary-gravity model to suppresses short wave length fluctuations. Here we shall instead simply forbid wavevectors with q≳π/ℓ in the Fourier expansion of the free energy,

    F_I   =   1/2∫ dx [b(∇⃗ f)^2+a f^2],

where ℓ is a microscopic length of order the interface thickness.  

On the other hand, neutrally buoyant membranes subject to an isotropic tension due to an area constraints, as in Fig. 1, can be described by setting a=0,

    F_M   =   1/2∫ dx [c(∇^2 f)^2+b(∇⃗ f)^2].
  
The probability of a given equilibrium configuration f(x) at absolute temperature T is given by P(f)∝exp[-F/k_BT], where k_B is Boltzmann's constant. Henceforth, we rescale energy units such that k_BT=1. It will be convenient to expand f(x) in Fourier modes,
f(x)=1/L∑_qf_qe^iqx, where L is the macroscopic interface or membrane size, and we assume periodic boundary conditions. 

We first calculate the 
height correlation function which describes fluctuations of interfaces and membranes in thermodynamic equilibrium.  
Upon passing to Fourier space, the correlation function for the out-of-plane fluctuation of the two dimensional interface is described, in the limit of small c, by 

    C(y)   =   ⟨ f(x_0)f(x_0+y)⟩≈1/2√(ab)e^-|y|/ξ_I

where the capillary-gravity interfacial correlation length is ξ_I=√(b/a). This limit sends the effective ultraviolet cutoff in Fourier space 1/ℓ=√(b/c), to large values, and allows us to use the contour integration for q in the complex plane with δ=b^2-4ac>0 (see Appendix A for details).
Similar methods show that, in contrast, the second derivative of this correlation function depends explicitly on the ultraviolet cutoff parameter c. When c becomes small, we have

    C^(2)(y)
    ≈-1/21/√(bc)e^-|y|/ξ_I,

where C^(2)(y)≡d^2C(y)/dy^2. To determine quantities such as the densities of zero crossings and specular points, we shall need (see below) the ratio of the correlation function to its second derivative at y=0. This quantity is governed by two distinct length scales,

    C^(2)(0)/C(0)   ≈   -√(a/c)=-1/ℓξ_I

where the last inequality holds, provided we regard the free energy Eq. <ref> with a, b and c>0 as model of an interface with an ultraviolet cutoff imposed by the parameter c and 
where the effective short distance cutoff for this model is ℓ=√(c/b). As discussed in Appendix A, we work in a regime such that ℓ≪ξ_M≪ L, where L is the macroscopic size of the interface. The first inequality is equivalent to the condition δ=b^2-4ac≫0.



 §.§ Capillary-gravity wave and fluctuating membrane models
 We now focus on the two limiting cases a=0 and c=0 in more detail. The free energy given by Eq. (<ref>) with finite a and b and vanishing c is a model of capillary-gravity waves at the liquid-air interface, as discussed above. In this case, a is proportional to the liquid density ρ and the gravitational constant g, a∝ρ g, and b describes the line tension at the one dimensional interface. It is instructive to 

recompute the correlation functions for the capillary-gravity model in the Fourier space with c=0, while imposing a hard upper cutoff π/ℓ on the allowed wavevectors, 


    C(0)   =1/2π∫^+π/ℓ_-π/ℓdq/a+b q^2=ξ_I/barctan[ξ_Iπ/ℓ]≈1/√(ab)=ξ_I/b
    
    
    C^(2)(0)   =1/2π∫^+π/ℓ_-π/ℓ-q^2 dq/a+b q^2≈-1/bℓ
    
    
    C^(4)(0)   =1/2π∫^+π/ℓ_-π/ℓq^4 dq/a+b q^2≈π^2/3bℓ^3


where ℓ, the interface thickness, sets the ultraviolet cutoff, and the capillary-gravity wave interfacial correlation length is ξ_I=√(b/a). In this hydrodynamic treatment, we again consider the regime ℓ≪ξ_I≪ L. Note that this alternative ultraviolet cutoff leads to C^(2)(0)/C(0)=-1/ℓξ_I, a result in agreement with Eq. (<ref>), where the short distance was imposed by c. Here, C^(4)(0)=lim_y→ 0 C^(4)(y)=lim_y→ 0d^4C(y)/dy^4=lim_y→ 0⟨ f”(x+y)f”(x)⟩, we shall need it in later sections. 

We now discuss the case a=0, i.e. a membrane with a tension in two dimensions, or equivalently a semi-flexible polymer in two dimensions with a force applied to the ends <cit.>.
In this case, c represents the bending rigidity and b, a force or tension.
The height correlation function for a 1-dimensional membrane (or stretched semi-flexible polymer) with a tension in two dimensions at zero separation is given by,


    C(0)   =1/π∫^∞_π/Ldq/b q^2+c q^4≈Lξ_M/√(bc)
    
    C^(2)(0)   =d^2C(y)/dy^2|_y=0=-1/π∫^∞_π/Ldq/b +c q^2≈-1/2√(1/bc)
    
    C^(4)(0)   =d^4C(y)/dy^4|_y=0=1/π∫^π/ℓ_π/Lq^2 dq/b +c q^2≈1/ℓ c


where we impose an infrared cutoff by forbidding wave vectors such that |q|<q_min=π/L where L is the system size.
We have also defined the membrane correlation length to be, ξ_M=√(c/b). In our long wavelength expansion, we consider the limit of large system size compared to the membrane correlation length. We also assume a short distance cutoff ℓ, of order the membrane thickness, such that ℓ≪ξ_M≪ L. Note that the fluctuation amplitude,
⟨ f^2(x)⟩=C(0), diverges with the system's size L, but L drops out of Eq. (<ref>) in this limit.




§ ZERO CROSSINGS OF INTERFACES AND MEMBRANES IN TWO DIMENSIONS

Before determining the density of specular points for interfaces and membranes, we first adapt ideas of O. Rice <cit.> to determine the density of zero crossings, a simpler problem illustrating functional integral techniques used elsewhere in this paper.  
Consider a single-valued random function f(x), which could be the two-dimensional profile of an interface or membrane which can assume the height z=f(x) at various points along the x-axis. We assume that the probability of different configurations is Gaussian and described by a hydrodynamic free energy such as 
Eq. (<ref>).
We again Fourier analyze the single valued height function f(x), into plane waves:

    f(x)=∑_q f_qe^iqx.

Upon assuming periodic boundary conditions, f(x+L)=f(x), the allowed wave vectors in Eq. (<ref>) are, q_m=2π/Lm for m=0,±1,±2,... . In the limit L→∞, we again replace the discrete sum, ∑_q_m by an integral over q with dq=2π/L. 
The requirement that f(x) be real leads to a condition on the complex Fourier amplitudes, namely f^*_q=f_-q with a Gaussian probability distribution function controlled by the negative exponential of the underlying free energy divided by k_B T.
The statistical properties of f(x) are fully controlled by the autocorrelation function C(x)=⟨ f(y)f(x+y) ⟩, which can be expressed as C(x)=∫ dq E(q)exp(iqx),
where E(q)=⟨ |f(q)|^2⟩ is the “power spectrum" in q space.

Consider the joint probability distribution p(f=y,f_1;x) such that p(f=y,f_1;x) df gives the chance that f(x) in the small interval (x+dx) takes a value in the range f∈ [y,y+df] with an undetermined slope f_1=df/dx. We are interested in the probability that f(x) passes through zero, f=0 (or in general f(x)=y), for any value of the slope df/dx in [x,x+dx]. This quantity is given by ∫^∞_-∞ df_1   p(f=y,f_1;x) df=dx∫^∞_-∞ p(f=y,f_1;x)|f_1|  df_1.
The probability density, p(f=y,f_1;x) can be expressed as a functional integral in terms of a weighted average on all possible configurations of
f(x),
 
    p[f(x)=g,f'(x)=g_1;x]   =   ∫𝒟f(y)δ[f(x)-g]δ[f'(x)-g_1]e^-F[f(y)]/∫𝒟f(y) e^-F[f(y)]
    
       =   1/2π[-C(0)C^(2)(0)]^1/2exp[-f^2/2C(0)+f^2_1/2C^(2)(0)],

where we define,


    C(x)   =⟨ f(y)f(x+y) ⟩
    
    C^(2)(x)   =d^2 C(x)/dx^2=-⟨df(y)/dydf(x+y)/dy⟩


The derrivation of Eq. (<ref>) is reported in Appendix B.

We can now calculate the average number of times the function f(x) in the interval -L/2<x<+L/2 passes through a specific value f=y (subject to periodic boundary conditions), 

    N_f=y   =   ∫^+L/2 _-L/2dx∫^∞_-∞ p(f=y,f_1;x) |f_1|  d f_1
    
    
       =   ∫^+L/2 _-L/2dx[-C^(2)(0)/π^2 C(0)]^1/2exp[-f^2/2C(0)]
       =   L[-C^(2)(0)/π^2C(0)]^1/2exp[-f^2/2C(0)]

Thus, the number density n_f=y=N_f=y/L of crossings is given by, 

    n_f=y=[-C^(2)(0)/π^2 C(0)]^1/2exp[-y^2/2C(0)]



In order to find the density of zero crossings, n_0, we need to simply substitute y=0 in Eq. (<ref>). 


We can now evaluate the density of crossings for interfaces and membranes from the correlation functions tabulated in Section II.
We can express the number of times per unit length the function f(x) crosses the value y in terms of the correlation lengths ξ_I or ξ_M and the microscopic system thickness ℓ and the macroscopic system size L, 

    n_I   ≈   1/π√(ℓξ_I)exp[-y^2/2⟨ f^2(x)⟩]           capillary-gravity wave interface
    
    
    n_M   ≈   1/π√(Lξ_M)exp[-y^2/2⟨ f^2(x)⟩]            membrane under tension

where ⟨ f^2(x)⟩=C(0) is given by Eqs. (<ref>) or (<ref>). Note that in the capillary-gravity interface model the geometric mean of ξ_I and the microscopic cutoff ℓ controls the density of zero crossings, whereas in the membrane model the geometric mean of ξ_M and the macroscopic system size L is the controlling factor. Thus there are many fewer zero crossings in the softer membrane where longer wavelength fluctuations dominate.  In both cases the number of crossings falls off rapidly where y≳√(⟨ f^2(x)⟩). Reasonable agreement with these predictions is found using Monte-Carlo simulations methods (see Sec. IV-C) in Figs 2 and 3 for interfaces and membranes, respectively. 
 




§ STATISTICAL GEOMETRY OF SPECULAR POINTS ON INTERFACES AND MEMBRANES IN THE PARAXIAL APPROXIMATION

In this section we use the approach sketched in the last section to determine the distribution of specular reflection points from thermally excited interfaces and membranes in two dimensions. We use simple ideas from the theory of focal singularities <cit.> to determine
the geometric condition for the formation of specular points.

In the following two subsections, we determine the density of specular points for two dimensional interfaces and membranes in two limits: large and small angles of incidence. 


 §.§ Nearly normal angles of incidence and reflection

We assume that the thermally fluctuating boundary, described by Eq. (<ref>), is a reflecting surface. A source of light emits wavelengths short compared to the undulations of the interface or membrane, hence we can use geometrical optics. The source resides above the surface, in the x-z plane, at S=(x_1,z_1) and we have the observer at O=(x_2,z_2). Provided  multiple reflections can be neglected, the condition for an arbitrary point, P=(x,f(x)) on the reflecting surface f(x), to be a specular point follows from a straightforward application of Fermat's principle <cit.>.
The ray from the source S needs to reach to the point P on the surface (see Fig.5.d) and then travel to the observer at O. The distances, l_1=OP and l_2=SP are given by,

    l_i=√((x-x_i)^2+(f(x)-z_i)^2)≈ z_i-f(x)+(x_i-x)^2/2z_i

which holds for i=1,2 and where we have taken the limit of positive z_i>>f, known as paraxial approximation <cit.>, in which both the observer and source reside far above a fluctuating surface with small deviations from flatness for an interface or membrane such that ⟨ f(x)⟩=0. We assume first that both source and observer have x_i≈0 as in Fig. 4a, i.e., they lie approximately above the origin of our coordinate system, although possibly with different heights z_1≠ z_2. In the next section we consider a more general case with both x_i 0. Upon assuming a constant propagation velocity in the region above the surface, the path of least time is given by <cit.>,

    d/dx[l_1(x)+l_2(x)]=0,

which leads via Eq. (<ref>) to,

    1/x d f(x)/dx   =   1/2[1/z_1+1/z_2]≡ K_1.

To determine the density of points satisfying this condition for fixed z_1 and z_2, it is useful to introduce two auxiliary functions, 


    g_1(x)   =1/xf'(x)=1/x d f(x)/dx,
    
    
     g_2(x)   =f”(x)=d^2 f(x)/dx^2,


and to define K_1≡1/2[1/z_1+1/z_2]. Similar to the last section we seek from Eq. (<ref>) the probability of a specular point with g_1=K_1, somewhere in the interval [x,x+dx],

    p_spec(x,K_1)   =   ∫^∞_-∞ dg_2 P(g_1=K_1,g_2;x)|g_2-K_1|dx/x,

where P(g_1,g_2)dg_1dg_2 gives the probability that the functions in Eq. (<ref>) and (<ref>) assume the specific values g_1 and g_2 in the intervals [g_1,g_1+dg_1] and  [g_2,g_2+dg_2]. 
To evaluate Eq. (<ref>), we used the relation |dg_1|=|dg_1/dx|dx=|g_1-g_2|/|x|dx. 
For the quadratic energy functionals used here, the probability distribution of g_1=1/xdf/dx and the curvature, g_2 can be determined by the methods of appendix B to be
 
    P(g_1,g_2;x)=(2π)^-1|M|^1/2exp[-1/2∑_i,j=1,2M_ijg_ig_j],

where the 2×2 matrix M_ij is given by,

    M_ij=
    [ -x^2/C^(2)(0)             0;                            ;             0    1/C^(4)(0) ]




where C^(2)(0)=-⟨(df(x)/dx)^2⟩ and C^(4)(0)=⟨(d^2f(x)/dx^2)^2⟩ are given in terms of the power spectrum E(q)=⟨ | f(q) |^2⟩ of f(x) by,


    C^(2)(0)   =-∫^∞_-∞ q^2 E(q)dq < 0, 
    
    C^(4)(0)   =∫^∞_-∞ q^4 E(q)dq.


The probability distribution entering Eq. (<ref>), is thus given by,

    P(g_1,g_2;x)=|x|/2π[-C^(2)(0)C^(4)(0)]^1/2exp[g_1^2x^2/2C^(2)(0)-g^2_2/2C^(4)(0)].

Upon inserting Eq. (<ref>) into Eq. (<ref>), we find the one-dimensional density of spectral points n_s≡ P_spec(x,K_1)
for near normal angles of incidence and reflection as,

    n_s   =   ∫^∞ _-∞ dg_2[|K_1-g_2|/|x|]  |x|/2π[-C^(2)(0)C^(4)(0)]^1/2exp[+K^2_1x^2/2C^(2)(0)-g^2_2/2C^(4)(0)]
    
    
    
    
    
    
       =   exp[+K_1^2x^2/2C^(2)(0)] K_1/[-2π C^(2)(0)]^1/2[(2/π)^1/2β^-1exp[-β^2/2]
    +erf(β/√(2))]
       =    N_sexp[K_1^2x^2/2C^(2)(0)]K_1/[-2π C^(2)(0)]^1/2

where N_s=∫^∞ _-∞n_sdx and we define β=K_1[C^(4)(0)]^-1/2.
Using Monte-Carlo simulation of the fluctuating interfaces and membranes, we test the validity of Eq. (<ref>) which predicts a Gaussian fall off of the density of specular points as a function of distance from the common location (x=0) of the source and the observer. Figure 4.b depicts the results of the numerical simulations (see Sec. IV) represented by points and the prediction of the continuum theory in Eq. (<ref>) which is given by the dashed curve. The small deviations from a Gaussian at small x and large x are likely due to higher order gradient couplings, neglected in our hydrodynamic approach.   
 

 










 §.§ Distribution of specular points in the large angle approximation:
 
We now consider a source and observer at more general coordinates (x_i,z_i), still making the paraxial approximation (z_j>> |x_j| and z_i>>f(x)) in Eq. (<ref>) but considering now the limit |x_i|>>x, i.e., source and observer reside far to the left and right of a specular point at x. Thus, we consider the configurations, in which the observation and incident angles measured form the normal and given by tanθ_i=|x_i|/z_i<<1  are finite and large compare to the angles tan^-1(|x|/z_j) that give rise to the highest density of specular points.
With these assumptions, vanishing of the first derivative of the total length l=l_1+l_2 defined in Eq. (<ref>) with respect to x 
leads to the condition

    d f(x)/dx   =   -1/2[x_1/z_1+x_2/z_2]≡ k_2,

where the source point and the observer reside at S=(x_1,z_1) and O=(x_2,z_2) respectively. Similar to the scheme we developed in the last section, we now define the auxiliary functions,

    g_1(x)   ≡   f'(x)= d f(x)/dx
    
     g_2(x)   ≡   f”(x)=d^2 f(x)/dx^2

where g_1(x)=k_2 is the condition for a specular reflection. To calculate the density of specular points, we need the joint probability distribution function P[g_2,g_1;x] which can be calculated explicitly via functional integral methods (see Appendix B),

    P(g_1,g_2;x)   =   ∫𝒟f(y)δ[f'(x)-g_1]δ[f”(x)-g_2]e^-F[f(y)]/∫𝒟f(y) e^-F[f(y)]
    
       =   1/2π[-C^(2)(0)C^(4)(0)]^1/2exp[g_1^2/2C^(2)(0)-g^2_2/2C^(4)(0)],

which is independent of x (provided |x|<<|x_j|) and where we have again assumed a quadratic free energy for interfaces and membranes such as Eq. (<ref>) and (<ref>). The probability distribution for specular reflections is now,

    p_spec(x,k_2)dx   =   ∫^∞ _-∞dg_2 P(g_1,g_2;x)|dg_1/dx| dx
       =   ∫^∞ _-∞dg_2 P(g_1,g_2;x)|g_2|dx,

which leads to an x-independent density of specular points,

    n_s=1/π(C^(4)(0)/- C^(2)(0))^1/2exp[-k^2_2/2C^(2)(0)],

where 

    C^(4)(0)=lim_y→0⟨ f”(x)f”(x+y)⟩.

Note that the density of specular points is proportional to [C^(4)(0)]^1/2=[⟨(∂^2_x f)^2⟩]^1/2, i.e. to the root mean square fluctuations in curvature of the reflecting curve. 
We can now use the fluctuating interface and membrane models introduced in Sec. II to study the statistical geometry of the specular points on a membrane in thermodynamic equilibrium. We express this density in terms of interface or membrane correlation length ξ_I, macroscopic size L and microscopic dimension  ℓ of order the thickness. We focus on the scaling of density of the specular points at infinity, so that k_2=[x_1/z_1+x_2/z_2]/2=0. In this extreme paraxial limit, the density n_s, is characterized only by derivatives of the correlation functions C(x)=⟨ f(x)f(x+y)⟩, as in Eq. (<ref>). Upon using the relations in Eqs. (<ref>) -(<ref>) for the capillary-gravity interface, in the limit ξ≫ℓ and c=0 we obtain,   

    lim _ξ→ 0 n_s=1/π(C^(4)(0)/- C^(2)(0))^1/2≈1/√(3)ℓ

Using Eqs. (<ref>)-(<ref>) for fluid membranes (or the equivalent model of fluctuating stretched two-dimensional polymer) yields,

    lim _ξ→ 0 n_s=1/π(C^(4)(0)/- C^(2)(0))^1/2≈1/√(ℓξ_M)

Eqs. (<ref>) and (<ref>) demonstrate that the density of specular points diverges for both models as the microscopic cutoff, ℓ→ 0. Although the specular density on a capillary-gravity wave interface is dominated by the short distance physics and is independent of the correlation length, the density of specular reflections for membranes, n^M_s is controlled by the reciprocal of the geometric mean of the membrane correlation length and microscopic cutoff. Fluid membranes generate large number of specular reflections only in the limit of short correlation lengths ξ_M > ℓ.


The densities of maxima, n_M and minima, n_m for a fluctuating random curve are closely related to the density of specular points in the extreme paraxial limit when both the observer and source reside at infinity. The specular point condition at infinity reduces to f'(x)=k_2=0 which of course is equivalent to being at a maximum or minimum. Since the average density of maxima and minima must be equal, they are related to the density of specular points calculated above by n_m=n_M=n_s/2. 
We can also investigate the ratio of n_0, the density of zero crossings to the density of maxima n_M or minima n_m for interfaces and membranes. 
With the aid of Eqs. (<ref>)-(<ref>) and (<ref>)-(<ref>) we find n_0/n_M=n_0/n_m≈√(ℓ/ξ_M)≪ 1 and n_0/n_m=n_0/n_M≈√(ℓ/L)≪ 1 for the capillary-gravity interface and fluid membrane model respectively. Thus the density of zero crossings is smaller than the density of maxima (or minima), which is plausible since for a highly fluctuating curve there can exist large number of maxima and minima between two consecutive zero crossings.














 §.§ A Monte Carlo simulations of fluctuating one dimensional chain model


To explore the validity of expressions for the density of zero crossings and specular points given for interfaces and membranes in terms of C(x)=⟨ f(y+x)f(x)⟩ and its derivatives at the origin (C(0), C^(2)(0) and C^(4)(0)), we carry out Monte-Carlo simulations. We simulate a fluctuating 1-d surface in two dimensions described by the discretized free energy,

    F_dis=∑_j=1^N[ a'/4(f_j^2+f_j+1^2)+b'/2(f_j-f_j+1)^2+c'/2(2f_j-f_j-1-f_j+1)^2 ],

where in the continuum limit our simulation parameters a',b' and c' are related to the parameters appearing in Eq. (<ref>) by a'=aℓ, b'=b/ℓ and c'=c/ℓ^3. 
The height function f(x) is defined in the interval (0,L) where L is the system size with periodic boundary conditions. It is discretized by considering the height at two neighboring points, f_j and f_j+1, where f_j=f(x=x_j) is defined in the small interval (x_j,x_j+1)=(x_j+ℓ), j=1,...,N with N=L/ℓ. The continuum energy expression in Eq. (<ref>) can be recovered via taking the limit ℓ/L→ 0 in eq. (<ref>).

Monte Carlo simulations are performed on this “chain model" at temperature k_BT (in units of bℓ) with the following move: a random vertex i residing at (x_i,y_i=f(x_i)) is chosen and displaced to (x_i+δ x,y_i+δ y) where δ x and δ y are small independent random real numbers. Moves are accepted or rejected based on the Metropolis criterion, using the Boltzmann factor exp[-F_dis/k_B T] associated with the free energy in Eq. (<ref>). By tracking the time-dependance of correlation functions, we found that system with N=150 equilibrates after ≈ 10^6 MC steps. One MC step is defined as N vertex moves.
Thermodynamic quantities like the average number of zero crossings or specular points, are measured after equilibration has been reached.   
We performed the simulations by starting from a configuration where all N points reside on a straight line, equally spaced with separation ℓ.
The results of the Monte Carlo simulations for the average density of zero crossings, n_0, of the
fluctuating membrane illustrated in Fig. 2 and Fig. 3. For the capillary-gravity interface model we take  c'=0 in Eq. (<ref>), while for the membrane or polymer-like model, we take a'=0. The zero crossings are shown as a function of the ratio of the lattice spacing to the correlation length √(ℓ/ξ) in Figs. 2 and 3 with the corresponding error bars. The dashed line represents the analytical result for the density of zero crossings in eq. (<ref>) with the correlation function calculated based on the free energy with c=0, given by Eqs. (<ref>)-(<ref>). In fig.2b we confirm the Gaussian fall-off of the expected number of f(x)=y passages, n(y), as a function of the scaled specific height y/√(C(0)), for √(ℓ/ξ)=0,0.5,0.9. A similar comparison has been done for the density of zero crossings for the membrane model in Fig 3b. 

Fig.4b shows the dependence of the average density of specular points on x in the small angle approximation using the Monte-Carlo simulations (points) and the predictions of the continuum theory (dashed line) given in Eq. (28). It confirms the Gaussian fall off of the average density of specular points in agreement with Eq. (28) for the thermally excited interfaces.  
Next we simulate the density of specular reflection points from a thermally excited interface dressed with the capillary-gravity waves (c=0) as a function of the correlation length, ξ and the distance of the observer from the 1d membrane in terms of k_2=-1/2[x_1/z_1+x_2/z_2] in the extreme paraxial limit (|x_j|≪ z_j). Figure 5a and 6a illustrate the density of specular points at infinity when k_2=0 as function of √(ℓ/ξ), where ξ is the relevant interface or membrane correlation length. We show the dependence of the density if specular points on the distance of the observer from the reflecting line which is encoded in k_2 and described by Eq. (<ref>). The density of specular points was evaluated in the simulation by measuring the slope of the fluctuating curves and computing the expected number of times the local slope satisfies the condition for specular reflection df/dx=k_2 given by Eq. (<ref>). 


The density of specular points reaches a constant value for large ξ. 

The schematic presented in Fig.5c shows a simulated reflecting interface and the reflected light rays of a light source at infinity. The incoming light rays parallel to z axis and has not been shown. The number of reflecting light rays reaching the point of observation determines the number of specular points   
Similarly Fig 6a-c, depicts the result of the simulations for the fluctuating curve described via a=0 which mimics the behavior of the fluctuating membrane or polymer under tensile stresses characterized by the coefficient b and the bending modulus c. Fig.6a shows the expected number of specular points as a function of √(ℓ/ξ_M) for k_2=0 (observer at infinity) governed by eq. (<ref>) and Fig.6b illustrates the k_2 dependence for different values of √(ℓ/ξ_M). We also show the dependence of the density of specular points on the distance from the reflection center (x=0) on the reflecting surface has been shown in Fig4.b for both membrane and interface models described by c=0 and a=0. The result of the simulations is compared to the analytical calculation presented in eq. (<ref>) in the small angle approximation.    









§ CONCLUSION

To explore the statistical geometry of specular reflections from fluctuating surfaces in thermal equilibrium, we developed a mapping which connects configurations of thermalized surfaces described by the relevant elastic constants to the statistics of zero crossings and specular reflection points. The density of the specular points is proportional to the root mean square curvature of the interface or membrane.  

This study can be generalized to the statistics of fold, cusp and higher order caustics to construct a more complex mapping between the geometry, mechanics and thermodynamic of a surface and the statistical geometry of reflected caustics patterns in higher dimensions. More generally, this study suggests new avenues for surface characterization based on the distribution of specular patterns forms in a plane with dimensions lower that the dimension of the embedding space. This framework can be used to infer material properties such as surface tension and correlation lengths of thermalized interfaces and membranes, by studying the lower dimensional pattern of reflections.    

Upon zooming in the blurred reflected pattern of coherent light even at smaller length scales, one expects to observe the fine structure of the diffraction pattern associated with the optical wave field, consisting of a complex pattern of dislocation lines, which are disruptions in the reflected wavefront <cit.>. Future studies would be needed to reveal how the geometry and statistical mechanics of the fluctuating interfaces control the fine structure of the diffraction pattern. 



§ ACKNOWLEDGMENTS

It is a pleasure to acknowledge stimulating interactions with Massimo Cencini, Daniel A. Beller and Yoav Lahini. We are grateful to Michael Berry for comments on the manuscript. This work was supported by the National Science Foundation, through grants DMR-1608501 and via the Harvard Materials Science Research and Engineering Center via grant  DMR-1435999. We would like to dedicate this article to the memory of Pierre C. Hohenberg, whose pioneering work on dynamic critical phenomena with B.I. Halperin and others had a profound influence early in the career of us.





10
vinci
Drawing of caustics and light reflections in curved mirrors by Leonardo da Vinci, British library, codex Arundel 263, folio 87 v. (1508).
nye
J. F. Nye. Natural focusing and fine structure of light, Institute of Physics Publishing; 1st edition (1999).
berry
M. V. Berry and C. Upstill, IV Catastrophe Optics: Morphologies of Caustics and Their Diffraction Patterns,
Chapter in Progress in Optics 18:257-346, (1980)
berry_2
M. V. Berry, Journal of the Optical Society of America A Vol. 4, Issue 3, pp. 561-569 (1987)

long_1
M. S. Longuet-Higgins, The statistical analysis of a random, moving surface Phil. Trans. R. Soc. A 249 321 87 (1957).
long_2
M. S. Longuet-Higgins, Statistical properties of an isotropic random surface Phil. Trans. R. Soc. A 250 157 74 (1957).
long_3
M. S. Longuet-Higgins, The statistical distribution of the curvature of a random Gaussian surface Proc. Camb.
Phil. Soc. 54 439 53 (1958)
ref1
I. M. Fuks and M.I. Charnotskii, Statistics of specular points at a randomly rough surface. J. Opt. Soc. Am. A 23, 73-80 (2006).
ref2
A.A. Maradudin and T. Michel, The transverse correlation length for randomly rough surfaces, J. Stat. Phys. 58, 485-501 (1990).
ref3
P. Swerling, Statistical properties of the contours of random surfaces, IRE Trans. On Information Theory 8, 315-321 (1962).
ref4
S. Simeonov, A.R. McGurn, and A. A. Maradudin, Proc. SPIE 3141, 152-
163 (1997).

Munk
C. Cox, and W. Munk, Measurement of the roughness of the sea
surface from photographs of the sun's glitter, Journal of the Optical Society
of America 44, 838 (1954).

rice
S. O. Rice, Mathematical analysis of random noise Selected Papers on Noise and Stochastic Processes edited by N. Wax, New York: Dover (1954).

halp
B. I. Halperin, Statistical mechanics of topological defects Les Houches Session XXV Physics of Defects
ed R Balian, M Kleman and J-P Poirier, Amsterdam: North-Holland, (1981).
speckle
J. W. Goodman, Laser Speckle and Related Phenomena, edited by J. C. Dainty, Springer, Berlin, (1975).
halp_2
A. Weinrib , and B. I. Halperin, Distribution of maxima, minima, and saddle points of the intensity of laser speckle patterns Phys. Rev. B, 26, 1362 (1982). 


lax
B. I. Halperin and Melvin Lax, Phys. Rev. 148, 722 (1966); 153, 802 (1967).
a
 V. E. Zakharov, V. S. L'vov and G. Falkovich,   Kolmogorov Spectra of Turbulence I: Wave Turbulence, Springer Science and Business Media (2012).
b
S. I. Badulin, A. N. Pushkarev, D. Resio and V. E. Zakharov,  Self-similarity of wind-driven seas, Nonlinear Processes in Geophysics, 12, 891 (2005).
c
V. E. Zakharov , P. Guyenne, A. N. Pushkarev, and F. Dias, Wave turbulence in one-dimensional models, Physica D: Nonlinear Phenomena 152, 573 (2001).
d
S. A. Orszag, Statistical theory of turbulence, in Fluid Mechanics, Les Houches, eds. R. Balian and J. L. Peube, Gordon and Breach, New York (1973).
e
U. Frisch, U., Turbulence: The Legacy of AN Kolmogorov, Cambridge University Press (1995).
f
L. D. Landau, L.D. and E. M. Lifshitz, Fluid mechanics, Course of Theoretical Physics, Elsevier, Boston (1987).
g
D. R. Nelson, D.R., Defects and geometry in condensed matter physics, Sec. 5.1, Cambridge University Press, Cambridge (2002).
h
F. C. MacKintosh J. Kas and P. A. Janmey, Elasticity of Semiflexible Biopolymer Networks, Phys. Rev. Lett. 75, 4425 (1995).  

min
G. E. Volovik, V. P. Mineev, Investigation of singularities in superfluid He3 in liquid crystals by the homotopic topology methods, Zh. Eksp. Theor. Fiz. 72, 2256 (1976).
fey
A. Hibbs and R. Feynman, Quantum Mechanics and Path Integrals, McGraw-Hill (1965).

Nelson_Seung
H. S. Seung and D. R. Nelson, Phys. Rev. A 38, 1005 (1988).

wave_1
I. Simonsen, A. A. Maradudin and T. A. Leskova, Phys. Rev. A 81, 13806 (2010)
wave_2
N. B. Baranova, N. B. Mamaev, A. V. Pilipetskii, et. al. J. Opt. Soc. Am. A 73 525 (1983)
wave_3
I. Freund, Phys. Rev. E 52, 2348 (1995)
wave_4
M. S. Soskin, V. N. Gorshkov, M.V. Vasnetsov, et. al. Phys. Rev. A 56, 4064 (1997) 






§ APPENDICES





§ CORRELATION FUNCTION CALCULATION



Here we provide details of the correlation function calculations reported in section II. A for the fluctuating interfaces or membranes in 1+1 dimensions. From Eq. (<ref>) we obtain, 

    C(y)=⟨ f(x_0)f(x_0+y)⟩=1/2π∫^+∞_-∞dq e^iqy/a+b q^2+ c q^4

Upon considering the complex q-plane, the denominator defines four poles at:

    q   =   ± i [b∓√(δ)/2c]^1/2
    δ   =   b^2-4ac
  
We assume that δ>0. Upon defining an interfacial capillary-gravity correlation length by ξ_I≡√(b/a) and the ultraviolet cutoff induced by the c-term as ℓ≡√(c/b), this condition simply means that ξ_I>2ℓ. Because all poles reside on the imaginary axis, we also have b>√(δ). By completing the contour in the upper half plane, which encloses the poles at q= i [b±√(δ)/2c]^1/2, we find,

    C(y)=1/2π∫^+∞_-∞dq e^iqy/a+b q^2+ c q^4   =    i (e^iqy/2bq+4cq^3)_q=i [b+√(δ)/2c]^1/2+ i (e^iqy/2bq+4cq^3)_q=i [b-√(δ)/2c]^1/2
       =   1/√(2 δ)[e^-√((b-√(δ)/2c))|y|/√(b-√(δ)/c)-e^-√((b+√(δ)/2c))|y|/√(b+√(δ)/c)]

In a similar fashion we find that the second derivative of the correlation function is given by,

    C^(2)(y)   =   -1/2π∫^+∞_-∞q^2 dq e^iqy/a+b q^2+ c q^4
       =   1/2√(2 δ)[√(b-√(δ)/c)e^-√((b-√(δ)/2c))|y|-√(b+√(δ)/c)e^-√((b+√(δ)/2c))|y|]

Note that C^(2)(0) diverges as c→0, reflecting the sensitivity of this quantity to the effective ultraviolet cutoff ℓ=√(c/b) as ℓ→0.



§ FUNCTIONAL INTEGRALS FOR THE PROBABILITY DENSITY OF A FLUCTUATING MEMBRANE

Here we calculate, via functional integral methods, the normalized probability of a function f(y) to have a specific value and slope at position x given by f(x)=g, f'(x)=g_1. The probability density for these quantities with a fluctuating interface described by our simple quadratic free energy functional F_I[f(y)], (we use the interface free energy Eq. (<ref>) for concreteness, but identical manipulation apply for the membrane free energy Eq. (<ref>)) is given by, 
 
    p[f(x)=g,f'(x)=g_1;x]   =   ∫𝒟f(y)δ[f(x)-g]δ[f'(x)-g_1]e^-F_I[f(y)]/∫𝒟f(y) e^-F_I[f(y)]
       =   ∫𝒟f(y)δ[f(x)-g]δ[f'(x)-g_1]e^-1/2∫ dy[a f(y)^2+b(df(y)/dy)^2]/∫𝒟f(y) e^-F_I[f(y)]
       =   ∫ ds∫ dt e^-2π i sge^-2π i t g_1⟨ e^2π i s f(x)+2π i t f'(x) ⟩
       =   ∫ ds∫ dt e^-2π i sge^-2π i t g_1e^-2π^2 ⟨ f^2(x)⟩ s^2 -2π^2 ⟨ f'^2(x)⟩ t^2-4π^2 ⟨ f(x) f'(x)⟩ st  
       =   ∫ ds∫ dt e^-2π i sge^-2π i t g_1e^-2π^2 C(0) s^2+2π^2 C^(2)(0) t^2

where ⟨ .⟩ represents a thermal average, and
⟨ f(x)f'(x)⟩=1/2⟨d/dxf^2(x)⟩=0 with our periodic boundary conditions. We also have C(0)=⟨ f^2(x)⟩ and C^(2)(x)=-⟨(df/dx)^2⟩.
Here we used the identity, valid for any Gaussian probability distribution, ⟨ e^h(x)⟩=e^1/2⟨h^2(x)⟩ and the integral representation of the delta function,

    δ[f(x)-g]=∫^∞_-∞ ds e^2π i [f(x)-g]s .

This probability density for a fluctuating membrane is normalized by construction, ∫^∞_-∞ p[g,g_1]dg_1,dg_2=1.
By the aid of standard Gaussian integrals we find from Eq. (<ref>) the probability density, p[f(x)=g, f'(x)=g_1;x] is in fact independent of x,

    p[f(x)=g,f'(x)=g_1;x]=1/2π[-C(0)C^(2)(0)]^1/2exp[-g^2/2C(0)+g^2_1/2C^(2)(0)].






Similarly, one can calculate the normalized probability of a function f(y) to have a given slope and curvature described by f'=g_1, f”=g_2 at the position x,

    p[f'(x)=g_1,f(x)”=g_2;x]   =   ∫𝒟f(y)δ[f'(x)-g_1]δ[f”(x)-g_2]e^-F[f(y)]/∫𝒟f(y) e^-F[f(y)]
       =   ∫ ds∫ dt e^-2π i sg_1e^-2π i t g_2⟨ e^2π i s f'(x)+2π i t f”(x) ⟩
       =   ∫ ds∫ dt e^-2π i sg_1e^-2π i t g_2e^2π^2 C^(2)(0) s^2-2π^2 C^(4)(0) t^2  
       =   1/2π[-C^(2)(0)C^(4)(0)]^1/2exp[g_1^2/2C^(2)(0)-g_2^2/2C^(4)(0)]

where C^(4)(0)=lim_y→0⟨ f”(x+y)f”(x)⟩=⟨(d^2 f/dx^2)^2⟩. In a similar fashion, we can find the probability density of a stochastic function with fixed higher order derivatives, at position x, with d^n f(x)/dx^n=g_n.

