Imaging grain microstructure in a model ceramic energy material with optically generated coherent acoustic phonons

Characterization of microstructure, chemistry and function of energy materials remains a challenge for instrumentation science. This active area of research is making considerable strides with methodologies that employ bright X-rays, electron microscopy, and optical spectroscopy. However, further development of instruments capable of multimodal measurements, is necessary to reveal complex microstructure evolution in realistic environments. In this regard, laser-based instruments have a unique advantage as multiple methodologies are easily combined into a single instrument. A pump-probe method that uses optically generated acoustic phonons is expanding standard optical characterization by providing depth resolved information. Here we report on an extension of this method to image grain microstructure in ceria. Rich information regarding the orientation of individual crystallites is obtained by noting how the polarization of the probe beam influences the detected signal amplitude. When paired with other optical microscopies, this methodology will provide new perspectives for characterization of ceramic materials.

O ptical spectroscopies have the demonstrated capability to noninvasively characterize several aspects of advanced energy materials in realistic environments. Examples include monitoring function of organic photovoltaics using photoluminescence spectroscopy, measuring in situ changes of chemistry of fuel cell anodes using Raman spectroscopy, monitoring thermal transport across individual interfaces, and investigation of hot carrier cooling mechanisms in photovoltaics using transient absorption spectroscopy [1][2][3][4] . Additional attributes include benchtop deployment and straightforward accommodation of multiple measurement methodologies using standard optics. Adding to this list of optics-based characterization tools is a technique that uses optically generated and detected acoustic phonons to image materials in the depth direction 5,6 . This method, termed time domain Brillouin scattering (TDBS), had been first applied to characterize homogeneous semiconductors and dielectrics 7 . Applications were extended to homogeneous liquids 8 , vegetal cells 9 , and polycrystalline materials [10][11][12] .
Applications were further extended to depth profiling of elastic inhomogeneities in nanoporous organosilicate 13 and ion irradiated GaAs 14 and to three-dimensional imaging of animal cells 15 . Recent work that is closely tied to the current study involves surface scanning using TDBS to image subsurface grain boundary orientation 16 and to resolve the position of different crystallites 17 . A review of other applications of TDBS can be found in reference 18 .
Although the velocity of longitudinal acoustic (LA) phonons are traditionally used for characterization, recent studies have demonstrated the potential of using transverse acoustic (TA) phonons to obtain additional information regarding grain microstructure 10 . Key examples include using TA and LA phonon velocities to obtain information about the first two Euler angles 17 and using surface acoustic wave velocities to obtain the orientation of large crystallites having dimensions of a few tens of millimeters 19 . However, measurement of bulk wave velocities alone cannot resolve the third Euler angle, which rotates the crystallite around the surface normal [20][21][22] . In addition to measuring phonon velocities, it is expected that detailed studies of the detection process will reveal additional information about grain microstructure 10 .
Here, we systematically study the detected amplitude of TA and LA phonons as a function of probe polarization. CeO 2 (ceria) was chosen for this study as it is a model fuel cell material owing to its mechanical stability and ionic conductivity 23,24 . Ceria is also an important surrogate material for oxide nuclear fuels 25,26 .
Applying TDBS to ceria presents an excellent test case for research on advanced ceramic materials with applications in the energy industry. As discussed in detail below, our experimental results show that the detected TA phonon amplitude is a strong function of probe polarization. The probe beam polarization dependence is explained by an optoacoustic model that considers the strain associated with each acoustic mode and the probe beam polarization. Based on our finding, a method is suggested for the determination of the orientation of individual crystallites. When coupled with subsurface imaging, this approach has the prospect to fully characterize the near surface grain microstructure in three-dimensions (3D). Moreover, combining TDBS with other optical microscopies has the potential to greatly extend opticsbased instrumentation science with broad application involving 3D characterization of polycrystalline ceramic materials.

