## Abstract

Recently, the ray tracing method has been used to derive the non-singular cylindrical invisibility cloaks for out-of-plane shear waves, which is impossible via the transformation method directly owing to the singular push-forward mapping. In this paper, the method is adopted to design a kind of non-singular acoustic cloak. Based on Hamilton's equations of motion, eikonal equation and pre-designed ray equations, we derive several constraint equations for bulk modulus and density tensor. On the premise that the perfect matching conditions are satisfied, a series of non-singular physical profiles can be obtained by arranging the singular terms reasonably. The physical profiles derived by the ray tracing method will degenerate to the transformation-based solutions when taking the transport equation into consideration. This illuminates the essence of the newly designed cloaks that they are actually the so-called eikonal cloaks that can accurately control the paths of energy flux but with small disturbance in energy distribution along the paths. The near-perfect invisible performance has been demonstrated by the numerical ray tracing results and the pressure distribution snapshots. Finally, a kind of reduced cloak is conceived, and the good invisible performance has been measured quantitatively by the normalized scattering width.

## 1. Introduction

It is interesting that the transformation method and conformal mapping method were proposed by Pendry *et al.* [1] and Leonhardt [2], respectively, almost at the same time to achieve the passive cloaking for electromagnetic waves. Owing to its intuitive physical picture and simple mathematical operation, the transformation method rapidly became the mainstream strategy for energy controlling to acoustic waves [3–7], linear surface water waves [8,9], bending waves [10–14], out-of-plane shear waves [15,16] and thermodynamics [17–19].

In the realm of elastodynamics, however, Milton *et al.* [20] demonstrated that Navier equations are not form invariant under a general mapping and will be transformed to Willis equations with a particular choice of gauge. Brun *et al.* [21] have designed a cylindrical cloak with recourse to the Cosserat-type materials for achieving control of in-plane elastic waves. In 2011, Norris & Shuvalov [22] developed the transformation theory in a more general elasticity framework and fuelled the interest in transformation elastodynamics. Nevertheless, the non-symmetric constitutive relations introduced in the aforementioned studies pose real challenges for the further development of elastodynamic cloaking. More recently, Chang *et al.* [23] showed that the elastic ray path can be controlled exactly with the help of transformation ray theory in which only a symmetric elastic tensor and density tensor are required. Hu *et al.* [24] also re-examined the transformation theory from the deformation perspective and obtained a set of symmetric transformed relations for controlling elastic waves approximately. Based on the former results, some functional devices with isotropic materials have been designed under local conformal mappings [25,26].

In this paper, we will limit our topics to the acoustic case. In 2007, Cummer *et al.* [3] derived the transformation-based solution for acoustic cloak in the two-dimensional geometry via a variable exchange owing to the isomorphism between the single polarization Maxwell equations and acoustic equations. Soon afterwards, the physical profile of a spherical acoustic cloak was obtained by an analogy mapping to the conductivity equation [4] and directly through the scattering theory [27], respectively. Although the transformation-based solutions were theoretically demonstrated to be perfect in the two- and three-dimensional cases, they are impractical to be fabricated in reality because of the infinite mass dilemma. To overcome this drawback, Norris [5] generalized the governing equation of acoustic waves from the inviscid ideal fluid to the metafluid and obtained the finite mass cloaks with anisotropic stiffness and tensor density. Another approach to eradicate the singular dilemma is to apply the so-called regularized transformation [16,28,29] which involves deforming a small region rather than a mere point to the cloaking region. The method solved the problem to some extent, however, it should be emphasized that with a small regularization parameter the material parameters near the inner boundaries are still extremely large and are impossible for practical fabrication. Inspired by the deformation description in the continuum mechanics, Chang *et al.* [6] regarded the geometry transformation as deformation in the language of finite elasticity and gained an intrinsic comprehension of the transformation method for electromagnetics, acoustics and elastodynamics. Based on the deformation perspective, Hu *et al.* [7] revealed the underlying mechanism of the non-uniqueness of transformation acoustics and derived several sets of transformation relations by considering energy conservation and form invariance of governing equations under the local affine transformation. Apart from the above-mentioned studies, some other approaches, such as impedance cloaking [30], scattering cancellation cloaking [31] and active exterior cloaking [32,33] were also proposed in the cloaking research.

