**200 mm×200 mm**aperture was chosen to test the feasibility of our proposed method. The results indicate that sparse stitching is suitable for measuring a large planar surface polished by ring polishing. The results of two chain lattice are closer to the fully covered lattice than a one chain lattice. However, more measuring time could be saved by the one chain lattice method. So, the usage of different lattices could be adapted for different periods of large aperture planar optical manufacturing.

## 1.

## Introduction

High-power solid-state laser equipment, such as National Ignition Facility (NIF),^{1}2.^{–}^{3} laser mégajoule (LMJ),^{4}^{,}^{5} requires a great number of large aperture planar optical elements. In NIF, more than 7000 planar optical elements with an aperture larger than $400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ are required.^{6} In the optical manufacturing, the surface shape measurement is a significant process to determine the precision of the optical elements.^{7}8.9.^{–}^{10} Thus, a measuring method with high precision and efficiency is required.

The surface shape measuring method of large aperture optical elements includes the Ritchey–Common test, interferometry, and subaperture stitching interferometry. Whether the Ritchey–Common test or large aperture interferometry is used, a reference optical element with an aperture equal to or larger than the measured ones is required.^{11}^{,}^{12} However, large reference optics is hard to manufacture and is very expensive. The basic principle of the subaperture stitching technique is synthesizing a full-aperture phase map from multiple subaperture phase shapes with motion capability and a stitching algorithm.^{13}14.15.^{–}^{16} So, the large aperture reference optics is omitted in subaperture stitching interferometry, because only a common small aperture interferometer is required. Meanwhile, the mid-frequency surface shape of the elements is also of concern in NIF.^{6}^{,}^{17} And higher spatial resolutions surface information could be obtained by subaperture stitching than by the other two methods. However, the efficiency of the subaperture stitching method is always criticized especially when the number of the subaperture is very large.

The ring polishing method is a primary method to polish large aperture planar optics and is widely used in the fabrication process in NIF and LMJ.^{17}^{,}^{18} The relative motion track is always rotationally symmetrical in ring polishing and the amount of material removal is related to the radius. Thus, the height errors are identical when the points are in the same radius,^{19}20.^{–}^{21} which can largely reduce the coverage area when the subaperture stitching method is used to measure the figure error of the surface. Some subapertures can be omitted from the fully-covered lattice, and some others are arranged as the sparse lattice to fit the rotationally symmetrical surface shape. It would be a time-saving method with a certain level of accuracy.

In this paper, the new stitching strategy called sparse subaperture stitching method is introduced first. The two chains lattice method is more accurate but less efficient than the one chain lattice method. So, the usage of different lattices could be adapted for the different steps of manufacturing. The accuracy should be the most concerned question when the number of the subapertures is reduced. After the algorithm is shown, the cumulative error is analyzed to certify that the accuracy is adapted for measuring. Finally, a planar optical element with $200\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 200\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ aperture polished by ring polishing is used for confirming the theory.

## 2.

## Principle of Sparse Subaperture Stitching and the Related Algorithm

## 2.1.

### Lattice of Sparse Subaperture Stitching

As mentioned in Sec. 1, the surface shape of large aperture planar optics polished by ring polishing is rotationally symmetrical. The lattice can be rearranged by utilizing the characteristics of rotational symmetry meanwhile, however, some rules should be obeyed. First, the lattice should satisfy the rotational symmetry to ensure that proper measuring results are obtained. In addition, the edge subapertures and the center subaperture are necessary to ensure coverage of surface minimum frequency-domain information. Overlapping areas between the subapertures must be available to calibrate the relative tilts and pistons.

Thus, the lattice of one chain structure could be as shown in Fig. 1(a). Measuring time could be largely saved with this lattice. However, in order to ensure the consistency of the two different directions, the lattice of double chains could also be adopted, as shown in Fig. 1(b). The two chains should be close to orthogonally oriented. This consistency detects the appearance of a singular surface and improves the accuracy of stitching results. In Figs. 1(a) and 1(b), the aperture of the subaperture is set for $\mathrm{\Phi}100\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ and the size of the square planar elements is set for $400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$.

