## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# A new regime of nanoscale thermal transport: Collective diffusion increases dissipation efficiency

Contributed by Henry C. Kapteyn, February 19, 2015 (sent for review February 3, 2015)

## Significance

A complete description of nanoscale thermal transport is a fundamental problem that has defied understanding for many decades. Here, we uncover a surprising new regime of nanoscale thermal transport where, counterintuitively, nanoscale heat sources cool more quickly when placed close together than when they are widely separated. This increased cooling efficiency is possible when the separation between nanoscale heat sources is comparable to the average mean free paths of the dominant heat-carrying phonons. This finding suggests new approaches for addressing the significant challenge of thermal management in nanosystems, with design implications for integrated circuits, thermoelectric devices, nanoparticle-mediated thermal therapies, and nanoenhanced photovoltaics for improving clean-energy technologies.

## Abstract

Understanding thermal transport from nanoscale heat sources is important for a fundamental description of energy flow in materials, as well as for many technological applications including thermal management in nanoelectronics and optoelectronics, thermoelectric devices, nanoenhanced photovoltaics, and nanoparticle-mediated thermal therapies. Thermal transport at the nanoscale is fundamentally different from that at the macroscale and is determined by the distribution of carrier mean free paths and energy dispersion in a material, the length scales of the heat sources, and the distance over which heat is transported. Past work has shown that Fourier’s law for heat conduction dramatically overpredicts the rate of heat dissipation from heat sources with dimensions smaller than the mean free path of the dominant heat-carrying phonons. In this work, we uncover a new regime of nanoscale thermal transport that dominates when the separation between nanoscale heat sources is small compared with the dominant phonon mean free paths. Surprisingly, the interaction of phonons originating from neighboring heat sources enables more efficient diffusive-like heat dissipation, even from nanoscale heat sources much smaller than the dominant phonon mean free paths. This finding suggests that thermal management in nanoscale systems including integrated circuits might not be as challenging as previously projected. Finally, we demonstrate a unique capability to extract differential conductivity as a function of phonon mean free path in materials, allowing the first (to our knowledge) experimental validation of predictions from the recently developed first-principles calculations.

- nanoscale thermal transport
- nondiffusive transport
- mean free path spectroscopy
- high harmonic generation
- ultrafast X-rays

Critical applications including thermoelectrics for energy harvesting, nanoparticle-mediated thermal therapy, nanoenhanced photovoltaics, and thermal management in integrated circuits require a comprehensive understanding of energy flow at the nanoscale. Recent work has shown that the rate of heat dissipation from a heat source is reduced significantly below that predicted by Fourier’s law for diffusive heat transfer when the characteristic dimension of the heat source is smaller than the mean free path (MFP) of the dominant heat carriers (phonons in dielectric and semiconductor materials) (1⇓⇓⇓⇓–6). However, a complete fundamental description of nanoscale thermal transport is still elusive, and current theoretical efforts are limited by a lack of experimental validation.

Diffusive heat transfer requires many collisions among heat carriers to establish a local thermal equilibrium and a continuous temperature gradient along which energy dissipates. However, when the dimension of a heat source is smaller than the phonon MFP, the diffusion equation is intrinsically invalid because phonons move ballistically without collisions. The rate of nanoscale heat dissipation is significantly lower than the diffusive prediction such that smaller heat sources appear increasingly inefficient in dissipating heat. Furthermore, heat-carrying phonons in real materials have a wide distribution of MFPs, from several nanometers to hundreds of microns (7). For a given nanoscale heat source size, phonons with MFPs shorter than the hot-spot dimension remain fully diffusive and contribute to efficient heat dissipation and a high thermal conductivity (or equivalently, a low thermal resistivity). In contrast, phonons with long MFPs travel ballistically far from the heat source before scattering, with an effective thermal resistivity far larger than the diffusive prediction. Phonons with intermediate MFPs fall in between; here, heat transfer is quasiballistic with varying degrees of reduced contributions to the conduction of heat away from the nanoscale source. Most work to date explored the reduction in heat transfer from functionally isolated microscale and nanoscale heat sources (1⇓⇓⇓–5). Indeed, characterizing heat transfer from micro/nanostructures with varying size can be used to experimentally measure cumulative phonon MFP spectra of materials (2, 8, 9), with the proof-of-principle demonstrated for long-MFP (>1-μm) phonons in silicon (3).

