nanomaterials Article Impact of the Regularization Parameter in the Mean Free Path Reconstruction Method: Nanoscale Heat Transport and Beyond Miguel-Ángel Sanchez-Martinez 1 , Francesc Alzina 1 , Juan Oyarzo 2 , Clivia M. Sotomayor Torres 1,3 and Emigdio Chavez-Angel 1, * 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and The Barcelona Institute of Science and Technology (BIST), Campus UAB, 08193 Bellaterra, Barcelona, Spain; miguel.sanchez.92@gmail.com (M.-Á.S.-M.); francesc.alzina@icn2.cat (F.A.); clivia.sotomayor@icn2.cat (C.M.S.T.) 2 Instituto de Química, Pontificia Universidad Católica de Valparaíso, Casilla 4059, Valparaíso, Chile; joyarzo89@gmail.com 3 ICREA-Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona, Spain  Received: 18 February 2019; Accepted: 6 March 2019; Published: 11 March 2019  Abstract: The understanding of the mean free path (MFP) distribution of the energy carriers in materials (e.g., electrons, phonons, magnons, etc.) provides a key physical insight into their transport properties. In this context, MFP spectroscopy has become an important tool to describe the contribution of carriers with different MFP to the total transport phenomenon. In this work, we revise the MFP reconstruction technique and present a study on the impact of the regularization parameter on the MFP distribution of the energy carriers. By using the L-curve criterion, we calculate the optimal mathematical value of the regularization parameter. The effect of the change from the optimal value in the MFP distribution is analyzed in three case studies of heat transport by phonons. These results demonstrate that the choice of the regularization parameter has a large impact on the physical information obtained from the reconstructed accumulation function, and thus cannot be chosen arbitrarily. The approach can be applied to various transport phenomena at the nanoscale involving carriers of different physical nature and behavior. Keywords: mean free path reconstruction; mean free path distribution of phonons; thermal conductivity distribution 1. Introduction In solid-state materials there is a variety of scattering mechanisms for energy carriers involved in different transport phenomena, such as impurities, boundaries, and collisions with other particles/quasi-particles. The average distance that a moving particle (photon, electron, etc.) or quasi particle (phonon, magnon, etc.) travels before being absorbed, attenuated, or scattered is defined as the mean free path (MFP). It is well known that energy carriers that propagate over different distances in a material, having different MFP contribute differently to the energy transport. Thus, the use of a single-averaged MFP may be inaccurate to describe the system [1,2]. It is possible to quantitatively describe how energy carriers with a specific MFP contribute to the total transport property by an MFP spectral function or MFP distribution [3], which contains the information of the specific transport property associated with the energy carriers with a certain MFP. By normalizing and integrating this spectral function we obtain the accumulation function, which describes the contribution of carriers with different MFPs below a certain MFP cut off to the total Nanomaterials 2019, 9, 414; doi:10.3390/nano9030414 www.mdpi.com/journal/nanomaterials Nanomaterials 2019, 9, 414 2 of 10 transport property, being very intuitive to identify which MFPs are the most relevant to the transport phenomenon under study by plotting this function. When studying transport at the nanoscale, boundary scattering becomes important as the characteristic size of the nanostructure approaches the MFP of the carriers involved. From the bulk MFP distribution, it is possible to predict how size reduction will affect a transport property in this material given that we know which MFPs are contributing the most. Inversely, it is possible to obtain the bulk MFP distribution from size-dependent experiments, where the critical size of each measurement acts as an MFP cut off due to boundary scattering. This relation between the transport property at the nanoscale and the bulk MFP distribution is given by an integral transform. A suppression function (SF), which accounts for the specific geometry of the experiment and depends upon the characteristic size of the structures and the MFP of the carriers, connects the bulk MFP distribution and the experimental measurements; acting in the kernel of the integral transformation. Using this integral relation, it is possible to recover the MFP distribution from experimental data. This is known as the MFP reconstruction method [4]. This mathematical procedure however is an ill-posed problem with, in principle, infinite solutions. To obtain a physically meaningful result from it, some restrictions must be imposed. These constraints are mainly related to the shape of the mean free path distribution. The distribution is a cumulative distribution function and it is subjected to some restrictions, e.g., the MFP distribution is unlikely to have abrupt steps because it is spread over such a wide range of MFPs. The distribution therefore must obey some type of smoothness restriction [4]. Now, the problem becomes a minimization problem, where the solution lies in the best balance between the smoothness of the reconstructed function (solution norm) and its proximity to the experimental data (residual norm). This balance is controlled by the choice of the regularization parameter. The role of this parameter has been widely overlooked in the literature, where the choice of its value has been poorly justified and, on some occasions, it remains unmentioned. In this work, we present a method to obtain the optimal value of this parameter using the L-curve criterion [5,6]. We apply this criterion to the thermal conductivity and the phonon-MFP distribution and we present a study of the impact that the choice of its value has in the reconstructed accumulation. We demonstrate that this methodology can be extended to several transport properties involving carriers of different physical nature and behavior. 2. Materials and Methods To perform the reconstruction of the MFP distribution of the energy carriers, the only input needed is a characteristic SF and a well distributed set of experimental data, i.e., a large amount the experimental measurements spread over the different characteristic sizes of the system. The suppression function strongly depends on the specific geometry of the sample and the experimental configuration. It relates the transport property of the nanostructure αnano and that of the bulk αbulk . The suppression has been derived for different experimental geometries from the Boltzmann transport equation [7–9]. The relation between the transport coefficient αnano (d) and the suppression function was originally derived for thermal conductivity [4,10], and more recently has been used to determine the MFP of magnons and the spin diffusion length distribution [11]. This relation can be expressed by means of a cumulative MFP distribution as Z ∞ αnano (d) dS(χ) dχ α= =− Facc (Λbulk ) dΛbulk (1) αbulk 0 dχ dΛbulk where S(χ) is the suppression function, χ is the Knudsen number χ = Λbulk /d, Λbulk is the bulk MFP and d the characteristic length of the sample and Facc is the accumulation function given by Z Λc 1 Facc = αΛ (Λbulk )dΛbulk (2) αbulk 0 Nanomaterials 2019, 9, 414 3 of 10 where αΛ is the contribution to the total transport property of carriers with a mean free path Λ. This function represents the contribution of carriers with MFPs up to an upper limit, Λc , to the total transport property, and is the object that will be recovered by applying the MFP reconstruction technique. From this definition it is easy to see that the accumulated function is subject to some physical restrictions: it cannot take values lower than zero for Λc = 0 or higher than one for Λc → ∞, and it must be monotonous [4]. We can recognize that Equation (2) is a Fredholm integral equation of the dS(χ) first kind that transforms the accumulation function Facc (Λbulk ) into α with K = dχ dΛdχ acting as a bulk kernel. As the inverse problem of reconstructing the accumulation function Facc is an ill-posed problem with infinite solutions [10]. Minnich demonstrated some restriction can be imposed on Facc to obtain a unique solution [4]. Furthermore, it is reasonable to require the smoothness conditions mentioned before on Facc , since it is unlikely to have abrupt behavior in all its domain. These requirements can be applied through the Tikhonov regularization method, where the criterion to obtain the best solution Facc is to find the following minimum: min{k A· Facc − αi k22 +µ2 k ∆2 Facc k22 } (3) where αi is the normalized i-th measurement, A = K (χi,j ) × β i,j is an m × n matrix, where m is the number of measurements and n is the number of discretization points, Ki,j is the value of the kernel at χi,j = Λi,j /di , and β i,j the weight of this point for the quadrature. The operators k k2 and ∆2 are the 2-norm and the (n − 2) × n trigonal Toeplitz matrix which represent a second order derivative operator, respectively. The first term of Equation (3) is related to how well our result fits to experimental data (residual norm), while the second term represents the smoothness of the accumulation function (solution norm). The balance between both is controlled by the regularization parameter µ. In other words, µ sets the equilibrium between how good the experimental data is fitted and how smooth is the fitting function. The choice of µ will thus have a huge impact on the final result of the accumulation function, and a criterion to obtain the optimal value must be established. The selection of the most adequate µ is still an open question in mathematics. Several heuristic methods are frequently used, such as the Morozov’s discrepancy principle, the Quasi-Optimally criterion, the generalized cross validation, L-curve criterion, the Gfrerer/Raun method, to name a few [5,6]. Among these methods, the L-curve criterion is one of the most popular methods due to its robustness, convergence speed and efficiency. This method establishes a balance between the size of the discrepancy of the fitting function and the experimental data (residual norm) with the size of the regularized solution (solution norm) for different values of µ. The curve has an L-like shape composed by a steep part where the solutions are dominated by perturbation errors and a flat part where the solution is dominated by regularization errors. The corner represents a compromise between a good fit of the experimental data and the smoothness of the solution. It has been found that the corresponding point lies in the corner of the L-curve in the residual norm-solution norm plane, which can be defined as the point of maximum curvature. The optimal value of µ can be found by locating the peak position of the curvature as a function of µ [5,6]. The method employed is depicted in Figure 1, and is common to all cases presented here. The experimental data in the case studies shown in this work was reproduced from the digitalization of the images and the corresponding uncertainties. A specific suppression function is selected for each case, depending on the particular geometry and experimental configuration. As explained above, the integral in Equation (2) is discretized, and an adequate integration interval depending on the span of experimental data is chosen. At this point, instead of introducing an arbitrary value of the regularization parameter, we determine the optimal value via the L-curve criterion. We have observed that the distribution of the curvature depending on log µ follows a Gaussian distribution, thus allowing us to obtain the peak, i.e., the highest curvature (corresponding to the corner of the L-curve), using a Gaussian fitting with a reduced number of computational points. With this optimal value of the Nanomaterials 2019, 9, 414 4 of 10 regularization parameter we can proceed to apply the Tikhonov regularization method and impose the conditions Nanomaterials 9, xFFOR 2019,on acc using a convex optimization package for MATLAB called CVX [12,13]. 4 of 10 PEER REVIEW Figure Figure 1. Flow-chart of 1. Flow-chart of the the reconstruction method used reconstruction method used to to obtain obtain the the accumulation accumulation function function from from experimental experimental data. data. 3. Results and Discussions 3. Results and Discussions The MFP reconstruction method, as well as the method here presented for the selection of µ, does The MFP reconstruction method, as well as the method here presented for the selection of μ, not require a priory any physical assumption about the carrier, such as band structure or velocity. does not require a priory any physical assumption about the carrier, such as band structure or In this section, we apply the method to heat transport by phonons in three examples cases where the velocity. In this section, we apply the method to heat transport by phonons in three examples cases characteristic size is given by (a) the thickness of the nanostructure in a layered system, (b) a length where the characteristic size is given by (a) the thickness of the nanostructure in a layered system, (b) scale of the measurement technique and (c) a combination of the former two. The extension of the a length scale of the measurement technique and (c) a combination of the former two. The extension method to magnon-mediated transport phenomena, namely, the spin Seebeck effect and spin-Hall of the method to magnon-mediated transport phenomena, namely, the spin Seebeck effect and spin- torque coefficient, is further examined in the supporting information. Hall torque coefficient, is further examined in the supporting information. 3.1. Phonons in Out-of-Plane Thermal Transport in Graphite from Molecular Dynamic Simulations 3.1. Phonons in Out-of-Plane Thermal Transport in Graphite from Molecular Dynamic Simulations Firstly, we will study the MFP reconstruction of phonons in cross-plane thermal transport along Firstly, the c-axis we will study of graphite. the MFP The thermal reconstruction conductivities (κ)of phonons were in cross-plane obtained from a set ofthermal transport simulations along performed the c-axis of graphite. The thermal conductivities (κ) were obtained by Wei et al. [14] From this numerical experiment, we have recovered the MFP distribution in bulk from a set of simulations performed graphite by the along Weic-axis et al. by [14] usingFrom this numerical the L-curve criterion basedexperiment, we have recovered MFP reconstruction method. the ThisMFP will distribution in bulk graphite along the c-axis by using the L-curve criterion illustrate the different steps and details which apply to all of the cases that we present here. based MFP reconstruction method. This will As shown illustrate in Figure the obtaining 1 after different steps and details which the “experimental data”apply to all the only of the input casesisthat needed the we SF. present here. As the thermal transport occurs in the c-direction, i.e., in the cross-plane direction, we will employ the As shown cross-plane in Figure 1 afterSF, Fuchs-Sondheimer obtaining described theby“experimental [15] data” the only input needed is the SF. As the thermal transport occurs in the c-direction, i.e., in the cross-plane direction, we will employ  Z 1  the cross-plane Fuchs-Sondheimer SF, described 1 by [15] 3 −1/(χy) S(χ) = 1 − 3χ − y e dy (4) 4 0 1 𝑆(𝜒) = 1 − 3𝜒 − 𝑦 e ⁄( ) d𝑦 (4) Equation (4) describes the modification of 4 the phonon transport along the perpendicular direction of a thin film due Equation (4) to the finitethe describes sizemodification effect. Once of thethe SF is fixed, the phonon next step transport is to the along findperpendicular the optimum direction of aby µ-parameter minimizing thin film due Equation to the finite (3) size for different effect. Once the SFThe µ-values. dependence is fixed, the nextofstep bothisthe to residual find the and the solution norms on µ can be visualized as a 3-D curve (see Figure optimum μ-parameter by minimizing Equation (3) for different μ-values. The dependence of both the 2a), and its projection on xz-plane is the L-curve, named after its characteristic shape. The residual and the solution norms on μ can be visualized as a 3-D curve (see Figure 2a), and itscurve was generated using the simulated projection on cross-plane xz-plane is of graphite µ the L-curve, as function named afterof itsthickness at T shape. characteristic = 300 K The [14]. curveFigure 2b shows was generated the curvature distribution of k A · using the simulated cross-plane μ ofaccgraphiteF − α k vs. k ∆ i 2 as functionacc 2F k 2 expressed of thickness in heat-like colour at T = 300 K [14]. Figure map2b representing the maximum shows the curvature (yellow) distribution of ‖𝐴 and· 𝐹minimum − 𝛼 ‖ (black) vs. ‖Δ 𝐹 values of the curvature. ‖ expressed We can in heat-like observe colour map that the maximum representing of the curvature the maximum (yellow) is andfound right at minimum the corner (black) valuesofofthe theL-curve curvature.for µWe ≈ can 0.95,observe which is extracted from the Gaussian distribution of the curvature. In this corner that the maximum of the curvature is found right at the corner of the L-curve for μ ≈ 0.95, which is the balance between a extracted from the Gaussian distribution of the curvature. In this corner the balance between a smooth function and a good fitting of the experimental data can be found. Figure 3a shows the Facc as a function of the MFP for different μ-values. For this particular example, we can observe that depending on the μ-values, the span of MFP can fluctuate from a very wide (5 nm–6.88 μm for μ = 10) to very narrow (60–100 nm for μ = 0.1) distribution. It is also important to notice that μ-values Nanomaterials 2019, 9, 414 5 of 10 smooth function and a good fitting of the experimental data can be found. Figure 3a shows the Facc as a function of the MFP for different µ-values. For this particular example, we can observe that depending Nanomaterials 2019, 9, x FOR PEER REVIEW 5 of 10 on the µ-values, Nanomaterials 2019, 9,the span x FOR ofREVIEW PEER MFP can fluctuate from a very wide (5 nm–6.88 µm for µ = 10) to5very of 10 narrow (60–100 nm for µ = 0.1) distribution. It is also important to notice that µ-values between 0.1–4.0 between 0.1–4.0 lead to quite good agreement between the simulated κ (“experimental data”) and lead betweento quite good 0.1–4.0 leadagreement to quitebetween the simulated good agreement κ (“experimental between the simulateddata”) and A·Facc product κ (“experimental data”) (see and A·F acc product (see Figure 3b), although these values yield to a completely different span of the MFP. Figure A·F 3b), although acc product (see these 3b), Figure values yield tothese although a completely values different yield to a span of the MFP. completely The span different MFP distribution of the MFP. The MFP distribution varies from very narrow (low μ-values) to very wide (high μ-values). This varies The MFP from very narrow distribution (lowfrom varies to very (low wideμ-values) (high µ-values). This widedifference in the MFP isvery narrow to verygiven (high eachμ-values). This µ-values) difference in the MFP distribution a direct consequence of the weight to component in distribution isthe a direct consequence of athe weight given to each component in thetoEquation (3). The low the Equation (3). The low μ-values give a large weight to the residual norm and, as a consequence, in difference in MFP distribution is direct consequence of the weight given each component it the Equationgive(3). a large The weight low to the give residual norm and, to as a consequence, it leads toaan overfitting of of μ-values a large weight At the the residual norm and,μ-values as consequence, it µ-values leads to an overfitting the reconstructed function. same time, large give higher the reconstructed leads to an to function. overfitting of theAt the same reconstructed time, large function. give At the same µ-values higher time, importance large μ-values to the solution give higher importance the solution norm, leading to and over smoothing of the reconstructed function. norm, leading importance to and to the over smoothing solution norm, leading of the to reconstructed function. and over smoothing of the reconstructed function. -2 x 10 x 10 -2 -2 5x 10 x 10 5 -2 (a) Δ ·F||acc||2 5 Δ ·F||acc||2 (a) 5 acc 2 acc 2 4 4 4 2 ||Δ ||·F 4 2 ||Δ||·F 2 Max 3 23 norm 278.0 norm [a.u.] Max56.92 3 3 norm 278.0 norm [a.u.] 11.65 2 Curvature 56.92 2 Solution 2.386 11.65 2 Curvature Solution 2 Solution 0.4884 2.386 1 z Solution 0.1000 -1 1 Min 0.4884 1 z x x 10 0.1000 -1 1 Min yx x 10 || 1.5 0 1.2 1.5 -α i| 2 0 2 y 0.9 1.2 ·F acc-α |i 2 0 2 4 6 0.6 0.9 ||Α ·F acc 0 4 0.3 0.6norm ||Α -1 μ 6 8 0.4 0.6 0.8 1.0 1.2 1.4 x 10 8 10 0.0 0.3dual norm -1 μ 0.4 Residual 0.6 norm ||A·F 0.8 1.0 acc-1.2 αi||2 1.4(b)x 10 10 0.0Resi dual Residual norm ||A·Facc-αi||2 (b) si Re Figure 2. (a) Three-dimensional visualization of the relation between the L-Curve (blue) and the Figure Figure 2. (a) Three-dimensional Three-dimensional visualization 2. values (a) visualization of of the the relation relation between between the L-Curve L-Curve (blue) and the different of μ for a cross-plane Fuchs-Sondheimer suppressionthe function at (blue) T = 300and K the for different different values valuesofofµ μ forfor a cross-plane Fuchs-Sondheimer a cross-plane Fuchs-Sondheimersuppression function suppression at T = 300 function at TKis =for graphite 300 K for graphite simulations. (b) L-curve or xz projection of 3D curve, the maximum curvature displayed simulations. graphite (b) L-curve or xz projection of 3D curve, the maximum curvature is displayed as heat-like simulations. as heat-like color bar. (b) L-curve or xz projection of 3D curve, the maximum curvature is displayed color bar. as heat-like color bar. 1.0 Zhang et al. 100 nm MD data Wei et al. 1.0 Zhang Wei al.al. 100 nm et et 0.60 μMD opt =data 0.95Wei et al. 0.60 function 0.8 Wei et al. μopt = 0.95 6.88 μm μμopt = = 4.0 0.95 function κbulk μ = 0.95 6.88 μm μ = 4.0 0.8 μ opt = 4.0 0.45 κ/κ/bulk 0.6 μ = 4.0 0.45 μ 0.6 Accumulation μ μ 10 oror μ 10 0.30 Accumulation 10 8.2 0.4 acc acc 10 8.2 0.30 8.2 0.4 A·F 8.2 5.5 A·F 5.5 5.5 0.2 5.5 0.15 2.8 0.2 5 nm 2.8 0.15 2.8 5 nm 2.8 0.10 0.0 60 nm 0.10 0.0 -1 60 nm 0.00 0.10 0.10 10 -1 10 0 0 10 1 1 10 2 10 3 2 3 4 10 4 0.00 0 30 60 90 120 150 180 10 10 10MFP10 [nm] 10 10 (a) 0 30 60 90 [nm] Thickness 120 150 180 (b) MFP [nm] (a) Thickness [nm] (b) Figure 3. Figure (a) Phonon 3. (a) Phonon meanmean free free path path distribution distribution of of the the accumulation accumulation function function reconstructed reconstructed for for Figure 3. different (a) values Phonon of µ mean (0.1–10) free using path the distribution FS cross-plane of the accumulation suppression function function and different values of μ (0.1–10) using the FS cross-plane suppression function and the simulated κ made the reconstructed simulated κ madefor different by Wei et values al. [14].of μ The (0.1–10) using black-dashed the and FS cross-plane black-dotted suppression lines represent function and the reconstructions by Wei et al. [14]. The black-dashed and black-dotted lines represent reconstructions obtained by simulated obtained κ by made using by optimum an Wei using anetoptimum al.(0.95) [14]. and The (0.95)black-dashed large and(4.0) andμ-value, µ-value, large (4.0) black-dotted respectively. lines The respectively. represent cyan solid The reconstructions and cyan grey and solid dotted obtained line grey by represent dotted line using the MFPan reconstruction represent optimum the MFP (0.95) and largeby obtained reconstruction (4.0) Zhang obtained μ-value, et byal. respectively. [16] and Zhang Wei et al. The cyan et al. [16] and[14]solid Wei and (b)etThermalgreyconductivity al. [14] dotted (b) line Thermal represent normalized the to MFP the reconstruction bulk value and/or A · obtained Facc by Zhang product as et al. function [16] of and Wei graphite et al. thickness. conductivity normalized to the bulk value and/or A·Facc product as function of graphite thickness. The [14] The(b) Thermal open dots conductivity represent open dotsthe normalized simulated represent tosimulated the bulk thethermal value conductivity thermaland/or A·Facc product (“experimental conductivity as function data”). of data”). graphite Heat-like (“experimental thickness.lines linesHeat-like correspond Theto open dots different A · represent F for the µ-values simulated in a range thermal of 0.1 < µconductivity < 10. correspond toaccdifferent A·Facc for μ-values in a range of 0.1 < μ < 10. (“experimental data”). Heat-like lines correspond to different A·Facc for μ-values in a range of 0.1 < μ < 10. Figure 3a also shows the MFP reconstruction generated by Zhang et al. [16] (cyan solid line) and Figure by Wei 3a(grey et al. also shows dottedthe MFP line). reconstruction The generated first reconstruction bycarried was Zhang et al.using out [16] (cyan solid line) κand the measured of graphite as function of the thickness and the suppression function given by Equation (4). We canκsee by Wei et al. (grey dotted line). The first reconstruction was carried out using the measured of graphite that as function for the of the of optimal value thickness μ = 0.95,and thethe suppression MFP function of the carriers given spans by Equation 36 nm < Λ < 270 (4). nm,Wein can see similar Nanomaterials 2019, 9, 414 6 of 10 Figure 3a also shows the MFP reconstruction generated by Zhang et al. [16] (cyan solid line) and by Wei et al. (grey dotted line). The first reconstruction was carried out using the measured κ of graphite as function of the thickness and the suppression function given by Equation (4). We can see that for the optimal value of µ = 0.95, the MFP of the carriers spans 36 nm < Λ < 270 nm, in similar range compared with the reconstructed values obtained by Zhang et al. 40 nm < Λ < 210 nm [16]. The second reconstruction was calculated using the simulated κ and a quasi-ballistic model to describe the suppression function (Equations (3) and (4) of ref. [14]). If the drastic change of the MFP-distribution is due to the reconstruction method being an ill-posed problem, then, a drastic change of SF can have a strong impact on its distribution. The desired Facc will be strongly affected by the “shape” of the selected SF. Therefore the choice of a different SF will lead to a different MFP distribution [11]. 3.2. In-Plane Thermal Transport in 400 nm Thick Si Membrane The second case of study is the in-plane thermal transport in a 400 nm Si membrane at room temperature probed by the Thermal Transient Grating (TTG) technique, obtained from Johnson et al. [17]. In the previous example, the characteristic length of the system was the thickness of the graphite sheets. Here the Si membrane has fixed thickness d = 400 nm and the variable length scale is the period of the thermal grating L. Since the measurement of in-plane thermal conductivity in the membrane will correspond to phonons with MFPs lower than L in each measurement, the grating period becomes the equivalent of the maximum MFP. In this case we have to consider two different effects: (i) the impact of the finite size of the membrane and the boundary scattering of the effective in-plane MFP and (ii) the crossover from non-diffusive to diffusive phonon transport given by the period of the thermal grating. Then, a combination of two suppression functions have to be introduced that takes both effects into account [8]. In the first place, we define an effective MFP (Λ’) to take into account the effect of the boundary scattering due to thickness of the membrane as follow: Λ0 = ΛS2 (d/Λ) (5) where Λ is the bulk MFP, d is the thickness of the membrane and S2 is Fuchs-Sondheimer SF for in-plane thermal transport given by [10] Z ∞ 3d 3d S2 (d/Λ) = 1 − + (y−3 − y−5 )e−yΛ/d dy (6) 8Λ 2Λ 1 With this effective MFP we proceed to perform the reconstruction using the suppression function for the specific geometry of the experiment given by [17] 3 arctan(qΛ0 ) S1 (qΛ0 ) = ( 1 − ) (7) q2 Λ 02 qΛ0 where q = 2Λ/L and L is the grating wavevector. If we define ζ = qΛ’, the kernel for the reconstruction is given by Z ∞ 3arctan(ζ ) 3 + 2ζ 2    dS(ζ ) dζ dΛ 2π 3 3 −3 −5 −tΛ 0 0 = − 1+ dt(t −t ) e d (1 − t ) (8) dζ dΛ dΛ L ζ3 ζ 1 + ζ2 2 1 The S1 function is unity in the limit qΛ’ << 1, in the diffusive limit, and goes like (qΛ’)−2 for qΛ’ >> 1, in the ballistic regime, thus describing the transition between both regimes necessary to adequately interpret the measured quantities in the experiment [17]. Similarly to the previous case, we apply a Gaussian fitting procedure to obtain a value µopt = 1.05 at room temperature. In Figure 4 we can observe the effect of the different values of µ. Note that, in this Nanomaterials 2019, 9, 414 7 of 10 case, the reduction of µ affects mainly the smoothness of the reconstructed function with large changes on the concavity and convexity of the accumulation function. However, there was not a significant change on the A·Facc product for small µ-values. The oscillations observed for low µ-values are a direct consequence of the overfitting of the reconstructed function, due to the large weight imposed to the minimization problem. The increase of the regularization parameter beyond the optimal value results in an increase of the span of the MFP of the carriers, as shown in Figure 4a, and a poor agreement with the experimental Nanomaterials 2019, 9, data, x FOR as canREVIEW PEER be seen in Figure 4b. 7 of 10 MFP distribution of 400 nm Si membrane 0.6 1.0 0.60 Accumulation function 0.8 0.48 κ/κbulk 0.6 0.36 κ/κbulk 0.5 0.4 0.24 μlow = 0.1 μlow = 0.1 μopt = 1.05 μopt = 1.05 0.2 0.12 μhigh = 4.0 μhigh = 4.0 Johnson 0.0 0.00 0.4 0.1 1 2 2 4 6 8 10 12 14 MFP [μm] (a) Transient grating period [μm] (b) Figure 4. (a) Figure 4. (a) Phonon Phonon meanmean free free path path distribution distribution of of the the accumulation accumulation function function for for aa 400 400 nm nm Si Si film film reconstructed for µ reconstructed for μopt (green line). The blue and red lines are the result obtained using a low and opt (green line). The blue and red lines are the result obtained using a low and high high value value of μ, respectively. of µ, respectively. (b)(b) Thermal Thermal conductivity conductivity of of 400 400 nm nm Si Si film film corresponding corresponding toto the the different different accumulation accumulation functions (red, green, and blue lines) for different transient grating periods in functions (red, green, and blue lines) for different transient grating periods in the the experimental technique experimental technique [17]. [17]. 3.3. In-Plane Thermal Transport in Si: Reconstruction by Changing the Thickness of the Membrane 3.3. In-Plane Thermal Transport in Si: Reconstruction by Changing the Thickness of the Membrane The dependence of the reconstructed accumulation function on the regularization parameter The dependence of the reconstructed accumulation function on the regularization parameter relies on the L-curve and the dependence of the residual norm k A· Facc − αi k2 and solution norm relies on the L-curve and the dependence of the residual norm ‖𝐴 · 𝐹 − 𝛼 ‖ and solution norm k ∆2 Facc k2 on µ. The following case illustrates an example of how the particular shapes of these ‖Δ 𝐹 ‖ on μ. The following case illustrates an example of how the particular shapes of these different curves can strongly affect the reconstruction. different curves can strongly affect the reconstruction. For this example, we used the thickness dependence of the in-plane κ measured by Cuffe et al. [10]. For this example, we used the thickness dependence of the in-plane κ measured by Cuffe et al. The experiment consisted in the measurement of the κ of free standing Si membranes ranging from [10]. The experiment consisted in the measurement of the κ of free standing Si membranes ranging 15–1518 nm using the transient thermal grating technique in the diffusive regime. As the thermal from 15–1518 nm using the transient thermal grating technique in the diffusive regime. As the transport occurs along the films, the Fuchs-Sondheimer for in-plane thermal transport (see Equation (6)) thermal transport occurs along the films, the Fuchs-Sondheimer for in-plane thermal transport (see was selected as SF. Equation (6)) was selected as SF. As we can see in Figure 5, the three-dimensional (3D) visualization of the minimization problem is As we can see in Figure 5, the three-dimensional (3D) visualization of the minimization problem manifestly different to that presented for graphite in Figure 2a. For graphite, we observed very smooth is manifestly different to that presented for graphite in Figure 2a. For graphite, we observed very behavior of each of the projections (xy, xz and yz planes) in the 3D curve. This leads a continuous smooth behavior of each of the projections (xy, xz and yz planes) in the 3D curve. This leads a variation of MFP from a narrow to a wide distribution as µ increases. On the other hand, Figure 5 continuous variation of MFP from a narrow to a wide distribution as μ increases. On the other hand, shows singularities and abrupt jumps of the different projections. For the case of low µ-values, this Figure 5 shows singularities and abrupt jumps of the different projections. For the case of low μ- leads to large oscillations in the MFP distribution function, as is shown in the blue solid line in Figure 6a. values, this leads to large oscillations in the MFP distribution function, as is shown in the blue solid However for large µ-values, we observe a linear distribution of the reconstructed curve as a function line in Figure 6a. However for large μ-values, we observe a linear distribution of the reconstructed of the logarithm of MFP (see red solid line in Figure 6a). For intermediate µ-values (0.8–3), we can curve as a function of the logarithm of MFP (see red solid line in Figure 6a). For intermediate μ-values observe no major differences either in the smoothness or in the MFP span of the accumulated function. (0.8–3), we can observe no major differences either in the smoothness or in the MFP span of the Similarly, we did not observe huge changes in the fitted function as it is displayed in Figure 6b. accumulated function. Similarly, we did not observe huge changes in the fitted function as it is Figure 6a also shows the MFP reconstruction generated by Cuffe et al. [10]. This reconstruction displayed in Figure 6b. was carried out using the same procedure than that shown in this section, including the same Figure 6a also shows the MFP reconstruction generated by Cuffe et al. [10]. This reconstruction suppression function and the minimization procedure. However, no information is given regarding was carried out using the same procedure than that shown in this section, including the same the regularization parameter or the procedure to estimate it. Despite this, we can observe that the MFP suppression function and the minimization procedure. However, no information is given regarding distribution estimated for the optimal value of µ = 0.8 and the one obtained by Cuffe et al. spans in the regularization parameter or the procedure to estimate it. Despite this, we can observe that the MFP distribution estimated for the optimal value of μ = 0.8 and the one obtained by Cuffe et al. spans in similar range with MFP ranging from 20 nm < Λ < 100 μm. Slight deviations can be observed for MFP around Λ = 1 μm, which are expected from small differences in the regularization parameter and the numerical calculation of the minimization problem. For comparison, the calculated and Nanomaterials 2019, 9, 414 8 of 10 similar range with MFP ranging from 20 nm < Λ < 100 µm. Slight deviations can be observed for MFP around Λ = 1 µm, which are expected from small differences in the regularization parameter and the numerical calculation of the minimization problem. For comparison, the calculated and measured MFP distribution of bulk Si are also included. The good agreement between our reconstruction distribution and the molecular dynamics simulations by Henry et al. [18] and the measurements performed by Regner et al. [19] is remarkable. In the latter case, it is important to keep in mind that the experiment is Nanomaterials 2019, 9, x FOR PEER REVIEW 8 of 10 completely different Nanomaterials 2019, from 9, x FOR the PEER one presented here. REVIEW 8 of 10 -2 xx 10 10-2 norm||Δ||Δ·F·Faccacc||2||2 33 2 2 22 Solutionnorm Solution 11 -1 xx 10 10-1 || 22 z α|| 33 i i z x x - α- y c cc 22 ·F·F y ac a 00 11 ||A||A 22 3 11 μ3 44 mrm onro μ 00 55 l nl udaua seisdi RRe Figure 5. Three-dimensional visualization of the relation between the L-Curve (blue) and the different Figure 5. 5.Three-dimensional Three-dimensionalvisualization of the visualization of relation between the relation the L-Curve between (blue) and the L-Curve the different (blue) and the values of μ for an in-plane Fuchs-Sondheimer suppression function at T = 300 K for silicon values ofvalues different μ for of an µin-plane Fuchs-Sondheimer for an in-plane suppression Fuchs-Sondheimer function function suppression at T = 300 at TK= for 300 silicon K for measurements. measurements silicon . measurements. 0.8 0.8 1.0 Henry (MD simulations) 1.0 Henry (MD simulations) Cuffe (reconstruction) 0.7 0.7 function Cuffe (reconstruction) Accumulationfunction 0.8 Exp. Regner 0.8 Exp. Regner 0.6 0.6 0.6 0.5 0.6 0.5 κκ/κ/κbulkbulk Accumulation 0.4 μlow = 0.1 0.4 0.4 μlow = 0.1 0.4 μlow = 0.1 μopt = 0.8 μlow = 0.1 0.3 μopt μopt = 0.8 0.3 = 0.8 μhigh = 3.0 0.2 μopt = 0.8 μhigh 0.2 μhigh = 3.0 0.2 = 3.0 μhighest = 5.0 μhigh = 3.0 0.2 μhighest = 5.0 0.0 μhighest = 5.0 Exp. Cuffe 0.0 μhighest = 5.0 0.1 Exp. Cuffe 0.1 100 0 101 1 102 2 103 3 104 4 105 5 10 100 1000 10 10 10MFP [nm] 10 10 10(a) 10 100 1000 Thickness [nm] (b) MFP [nm] (a) Thickness [nm] (b) Figure 6.6. (a) Figure (a) Phonon Phonon mean mean free free path path distribution distribution of of bulk bulk Silicon Silicon reconstructed reconstructed from from thethe thickness thickness Figure 6. (a) Phonon mean free path distribution of bulk Silicon reconstructed from the thickness dependence dependence of of thermal thermal conductivity. conductivity. The The dark dark green, green, red, red, dotted dotted grey grey and and dotted dotted dark dark cyan cyan lines lines dependence of thermal conductivity. The dark green, red, dotted grey and dotted dark cyan lines represent reconstructions for different µ-values. The black stared and blue dotted represent reconstructions for different μ-values. The black stared and blue dotted lines represent the lines represent represent reconstructions for different μ-values. The black stared and blue dotted lines represent the the MFP MFP reconstruction reconstruction and molecular and molecular dynamic dynamic calculation calculation performedperformed by CuffebyetCuffe et and al. [10] al. [10] Henryand et MFP reconstruction Henry et al. [18], and molecular respectively. dynamic The calculation orange-solid dots performed represent by theCuffe et al. [10] phonon-MFP experimental and Henry et al. [18], respectively. The orange-solid dots represent the experimental phonon-MFP distribution of al. [18], respectively. distribution The orange-solid bydots represent the experimental phonon-MFP distribution of bulk siliconof bulk silicon measured measured by Regner et al. Regner [19]. (b) et al. [19]. Thermal (b) Thermal conductivity conductivity corresponding corresponding to the differentto bulk the silicon different measured accumulationby Regner functionset al. for [19]. (b) different Thermal samples conductivity with corresponding different thickness to [10]. the different accumulation functions for different samples with different thickness [10]. accumulation functions for different samples with different thickness [10]. The L-curve criterion is efficient for obtaining the most adequate regularization parameter, but The L-curve criterion is efficient for obtaining the most adequate regularization parameter, but in this case the reconstructed accumulation function is very robust against changes in μ due to the in this case the reconstructed accumulation function is very robust against changes in μ due to the particular flatness of the residual and solution norm depending on the values of μ. It is also important particular flatness of the residual and solution norm depending on the values of μ. It is also important to remark that the estimation of the optimum μ-value has to be carried out for each set of to remark that the estimation of the optimum μ-value has to be carried out for each set of measurements. As we showed in the supporting information, for the case of temperature dependence measurements. As we showed in the supporting information, for the case of temperature dependence Nanomaterials 2019, 9, 414 9 of 10 The L-curve criterion is efficient for obtaining the most adequate regularization parameter, but in this case the reconstructed accumulation function is very robust against changes in µ due to the particular flatness of the residual and solution norm depending on the values of µ. It is also important to remark that the estimation of the optimum µ-value has to be carried out for each set of measurements. As we showed in the supporting information, for the case of temperature dependence measurements of the Spin Seebeck coefficient, the optimal µ-value varies among the measurements (See Figure S2). This effect comes from the differences in the spread of the experimental data in each temperature regime and the importance of carriers with longer MFPs as the temperature decreases [11]. 4. Conclusions We have demonstrated that the choice of the regularization parameter µ has a large impact on the physical information obtained from the reconstructed accumulation function, and thus cannot be chosen arbitrarily. On one hand, small µ-values lead to a very narrow MFP distribution but also very good agreement between the fitting function and the experimental data. However, they also introduce artefacts (oscillations) to the reconstructed MFP distributions. These oscillations are the consequence of an overfitting of the reconstructed function due to the low weight given to the solution norm. On the other hand, large µ-values lead to a wider MFP distribution but worse agreement between the fitting function and the experimental data. This is a direct consequence of the larger weight to the solution norm and, as consequence, an over smoothing of the fitted function. It is also important to notice that the estimation of the optimum µ-value has to be carried out for each set of measurements. It cannot be assumed the same value for temperature dependence measurements will occur, as we show in the supporting information. The method presented here is not limited only to phonon transport and it can be applied to many different cases. The only input needed to reconstruct the mean free path distribution of the carriers is a well-known suppression function and a well distributed set of experimental data points. Regardless of the carrier physical behavior or the span of MFP, we have proven the method to be applicable, and we have seen that the robustness of the reconstruction against deviations from the estimated value of µ will depend on the particular variation of the solution norm and reduced norm with µ. Supplementary Materials: The following are available online at http://www.mdpi.com/2079-4991/9/3/414/s1, Figure S1: (a) Magnon mean free path distribution of the accumulation function for spin Seebeck effect experiments. (b) Normalized Longitudinal Spin-Seebeck coefficient for different thickness of the YIG. Figure S2: Optimal values of µ for the Magnon-MFP reconstruction. Figure S3: (a) Spin diffusion length distribution of the accumulation function reconstructed for spin-Hall torque experiments. The blue and red lines are the result obtained using a low and high value of µ, respectively. The green line represent the MFP distribution reconstructed using µopt . (b) Normalized Longitudinal spin-Hall torque coefficient for different thickness of the Pt filmvalues of µ for the Magnon-MFP reconstruction at different temperatures for LSSE experiments. Author Contributions: M.-Á.S.-M., F.A., E.C.-A. wrote the article. M.-Á.S.-M., J.O. and E.C.-A. designed the different numerical programs. C.M.S.T. and F.A. supervised the work and discussed the results. All authors discussed the results and commented on the manuscript. Funding: This research was funded by the EU FP7 project QUANTIHEAT (Grant No 604668) and the Spanish MINECO-Feder project PHENTOM (FIS2015-70862-P). M.-Á.S.-M. acknowledges support from SO-FPI fellowship BES 2015-075920. ICN2 acknowledges support from Severo Ochoa Program (MINECO, Grant SEV-2017-0706) and funding from the CERCA Programme/Generalitat de Catalunya. Acknowledgments: We thank G. Whitworth for a critical reading of the manuscript. Conflicts of Interest: The authors declare no conflict of interest. References 1. Ashcroft, N.W.; Mermin, N.D. Solid State Physics; Garbose Crane, D., Ed.; Saunders College: Rochester, NY, USA, 1976; ISBN 0-03-083993-9. 2. Esfarjani, K.; Chen, G.; Stokes, H.T. Heat transport in silicon from first-principles calculations. Phys. Rev. B 2011, 84, 085204. [CrossRef] Nanomaterials 2019, 9, 414 10 of 10 3. Yang, F.; Dames, C. Mean free path spectra as a tool to understand thermal conductivity in bulk and nanostructures. Phys. Rev. B 2013, 87, 35437. [CrossRef] 4. Minnich, A.J. Determining phonon mean free paths from observations of quasiballistic thermal transport. Phys. Rev. Lett. 2012, 109, 205901. [CrossRef] [PubMed] 5. Hansen, P.C.; O’Leary, D.P. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM J. Sci. Comput. 1993, 14, 1487–1503. [CrossRef] 6. Hansen, P.C. The L-curve and its use in the numerical treatment of inverse problems. In Computational Inverse Problems in Electrocardiology; Johnston, P., Ed.; WIT press: Ashurst, UK, 2000; pp. 119–142. 7. Chiloyan, V.; Zeng, L.; Huberman, S.; Maznev, A.A.; Nelson, K.A.; Chen, G. Variational approach to extracting the phonon mean free path distribution from the spectral Boltzmann transport equation. Phys. Rev. B 2016, 93, 155201. [CrossRef] 8. Hua, C.; Minnich, A.J. Transport regimes in quasiballistic heat conduction. Phys. Rev. B 2014, 89, 94302. [CrossRef] 9. Zeng, L.; Collins, K.C.; Hu, Y.; Luckyanova, M.N.; Maznev, A.A.; Huberman, S.; Chiloyan, V.; Zhou, J.; Huang, X.; Nelson, K.A.; et al. Measuring Phonon Mean Free Path Distributions by Probing Quasiballistic Phonon Transport in Grating Nanostructures. Sci. Rep. 2015, 5, 17131. [CrossRef] [PubMed] 10. Cuffe, J.; Eliason, J.K.; Maznev, A.A.; Collins, K.C.; Johnson, J.A.; Shchepetov, A.; Prunnila, M.; Ahopelto, J.; Sotomayor Torres, C.M.; Chen, G.; et al. Reconstructing phonon mean-free-path contributions to thermal conductivity using nanoscale membranes. Phys. Rev. B 2015, 91, 245423. [CrossRef] 11. Chavez-Angel, E.; Zarate, R.A.; Fuentes, S.; Guo, E.J.; Kläui, M.; Jakob, G. Reconstruction of an effective magnon mean free path distribution from spin Seebeck measurements in thin films. New J. Phys. 2017, 19, 13011. [CrossRef] 12. Grant, M.; Boyd, S. Graph implementations for nonsmooth convex programs. In Recent Advances in Learning and Control; Blondel, V., Boyd, S., Kimura, H., Eds.; Lecture Notes in Control and Information Sciences; Springer-Verlag Limited: New York, NY, USA, 2008; pp. 95–110. 13. Grant, M.; Boyd, S. CVX: Matlab Software for Disciplined Convex Programming, Version 2.1. Available online: http://cvxr.com/cvx (accessed on 29 July 2017). 14. Wei, Z.; Yang, J.; Chen, W.; Bi, K.; Li, D.; Chen, Y. Phonon mean free path of graphite along the c-axis. Appl. Phys. Lett. 2014, 104, 81903. [CrossRef] 15. Hua, C.; Minnich, A.J. Semi-analytical solution to the frequency-dependent Boltzmann transport equation for cross-plane heat conduction in thin films. J. Appl. Phys. 2015, 117, 175306. [CrossRef] 16. Zhang, H.; Chen, X.; Jho, Y.-D.; Minnich, A.J. Temperature-Dependent Mean Free Path Spectra of Thermal Phonons Along the c-Axis of Graphite. Nano Lett. 2016, 16, 1643–1649. [CrossRef] [PubMed] 17. Johnson, J.A.; Maznev, A.A.; Cuffe, J.; Eliason, J.K.; Minnich, A.J.; Kehoe, T.; Sotomayor Torres, C.M.; Chen, G.; Nelson, K.A. Direct Measurement of Room-Temperature Nondiffusive Thermal Transport Over Micron Distances in a Silicon Membrane. Phys. Rev. Lett. 2013, 110, 25901. [CrossRef] [PubMed] 18. Henry, A.S.; Chen, G. Spectral Phonon Transport Properties of Silicon Based on Molecular Dynamics Simulations and Lattice Dynamics. J. Comput. Theor. Nanosci. 2008, 5, 141–152. [CrossRef] 19. Regner, K.T.; Sellan, D.P.; Su, Z.; Amon, C.H.; McGaughey, A.J.H.; Malen, J.A. Broadband phonon mean free path contributions to thermal conductivity measured using frequency domain thermoreflectance. Nat. Commun. 2013, 4, 1640. [CrossRef] [PubMed] © 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).