Results
Crystallite mapping and orientation measurement. The experimental geometry is shown in Fig. 1a. A thin gold film is coated on the sample surface to ensure strong optical absorption of the laser pump beam. A detailed description of experimental and material parameters is given in the methods section. The crystallite orientation at the surface of the sample, recorded using electron backscatter diffraction (EBSD), is shown in Fig. 1b. A platinum fiducial square, indicated by the white arrows in the EBSD image, is used to optically locate specific crystallites and to define the coordinate system. As indicated, x 1 and x 2 are aligned with the sides of the fiducial square and x 3 is antiparallel to the surface normal. The skewed appearance of the fiducial square is related to projection issues associated with EBSD. Several representative crystallites (A-F) were identified for this measurement. The orientation of the crystallites investigated is listed in Table 1. Using metallurgical terminology, we define the direction normal to the plane (hkl) as the normal direction (ND), and a direction parallel to [uvw] as the rolling direction (RD). The relation between ND/RD and the three Euler angles is given in Fig. 1c.
Time domain Brillouin scattering spectroscopy. An example of a reflectivity transient exhibiting strong Brillouin oscillations is shown in Fig. 2a, corresponding to crystallite A. To highlight the oscillatory features, the thermal background has been subtracted. The data exhibit two oscillatory components having different frequencies and amplitudes.
[010] 1 2 x 1 x 1 x 1 x 3 x 2   Table 1. The platinum fiducial square, delineated by white arrows, is used to define the laboratory coordinate system. c The three Euler angles are used to define the crystallite orientation relative to the laboratory coordinate system. The first rotation, φ 1 , is about [001] and the last rotation, φ 2, is about ND.
The observed Brillouin oscillations stem from optical interference between a portion of the probe that reflects from the gold film and a portion that reflects from thermoelastically generated phonons propagating in the ceria substrate (see Methods). The frequency of these oscillations is governed by f = 2nν m /λ, where n is the refractive index, ν m is the mode specific phonon velocity and λ is the photon wavelength of the probe pulse. Ceria is a cubic crystal and optically isotropic (i.e., only one electric permittivity tensor component). However, for birefringent materials, multiple frequency would be recorded for each acoustic mode 10,12 . The Fourier amplitude spectrum shown in Fig. 2b, reveals two phonon modes, located at 93.4 GHz and 48.3 GHz that are associated with the quasi longitudinal (LA) and fast quasi shear acoustic (FTA) modes, respectively. The slow quasi shear mode (STA) is not detected.
The acoustic velocities are obtained by calculating eigenvalues of the Christoffel equation using the rotated elastic stiffness tensor 16 . The Euler angles used to derive the rotated elastic stiffness tensor are φ 1 = −63°, Φ = 25°, and φ 2 = 86°. The refractive index of ceria at the probe wavelength (408 nm) was measured to be 2.57 using optical ellipsometry. The imaginary part of the index of refraction is negligible giving rise to long-lived Brillouin oscillations. Measurement of the single crystal elastic constants either requires large single crystals or methods applied to polycrystals that can measure elastic properties on length scales below the crystallite diameter. It is well known that single crystal elastic constants can be obtained by frequency-domain Brillouin scattering measurement (see, for example, 27,28 and the references therein). The elastic constants of ceria used in this study were reported in previous work based on time domain Brillouin measurements made on several crystallites with known orientations 17 . Ceria has three independent elastic tensor components represented using engineering notation as C 11 , C 12 , and C 44 . We find good agreement between the experimental and calculated Brillouin frequencies for each detected mode and for each crystallite investigated as summarized in Table 1. An expression for the Brillouin frequency as a function of phonon velocity and index of refraction is given in the methods section.
Probe polarization dependent measurements. Next, we investigate the influence of the probe beam optical polarization on the measured signal. A half-wave plate was used to rotate the probe polarization counterclockwise in steps of 10°. The power of the pump and probe beams were kept constant. The polarization angle, Ω, is defined as the angle between the optical polarization vector and the x 1 axis (see Fig. 1a). The amplitude profiles of the detected acoustic modes for crystallite A as a function of Ω are shown in Fig. 3a. The major axis (connecting points of maximum amplitude) is perpendicular to the minor axis. The major/minor axes coincide with axes of mirror symmetry (represented as double-arrow lines in Fig. 3). The FTA amplitude is normalized using the largest LA amplitude. The remaining discussion will only consider the amplitude profiles of the TA modes as they exhibit a stronger dependence on probe polarization.
The normalized TA amplitude profiles for other crystallites is shown in Fig. 3b-d. The amplitude profiles indicate twofold rotation symmetry, which is consistent with previous observations of TA modes in BiFeO 3 12 . However, it is noted that there is a slight left/right asymmetry in the experimental data. This asymmetry is not reproducible and we thus attribute it to an experimental artifact such as beam walking of the pump relative   NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15360-3 ARTICLE to the probe caused by rotation of the λ/2 plate. Using other data sets, not presented in this manuscript, we find that these asymmetries minimally influence the extracted value of φ 2 . For most crystallites, only the FTA mode was detected. However, in crystallite D, both the FTA and STA modes were observed. No shear modes were detected in crystallites E and F as both correspond to high symmetry directions.

