Quantitative analysis of observed signatures in phonon reflection experiment

For phonons propagating in a crystal at low temperature, we develop a phonon reflection model and an analysis framework to fully understand reflection peaks from underlying physics principles. Our study shows that the specular reflection structure of peaks contains detailed information about properties of the crystal material, and provides a solid quantitative link between phonon theory and experimental signals. For both kinds of experiments that either use crystal as a particle detector or as a target material for study, this analysis provides a novel and essential method to understand the experimental signatures.


Introduction
In the particle detection experiments of low noise environment, a large number of experiments use cryogenic bolometers at low temperature [1]. The detection principle is based on collection of phonons, the acoustic wave quanta. Especially in rare event detection of dark matter or double beta decay, bolometers using crystal target have been widely used due to the excellent resolution and sensitivity [2]. When a nuclear or electron recoil occurs due to energy deposit, phonons and electron-hole pairs are created. The electron-hole pairs will recombine to create additional phonons or scintillation lights. The net effect for the energy which goes to heat is the creation of high frequency phonons at the interaction point. One or more phonon collectors which are usually attached to the crystal surface, absorb hitting phonons. The detected signals depend on heat transfer to the collector through phonon travel.
At present, a full understanding of the experimental signal shape has not been achieved, and empirical data are used as templates in signal processing [3,4]. The detailed phonon scattering and down-conversion at crystal medium has been Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. studied and those models agree to experimental data to some extent [5]. However, for practical use of cryogenic crystal bolometers, most of the athermal phonons reflects many times at crystal surfaces before propagating to the phonon collector [3]. The crystalline material is always anisotropic due to the cell structure. This leads to three non-trivial vibration modes of a phonon, for each of which, there are three different reflection transitions. So a full 3 × 3 matrix is required to describe the reflection properties at the interface, and the matrix itself is a complicated function of both incident and interface directions. Thus, it is necessary to develop a systematic method to obtain reflection and transmission coefficients using boundary conditions. However, the complexity of experimental setups usually prevents direct testing and validation of such methods. In fact, only a simple mirror reflection model is used in reference [5].
On the other hand, in the field of solid state physics, experiments have attempted to study the propagation [6,7] and reflection [8,9] of thermal phonons through a pure crystal at low temperature. Quantitative investigations has been done in phonon transmission [10] and reflection [11] experiments using analytic expressions. However, due to strong phonon-phonon interaction, poor time resolution and high scattering background, quantitative extraction of reflection peak sizes from experiment data was not possible in reference [11]. Recently, practical computation methods have been tried to determine phonon reflection probabilities in reference [12], but the results are questionable because the evanescent waves and corresponding wave equations were just ignored, possibly causing the large jumps in the reflection probability plots.
Here, we focus on experiments with high resolution of time of flight, that have shown distinct peaks [13,14] at expected arriving times from a single reflection. They provide a useful tool to study the specular reflection phenomena and the anisotropic propagation of phonons. The simple and dedicated experimental geometry leaves direct comparison between results from theory and experiment possible, because the peaks correspond to short travel path and the diffusion effect is small. Only preliminary qualitative analysis has been done for the reflection peak sizes in the original experiment [13].
A full simulation provides a natural method to obtain the experimental spectra by implementing the phonon propagation principles. Reference [15] applies the full simulation to the geometry of the reflection experiment [13], with the scattering process included. The result shows that the simulation can well explain the experimental data.
However, it remains to explain the structure of reflection peaks by directly using the reflection coefficients and related properties. In this paper, we will discuss the development of a phonon reflection model that can be used in a framework to directly compute the experimental energy deposit spectra from crystal properties.