In this work, we show through both experiment and theory that the size of the heat source is not the only important scale that determines nanoscale heat dissipation. We identify a new regime of thermal transport that occurs when the separation between nanoscale heat sources is smaller than the average phonon MFP. Surprisingly, we find that when phonons from neighboring heat sources interact, more of them dissipate heat in a diffusive regime, thus counteracting the inefficiency caused by ballistic effects in the case of isolated heat sources. This collective behavior can increase heat transfer to near the diffusive limit. Most importantly, the appearance of this “collectively diffusive” regime mitigates the scaling problems for thermal management in nanoelectronics, which may not be as serious as projected (2, 10, 11). Finally, we use this previously unobserved phenomenon to extract the contribution to thermal transport from specific regions of the phonon MFP spectrum, opening up a new approach for thermal transport metrology and MFP spectroscopy. This is because by varying both nanostructure size and separation, an effective phonon filter is introduced that suppresses specific sections of MFP contributions to thermal conductivity. We compare our extracted phonon MFP spectra with predictions from first-principles calculations and find excellent agreement between experiment and theory. Looking forward, we have a unique capability for characterizing phonon transport in novel materials where predictions do not yet exist.

Fig. 1 illustrates the differences between the three regimes of heat transport from nanoscale heat sources—purely diffusive, quasiballistic, and collectively diffusive. Quasiballistic transport (Fig. 1*B*) dominates when the size of isolated nanoscale heat sources is smaller than dominant phonon MFPs. In the new collectively diffusive regime we uncovered (Fig. 1*C*), the separation between heat sources is small enough that long-MFP phonons, whose contribution to heat dissipation would normally be limited by the small size of nano-heat sources, can once again play a significant role and restore heat transfer efficiency to near the diffusive limit. Although these phonons travel ballistically away from each individual heat source, they can scatter with phonons originating from a neighboring heat source, thus creating an effectively larger heat source size. In the limiting case, the spacing between heat sources vanishes and this regime approaches heat dissipation from a uniformly heated layer.

## Experiment

In our experiment, arrays of nickel nanowires were fabricated by e-beam lithography and lift-off techniques to form periodic gratings on the surface of sapphire and silicon substrates. The nanowire line widths *L* range from 750 nm down to 30 nm, with period *P* = 4*L* and a rectangular profile height of ∼13.5 nm. The use of nanopatterned structures rather than optical absorption allows us to explore the dynamics of heat sources much smaller than the diffraction limit of visible light. The metallic nanowires are heated by a 25-fs pump pulse centered at a wavelength of 800 nm. The sapphire substrate is transparent at this wavelength, whereas the silicon substrate has such a long absorption depth that any small, uniform heating of the substrate can be neglected. Laser excitation thus creates an array of nanoscale hot spots (lines) on the surface of a cold substrate. Because all nanostructures are fabricated on the same substrate at the same time, the intrinsic thermal boundary resistivity at the interface between the metallic line and the substrate will be constant across all samples: any variation in efficiency of heat dissipation as the hot-spot size or spacing is varied can thus be attributed to different regimes of thermal transport.

The laser-induced thermal expansion and subsequent cooling of the nanogratings is probed using coherent extreme UV (EUV) light centered at a wavelength of 29 nm, created by high harmonic up-conversion of an 800-nm Ti:sapphire laser (12). The time delay between the EUV probe pulse and the laser pump pulse is adjusted using a mechanical delay stage between −400 and +8,000 ps, with step size as small as 1 ps. As the EUV light diffracts from the periodic array of Ni nanowires, expansion and cooling of the nanogratings changes the diffraction efficiency, and this signal is recorded by a CCD camera as a function of delay time between pump and probe pulses. Examples of this dynamic signal are shown in Fig. 2*A*. [Note that more details of the experimental setup are described elsewhere (13, 14) and in *SI Text*, section S1.] Because the reflectivities of these materials do not change with temperature at EUV wavelengths (15), the change in the diffraction signal can be uniquely attributed to physical deformations in the surface profile. Moreover, because the thermally expanded nanowires will change the diffraction signal as long as their temperature is higher than their initial state, this serves as a measure of the cooling rate of the nanowires rather than the rate of heat transfer across a specific distance as in optical transient grating (6) or buried-layer heating experiments (16). Thus, the data of Fig. 2*A* can be used to directly extract the average thermal expansion and relaxation of each individual nanowire induced by laser heating and subsequent heat dissipation into the substrate, in addition to the surface deformations caused by acoustic waves launched by the initial impulsive expansion (13, 14, 17).