The generalized transformation acoustics results in metamaterials with tensor density or anisotropic stiffness which cannot be found in nature [5,7]. In order to verify the invisible performance experimentally, some studies about acoustic metamaterials and fabrication of invisibility cloaks have been implemented recently. Some approaches, such as the alternating layered system [34–37], the composite materials with periodically arranged sound scatterers embedded in a fluid or gas matrix [38–40] and metamaterials made of perforated plastic plates [41], have been reported to realize the anisotropic density. It should be emphasized that Norris *et al.* [36] and Urzhumov *et al.* [37] tried to devise cloaks by choosing the cloaking transformation that fits the multilayered material model with attainable materials, instead of seeking challenging metamaterials to meet the demand of a particular cloaking transformation. The pentamode materials, which were proposed by Milton & Cherkaev [42] in 1995, are another generic metafluid for achieving anisotropic stiffness promisingly. A simple sample has been fabricated with the current state-of-the-art lithography recently [43]. On the premise of the above studies, some invisibility cloaks have been fabricated and verified experimentally. For instance, Popa *et al.* [44] and Zigoneanu *et al.* [41] have demonstrated the effectiveness of ground cloaks experimentally in the two- and three-dimensional cases, respectively. An acoustic cloak comprised perforated plates for generating virtual images was validated by Hu *et al.* [45]. Apart from the above-mentioned homogeneous cloaks, an annulus cloak constructed with a network of acoustic circuit elements [46], and a directional three-dimensional cloak designed on the basis of the scattering cancellation technique [47] were also validated experimentally.

Norris [5] proved that the singular mapping will inevitably result in the infinite mass dilemma for inertial cloaks, a generic class of cloaks defined by scalar bulk modulus combined with anisotropic density. However, avoiding the singular profile and obtaining a near-perfect cloak is possible through the ray tracing method which is based on the high-frequency approximation theory, as manipulated in the recent paper about cloaking for out-of-plane shear waves [48]. Thanks to the isomorphism in the two-dimensional case, the results reported in [48] can be extended to the acoustic case directly via variable exchanges. Nevertheless, similar analogy does not exist in the three-dimensional case. In this paper, we try to devise a broad family of non-singular acoustic cloaks by examining push-forward mappings from a ray trajectory perspective. It is shown that cloaks constructed by the present method are actually the so-called eikonal cloaks [5,37], which will not cause any scattering at the cloaks' outer boundary and will return to perfect cloaks when the transport equation is taken into consideration. The results contribute to illuminate the difference between eikonal cloaks and perfect cloaks essentially.

This paper is structured as follows. In §2, several constraint equations governing the variation of material parameters of eikonal cloaks are derived through Hamilton's equations of motion, the eikonal equation and the pre-designed ray equations. The perfect matching conditions are considered in order to reduce radiation strength caused by the impedance mismatch. As an illustration, the classical linear transformation mapping [1,4] is adopted to design a spherical cloak without singular parameters for both density and bulk modulus. Actually, the kind of cloak still has singularity, because infinite speed of sound is required at the inner boundary. To thoroughly overcome this problem, a ‘reduction’ manipulation is implemented to a cylindrical cloak for later numerical demonstration. In §3, some numerical simulations are implemented to demonstrate the cloaking performance for both eikonal cloak and ‘reduced’ cloak. Finally, some concluding remarks from this study are summarized in §4. Here, the summation convention over repeated lowercase Latin indices is adopted over the range 1–3.

## 2. Design of the spherical cloak

There is no doubt that a push-forward mapping, which involves deforming a region such that a point is mapped to the inner boundary of a cloak, will inevitably result in a singular inertial cloak. This paper adopts the mapping, but from a ray trajectory perspective, to derive near-perfect spherical cloaks (PSCs) with adjustable non-singular material parameters. Figure 1 illustrates the analytical model in which a spherical cloak, with the inner radius *a* and outer radius *b*, is presented with some rays winding around it. In the spherical coordinates, the cloak is orthotropic and should be characterized by two physical quantities, namely the scalar bulk modulus λ and tensor density ** ρ** with principal components

*ρ*

_{r}(

*r*),

*ρ*

_{θ}(

*r*) and

*ρ*

_{φ}(

*r*). For convenience but without loss of generality, a time harmonic plane wave

*p*

_{0}=e

^{i(k0x2−ωt)}which propagates along the positive direction of the