Phonon reflection model
In anisotropic material of a pure crystal, the acoustic wave travels along its group wave velocity w, which is not the same as the wavevector q. The three-dimensional wave equation is [16]: where ρ is the mass density of the crystal, ω is the wave frequency, and C i jlm is the elastic constant tensor. The normalized matrix (1/ρ)C i jlm q n, j q n,m (q n = q/q) gives three eigenvalues as phase velocity c = ω/q squared, with three eigenvectors as polarization vectors (e). The slowness vector is defined as the wave vector divided by the wave frequency s = q/ω, whose magnitude is the reversal of the phase velocity c. The three states are usually labeled as longitudinal (L), fast transverse (FT) and slow transverse (ST) modes, ranked as ascending order of slowness magnitude. The group velocity w can be found by differentiating ω with respect to the wavevector q, to give: The unit energy flow density vector S is closely related to the group velocity w: S = ρωw. The energy flow per unit area of surface with unit amplitude is then S · n, the projection of S to the surface unit normal vector n.
Aside from a common space-time dependent term, the second-rank strain tensor for spacial derivatives of displacement is lm = 1 2 (e l q m + e m q l ) with unit amplitude. And the stress tensor σ is related to the strain tensor as: σ i j = C i jlm lm . The total traction force to a surface is the projection of the stress tensor on the normal vector n, as σ i j n j in component form. During a specular reflection at a vacuum interface, there are three additional reflected modes labeled as (τ = 1, 2, 3), and the total traction force across the interface should be zero: Here the Γ τ are relative reflection amplitudes for three reflected modes, to be solved with three equations for i = 1, 2, 3. Since number of unknown variables and independent equations are same, we obtain three Γ τ values, which could be complex if there are complex wavevectors involved. To construct equation (3), we need the slownesses s τ of the three reflected modes, which can be expressed as a sum of a component that is parallel to the surface s , and a vertical component with the direction of surface normal n: s τ = s + hn. The s is fixed to the incident wave's value, due to the Snell's law at interface. The magnitude h of the vertical component is determined by plugging s τ into equation (1), where the wavevector becomes slowness if divided by ω 2 : The above equation is a general eigenvalue problem for h and can be solved as in reference [17]. We developed software to obtain six solutions of complex h values in three pairs, for any general direction of s and n. A real pair of h corresponds to two solutions with opposite energy flow onto the interface: one with w · n > 0 for reflection, and another with w · n < 0 for transmission to the next medium. Note that the direction of energy follows the direction of group velocity w. A complex pair of h makes two solutions with non-zero imaginary parts Im(h) of opposite sign, that represent the evanescent waves. In the convention of this paper, the wave is defined with phase i(−q · r + ωt). Thus lm(h) < 0 (> 0) correspond to the reflection (transmission) modes, since the wave should be evanescent as r runs further from the interface in both cases.
There will always be three solutions of h for the reflected wave, whether it includes the evanescent wave or not. The corresponding three polarization vectors which could be complex, are then the eigenvectors in equation (4).
The reflection coefficient R τ for reflected mode τ is defined as the ratio of two energy flow per unit surface area S · n for reflected and incident waves [17]. Here the amplitude squared |Γ τ | 2 should be multiplied to the unit energy flow vector for the reflected waves, where we use a minus sign to balance the in-and-out opposite sign of two surface energy flows:

Analysis of reflection peak structures
Though the reflection model is elegant theoretically, its validity should be tested in application to real experimental data.
In setup A of reference [15], phonons travel from an aluminum radiator to a tin absorber through a sapphire crystal to give energy deposit signals, using the experimental setup in reference [13]. A large portion of phonons travel through a reflection at the opposite surface with normal n, as shown schematically in figure 1. The circular shaped radiator (0.4 mm 2 ) and absorber (0.5 mm 2 ) are point-like, when compared with the 2 mm distance between them and the crystal thickness of 9.53 mm. Consequently, phonon energy spectrum shows sharp peaks at arrival times for several reflection transition modes. We define the nominal paths as a point to point phonon travel path that starts from radiator center to absorber center with one reflection, whose arrival times are representative among the peaks. As from figure 2 of reference [15], four distinct reflection peaks are clearly seen above the broad shape from scattering, for both the experimental and simulation spectra. Figure 2 explicitly singles out the four reflection peaks selected from the information of the full simulation.
In this paper, the reason for different size of reflection peaks will be explained using the reflection coefficient and other propagation characteristics for phonons in anisotropic media. The elastic tensor in Voigt's notation and density used for the sapphire crystal oriented in figure 1 are: C 11 = 4.952, C 12 = 1.654, C 13 = 1.130, C 14 = −0.232, C 33 = 4.932, C 44 = 1.488 in units of 10 11 Pa, and ρ = 3.986 g cm −3 from figures 6(c) and 10 of reference [6].
In the YOZ plane, the slowness and group velocity of incident and reflected phonon lie in this same plane due to the symmetry of elastic tensor. Since the center of the radiator and absorber align along the z axis, the nominal path of the travel  [15]. The labels indicate different phonon reflection transitions. Here the two transverse modes are classified as SV and SH using polarization directions, which are discussed in the text. lies in this YOZ plane too. In figure 3, three-dimensional slowness surfaces reduce to slowness curves on YOZ plane, while s and direction of w can be plotted without projection effect.
Here, it is more convenient to classify two transverse modes according to their polarization vector directions, that either lies in or is orthogonal to the YOZ plane. We label the former as shear vertical (SV) and the latter shear horizontal (SH) modes. The SH mode is decoupled from the other two modes, because the polarization vector is always perpendicular to them. For SV and L modes, they can reflect to both of themselves. So that we have a total five possible reflection modes out of the general nine reflection transitions. In figure 3, exact slowness positions and group velocity directions are plotted as lines and arrows for both incident and reflected phonons in five nominal paths. The Snell's law requires the parallel component of the slowness to be the same for the incident and reflected waves, which can be confirmed visually in the plot. The five travel paths and group velocities are all different, even with common starting and ending point.
When the group velocity points toward the vacuum, the slowness curve is drawn thinner, otherwise it is thicker in figure 3. For any incident wave, the energy flow should always hit the surface to the vacuum, i.e., the slowness should lie in the thin part of the curves. Now let us take a close look at the very small θ region for the SV slowness curve, marked with 'SV reversed' in figure 3. Here the wavevector points toward the vacuum, but the group velocity points toward the crystal, as represented by the thick line segment. When calculating the reflection coefficient in this small region, the incident wavevector direction is reversed. In the SH slowness curve, similar region also exists.
Of course, the peak sizes in figure 2 are directly proportional to the reflection coefficients. Figure 4 shows four reflection coefficients from equation (5) involving SV and L modes, as a function of incident wavevector angle θ. Note that the coefficient sum of all reflected modes for a certain incident mode is unity, as required by the law of energy conservation.  Besides the obvious reflection coefficient, the size of a single reflection signal depends on the following additional factors: • Density of states (DOS) factor for initial phonon mode σ from radiator emission. This is expressed as (−n · w σ (q)) s 3 σ (q). Here the factor in the bracket is the cosine law, applied as the projection of group velocity on the normal of the radiation surface −n pointing toward the crystal when viewing in the radiator. The second factor s 3 σ accounts for the proportionality of the DOS to cubic power of slowness magnitude for a certain wavevector direction q.
• Transmissivity for the initial phonon. The initial phonon intensity also depends on the total energy transmissivity through the metal radiator to the crystal through the interface. Due to detailed balancing [11], it is the sum of reverse transmission coefficients for a certain phonon mode σ in the crystal, to all modes σ in the radiator metal, expressed as 3 σ =1 t (01) σσ (−q). Here t represents the transmission coefficient for transmission from mode σ in crystal (0) to mode σ in radiator (1), assuming the initial phonon is traveling backward with a reverse wavevector −q. It is the ratio of transmitted to incident energy flow, where the amplitude ratio is determined from the equations for continuity of displacement and stress vector across the interface 1 .
• Geometric factor (G). Due to the finite size of absorber, incident phonons' wavevectors have to be inside a geometric solid angle Ω q to hit the absorber by one reflection.
Here the solid angle should be evaluated in the initially isotropic wavevector q space instead of the group velocity space. Total number of phonon the absorber can receive is proportional to the area of the radiator and the size of Ω q . A given infinitesimal solid angle of phonon wavevector ΔΩ q in radiation has a corresponding hit area ΔA in the absorber for the travel path with a single specular reflection, as illustrated in figure 1. The solid angle per hit area value Ω A = dΩ q /dA then accounts for the geometric effect from the size of Ω q , since the total area of absorber is fixed. For an isotropic medium of thickness h and a distance d between radiator and absorber, a samemode transition goes through a simple mirror reflection with equal incident and reflected angles. This value is then However, in case of a general anisotropic medium, the value Ω A is affected first by the phonon focusing effect [16] of group velocity, and then by the complicated relationship between the incident and reflected wavevectors which depends on two slowness surfaces and the interface's direction. We shall proceed the calculation numerically. By applying a change of ±0.01 rad in the θ and φ coordinate around the nominal incident direction, we obtain the corresponding change of the hitting point Δr θ and Δr φ in the absorber surface. Then the infinitesimal solid angle and hit area are calculated as ΔΩ q = sin θΔθΔφ and ΔA = |Δr θ × Δr φ |. Finally Ω A value is then ΔΩ q /ΔA. We define a geometric factor G = Ω A /Ω A,iso , that measures Ω A in terms of the isotropic value Ω A,iso for the same geometry. If G > 1, the energy flow is focused, otherwise it is defocused. 1 In case of transmission from medium (0) to (1), the relative amplitudes of mode τ for reflection Γ (0) τ and transmission Γ (1) τ are determined by the following two equations.
Here the first and second equation are continuity of displacement and stress vector, respectively. a. This is evaluated at the nominal path at φ = 0 in figure 6. b. Since there are multiple paths, geometry factor at nominal path is meaningless. c. Total factor value for SV → SV mode is computed in section 4, instead of simple product of the individual values.
The above procedure for calculation of the factor G works for most cases. However, it breaks down when the relationship between the hit point and the incident wavevector direction can not be approximated to be linear in the region of travel, such as for the SV → SV mode discussed in the next section.
• Absorption factor. A phonon that hits the absorber can transmit into the absorber metal with a certain modedependent probability, and then goes through attenuation and multiple reflection to be absorbed or re-transmitted back. The total absorption factor is discussed in reference [15]. Note that reference [11] did not consider this factor in calculation. Our analysis provide the key clue of signal structure, compared to the 'black-box' type output from a full simulation. Here, a direct link between the observed peak structure to the phonon reflection and transportation theory is made. This study also provides a novel method to directly use the energy spectrum structure to understand the underlying physical properties of the material.