## Theory

To understand the different regimes of thermal transport illustrated in Fig. 1, we consider three models: (*i*) the model described in our previous work that assumes isolated heat sources (1); (*ii*) an analytical model we develop here to account for interactions of phonons originating from neighboring heat sources using a gray, single-phonon–MFP approximation; and (*iii*) a more advanced interacting model that includes a distribution of phonon MFPs. As discussed in detail below, this interacting multi-MFP model allows us to extract MFP-dependent contributions to thermal conductivity for MFPs as short as 14 nm for the first time (to our knowledge), providing data that can be directly compared with predictions from first-principles density functional theory (DFT).

To quantify the deviations from diffusive heat transport, we first build upon methods similar to those described by Siemens et al. (1), but including more comprehensive finite-element physical modeling (18) to improve data reduction accuracy. We model our system using diffusive heat conduction theory, while allowing the effective thermal boundary resistivity (which sets the temperature discontinuity across the boundary between the nickel nanowires and the substrate) to vary as a function of line width to account for nondiffusive effects in the substrate near the heat source.

We use accurate sample dimensions (height, line width, and period) characterized by atomic force microscopy (*SI Text*, section S1). Fresnel optical propagation is then used to calculate the diffraction signal from the simulated surface deformations. The effective boundary resistivity, *r*_{eff}, that provides the best fit to the experimental data of Fig. 2*A* represents the sum of two terms: first, the constant intrinsic thermal boundary resistivity, *r*_{TBR}, that originates from the material difference between nickel and substrate; and second, a correction term, *r*_{Corr}, due to nondiffusive transport in the substrate close to the heat source when either *L* or *P* is smaller than MFPs. Although this is a similar effective correction for nondiffusive transport near isolated nanoscale heat sources as that used in previous works (2⇓⇓⇓–6), by assigning the nondiffusive contribution to the thermal boundary resistivity rather than to changes in the thermal conductivity of the substrate, we maintain a simple modeling geometry and avoid the need to assume a particular region of the substrate in which an effective conductivity change should apply. Further discussion can be found in *SI Text*, section S2.

The effective resistivity results are plotted in Fig. 2*B*. For large line widths on both sapphire and silicon substrates, the effective resistivity converges toward a constant value—the intrinsic thermal boundary resistivity. As the line width approaches the dominant phonon MFPs in the substrate, the effective resistivity rises as thermal transport becomes quasiballistic and the contribution to heat dissipation of long-MFP phonon modes is suppressed (1, 3, 8). This behavior was successfully described in past work using a simple gray model for sapphire and fused silica, which assumes a single-phonon MFP to loosely describe a weighted average of the MFPs from all of the phonon modes contributing to thermal transport in a given material. According to this model, a ballistic correction term proportional to Λ_{gray}/(*L*/2) can be added to the intrinsic thermal boundary resistivity (1, 19); this prediction is plotted in dotted red in Fig. 2*B*.

However, as the line width (and period) shrinks further, Fig. 2*B* shows that the effective resistivity starts to decrease rather than continuing to increase. The constant grating duty cycle for our series of samples means that the smallest-line width nanowires are also those with the smallest separation between neighboring heat sources. Thus, for small line widths, the separation becomes comparable to dominant phonon MFPs. For silicon, this peak in *r*_{eff} is shifted toward longer line widths/periods compared with sapphire because the phonon MFP distribution in silicon is also shifted toward longer MFPs, i.e., silicon has a longer average MFP than sapphire (1, 6). As illustrated in Fig. 1*C*, in this collectively diffusive regime, longer-MFP phonons from neighboring heat sources interact with each other as they would if they originated from a single, large heat source, leading to diffusive-like heat dissipation and decreasing the effective resistivity. The quasiballistic model for isolated heat sources clearly fails to capture this experimental observation, and a new model for *r*_{Corr} is required to account for the transition to this previously unpredicted and unobserved collectively diffusive regime.