Discussion
LA acoustic pulses are generated thermoelastically in the elastically isotropic gold film (see methods). TA acoustic pulses are generated owing to mode conversion at the film/substrate interface [29][30][31] . The detected amplitude profile of each acoustic mode is a product of components of the strain and photoelastic tensors (see methods), the latter being a strong function of probe wavelength. A thorough literature search revealed no data for the photoelastic tensor components for ceria. Nevertheless, a set of measured profiles over various crystallites (Fig. 3) allowed us to determine the relative values of photoelastic coefficients for ceria at 408 nm. We compare the theoretical TA and LA amplitudes with the experimental results by adjusting the photoelastic ratio ([P 11 −P 12 −2P 44 ]/P 12 ). We obtain an optimal fit by adopting the following ratio: (P 11 −P 12 −2P 44 )/P 12 =1.2. We note that calcium fluoride, with a similar fluorite crystal structure, has (P 11 −P 12 −2P 44 )/P 12 = 1.1 32 . In the following discussion, we show that the dependence of the amplitude profiles of the TA modes on the photoelastic constants is only a function of this ratio. Theoretical and experimental results for several crystallites are compared in Fig. 3. The theoretical amplitude profiles (see methods) are presented as solid lines. Overall, we find a good agreement between experimental and modelling results. The modelling shapes all have twofold rotation symmetry and exhibit a dipole-like feature. This is in agreement with the structural invariance after rotating the third Euler angle by 180°.
In most cases, involving crystallites aligned along low symmetry directions, only the LA and FTA modes are detected. This is confirmed using our theoretical model (see methods) that calculates relative shear strain displacement for all possible crystallite orientations (shown in Fig. 3e, f). Generally, the FTA mode has a much larger strain than the STA mode and can be observed for most orientations. Moreover, to efficiently mode convert the LA mode in the film into TA modes in the substrate require that the TA modes have a sizeable longitudinal component. For most orientations, the FTA mode has a relatively large longitudinal component as compared to the STA mode. The STA mode can only be detected in the vicinity of the (111) orientation (red region in Fig. 3e). Crystallite D is located in this region and as expected we detect both shear phonon modes. For a range of probe polarizations satisfying 40°< Ω < 90°and 220°< Ω < 270°, the amplitude of the STA mode is larger than the FTA mode. The STA also exhibits additional lobes around the minor axis.
The photoelastic tensor governs the amplitude profile of the FTA mode. Using engineering notation, the rotated photoelastic tensor is represented by a 6 × 6 matrix, P ij . The components of interest for x 1 polarized light are P 13 ; P 14 ; and P 15 and P 23 ; P 24 ; and P 25 for x 2 polarized light. Examination of these terms reveals that they have the following form: where Δ i (φ 1 , Φ, φ 2 ) are constants defined by the crystallite orientation. P 23 ; P 24 ; and P 25 share the same form, although Δ i terms are different. The fact that all terms in Eq. (1) have a common factor P 11 −P 12 −2P 44 shows its significance in the detection process. Its value physically quantifies the coupling strength between the shear strains of the acoustic phonon modes and the electric field of the polarized probe beam. Therefore, the detectability of TA phonon is strongly related to the amplitude of this factor. In reality, by rotating the probe polarization, the second term on the right-hand side of P 13 has at most a comparable magnitude to that of the first term. Therefore, rather than giving the values of P ij , we fit the ratio (P 11 −P 12 −2P 44 )/P 12 . From this observation, it is possible to study P 12 alone if the crystallite has an orientation where the Δ i terms are very small. A similar experiment has been undertaken to study the evolution of P 12 in proton irradiated GaAs 33 . Here, we define the anisotropy of the photoelastic tensor α r ¼ 2P 44 =ðP 11 À P 12 Þ, in analogy to the Zener ratio for the elastic tensor. In the case of photoelastic isotropy, α r = 1, only P 13 ¼ P 12 is nonzero and the detected amplitude does not depend on probe polarization. Photoelastic isotropy has been observed in several alkali halides at certain wavelength where optical birefringence does not depend on the strain 34 . Thus, the dependence on the optical polarization of the probe beam, shown in Fig. 3, is a consequence of photoelastic anisotropy. It is expected that as the photoelastic tensor becomes more anisotropic, or as α r becomes more distinct from unity, the optical polarization dependence of the detected LA and TA phonon amplitudes becomes more pronounced. The shape of mth mode amplitude profile can be written as a function of the strain induced change in the dielectric permittivity. Using engineering notation, the relative reflectivity change can be expressed as 35 δR R / Δε 1 cos 2 Ω þ Δε 2 sin 2 Ω þ Δε 6 sin2Ω ð2Þ Modulation of the dielectric permittivity, Δε 1 , Δε 2 , and Δε 6 , are calculated for a given crystallite orientation. Equation (2) From Eq. (3), it is straightforward to see that the polarization dependence has a dipole shape and that it has the largest amplitude (N+M) at Ω max ¼ πÀ2Ω 0  Fig. 3a-d. Although the shape of the theoretical amplitude profiles closely follows the experimental data it is noted that the magnitude of the amplitude profile in some cases is over/under estimated. This is due to the fact that we used a global fit, using all of the data acquired, to extract the photoelastic ratio. In the case where N < M (STA mode in crystallite D), the relative reflectivity given in Eq. (3) will experience a crossover from positive to negative values during rotation, and the absolute value will reach a minimum value at angles defined by N þ M sin 2Ω þ Ω 0 ð Þ¼0. As a result, an additional set of lobes are observed near the minor axis of crystallite D shown in Fig. 3d. In both cases (N > M and N < M) the responses are mirror symmetric along the major axis defined by Ω 0 and have twofold rotation symmetry.
It is important to note that the parameter Ω 0 is a strong function of crystallite orientation (φ 1 , Φ, φ 2 ). Consider the FTA mode in crystallite A as an example. From Eq. (4) we calculate Ω 0,2 = 148°, indicating min/max amplitudes of the FTA mode to be along 61°/−29°directions. This agrees very well with the result presented in Fig. 3a. The same calculation was also conducted on crystallites B, C, D and is summarized in Table 1 (negative angles are located in the third and fourth quadrants). Experimental results and model prediction show close agreement (all cases are within 7°).
Next, we discuss how the Euler angles affect the amplitude profiles of the FTA mode. The shape is determined by the first two Euler angles, φ 1 and Φ. We categorize the shapes into three groups: capsule shape, dipole shape, and quadrupole shape. These shapes are presented in Fig. 4a-c, respectively. In the first case, the shape does not have an eminent dip near the minimum, which corresponds to N ≫ M. The second group corresponding to N > M, is the most common. The last group is the least common and corresponds to N < M. In Fig. 4d, we plot the M and N ratios and their corresponding shape onto a color-coded inverse pole figure for FTA mode. The majority of the orientations produce a dipole shape, in agreement with our experimental results. It is also noted that the third Euler angle φ 2 by definition rotates the crystal around x 3 axis (ND) and thus does not change the dipole shape, only its orientation. This point is illustrated in Fig. 4e, where two profiles that have the same values of φ 1 and Φ but different values for φ 2 (solid and dashed line) are compared.
The result of these findings suggests a method for the complete crystallite orientation using TDBS. The modulated optical permittivity components, Δε k;i , implicitly depend on three angles (φ 1 , Φ, φ 2 ) and elasticity (C ij ) as well as photoelasticity P 12 and P 11 À P 12 À 2P 44 ð Þ . In principle, all three angles can be obtained from measurement of Brillouin frequencies and the FTA amplitude profile. Applying this approach to grain C, we obtain φ 1 = −67°, and Φ = −23°from the measured Brillouin frequencies. In the next step we calculate the sum of the square of the differences, χ 2 , between the experimental and theoretical amplitude profiles versus φ 2 (Fig. 4f). The minimum value of χ 2 corresponds to φ 2 = 91. These angles agree closely with the values obtained from EBSD, φ 1 = −64°, and Φ = −25°and φ 2 = 85°. The shape of χ 2 is fairly steep and exhibits a narrow minimum, suggesting high sensitivity to φ 2 . It is noted the procedure sketched above requires knowledge of the photoelastic ratio, P 12 /(P 11 −2P 44. ). Because information on the photoelastic constants for most materials is scarce, the photoelastic ratio can be obtained by making measurements on several crystallites with known orientation and fitting predicted polarization plots to the measured ones by varying the photoelastic ratio. In addition, although the STA mode is not required to uniquely determine all three Euler angles, it is expected that the accuracy of this approach would be even higher for orientations for which the STA mode is generated and detected.
To place this new application of TDBS in the proper perspective it is important to understand current limitations. First, this technique is only appropriate for materials that are transparent or partially transparent at the probe wavelength. In addition, for materials that are transparent at the pump wavelength, a thin metallic film must be applied to ensure strong optical absorption. The second limitation is related to spatial resolution. Currently the lateral spatial resolution of the technique is limited by the overlap of the foci of the pump and the probe laser pulses. Using far-field focusing techniques, a spot size of 1 µm is easily achievable, and spot sizes of~100 nm can be obtained by employing near-field focusing techniques 36,37 .
Although these values are at least an order of magnitude larger than that of scanning electron microscope-based methods, this approach does not require serial sectioning to obtain depth resolved information and provides volumetric information over a larger depth than accessible using focused ion beam liftout techniques. Lastly, this technique uses anisotropy in material's property to extract crystallite orientation. As a consequence, the sensitivity to orientation increases with increasing elastic and photoelastic anisotropy.
In conclusion, we applied TDBS to characterize grain microstructure in a polycrystalline ceria sample. It was shown that for crystallites oriented along low symmetry directions, both LA and TA modes are generated and detected. The generation of the TA modes is due to mode conversion of the LA mode at the interface between the elastically isotropic gold transducer film and the elastically anisotropic substrate. The TA amplitude profiles exhibit a strong dependence on the direction of the probe polarization vector. Using an optoacoustic model that accounts for elastic and photoelastic anisotropy, it is shown that the shape of the amplitude profile for the FTA mode falls into one of three categories (capsule, dipole, quadrupole). A method was suggested to extract all three Euler angles by using the orientation of the major axis of the FTA amplitude profile and the velocities of the detected phonon modes. While outside the scope of the current manuscript, future work will involve refinement of this approach by considering the shape, orientation and relative magnitude of the amplitude profiles for all of the detected phonon modes.
The results presented in this manuscript are important for a broad class of applications involving 3D imaging of polycrystalline ceramic materials relevant to the energy industry. Energy materials are routinely exposed to extremes in temperature, pressure and radiation. As a consequence, further materials development will require instruments that can record multiple aspects of material performance in realistic environments. The current work involving the development of non-invasive imaging capability of 3D microstructure is an important step in this direction. Future work aimed at combining 3D imaging using TDBS with other optical microscopies will provide synthesis scientists with essential information regarding in-situ and/or inoperando changes in material structure and function.