*x*

_{2}-axis with unit amplitude is considered, where

*k*

_{0}is the wave number in the homogeneous surrounding medium

*Ω*

_{0}, and

*ω*is the circular frequency of the incident wave. Correspondingly, the transmitted wave in the cloak

*Ω*

_{c}can be represented as

*p*

_{c}=e

^{i(kc⋅x−ωt)}, where

**k**

_{c}is the wavevector varying from position to position.

### (a) Ray trajectory formulation

The ray tracing method is a straightforward manner to investigate propagation characteristics for various kinds of wave phenomena. Its theoretical foundation is the high-frequency asymptotic method that uses an approximate time harmonic high-frequency solution of the governing equation to deduce the eikonal equation [49]. Subsequently, the eikonal equation is applied to formularize the Hamiltonian for tracing the energy flux. The governing equation of acoustic waves in metafluid with tensor density reads as
*α*=1/λ is the compressibility, and ** ρ**=

*ρ*

_{0}

**I**corresponds to the ideal fluid. Inserting the ansatz solution

*ω*, we have

*P*(

**x**) and

*T*(

**x**) are the smooth scalar functions of coordinates, equation

*t*=

*T*(

**x**) represents the moving wavefront and ∇

*T*=

**k**/

*ω*is the slowness vector of the wave under consideration. Because (2.2) should be satisfied for any frequency, the expressions with

*ω*

^{2},

*ω*

^{1}and

*ω*

^{0}must vanish and we obtain the following equations