We propose to use the concept of a notch filter in the MFP spectrum to describe the effects of grating line width and separation, shown schematically in Fig. 3*A*. The notch filter suppresses the contribution of phonon modes with MFPs that fall between the line width *L* and period *P* of the nanogratings. Thus, if the grating period (separation) remains large while the line width is decreased, one would expect the effective boundary resistivity to continue to rise, as shown in the red dotted lines of Fig. 2*B*. This is because the contributions of all phonon modes with MFPs longer than the line width *L* are suppressed in the quasiballistic regime of isolated heat sources. On the other hand, if the grating period shrinks, long-MFP phonon modes start to contribute again because phonons originating from neighboring heat sources interact with each other as they would in a bulk system, so the effective boundary resistivity should recover toward the bulk value, as seen experimentally in Fig. 2*B*. If this picture is true, then a nanograting with a small line width and large period (i.e., wide filter) should exhibit a slower thermal decay than the same line width grating with shorter period (i.e., narrower filter). Using a set of periodic nickel nanowires on silicon with varied duty cycle, we directly observed this pronounced difference as shown in Fig. 3*B*.

To build an analytical expression for *r*_{Corr} based on this idea, we use the concept of a phonon conductivity suppression function, *S*(*L*, *P*, Λ), similar to those described by others (8, 20). This suppression function is applied to a bulk differential conductivity spectrum versus phonon MFP of the substrate, *k*(Λ_{i}), to calculate an effective nanoscale conductivity *K*_{nano}:*r*_{Corr} to the change in conductivity represented by this suppression:*A* collects geometrical constants (discussed further in *SI Text*, section S3) and *K*_{bulk} is the bulk conductivity of the substrate, simply given by _{i}, *S* must approach unity in the diffusive regime when both *L* and *P* are large and at the limit of uniform heating when *L* = *P*. For the limit of small, isolated heat sources when *L* → 0 but *P* is large, *S* → *L* /(2 Λ) to reproduce the behavior of the previously published quasiballistic model (1, 21). Finally, the effects of *L* and *P* should be the same but opposite to each other so that *L* suppresses phonon mode contributions in the same way as *P* reactivates them. To capture these behaviors along with smooth transitions among regimes, the two effects are represented by a special case of the generic family of logistic functions, where the total suppression function is written as follows:*SI Text*, section S3. Although more rigorous methods of deriving suppression functions for various experimental geometries are currently being explored (9, 20, 22), none has yet sought to account for closely spaced heat sources or an accompanying reintroduction of phonon modes. Eq. **3** represents the first attempt (to our knowledge) to include the contribution of heat source spacing and offers a model that is simple enough for fast integration into existing models of heat transfer in nanoscale devices, for example, but complex enough to capture the previously unobserved behavior and make successful predictions.

To test this new model for *r*_{Corr}, we first assume the simple single-MFP (gray) model (where the MFP distribution is a delta function). The resulting predictions are shown in the blue curves in Fig. 2*B*. Specifically, *r*_{Corr} in this case is given by the following:*r*_{eff} data for sapphire, we extract values for *r*_{TBR} and Λ_{gray} that are consistent with previous results (1): Λ_{gray,} _{int} = 131 ± 11 nm, *r*_{TBR,} _{int} = 2.58 ± 0.19 × 10^{−9} m^{2}⋅K/W. This good agreement with the previous larger-line width data and the accurate fit for the full range of our data validate our interacting *r*_{Corr} model as an improved method to account for nanoscale size effects in heat transport—for both quasiballistic and collectively diffusive regimes. Interestingly, this single-MFP model provides a reasonable approximation for the entire range of heat transport in sapphire.