Obviously, many more subapertures could be omitted when the scale of the measured surface is larger by the sparse subaperture method. For the example of the surface with the $400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ scale, the number of subapertures is 7 or 15 with one or two chains while the number is 49 with full covered subaperture measuring. Almost $6/7$ or $7/10$ of the time could be saved in the measuring process.

## 2.2.

### Principle of the Stitching Algorithm

An example of achieving high-precision stitching requirements is presented: Two subaperture surfaces, ${W}_{m}(x,y)$ and ${W}_{n}(x,y)$, have a common area but reference wavefront errors are overlying them. Thus, the errors ${W}_{r}(x,y)$ should be removed from ${W}_{m}(x,y)$ and ${W}_{n}(x,y)$ and the tilt and piston should be calibrated. The relation of the phase distributions in the overlapping area can be determined according to Eq. (1):

## (1)

$${W}_{m}(x-{x}_{m},y-{y}_{m})={W}_{n}(x-{x}_{n},y-{y}_{n})+{a}_{m-n,0}+{a}_{m-n,1}x+{a}_{m-n,2}y.$$Equation (1) indicates the ideal scenario of two subapertures in the overlapping area. However, to allow the phase distributions of these overlapping areas to approach each other, the global residual sum of the squared differences of all overlapping areas should be minimized, i.e.,

## (2)

$$f({a}_{\mathrm{0,0}},{a}_{\mathrm{0,1}},{a}_{\mathrm{0,2}},{a}_{\mathrm{1,0}},{a}_{\mathrm{1,1}},{a}_{\mathrm{1,2}},\dots ,{a}_{N,0},{a}_{N,1},{a}_{N,2})=\sum _{j=2}^{N}\sum _{i=0}^{j-1}\left\{\sum _{i\cap j}{[{W}_{i}^{*}(x,y)-{W}_{j}^{*}(x,y)]}^{2}\right\}\to \mathrm{min}.$$Differentiating Eq. (2), ($3N$) equations could be obtained and these may be expressed as a simple matrix $QR=P$, where $R$ refers to the sequence of the tilts, pistons, and reference errors of the Zernike coefficients. $R$ could be expressed as

## (3)

$$R={({a}_{0,0,}{a}_{0,1},{a}_{0,2},{a}_{1,0},{a}_{0,1},{a}_{0,2},\dots ,{a}_{\mathrm{Ns},0},{a}_{\mathrm{Ns},1},{a}_{\mathrm{Ns},2})}^{T},$$## (4)