Multiple travel path for the SV → SV transition
The calculation for the Ω A and all other factors in the previous section works since they can be assumed as a constant for a focused direction of travel due to the small size of radiator and absorber. This section deals with the SV → SV mode where the above assumption breaks down. Figure 5 shows the x travel distance from emission to absorption as a function of angle φ of the initial emitting wavevector, with θ fixed at the nominal value. We see that the hit point makes an 'S' shape as φ changes. It indicates that there are three φ values associated with a same hit point at small x distance.
As long as the x distance of phonon travel is less than 0.756 mm, the radius sum for radiator and absorber, it will have a chance to travel from radiator to absorber and contribute to the signal. The shaded region of Δx ∈ [−0.756, 0.756] mm  in figure 5 shows that for a wide range of φ ∈ [−0.2, 0.2] rad, the phonon from radiator is able to hit the absorber. As a result, we need to consider this wide range of φ in calculations of the geometry effects. It then turns out that the dependence of other four factors in table 1 on φ should be studied too, since all of them can have a substantial change due to large change of φ. Among those factors, we found that the reflection coefficient has significant dependence on φ, as shown in figure 6. While other factors can be considered as a constant even in the wide φ range.
To understand the strange behavior in figure 5, we plot the projection of slowness s and group velocity w on the XOY plane for the incident phonon, as shown on figure 7. Here the two vectors do not lie inside the YOZ plane for φ = 0. We notice the slowness curve has a saddle-like shape with a dip at nominal direction (φ = 0). As |φ| increases from zero, x components of the group velocity and the slowness have opposite sign at small value of |φ| (such as 0.08 rad), before returning to normal behavior of same sign for sufficient large |φ| value. The same argument holds for the reflected waves. Thus, this leads to the non-monotonic change of the hitting point at x direction.
Since there exists strong dependence of the x distance and reflection coefficient over φ in the SV → SV transition, a simple calculation in the previous section is not appropriate. First, we point out that all the factors in table 1 only depend on the relative distance between the source point (x 1 , z 1 ) and the destination point (x 2 , z 2 ) in x and z direction: Δx = x 2 − x 1 and Δz = z 2 − z 1 . The most general formula for the total factor f tot in table 1 is the average of all the factors for both the radiator (1) and absorber (2), as an integration over both active areas divided by the total areas A 1 and A 2 of them: Note that the integration range is for both the source point (x 1 , z 1 ) and the destination point (x 2 , z 2 ). The term f other represents the product of all the factors appearing in table 1, except the geometric factor G = Ω A /Ω A,iso , which is explicitly shown. In normal cases for small area of A 1 and A 2 , the integration is not necessary, since an evaluation of the factors at the nominal path is sufficient to represent the averaged f tot factor. To proceed the calculation of equation (6), we first simplify it by separating the dependence of (Δx, Δz) to one of φ or θ. Due to the symmetry on both side of the YOZ plane at small |Δx|, we have a simplified dependence between the polar angles of the incident wavevector and (Δx, Δz): Δz only depends on θ and Δx only depends on φ. So we have Ω A = dΩ q /dΔxdΔz = sin θ · dθ/dΔz · dφ/dΔx, where the term dφ corresponds to the sum of all possible dφ values that can give a certain Δx value.
Unlike the strange behavior of Δx − φ dependence, the dependence of Δz on θ is normal. So we can fix θ and Δz to the value for the nominal path θ 0 and Δz 0 , for the small source and destination area. Since now there is no dependence on Δz = z 2 − z 1 for all the terms, the variable z 1 can be integrated out to give a factor as 2 R 2 1 − x 2 1 , and f other (Δx, Δz) becomes f other (Δx). Then, we expand the G factor to have: Due to the multiple φ to single Δx correspondence, a φ instead of a Δx value uniquely defines the travel path. So instead of binning in Δx, we scan and divide the azimuthal angle φ into bins to proceed the integration as a summation: As usual, the term dφ/dΔx| i is the ratio of the range for φ and Δx in ith bin, serving as an approximate value of dφ/dΔx for the bin. For a change in φ value from low to high bin boundary for ith bin, the change in Δx value is recorded in the range [B 0i , B 1i ]. The term Δx in ith bin dx 2 dz 2 now evaluates to the geometric area of a slice in a circle of radius R 2 with the x coordinate constrained between x 1 + B 0i to x 1 + B 1i . For the SV → SV mode in table 1, the total factor is directly computed using equation (8), while all other factors are only reference values evaluated at the nominal path. The bin size of φ used in calculation is 0.01 rad and the corresponding bin ranges of Δx can be viewed from figure 5.