For silicon, the interacting *r*_{Corr} follows the general shape of the data and yields Λ_{gray,} _{int} = 306 ± 17 nm, which is consistent with previously reported values (6). However, the interacting gray-model approach, although more successful than the isolated model, fails to capture the additional structured tail in effective resistivity that appears for very small line widths and periods, below *L* = 100 nm. The failure of this approach is not surprising, because the single-MFP model is known to be a poor approximation for silicon with its broad distribution of phonon MFPs (8, 20).

## Discussion of Collectively Diffusive Regime and Extension to Phonon Mode Conductivity Spectra

Having developed and validated the new model to capture the transitions from diffusive, to quasiballistic, to collectively diffusive regimes, we can now extend our calculations beyond the simple single-MFP model and use our data to extract the MFP-dependent contributions to thermal conductivity in different materials down to MFPs as small as 14 nm for the first time (to our knowledge). Because line width and period set the location and width of the effective notch filter in the phonon MFP spectrum, each configuration uniquely samples the contribution to thermal conductivity of different MFP ranges of phonon modes with a resolution controlled primarily by the number of configurations tested. The larger the resistivity correction needed for a given nanograting, the stronger the conductivity contribution of phonon modes that were suppressed.

To extract information about the differential conductivity spectrum, we use the full multi-MFP form of *r*_{Corr}, given by the following:*SI Text*, section S3). Then by fitting our set of *r*_{eff} data as *r*_{eff} *= r*_{TBR} + *r*_{Corr}, we extract the average *k*(Λ_{i}) that is associated with all modes Λ_{i} within each bin, thus assessing the relative contributions to thermal conductivity of each region of the phonon MFP spectrum (plotted in Fig. 4). By limiting the number of bins to be no more than one-half of our number of data points, we ensure a conservative, well-conditioned fit, although we note that changing the bin number does not substantially alter the trends we observe. As shown by the purple curves in Fig. 2*B*, for sapphire, this procedure matches the experimentally measured *r*_{eff} as well as the gray model for interacting heat sources. For silicon, this more complete multiple-MFP interacting model is able to match our experimental measurements of thermal boundary resistivity over diffusive, quasiballistic and collectively diffusive regimes, including the exceptionally slow drop in resistivity for small line widths below 100 nm.

Although the number of experimental data points limits the number of regions we can reasonably consider in this first demonstration, this approach still offers unprecedented experimental access to the differential thermal conductivity contributions of phonons with different MFPs and for benchmarking theoretical predictions, including those from first-principles DFT calculations shown in Figs. 3*A* and 4. In particular, our experimental data across all MFP ranges measured in silicon are in very good agreement with our DFT calculations [which also agree with those in the literature (8)]. However, some small discrepancies appear for phonon MFPs around ∼100 nm, where experimental verification was not possible before. Differences between the experimental and theoretical spectra in this region may also be exaggerated by our limited set of small-line width gratings; increased resolution with a larger sample set can address this issue. Our data are also consistent with observations by others that emphasize significant contributions from long-MFP (>1-μm) modes (2, 3, 6), but the limited number of data points we have for structures much larger than the average phonon MFP results in a relatively large uncertainty in this region. For the purpose of comparison in Fig. 4, we normalized the experimental spectra in silicon by assuming the integrated conductivity up to 1 μm should match that calculated by DFT.

For sapphire, both calculation and experimental data imply that phonons with MFPs shorter than 1 μm are responsible for >95% of the thermal conductivity. The discrepancy below 300 nm between experimental and theoretical curves (apparent in the cumulative distribution) is due to two factors. First, the sharper rise in the DFT cumulative curve arises due to the very strong short-MFP peak in the conductivity spectrum—a peak that lies at ∼5 nm, below the lower bound of our experimental sensitivity (14 nm) using 30-nm structures. Thus, the experimental data simply do not include the shortest phonon MFP peak. Second, because of the complex crystal structure of sapphire, the DFT calculations required the use of relatively small interaction-distance cutoffs for determining the harmonic and anharmonic force constants, which may cause a larger error in the predictions than for silicon.