$${q}_{ij}=\{\begin{array}{l}(2{\delta}_{mn}-1){(\sum {x}^{2})}_{mn},1\le m\le n\le {N}_{S},i=3m-2,j=3n-2\\ (2{\delta}_{mn}-1){(\sum xy)}_{mn},1\le m\le n\le {N}_{S},i=3m-2,j=3n-1\\ (2{\delta}_{mn}-1){\left(\sum x\right)}_{mn},1\le m\le n\le {N}_{S},i=3m-2,j=3n\\ (2{\delta}_{mn}-1){(\sum {y}^{2})}_{mn},1\le m\le n\le {N}_{S},i=3m-1,j=3n-1\\ (2{\delta}_{mn}-1){\left(\sum y\right)}_{mn},1\le m\le n\le {N}_{S},i=3m-1,j=3n\\ (2{\delta}_{mn}-1){({N}_{P})}_{mn},1\le m\le n\le {N}_{S},i=3m,j=3n\\ {q}_{ji},i>j\end{array},$$## (5)

$${p}_{i}=\{\begin{array}{l}\sum x({W}_{m*}-{W}_{m}),1\le m\le {N}_{S},i=3m-2\\ \sum y({W}_{m*}-{W}_{m}),1\le m\le {N}_{S},i=3m-1\\ \sum ({W}_{m*}-{W}_{m}),1\le m\le {N}_{S},i=3m\end{array}.$$Then using Doolittle decomposition to solve the $Q$ coefficients:

1. To $k=\mathrm{1,2},\cdots ,n$, do:

2. ${u}_{kj}={q}_{kj}-\sum _{s=1}^{k-1}{l}_{ks}{u}_{sj},j=k,k+1,\cdots ,n$,

3. ${l}_{ik}=({q}_{ik}-\sum _{s=1}^{k-1}{l}_{is}{u}_{sk})/{u}_{kk},i=k+1,\cdots ,n$.

After $Q$ matrix decomposition, the coefficients of $R$ may be obtained by the equation:

## (6)

$$\{\begin{array}{l}{y}_{i}={p}_{i}-\sum _{j=i}^{i-1}{l}_{ij}{y}_{j},i=\mathrm{1,2},\cdots ,n,\\ {r}_{i}=({y}_{i}-\sum _{j=i+1}^{n}{u}_{ij}{x}_{j})/{u}_{ii},i=n,n-1,\cdots ,1\end{array}.$$Then the tilts, pistons, and Zernike coefficients could be obtained using ${r}_{i}$ as a reference and Eq. (6). The wavefront of each subaperture ${W}_{m}^{*}(x,y)$ could be stitched directly thereafter.

## 3.

## Analysis of Precision

## 3.1.

### Uncertainty of the One-Dimensional Coefficient

The measurement errors of the subaperture stitching method can be classified into human errors, instrument errors, methodological errors, and environmental errors. When measuring the same surface several times under the same conditions, the results are likely to differ because of random errors. The measured value of a single instance is not static. However, random errors obey statistical rules if large amounts of measured values are observed and analyzed. Multiple error sources impact random errors following a normal distribution.^{22}

Errors can be calibrated through the standard deviation or variance of the relative tilt and shift. For simplicity, the one-dimensional (1-D) scenario is first studied. Consider a case where the common area of two tested 1-D areas $A$ and ${A}^{\prime}$ has $n$ pixels and the rms error of each measured phase distribution is $\sigma $. The tilt and shift are the two directions of the degrees of freedom in the single axis stitching test. When the two areas match perfectly, the relation between the phase distributions $z(x)$ and ${z}^{\prime}(x)$ in the overlapping area is

According to the calculation equation of the 1-D least-squares method, the tilt coefficient $a$ and shift coefficient $c$ can be expressed as

## (8)

$$a=\frac{n\sum x\mathrm{\Delta}(x)-\sum x\sum \mathrm{\Delta}(x)}{n\sum {x}^{2}-{(\sum x)}^{2}},$$## (9)

$$c=\frac{\sum {x}^{2}\sum \mathrm{\Delta}(x)-\sum x\sum x\mathrm{\Delta}(x)}{n\sum {x}^{2}-{(\sum x)}^{2}},$$## (10)

$${\sigma}_{a}^{2}=2\frac{{n}^{2}\sum {x}^{2}-n{\left(\sum x\right)}^{2}}{{[n\sum {x}^{2}-{\left(\sum x\right)}^{2}]}^{2}}{\sigma}^{2}=\frac{2n}{n\sum {x}^{2}-{\left(\sum x\right)}^{2}}{\sigma}^{2},$$## (11)

$${\sigma}_{c}^{2}=2\frac{n{(\sum {x}^{2})}^{2}-{(\sum x)}^{2}\sum {x}^{2}}{{[n\sum {x}^{2}-{(\sum x)}^{2}]}^{2}}{\sigma}^{2}=\frac{2\sum {x}^{2}}{n\sum {x}^{2}-{(\sum x)}^{2}}{\sigma}^{2},$$The variances of $c$ differ according to the zero point of $x$, namely taking the variance of point $A$ when the point is a zero point. By the calculations, the variance in the middle is the minimum, and its quadratic gradually increases to both sides. So, the variance of the overlapping edge is the largest.

If the zero point of $x$ is in the center position, the values of $x$ are $-(n-1)/2$, $-(n-1)/2+1,\dots ,(n-1)/2$ and

If the zero point of $x$ is in the edge position, the values of $x$ are 0, 1, $2,\dots ,N-1$. Then Eq. (11) becomes

## (14)

$${\sigma}_{c}^{\prime 2}=2\xb7[\frac{1}{6}n(n-1)(2n-1)]\xb7\frac{1}{\frac{{n}^{4}}{12}-\frac{{n}^{2}}{12}}{\sigma}^{2}=\frac{8n-4}{n(n+1)}{\sigma}^{2}.$$The point close to the edge increases because of the uncertainty of the calculated relative tilt. More specifically, the existence of ${\sigma}_{a}^{2}$ will directly increase the value of ${\sigma}_{c}^{2}$. Equation (14) can be obtained by the variances of tilt ${\sigma}_{a}^{2}$ and the middle point ${\sigma}_{c}^{2}$, which can be calculated by the equation

## (15)

$${\sigma}_{c}^{\prime 2}={\sigma}_{a}^{2}{(\mathrm{\Delta}x)}^{2}+{\sigma}_{c}^{2}=\frac{24{\sigma}^{2}}{{n}^{3}-n}\xb7{\left(\frac{n-1}{2}\right)}^{2}+\frac{2{\sigma}^{2}}{n}=\frac{8n-4}{n(n+1)}{\sigma}^{2}.$$Thus, when the distance to the center point and the variance of the test error in the overlapping area are known, the variance of the point can be calculated by Eq. (16):

## (16)

$${\sigma}_{c}^{2}(\mathrm{\Delta}x)=\frac{24{\sigma}^{2}}{{n}^{3}-n}\xb7{(\mathrm{\Delta}x)}^{2}+\frac{2{\sigma}^{2}}{n}.$$Suppose two tested areas $z(x)$ and ${z}^{\prime}(x)$ have a common area $O(x)$ and that 21 overlapping points may be found in this area. Through calculation by the least-squares method, the relative tilt $a$ is 0.1, the relative shift $b$ is 0, and the unbiased estimation variance ${\sigma}^{2}$ is 1. The fitting line $|z(x)-z\prime (x)|$ and standard deviation ranges of the overlapping area are shown in Fig. 2.

In Fig. 2, the black line indicates the result of 1-D linear fitting, and the red and blue lines, respectively, represent the upper- and lower-limit standard deviations of the fitting line $|z(x)-{z}^{\prime}(x)|$. The region between the red and blue lines is the range of the fitting standard deviation, which indicates that the probability of the actual stitching data of the overlapping area in ($|z(x)-{z}^{\prime}(x)|-{\sigma}_{c}^{2}(x)$, $|z(x)-{z}^{\prime}(x)|+{\sigma}_{c}^{2}(x)$) is 68.26%.

## 3.2.

### Analysis of One-Dimensional Cumulative Error

The previous section discussed the variance distribution in overlapping areas, and the variance was found to increase with increasing distance from the center. During actual stitching, the variance of tilt causes increases in the variance of the phase distribution outside the overlapping area and also presents a trend of the quadratic amplification.

As the number of sparse subapertures on the chains increases, the error also increases. This type of error is called cumulative error. In the 1-D linear chained stitching model, assuming that the number of the subapertures is $n$, the number of sampling points of each overlapping area is $c$, and the number of sampling points of a single subaperture is $m$, the error variance of the edge point of the end subaperture is

## (17)

$${\sigma}_{c}^{2}(x)=\sum _{i=1}^{n-1}[{\sigma}_{a}^{2}{(x-\frac{c}{2}-im+ic)}^{2}+{\sigma}_{c}^{2}]\phantom{\rule{0ex}{0ex}}={\sigma}_{a}^{2}(n-1)[{(x-\frac{c}{2})}^{2}-(m-c)n(x-\frac{c}{2})+\frac{1}{6}n(2n-1){(m-c)}^{2}]+(n-1){\sigma}_{c}^{2},$$The cumulative errors show an increasing trend as the number of subapertures increases. Thus, the number of subapertures should be kept low. However, overlapping areas decrease when the number of subapertures is low, which directly affects the stitching accuracy. Thus, the number of subapertures should be balanced to achieve high levels of accuracy.

## 3.3.

### Analysis of Two-Dimensional Cumulative Error

As the number of subapertures increases, the standard deviation gradually increases as a trend of the quadratic amplification. This amplification can cause increases in deviation between the real values of various parameters, such as the peak-to-valley (PV) value, rms value, Zernike polynomial coefficients, and the calculated results. Generally, the transverse distance of the PV points is high, and stitching must cross many subapertures, so the effect of the PV value is large. Thus, the subaperture in a diameter cannot be very much. Otherwise, the cumulative error will be too large and the PV value and the real value will deviate.

If the subaperture and the overlapping areas are rectangular, the variance of the tilt and shift should be

## (18)

$${\sigma}_{a}^{2}=\frac{24{\sigma}^{2}}{{c}_{y}({c}_{x}^{3}-{c}_{x})},\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}{\sigma}_{b}^{2}=\frac{24{\sigma}^{2}}{{c}_{x}({c}_{y}^{3}-{c}_{y})},\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}{\sigma}_{c}^{2}=\frac{2{\sigma}^{2}}{{c}_{x}{c}_{y}}.$$The term ${\sigma}_{c}^{2}$ is the variance at the center position of the overlapping area. The variance of the four points in the rectangular vertex angle is

## (19)

$${\sigma}_{c}^{\prime 2}={\sigma}_{c}^{2}+{\sigma}_{a}^{2}{\left(\frac{{c}_{x}-1}{2}\right)}^{2}+{\sigma}_{b}^{2}{\left(\frac{{c}_{y}-1}{2}\right)}^{2}\phantom{\rule{0ex}{0ex}}=\frac{2{\sigma}^{2}}{{c}_{x}{c}_{y}}+\frac{24{\sigma}^{2}}{{c}_{y}({c}_{x}^{3}-{c}_{x})}{\left(\frac{{c}_{x}-1}{2}\right)}^{2}+\frac{24{\sigma}^{2}}{{c}_{x}({c}_{y}^{3}-{c}_{y})}{\left(\frac{{c}_{y}-1}{2}\right)}^{2}\phantom{\rule{0ex}{0ex}}=\frac{2(7{c}_{x}{c}_{y}+{c}_{x}+{c}_{y}+5)}{{c}_{x}{c}_{y}({c}_{x}+1)({c}_{y}+1)}{\sigma}^{2}.$$Supposing the number of chain subapertures is $n$, the sampling point numbers of each two-dimensional (2-D) overlapping area are ${c}_{x}$ and ${c}_{y}$, and the sampling point numbers of single subapertures are ${m}_{x}$ and ${m}_{y}$. The sampling point numbers of the full stitching aperture are ${p}_{x}$ and ${p}_{y}$. As the distance from the center increases, the variance ${\sigma}_{c}^{2}$ also gradually increases. In sparse subaperture stitching, the variance of the cumulative errors is maximum between the edge points of the edge subaperture, and its value is

## (20)

$${\sigma}_{\mathrm{max}}^{2}(x)=\sum _{i=1}^{n-1}[{\sigma}_{a}^{2}{({p}_{x}-\frac{{c}_{x}}{2}-i{m}_{x}+i{c}_{x})}^{2}+{\sigma}_{b}^{2}{({p}_{y}-\frac{{c}_{y}}{2}-i{m}_{y}+i{c}_{y})}^{2}+{\sigma}_{c}^{2}].$$According to the geometrical relationship:

## (21)

$${c}_{x}=\frac{(2n-1){m}_{x}-{p}_{x}}{2n-2},\phantom{\rule{0ex}{0ex}}{c}_{y}=\frac{(2n-1){m}_{y}-{p}_{y}}{2n-2}.$$Equation (14) provides a quantitative description of the influence of cumulative error. The influence of the cumulative error can be calculated from the results obtained during simulated sparse subaperture stitching measurement.

Figure 4 shows the 2-D cumulative error when $n$ is 4. The cumulative error between the two red points is the maximum. Equation (11) provides a quantitative description of the influence of cumulative error. The influence of the cumulative error can be calculated from the results obtained during simulated sparse subaperture stitching measurement.

To explain the influence of the cumulative error, a simulation experiment is designed. A piece of plano-optics with an aperture of $400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 400\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ is simulated to evaluate the precision of the sparse subaperture stitching method. The size of the subaperture is $100\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 100\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$, and the sampling resolution is set to $0.2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}/\mathrm{pix}$ such that ${p}_{x}=2000$, ${p}_{y}=2000$, ${m}_{x}=500$, ${m}_{y}=500$. The subaperture number of a single chain is $n=3$. The parameters simulated above are common in actual measurements. A maximum cumulative error of ${\sigma}_{\mathrm{max}}^{2}=0.758{\sigma}^{2}$ is obtained through calculation. The cumulative error is very low, which shows that the $X$-type lattice can be used in subaperture stitching with adequate precision.