Methods
TDBS measurement. The sample investigated was a 3 × 4 × 5 mm sintered ceria pellet, purchased from AlfaAesar. The as-received sample was annealed at 1600°C for 4 h in an oxygen rich environment to prevent reduction and promote grain growth. The average grain size of the annealed sample was 50 µm. The sample was polished to be optically flat and a fiducial platinum square measuring 500 × 500 µm 2 was deposited using a focused ion beam. The last preparation step involved sputter coating a thin (~7 nm) gold transducer film to ensure strong optical absorption and to promote acoustic mode conversion at the film/substrate interface. The orientation of the crystallites at the sample surface were imaged using EBSD.
An optical pump-probe setup was used to measure the propagation of acoustic phonons in the ceria substrate. A mode-locked Ti-sapphire laser pump beam of pulse duration~1 ps, repetition rate 76 MHz and wavelength 816 nm is focused onto the sample with a spot size of~2 µm and a fluence~0.01 mJ cm −2 . The pump beam is amplitude modulated at 1 MHz using an acousto-optic modulator. A timedelayed, normal incident 408 nm probe with similar fluence and beam waist as the pump beam is focused onto the sample in the same location as the pump beam. The time delay was imposed by a mechanical shaker with a maximum scan time of 350 ps. Small changes in the signal intensity governed by the photoelastic effect are demodulated using a lock-in amplifier. To ensure strong photoelastic coupling, the output wavelength of the laser is selected such that the probe photon energy is close to the ceria band gap 38 . Optical attenuators are used to control the incident powers of the pump and probe beams. A half-wave plate is placed in the probe leg to rotate the probe polarization. The sample is placed on a 2D translation stage, with micron resolution, to precisely locate the grains.
Optoacoustic model. The acoustic pulses propagating in ceria are thermoelastically generated in the gold film. Absorption of femtosecond laser pulses in gold is known to create a non-equilibrium distribution of overheated charge carriers, which transfer energy to the lattice (to the thermal phonons) in a characteristic time τ E ≈ 0.5 ps 39 . Because the thickness of the gold film (H =7 nm) is much shorter than the diffusion length of the non-equilibrium electrons, the gold film is heated uniformly in the depth direction 39 . The strain pulses in ceria have a characteristic duration equal to τ a ¼ 2H v Au % 4 ps, and are composed of two segments of the same shape but opposite polarity. The acoustic waves propagating in the direction of ceria substrate contribute to the leading edge of the strain pulses, whereas waves that are reflected from the mechanically free surface of the gold film contribute to the trailing edge. The characteristic duration of the emitted strain pulse τ a is much longer than both the pump laser pulse duration and τ E . Thus, the heating and the accompanying stress initiation can be considered as instantaneous. Cooling of the gold takes place on a time scale exceeding 50 ps, thus thermo-elastic stresses associated with cooling provide negligible contribution to the high frequency acoustic waves monitored in our experiments (see Fig. 2b). The strain pulses emitted in ceria (x 3 ≥ 0) can thus be modeled accurately using the following spatial profiles: In Eq. (6), θ represents the sign function, v m is the velocity of the acoustic mode and m is a summation index accounting for the three acoustic modes. The shear polarized modes are launched owing to mode conversion at the film/ceria interface (x 3 = 0) of the purely LA pulse generated in the gold film. The index i in Eq. (6) indicate one of the three nonzero strain components in the emitted plane acoustic waves. The amplitudes A im of the strain components are all proportional to the photo-induced temperature rise in the gold film, its thermal expansion coefficient and mode conversion coefficients. The mode conversion coefficients are a function of the rotated elastic stiffness tensor for each crystallite and are obtained from the elastic boundary conditions at the film/substrate interface 29,35 . Owing to good acoustic matching between gold and ceria, we have neglected the emission of the acoustic waves in ceria at times t > τ a caused by multiple reflections within the gold film. We estimated that the displacement amplitude of acoustic waves emitted at t > τ a is < 15% of the incident wave amplitude for all possible orientations of the ceria crystallites.
Governed by the photoelastic effect, the strain field in ceria (described by Eq. (6)) induces a local change in permittivity, Δϵ ij ¼ p ijkl η kl . Similar to the elastic stiffness tensor, the rotated photoelastic tensor, p ijkl , is obtained using the Euler angles for each crystallite. Owing to the small amplitude of the acoustic pulses, the optical reflection coefficients associated with this change in permittivity can be evaluated using the single-scattering approximation 40 . The optical reflection coefficients r ij for the jth component of the incident electric field into the ith component of the reflected electric field is given by: where k p and n denote the wavenumber of the probe light in vacuum and the refractive index, respectively. For the 408 nm probe light, the imaginary portion of the refractive index is set to zero. At times t > τ a , when the acoustic field in Eq. (6) is completely launched/located inside ceria, integration of Eq. (7) can be done explicitly, demonstrating that light scattering, as expected, is induced not by a homogeneous strain but by strain gradients: r ij ¼ 2 p ijk3 ik p n e 2ik p nH sin k p nH h i 2 A km e À2ik p nv m t ð8Þ Here, the phase shift of the probe light scattered by each acoustic mode grows proportionally with penetration depth v m t in ceria ðφ m ¼ 2Reðk p Þnv m tÞ. Equation (8) also illustrates that only phonons with wave-vectors equal to twice the wavevector of the optical probe light in ceria are detected. Two components of the probe beam, one reflected from the film surface and one reflected from the phonon propagating in the ceria sample, mix at the photodetector. This mixing sets up Brillouin oscillations with a characteristic frequency ω B;m ¼ ∂φ m ∂t ¼ 2Reðk p Þnv m .

Data availability
All data generated in this study are available from the corresponding authors upon request.
Received: 15 August 2019; Accepted: 28 February 2020; Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.