As for this study, the eikonal equation along the ray trajectory can be written in the cloak and surrounding medium, respectively, as
*ρ*_{ij} and *α*_{c} are the density component and compressibility of the spherical cloak, respectively, *ρ*_{0} and *α*_{0} are the density and compressibility of the homogeneous medium in *Ω*_{0}, respectively, and **1**=**k**_{c}/*k*_{c}, where **1** (with component *l*_{i}) represents the unit direction vector. Eliminating *ω*^{2} from (2.6) and (2.7) and notating the normalized inverse density tensor, compressibility and wavenumber, respectively, as *α*=*α*_{c}/*α*_{0} and *k*=*k*_{c}/*k*_{0}, we obtain the Hamiltonian with the identical form as that in [48]
*n* is the refractive index varying with position **x** and orientation **I** and its square can be expressed in the Cartesian coordinates as *n*^{2}=*α*/*c*_{ij}*l*_{i}*l*_{j}. With regard to ray tracing for wave propagation, Hamilton's equations of motion have the form [48–51]
*τ* parametrizes the path and **s** denotes the ray vector which is tangential to the ray trajectory and indicates the direction of energy flux. Following (2.9), it can be seen that *τ* and **s** are connected by the arc length |d**x**|=|**s**|d*τ*. The ray vector **s** can be written in the Cartesian coordinates as *s*_{i}=*c*_{ij}*k*_{j}/*c*_{pq}*l*_{p}*l*_{q}.

The transformation method interprets a push-forward mapping as a geometry transformation from a virtual space to the physical space. Actually, an alternative perspective, regarding the space transformation as the winding rays of energy flux, is more straightforward, as shown in Figure 1*b*. Here, we examine the following push-forward mapping
*R*,*Θ*,*Φ*) and lowercase letters (*r*,*θ*,*φ*) are the spherical coordinates in the virtual space and physical space, respectively. For the sake of simplicity, a line which is parallel to the *x*_{2}-axis in the virtual space is considered, and it can be formulated in the spherical coordinates as *b*,*θ*_{0},*φ*_{0}) are the spherical coordinates of the incident point. According to (2.11), the ray equations in the physical space are formulated as

### (b) Constraint equations of material parameters

The components of a tensor in the Cartesian coordinates (*x*_{1},*x*_{2},*x*_{3}) and those in the spherical coordinates (*r*,*θ*,*φ*) are related through the transformation of coordinates. In regard to the vector and second-order tensor, the transformation relations are formulated in the following form, respectively
*h*_{r}=1, *h*_{θ}=*r*, *J*_{mm′} is the Jacobian transformation matrix and the quantities with and without prime represent the components in the spherical and Cartesian coordinates, respectively. The underline for repeated indices means equal values only, and the summation convention is not used here.

After some mathematical manipulation, we have the following differential equations, which link cloaks' material parameters to a specific transformation mapping, from the first Hamilton's equation of motion.
*c*_{r}, *c*_{θ} and *c*_{φ} are the three non-vanishing principal components of the normalized inverse density tensor. The Hamiltonian equals to zero along the ray trajectory, which means that the eikonal equation is satisfied if the ray propagates along its real path. Thus, the combination of (2.8), (2.9) and (2.13) yields the following equation
*k*_{r},*k*_{θ},*kφ*)). Therefore, the specific expressions of *k*_{r}, *k*_{θ} and *k*_{φ} can be obtained tactfully by using the ellipsoid parametric equations as long as the form of (2.14) and (2.15) holds.
*b*, it is evident that *k*_{φ} should take a positive sign. Subsequently, the signs of *k*_{r} and *k*_{θ} can be easily defined according to (2.15) and (2.14), respectively. Finally, the combination of (2.14), (2.15) and (2.17) yields the following constraint equations
*τ* and *φ* are connected by the differential relation |d**x**|=|**s**|d*τ*, then we have
*τ* increases with the azimuthal coordinate *φ*, so the positive sign is adopted in the expression. Substituting (2.14), (2.15) and (2.22) into (2.23), we have after some manipulation
*n*^{2} with respect to the polar angle *θ* and using transformation relations in (2.13), the right-hand side of (2.20) is expressed, after some manipulation, as

### (c) Perfect matching conditions

Some boundary conditions are necessary in order to solve the differential equation in (2.26). In regard to cloaking, the interface should be non-reflecting, which means the coexistence of incident wave and transmitted wave only is capable of keeping the continuity conditions. Following the time harmonic plane wave assumption and (2.13) in [5], the continuity conditions at *r* = *b* can be written as

The continuity of the tangential component of normalized wavevector is attained owing to the continuity of pressure, thus
**k**_{0}=(0,1,0) is the normalized wavevector in the surrounding medium. Following (2.28) and the continuity of normal acceleration as expressed in (2.27), we have the following two boundary conditions
*b*,*θ*_{0},*φ*_{0}) to (*r*,*θ*,*φ*) and according to (2.29), we have

### (d) Physical profiles of spherical cloak

In order to clarify the similarity and difference between the transformation method and the ray tracing method, we re-examine the classical linear transformation mapping [4]
*p*_{c}=e^{i(kc⋅x−ωt)} in the cloak. It is clear, as a consequence of the constant amplitude of transmitted wave, that the identity (2.5) gets satisfied naturally and (2.4) degenerates to the following expression
**k**_{c}=*k*_{0}**k** and inserting (2.17) and (2.33) into (2.34), we have after some manipulation

In this paper, our main concerns are not paid on achieving the perfect invisibility, but rather on eliminating the singularity. Thus, (2.4) and (2.5) which correspond to the lower-order terms in (2.2) will be neglected, and the eikonal equation is sufficiently accurate, as demonstrated in the out-of-plane shear wave case [48], for controlling high-frequency acoustic waves. With four material parameters controlled by three constraint equations, there are numerous schemes to formulate physical profiles as long as the continuity of radial impedance is satisfied, which would make it possible to eliminate the singularity of material parameters. As an example, the following non-singular profile will be considered for numerical illustration.
*et al.* [52] in the area of acoustics, it is found that the same profile as (2.36) can be obtained by introducing a virtual principal stretch λ_{v}=[*b*/(*b*−*a*)]^{2}[(*r*−*a*/*r*)]^{2} in a fictitious four-dimensional space. Although the mathematical manipulation in [52] is easier than that of the ray tracing method, it cannot be ignored that the deformation perspective is too abstract to have a fully understanding of the physical meanings underlying the non-singular profile. For the narrative convenience, cloaks with physical profiles as shown in (2.35) and (2.36) will be called PSC and ESC, respectively.

### (e) Reduction manipulation of cylindrical cloak

Actually, the aforementioned process has not solved the singularity problem completely, because the infinite azimuthal speed of sound still exists for the profile listed in (2.36). It has been demonstrated that the maximum azimuthal speed of sound in the metamaterial determines the smallest possible cloaking deficiency [37]. Thus, a ‘reduction’ process is considered based on some results shown in [48] for out-of-plane shear waves. It is shown that the cloak derived by the ray tracing method tends to have constant cloaking performance to a change in boundary condition. The underlying reason is that, with the ray trajectories bulging outwards, the area near the inner boundary tends to have little influence on the whole cloaking performance, which is comprehensible from the point of high-frequency ray theory. Consequently, we expect to construct a kind of cloak in which the area with singular parameters could be removed without sacrificing much invisibility. For the sake of simplicity, a cylindrical cloak is considered to illustrate some concepts during the reduction process. The following power law mapping is efficient to bulge ray trajectories outwards
*η*. The adjustable parameter *η* provides us much freedom in the reduction operation. As discussed above, constraint equations are easily available through direct variable exchanges. As an example, we choose the following set of material parameters.
*a* and *b* have the same meanings as their counterparts in the spherical cloak. For convenience, cylindrical cloaks with and without reduction operation will be called ‘reduced’ cylindrical cloak (RCC) and eikonal cylindrical cloak (ECC), respectively. The parameter *g*=(*r*_{0}−*a*)/*b*−*a*) is defined as the reduction factor which means that the annulus *a*<*r*<*r*_{0} is removed in practical fabrication. As will be shown in §3, an appropriate reduction factor is necessary in order to keep a balance between material availability and cloaking performance.

## 3. Results and discussions

In §2, parallel incident rays which correspond to the energy flux of plane wave are considered for deriving the physical profiles of spherical cloak. However, it should be shown that the effectiveness of the derived cloak, owing to its geometrical symmetry, is waveform independent. In this section, several numerical simulations about ray tracing and pressure distribution were implemented by the COMSOL Multiphysics solver to validate the eikonal cloak's invisible performance. To be specific, the ray tracing simulations for plane wave and spherical wave were carried out by the mathematical particle tracing model with calculation model shown in Figure 1*a*, whereas the pressure distribution snapshots with and without cloaks were simulated by the PDE solver with perfectly matched layers [53] enclosing the computation domain. The spherical cloak is located at the centre of the simulation domain and a point source is applied at (−1.1 m, 1.1 m, 1.1 m). The sound-hard scatterer is achieved by setting the inner boundary *r* = *a* to be Neumann boundary. To facilitate the computation, non-dimensional material parameters are taken as *ρ*_{0}=λ_{0}=1 for homogeneous surrounding medium, and geometry parameters of cloak are taken as *a*=0.5 m, *b*=1 m.

Figure 2 shows the ray tracing results for plane wave (part (*a*)) and spherical wave (part (*b*)). For clarity, rays emitted from the point source are traced in a plane passing through the point source and the centre of the spherical cloak, as illustrated in Figure 2*b*. It can be seen that the incident rays smoothly wind around the scatterer and perfectly return to their original paths as if nothing were there. The colour scales depict the energy flux velocity along the rays. Intuitively, the rays have to propagate faster when they round the obstacle to keep the phase consistent. Because the material parameters given by (2.36) make the eikonal equation get fully satisfied, the energy flowing trajectories should be precisely controlled, which is in good agreement with the simulation results shown in Figure 2.

Figure 3 shows the simulation model (part (*a*)) and pressure distribution snapshots for uncloaked field (part (*b*)) and cloaked fields with PSC (part (*c*)) and ESC (part (*d*)). The snapshots in Figure 3*b*–*d* are taken in a plane, as shown in Figure 3*a*, and all the figures share the same colour scale. It can be found that both the PSC and ESC reduce the scattering field significantly, which verifies the validity of the ESC sufficiently. As mentioned in the previous discussions, the non-singular physical profile makes near-perfect cloaking for high-frequency acoustic waves, because we neglect the lower-order terms in (2.2). Thus, it is necessary to demonstrate the difference between the PSC and ESC. Figure 4 illustrates the pressure distribution for several cases along a line passing through the point source and the centre of the cloak, as the line shown in the lower right corner. The dashed line, dash-dot line, dotted line and solid line denote the ideal field, uncloaked field with sound-hard scatterer, cloaked fields with PSC and ESC, respectively. Compared with the ideal field, it is evident that the sound-hard scatterer causes not only serious amplitude disturbance, but also significant phase shift, whereas the curves for PSC and ESC are nearly overlapping with the curve for the ideal field, which demonstrates the effectiveness of the ESC clearly. It should be pointed out that the small difference between the curves for the ideal field and cloaked field with PSC results from numerical error. Compared with the curve for PSC, the near overlapping curve for ESC is phase consistent but with small difference in amplitude distribution, which illuminates the eikonal cloak's near-perfect invisible performance essentially. Namely, cloaks governed by (2.18) and (2.30) can control the paths of energy flux accurately but with small disturbance in energy distribution in that the eikonal equation is fully satisfied, whereas the transport equation is not the case.

Scattering cross-section and scattering width, which denote the ratio of scattering energy flux to incident energy flux, are often used to quantify the far field scattering strength for electromagnetic waves and acoustic waves [54–56]. Here, the normalized scattering width *σ*_{sw} (normalized to a bare scatterer with radius *r*=*a*) is adopted to measure the cloaking performance quantitatively. Numerical simulations are implemented with built-in program algorithm in the finite-element software COMSOL Multiphysics. The scatterer is achieved by setting inner boundary as Dirichlet boundary. The material parameters of surrounding medium and geometry parameters of cylindrical cloak are the same as those in the previous simulations, and some other parameters are set as: *η*=3 and *g*=0.2. Figure 5 shows the variation of normalized scattering width *σ*_{sw} for a cloaked object with angular frequency *ω*. It can be observed that *σ*_{sw} increases with increasing *ω* for three kinds of cloaks, whereas their difference is decreasing. In addition, the normalized scattering width of RCC almost drops 10 dB when *ω*<13 rad s^{−1}. These findings imply that the reduced cloak is efficient enough to reduce the far field scattering caused by the enclosed scatterer. Figure 6 shows the distribution of material parameters of the cylindrical cloak along the radial coordinate, and the solid lines marked with square, triangle and circle symbols represent *ρ*_{r}, *ρ*_{θ} and λ, respectively. The purple dash-dot line denotes the position of the RCC's inner boundary. Evidently, the annulus area with extreme material parameters has been removed, and the parameters required in the remaining part of the cloak are within the range of natural available materials. However, achieving independent control of density and bulk modulus of a material is really a challenging task owing to their correlation which is often expressed in the form of an Ashby chart [37]. Although some novel ideas, such as seeking desired cloaking transformation to fit a pre-designed material model [36,37], were proposed to construct cloaks with attainable materials; the fact that a small filling fraction of the gaseous fluid is required near a cloak's inner boundary poses a real challenge to practical fabrication. Obviously, there is still a lot of work ahead.

## 4. Concluding remarks

We extended the ray tracing method to the acoustic field, which makes it possible to construct a generic acoustic cloak without singular material parameters. The eikonal equation which corresponds to the higher-order term in (2.2) was adopted to formulate the Hamiltonian. Following Hamilton's equations of motion, three differential equations (2.14), (2.15) and (2.26) which determine the physical profiles of spherical cloak were derived. According to the pre-designed ray equations and perfect matching conditions, three constraint equations which specify the near-perfect eikonal cloaks were expressed in equations (2.18) and (2.30). The number of constraint equations is less than that of material parameters, thus a series of non-singular physical profiles can be obtained by arranging the singular terms *r*−*a* tactfully. This facilitates the experimental fabrication and verification tremendously. Cloaks governed by equations (2.18) and (2.30) are not perfect, because we omitted the lower-order terms in equation (2.2). A set of non-singular material parameters, as expressed in equation (2.36), was employed to validate the invisible performance. The ray tracing results and pressure distribution snapshots showed that the eikonal cloak was accurate enough for controlling the paths of energy flux but with neglectable disturbance in energy distribution. Finally, we conceived a kind of reduced cloak to resolve the singular problem thoroughly. The numerical analysis of normalized scattering width shows that the reduction manipulation will not sacrifice much invisibility if an appropriate reduction factor is chosen under the power law mapping.

## Ethics

This work does not involve human subjects or animals.

## Data accessibility

This work does not contain any experimental data.

## Authors' contributions

P.L.G. carried out the mathematical work and numerical simulations, and wrote the manuscript with the assistance and supervision of L.Z.W.

## Competing interests

We have no competing interests.

## Funding

This work was supported by the National Science Foundation of China under grant no. 11072068.

## Acknowledgements

The authors appreciate the constructive comments from anonymous reviewers, and P.L.G. acknowledges X. He for her helpful discussion.

- Received May 28, 2015.
- Accepted January 14, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.