## 4.

## Experiment and Discussion

An experiment was designed to evaluate the precision of the sparse subaperture stitching method. A planar element with the aperture of $200\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 200\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ polished by ring polishing was chosen to be measured by the sparse and fully covered subaperture methods. The lattice of one chain with two orientations, two chains and fully covered is designed. The one chain lattice includes 4 subapertures while the full covered lattice includes 12 subapertures. The lattices of four situations are shown in Fig. 5. The shape figures of the four situations and the large aperture interferometry are shown in Fig. 6.

There is a slightly deviation for the rms values between sparse stitching and the fully covered stitching because the numerical distribution of the blank is deviated from the overall distribution. The rms value is closer to the real one with the increasing of chains. The PV value of the two chains lattice method is also closer than that of the one chain lattice because the astigmatism of the surface shape always exists in actual manufacturing, and the surface shape is not absolute symmetrical. However, there is not much difference between the values of the one chain method and the fully covered method. But the one chain method could save half of the time of the two chains lattice method. In actual manufacturing, the one chain lattice could be used for a fast inspection, and the two chains lattice could be used in concrete judging. In the final inspection, the fully covered lattice or large aperture interferometry should be used.

Although the number of subapertures is just eight or four less in this case, more subapertures will be omitted to measure in optics with a larger aperture.