There are two significant advantages of our EUV-based technique compared with previously reported MFP spectroscopy techniques. First, our approach that combines nanoheaters with the phase sensitivity of short-wavelength probes is the only way to experimentally access dimensions far below 100 nm to directly resolve the contributions of phonons with MFP down to 14 nm. Other approaches require numerical extrapolation techniques and interpretation, which are still being developed (9). Second, characterizing phonon MFP spectra by harnessing interacting nanoscale heat sources allows us to probe arbitrary segments of the MFP spectrum for any novel material, enabling direct access to the differential, rather than only the cumulative, MFP distribution.

The ability to experimentally extract a phonon MFP distribution down to such small MFPs offers an exceptional, useful method for validating existing first-principles predictions across the whole range of phonon MFPs significant for heat conduction, as well as the first access to such information for more complex materials where calculations have not yet been performed. Furthermore, combined knowledge of both the differential and cumulative spectra may offer intriguing insight into the full MFP spectrum with the detail necessary for accurate prediction of heat transfer in nanostructured systems.

## Conclusion

In summary, we uncover a new regime of nanoscale thermal transport that dominates when the separation between nanoscale heat sources is small compared with the dominant phonon MFPs. We also present a new approach for characterizing the relative contributions of phonons with different MFPs that participate strongly in heat conduction. In particular, we can probe the small-MFP region, which has been previously inaccessible to experiment. This unique capability is important as the need for precise phonon MFP distributions in complex nanostructures becomes more pressing—for both fundamental understanding and to harness systems where modeling does not yet exist. In the future, because bright soft X-ray high harmonic sources can now reach wavelengths below 1 nm (23), this approach can be extended even further into the deep nanoregime. In addition, more comprehensive and fundamental insight into nanoscale thermal transport may be possible by adopting the framework created to bridge all types of anomalous diffusion (24⇓⇓–27), which does not rely on effective diffusion models for inherently nondiffusive heat transfer. Finally, the efficient, collectively diffusive regime of thermal transport that we observe for the first time can potentially mitigate projected problems for thermal management in nanoelectronics, where the power density is likely to increase as the individual nanostructures shrink in size (10, 11). It also highlights important design implications for nanostructured materials and devices for energy and biomedical applications.

## Acknowledgments

We gratefully acknowledge support from the US Department of Energy Basic Energy Sciences and the Semiconductor Research Corporation, and used facilities provided by the National Science Foundation (NSF) Engineering Research Center for EUV Science and Technology and a National Security Science and Engineering Faculty Fellowship award. K.M.H.-P. acknowledges support from the NSF under Award DGE 1144083. X.G. and R.Y. acknowledge the NSF CAREER award and Air Force Office of Scientific Research support.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: henry.kapteyn{at}colorado.edu.

This contribution is part of the special series of Inaugural Articles by members of the National Academy of Sciences elected in 2013.

Author contributions: K.M.H.-P., J.N.H.-C., E.H.A., W.C., R.W.F., R.Y., M.M.M., H.C.K., and D.N. designed research; K.M.H.-P., J.N.H.-C., T.D.F., and D.N. performed research; K.M.H.-P., J.N.H.-C., X.G., R.Y., and D.N. contributed new reagents/analytic tools; E.H.A. and W.C. fabricated samples; K.M.H.-P., J.N.H.-C., T.D.F., and D.N. analyzed data; and K.M.H.-P., J.N.H.-C., X.G., T.D.F., R.Y., M.M.M., H.C.K., and D.N. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503449112/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- King SW,
- Simka H,
- Herr D,
- Akinaga H,
- Garner M

- ↵
- ↵.
- Rundquist A, et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵COMSOL, Inc. (2013) COMSOL Multiphysics, Version 4.3b (COMSOL, Inc., Burlington, MA).
- ↵.
- Chen G,
- Borca-Tasciuc D,
- Yang RG

*Encyclopedia of Nanoscience and Nanotechnology*, ed Nalwa HS (American Scientific Publishers, Valencia, CA), Vol 7, pp 429–459 - ↵
- ↵
- ↵
- ↵.
- Popmintchev T, et al.

- ↵
- ↵
- ↵
- ↵