Conclusion and outlook
The theory of phonon reflection on crystal surface at low temperature is known. However, full explanation of the real experimental data has not been accomplished in a systematic and quantitative way.
We developed a model on the specular phonon reflection at a general surface of anisotropic medium, by applying acoustic mismatch principle. And the strength of reflection is demonstrated to be a combined effect of reflection coefficients and other phonon travel properties. In the previous approach to explain the reflection peak sizes on reference [13], the geometry factor was not considered, nor was the reflection efficiency correctly calculated. The simple scalar product of incoming and reflected polarization vectorsê in ·ê ref used as reflection efficiency in [13] deviates from the correct value; for instance, in SV → L and L → SV reflections,ê in ·ê ref evaluates to 0.007 005 and 0.254, while the correct values from column 2 of table 1 should be 7.354 × 10 −5 and 0.0928. Thus, we provide a novel framework to analyze the peaks in a solid quantitative way from phonon theories.
This model also provides an efficient calculation of precise reflection coefficients and group velocity in real time, and eliminates the cumbersome procedure to be obtained as interpolated approximate values from a preloaded huge lookup tables or maps, as in some other applications [12,18]. It recorded nice performance in the application of full phonon transport simulation [15].
A comprehensive phonon propagation simulation using this generic reflection model and other transport processes, such as phonon scattering and down-conversion, will enable the community to quantitatively model observed signal shape for any crystal and geometry. For example, phonon arrival times as a function of calibration source position [19] at the SuperCDMS detector, and a strong position dependence in the CaF 2 (Eu) bolometer [20], which should be contributed from the event location of initial energy deposit, can be studied in detail using the simulation.
The development of this model comes from the effort to understand the mechanism of cryogenic particle detectors, through the dedicated phonon experiments of solid state physics. As a consequence, this model will have application in rare event search experiments in particle physics, as well as in general material science related problems.