The surface shapes of the fully covered stitching and the full aperture test with a large aperture interferometer are similar but still exist slightly deviation. The slight trail of the surface shape of full aperture stitching can be seen in Fig. 6(d). These situations may be caused by the reference surface error mingled in each subaperture. If higher precision stitching results are required, the reference surface error should be calibrated. Further analysis and experimentation on the influence of the reference surface error can be done in our future work. But the error has little influence on the sparse stitching method. The characteristic of sparse stitching is embodied in the comparison to the fully covered lattice stitching, and the superiority is shown by its high efficiency in the test process.

## 5.

## Conclusion

Planar large aperture optical surfaces manufactured by ring polishing usually have rotational symmetry; thus, they can be tested by incompletely covered subapertures. The method adopting the stitching strategy of the one-chain lattice and two-chain lattice could increase the measuring efficiency. The cumulative errors were deduced and simulating data were obtained to certify that the precision is adequate for sparse subaperture stitching.

The experiment results indicate that sparse stitching is suitable for measuring the large planar surface polished by ring polishing. The results of the two-chain lattice are closer to that of the fully covered lattice than those of the one-chain method. Meanwhile, the one-chain method could largely save time compared to the two-chain lattice. So, the usage of different lattices could be adapted for different periods of large aperture planar optical manufacturing.

In addition, the most important meaning is that the measurements of eight or four subapertures could be omitted, which could save time in the measuring process. More subapertures could be omitted when the relative size of a planar element is larger. This method has been used in the actual process of planar optical elements with a large aperture polished by ring polishing. The measuring results guide the manufacturing and save the manufacturing and measuring cycle time. Further analysis and experimentation on this subject are required.

## Acknowledgments

This work is sponsored by the National Natural Science Foundation of China (Grant Nos. 11427804, 61205124, and 61235011), National Science Instrument and Equipment Development Major Project of Ministry of Science and Technology of China (Grant No. 2012YQ04016403) and the Fundamental Research Funds for the Central Universities (Grant No. 1370219227).