Valley-resolved edge plasmons are relevant to nano-optics at subwavelength scales. However, less attention has been paid to their tunable properties in time domain. In this work we investigate edge pseudomagnetoplasmons in a strained graphene modulated by multiple harmonics with frequency in the THz regime. The edge plasmon is described by a set of nonlinear hydrodynamic equations, which are self-consistently solved by the flux-corrected transport method. Without the applied voltage, there exist two unidirectional-propagating edge-plasmon modes with weak valley polarization . It is demonstrated that by varying the amplitude of multiple harmonics one can alter both the amplitude and the polarity of the valley polarization in the edge plasmon. One can achieve a full valley polarization at the instant of half cycle of the multiple harmonics and at the instant of one cycle. The edge-plasmon density and the transverse velocity vanish for the frozen valley.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
18 January 2019: A typographical correction was made to the title.
The idea of compressing electromagnetic waves laterally to subwavelength scales  has attracted considerable interest. Edge plasmons combine the confinement of electromagnetic waves and collective electron motions near the boundary of two-dimensional (2D) materials . In comparison to surface plasmons, edge plasmons can bring a stronger confinement of electromagnetic waves [3–6]. Fetter and Mast [3, 4] have provided the first experimental evidence for the presence of edge plasmons in a bounded two-dimensional electron gas (2DEG) under perpendicular magnetic fields. Aleiner and Glazman have discovered a new acoustic edge magnetoplasmons in a 2D electron liquid . For a 2DEG system modulated by opposite magnetic domains, a new type of one-way edge magnetoplasmons has been predicted .
Edge plasmons have also been predicted in graphene systems without external magnetic field. This kind of edge plasmons may elucidate different paradigms for realizing optical nonreciprocal devices at subwavelength scales and circumventing the requirement of strong magnetic fields. In a gapped graphene, a finite Berry flux provides a pseudomagnetic field in momentum space [9, 10]. Based on a hydrodynamic model, Song et al.  has predicted counterpropagating acoustic edge-plasmon modes (termed as “chiral Berry plasmons”) occurring on the boundary of a gapped graphene. These edge plasmons have splitting energy dispersions and could be useful for realizing optical nonreciprocity in the terahertz (THz) range [11, 12]. The quantum statistical effect and quantum diffraction effect render the chiral Berry plasmons unidirectionally-propagating . In a strained grapheme [13, 14], electrons in and valley will suffer different pseudomagnetic fields with opposite direction but identical strength. The strain-induced pseudomagnetic field can be uniform when the strain is designed to align along three main crystallographic directions . Under such a uniform pseudomagnetic field, the superposition of density oscillations in both valleys can lead to two counterpropagating acoustic modes localized on the boundary of a strained grapheme [15, 16]. Considering the nonlinear effect in the hydrodynamic model, we have demonstrated numerically that the edge plasmons can reach a full valley polarization under a pseudomagnetic field of tens of Tesla .
For the applications of valley-polarized edge plasmons, it is desirable to control dynamically the amplitude and/or polarity of the valley polarization. The dynamical generation of valley polarization in MoS monolayers has been reported in many works. In a pristine monolayer MoS Zeng et al.  achieved a valley polarization of by means of optical pumping. It has been predicted that in MoS monolayers under circularly-polarized light field one can obtain a complete dynamic valley polarization . Such a dynamical valley polarization can be suppressed/improved by tensile/compressive strain . For an epitaxial single-layer MoS on lattice-matched GaN substrates, enhanced valley helicity has been observed . These results demonstrated the feasibility of optical valley control in MoS monolayers [17–20]. All-optical control of the valley coherence has been demonstrated by means of the pseudomagnetic field associated with the optical Stark effect by using below-band gap circularly polarized light in monolayer WSe, which was manipulated optically by tuning the dynamic phase of excitons in opposite valleys . Wang et al.  probed the valley dynamics in monolayer WSe by monitoring the emission and polarization dynamics of the well-separated neutral excitons (bound electron-hole pairs) and charged excitons (trions) in photoluminescence, where a typical valley polarization decay time of the order of 1 ns was inferred. Valley-resolved pump-probe experiments demonstrated inverted valley polarization in optically-excited transition metal dichalcogenides . In a graphene mechanical resonator, quantum pumping has been proposed to generate valley-polarized bulk currents . For edge plasmons in graphene, the dynamical control of valley polarization has not been studied so far.
In this work, we study edge pseudomagnetoplasmons in a strained graphene modulated by a time-dependent voltage consisting of multiple harmonics with frequency in the THz regime. Our aim is to explore whether full valley polarization with a selected polarity can be achieved for an edge plasmon under a weak pseudomagnetic field by means of varying harmonic parameters (voltage amplitude, number of harmonics, fundamental frequency, and relative phase). The underlying mechanism of valley selection induced by multiple harmonics will be revealed. The outline of the paper is as follows. In Section 2, we present the two-component nonlinear hydrodynamic equations coupled with the Poisson equation to describe the considered plasmon motion. The flux-corrected transport algorithm  is adopted to self-consistently solve the nonlinear coupled equations. In Section 3, we show and discuss the voltage-amplitude dependence or time dependence of the valley-resolved electron densities and plasmon velocities as well as valley polarization for the edge plasmons. Finally, a brief conclusion is drawn in Section 4.
2. Theoretical model and numerical method
We consider a strained graphene sheet in the plane with a boundary at , which is schematically illustrated in Fig. 1. As in Refs. [15, 16], we assume that the strain-induced pseudomagnetic field is uniform in the graphene sheet. Electrons in the and valleys experience the pseudomagnetic field and for the edge plasmon propagating along the axis [15, 16]. For this "right-propagating edge plasmon", the propagating direction (), the direction of pseudomagnetic field , and the normal () of the boundary form a right-hand coordinate system [see Fig. 1(a)]. The edge is powered by a time-dependent voltage , which consists of multiple harmonics
Here is the voltage amplitude, is the number of the harmonics, is the fundamental frequency, and is the initial phase of the Ith harmonic.
The fundamental frequency is required to be larger than both the classical cyclotron frequency and the 2D bulk plasmon frequency . Here is the basic charge, is the plasmon mass, is the light velocity, is the reduced Planck’s constant, is the Fermi wave number, cm/s is the Fermi velocity, is the electron density in equilibrium, and is the propagating wave vector of the bulk plasmon. For cm, and Tesla, one has THz and THz. This indicates that the required voltage frequency is in the THz region. THz voltage waveforms have already been realized with widely tunable femtosecond pulses from a erbium-doped fiber source over ten years . The electron response to the THz fields has been reported by time-resolved THz spectroscopy . The conditions and accomplish the fastest response of graphene electrons to the applied voltage waveform.
The considered 2D electron system is treated as a two-component ( and ) Fermi liquid and described by the two-component hydrodynamic model . Microscopically, the electron motions are governed by the Boltzmann equation [28–30]. The macroscopic quantities of the Fermi liquid, such as the number density and mean velocity with representing the / valley, are described by the moments of the Boltzmann equation . Outside the graphene sheet () both and vanish. The basic equations for the and component under absolute zero temperature  include the continuity equations
Here , is the unit vector perpendicular to the graphene sheet, and is the self-induced electrostatic potential by the collective excitations of the electron system. satisfies the Poisson equation with the background of immobile positive carbon ions,Eq. (3), the first term is the electrostatic force including the external voltage and the induced potential, the second term is regarded as the Lorentz force due to the pseudomagnetic field, and the third term is the Fermi pressure due to the Pauli exclusion principle . For an electron in valley and another electron in valley with the same velocity, the Lorentz forces acting on them have opposite directions. As a consequence, the presence of pseudomagnetic fields may induce a valley-polarized edge plasmon.
The validity of Eqs. (2)–(3) is guaranteed by following principle conditions. We assume the edge as a straight line and the inter-valley scattering process is neglected . The effects due to the off-diagonal components of the Fermi pressure tensor are not included in the model . The inherent limitations of a fluid model, i.e. resonance effects of electron velocities matching the phase velocity or group velocity are neglected . The damping effect is neglected at zero temperature owing to the Pauli blocking . The applicable range of the hydrodynamic model  is , where for Tesla is the cyclotron radius, is the width of the bounded graphene system, is the Fermi wave length, and is the electron momentum-relaxation length.
The valley polarization of the acoustic pseudomagnetoplasmon at the edge and time is characterized by the relative difference between electron densities and in and valleys ,
In Ref. , the nonlinear hydrodynamic equations are linearized under the approximation so that an analytical solution can be obtained by means of the Wiener-Hopf technique. The linearized model neglects the mean velocity gradient that may contribute through the convective derivative in the momentum balance equation. It is thus not valid any more understrong modulations where the plasmon density may exceed the limitation . As pointed out in our previous work , the linearized hydrodynamic equations are not balance equations and cannot be solved numerically by a stable algorithm. It is evident that equations (2)-(4) are nonlinearly coupled each other. We have to solve self-consistently the plasmon density () and transverse and longitudinal velocities ( and ) for under the coexistence of the pseudomagnetic field and time-dependent voltage. It is hard to find analytical or numericalsolutions under linearized hydrodynamic equations. We adopt a 2D flux corrected transport (FCT) method  to solve the nonlinear hydrodynamic equations (2)-(3). The FCT method has been justified and described in detail in Refs.  and . It has the flexibility of adding source terms (such as the external voltage ) in the momentum balance equation. A two-stage Runge-Kutta time integration is performed to implement the 2D FCT algorithm with time step . The solutions are second-order accurate in the global truncation error. The Poisson equation (4) is discretized with the space step to form a set of linear equations. The latter have a symmetric and positive-definite coefficient matrix and can be solved by the successive over relaxation method.
3. Results and discussions
In the following numerical calculations, the fundamental frequency, pseudomagnetic field strength, and time step are fixed at THz, Tesla, and s. The space step satisfies the Courant limit . The initial phases are chosen as and for . In the initial time the system is in equilibrium and unstrained where the electron density in either valley is cm. The initial velocity of edge plasmon is chosen as , which is the same as in Ref. .
The pseudomagnetic field turns on from . A steady edge plasmon is launched after where is the period of the voltage waveform. The external voltages are applied after time to modulate the edge plasmon. In the following we express relevant quantities in dimensionless units where the density, velocity, length, and voltage amplitude are respectively in unit of , , , and .
In Fig. 2 we inspect the modulation effects on the edge pseudomagnetoplasmon at the end of one cycle () of voltage waveform. The values of the following physical quantities at this instant are plotted as a function of the voltage amplitude : the valley-resolved edge-plasmon density (), valley polarization , longitudinal velocity and transverse velocity at the edge . The waveform parameters are and . Without the voltage (), the edge plasmon density differs slightly from [Fig. 2(a)]. A weak pseudomagnetic field such as T brings only a small valley separation in the edge-plasmon density. This result agrees with that in . When applying the time-dependent voltage with a small amplitude , one can distinguish significantly and after one cycle. The reason for this modulation effect is that the Coulomb force at the edge enhances the electron collective motion at and near the sample edge. With increasing , decreases nonlinearly, while a nonmonotous variation is observed for . This characteristic variation shows some similarities with that under the change of pseudomagnetic field strength . The change of pseudomagnetic field strength requires keeping an appropriate geometry of the applied strain [35, 36] and the elastic deformation. The results in Fig. 2(a) offer a new method to modulate edge plasmons under a fixed (weak) pseudomagnetic field. This method establishes a link between real-space edge pseudomagnetoplasmons in graphene and multiple harmonics in the THz range.
Another feature in Fig. 2 is that after one cycle of the multiple harmonics with the edge-plasmon densities in both valleys are generally lower than the initial equilibrium density . This density reduction is caused by the negative part of the voltage waveform which spans a large fraction of the period [ for , see Fig. 1(b)]. It agrees well with the results in infrared nano-imaging experiments  in graphene. At the edge-plasmon density and , which violates the condition used in the linearized model. The linear-response theory is reasonable in the perturbative regime with weak modulations, leading to a small perturbation of equilibrium electron density. In the nonperturbative regime the nonlinear effect becomes important [16, 38]. This is because the density gradient in the continuity equation and the convective velocity term in the momentum-balance equation take steep changes at the edge, resulting in stronger nonlinearity and vanishing edge plasmon density. The FCT method can deal with the problem with strong nonlinearity under general initial and boundary conditions. In graphene, strong nonlinear response ofelectrons has been verified by time-resolved THz spectroscopy under a strong external THz electric field . Our results open up the possibility of investigating nonlinear behaviors of edge plasmons in strong modulation and ultrafast electrodynamics.
The relative difference between and results in the valley polarization of the edge plasmon, which is plotted in Fig. 2(b). It can be seen that is always valley-polarized () because . The amplitude of increases monotonously in a nonlinear way and finally a full valley polarization appears for . The valley polarization can be understood from the difference in the transverse and longitudinal components of the valley-resolved edge-plasmon velocities. From Fig. 2(c) one observes that the transverse velocities and have the same direction and increases with for . They are slightly separated at . The acceleration under the negative part of the voltage waveform distinguishes from . When is in the interval , the transverse velocity decreases monotonously to zero, while the corresponding edge-plasmon density decreases and eventually vanishes. The reduction of arises from the balance between the acceleration caused by the negative part of the waveform and the deceleration caused by the positive parts of the waveform. In contrast, in the same region of , decreases from to , corresponding to a finite edge-plasmon density in the valley. As shown in Fig. 2(d), under the multiple harmonics with small () the longitudinal plasmon velocities and decease almost linearly with and are well separated. With further increasing of , decreases while increases. They are equal at and have opposite directions for in the region and . Note that for near or , although is comparable to , the difference between and is remarkable. In this case the density oscillation in one valley is almost frozen (with a low or zero propagating speed). Accordingly, only one valley component of the density oscillation is propagating, indicating the formation of a new valley-polarized chiral edge-plasmon mode. The valley population of this mode can be tuned by ( valley for and valley for ), which is similar to the valley-selective circular dichroism of transition metal dichalcogenides. This dynamical valley-polarized edge-plasmon mode could be useful in the field of valleytronics [17, 39].
The polarity tunability of valley polarization is required for some applications of valley-polarized edge plasmons. This aim can be achieved if we focus on signals at the end of half cycle () of the voltage waveform, as shown in Fig. 3. At this instant, the edge-plasmon density decreases monotonously from to as varies from to [see Fig. 3(a)]. The edge-plasmon density increases firstly to the maximum (at ) and then decrease quickly tozero. At , the two density components are equal. The corresponding valley polarization is -polarized for and -polarized for [see Fig. 3(b)]. The amplitude of with negative polarity can approach . A full valley polarization appears at . Note that under the same voltage waveform with , full valley polarization with negative polarity appears at the end of one cycle [Fig. 2(b)]. Thus we can switch the polarity of valley polarization either by the voltage amplitude or the instant for output.
Figures 3(c) and 3(d) show the valley-resolved edge-plasmon velocities at the end of half cycle () of the voltage waveform. The transverse velocity increases nonlinearly to the maximum (at ) and then decreases, while shows an oscillatory variation for and then becomes negligible. Within half cycle, the deceleration caused by the positive parts of the waveform can be comparable to the acceleration caused by the negative part of the waveform for all values, leading to a shoulder feature in the curve and oscillation behavior in the curve. The longitudinal edge-plasmon velocities and are always positive, which is in contrast to the results in Fig. 2(d). In this situation there exist unidirectional edge-plasmon modes for both and valleys with a -tunable velocity and population. The branch disappears at where vanishes.
Next we examine the time evolution of edge-plasmon densities and under a THz optical pump pulse with several cycles [26, 40]. The result is shown in Fig. 4 where the parameters for multiple harmonics are , , and . The voltage waveform with three cycles is plotted in Fig. 1(a). It can be seen that globally decreases while increases with time. Both and peak at the beginning of the negative part of each cycle. The reason is that with increasing time, more electrons away from the edge can contribute to the edge plasmon excitation and move near or further to the edge during the positive or negative part of the voltage waveform. For a small voltage amplitude the difference between the edge-plasmon densities and as well as the valley polarization increases with time and becomes significant after several cycles. The infrared THz optical absorption measurements  and time-resolved terahertz THz spectroscopy  could be used to detect the time-dependent excitation of edge plasmons.
4. Conclusions and remarks
In conclusion, we have employed a two-component nonlinear hydrodynamic model to investigate the valley-resolved edge pseudomagnetoplasmons in a strained graphene sheet. The electrons in graphene is modulated by multiple harmonics with frequency in the THz regime. The nonlinear hydrodynamic model is solved by means of the FCT method. We have calculated the edge-plasmon density and velocity for the two valley components. We start from an equilibrium state where there exist two unidirectional-propagating edge-plasmon modes with weak valley polarization . After the application of multiple harmonics, one can achieve a full valley polarization at the instant of half cycle and at the instant of one cycle. In the case of full valley polarization, the edge-plasmon density and the transverse velocity vanish for one valley. It is demonstrated that by varying the amplitude of multiple harmonics, one can alter both the amplitude and the polarity of the valley polarization in the edge plasmon for a given instant. Under a fixed small amplitude of multiple harmonics, the valley polarization increases with time and becomes significant after several cycles. Our results demonstrate the nonlinear response of edge plasmons in strained graphene to strong modulation of THz pulse, which is relevant to dynamical control of edge plasmons.
For practical applications, the propagation loss of edge plasmon in graphene should be considered. Scattering processes that relax momentum conservation contribute to propagation losses. In addition, intervalley scattering introduces a mechanism of nonconservation of the valley density, yielding propagation loss of edge plasmon in a specific valley. In the FCT numerical code, both intravalley and intervalley scattering terms can be treated as a source term and can be added flexibly to the momentum balance equation.
National Natural Science Foundation of China (NSFC) (11775090, 11775164, 11774314, 11575135); Fundamental Research Funds for the Central Universities (WUT: 2017IVA79, 2018IB011).
2. A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-González, M. Carrega, K. Watanabe, T. Taniguchi, G. Vignale, and M. Polini et al., “Highly confined low-loss plasmons in graphene–boron nitride heterostructures,” Nat. Mater. 14, 421 (2015). [CrossRef]
4. A. L. Fetter, “Edge magnetoplasmons in a bounded two-dimensional electron fluid,” Phys. Rev. B 32, 7676 (1985). [CrossRef]
5. J. C. Song and M. S. Rudner, “Chiral plasmons without magnetic field,” Proc. Nat. Acad. Sci. 1134658 (2016).
6. Y. Zhang, F. Zhai, B. Guo, L. Yi, and W. Jiang, “Quantum hydrodynamic modeling of edge modes in chiral berry plasmons,” Phys. Rev. B 96, 045104 (2017). [CrossRef]
8. D. Jin, L. Lu, Z. Wang, C. Fang, J. D. Joannopoulos, M. Soljačić, L. Fu, and N. X. Fang, “Topological magnetoplasmon,” Nat. Commun. 7, 13486 (2016). [CrossRef]
9. D. Xiao, W. Yao, and Q. Niu, “Valley-contrasting physics in graphene: magnetic moment and topological transport,” Phys. Rev. Lett. 99, 236809 (2007). [CrossRef]
10. D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, “Coupled spin and valley physics in monolayers of mos 2 and other group-vi dichalcogenides,” Phys. Rev. Lett. 108, 196802 (2012).
12. D. Jalas, A. Petrov, M. Eich, W. Freude, S. Fan, Z. Yu, R. Baets, M. Popovic, A. Melloni, and J. D. Joannopoulos, “What is–and what is not–an optical isolator,” Nat. Photonics 7, 579 (2013). [CrossRef]
13. F. Guinea, M. I. Katsnelson, and A. K. Geim, “Energy gaps and a zero-field quantum hall effect in graphene by strain engineering,” Nat. Phys. 6, 30 (2010). [CrossRef]
14. M. A. Vozmediano, M. Katsnelson, and F. Guinea, “Gauge fields in graphene,” Phys. Rep. 496, 109–148 (2010). [CrossRef]
15. A. Principi, M. I. Katsnelson, and G. Vignale, “Edge plasmons in two-component electron liquids in the presence of pseudomagnetic fields,” Phys. Rev. Lett. 117, 196803 (2016).
16. Y. Zhang, B. Guo, F. Zhai, and W. Jiang, “Valley-polarized edge pseudomagnetoplasmons in graphene: A two-component hydrodynamic model,” Phys. Rev. B 97, 115455 (2018). [CrossRef]
18. K. F. Mak, K. He, J. Shan, and T. F. Heinz, “Control of valley polarization in monolayer mos 2 by optical helicity,” Nat. Nanotechnol. 7, 494 (2012). [CrossRef]
19. S. Aas and C. Bulutay, “Strain dependence of photoluminescence and circular dichroism in transition metal dichalcogenides: ak· p analysis,” Opt. Express 26, 28672–28681 (2018). [CrossRef]
20. Y. Wan, J. Xiao, J. Li, X. Fang, K. Zhang, L. Fu, P. Li, Z. Song, H. Zhang, and Y. Wang et al., “Epitaxial single-layer mos2 on gan with enhanced valley helicity,” Adv. Mater. 30, 1703888 (2018). [CrossRef]
21. Z. Ye, D. Sun, and T. F. Heinz, “Optical manipulation of valley pseudospin,” Nat. Phys. 13, 26 (2017). [CrossRef]
22. G. Wang, L. Bouet, D. Lagarde, M. Vidal, A. Balocchi, T. Amand, X. Marie, and B. Urbaszek, “Valley dynamics probed through charged and neutral exciton emission in monolayer wse 2,” Phys. Rev. B 90, 075413 (2014). [CrossRef]
23. G. Berghäuser, I. Bernal-Villamil, R. Schmidt, R. Schneider, I. Niehues, P. Erhart, S. M. de Vasconcellos, R. Bratschitsch, A. Knorr, and E. Malic, “Inverted valley polarization in optically excited transition metal dichalcogenides,” Nat. Commun. 9, 971 (2018). [CrossRef] [PubMed]
25. J. P. Boris, A. M. Landsberg, E. S. Oran, and J. H. Gardner, “Lcpfct-a flux-corrected transport algorithm for solving generalized continuity equations,” Tech. Rep., NRL Memorandom Report No. 6410 (Naval Research Laboratory, Washington, D.C), pp. 20375-5320 (1993).
27. P. Bowlan, E. Martinez-Moreno, K. Reimann, T. Elsaesser, and M. Woerner, “Ultrafast terahertz response of multilayer graphene in the nonperturbative regime,” Phys. Rev. B 89, 041408 (2014). [CrossRef]
28. L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics: Fluid Mechanics(Pergamon Press, Oxford,1987).
30. U. Briskot, M. Schütt, I. Gornyi, M. Titov, B. Narozhny, and A. Mirlin, “Collision-dominated nonlinear hydrodynamics in graphene,” Phys. Rev. B 92, 115426 (2015). [CrossRef]
31. F. Haas, Quantum Plasmas An Hydrodynamic Approach(Springer, 2011). [CrossRef]
33. Z. Li, E. A. Henriksen, Z. Jiang, Z. Hao, M. C. Martin, P. Kim, H. Stormer, and D. N. Basov, “Dirac charge dynamics in graphene by infrared spectroscopy,” Nat. Phys. 4, 532 (2008). [CrossRef]
35. F. Guinea, M. Katsnelson, and M. Vozmediano, “Midgap states and charge inhomogeneities in corrugated graphene,” Phys. Rev. B 77, 075422 (2008). [CrossRef]
36. A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. 81, 109 (2009). [CrossRef]
37. Z. Fei, A. Rodin, G. Andreev, W. Bao, A. McLeod, M. Wagner, L. Zhang, Z. Zhao, M. Thiemens, and G. Dominguez et al., “Gate-tuning of graphene plasmons revealed by infrared nano-imaging,” Nature 487, 82 (2012). [CrossRef] [PubMed]
38. Y. Zhang, F. Zhai, and L. Yi, “Study of spin-polarized plasma driven by spin force in a two-dimensional quantum electron gas,” Phys. Lett. A 380, 3908–3913 (2016). [CrossRef]
39. J. Kim, X. Hong, C. Jin, S.-F. Shi, C.-Y. S. Chang, M.-H. Chiu, L.-J. Li, and F. Wang, “Ultrafast generation of pseudo-magnetic field for valley excitons in wse2 monolayers,” Science 346, 1205–1208 (2014). [CrossRef] [PubMed]
40. J. Sotor and G. Sobon, “24 fs and 3 nj pulse generation from a simple, all polarization maintaining er-doped fiber laser,” Laser Phys. Lett. 13, 125102 (2016). [CrossRef]