Introduction {#Sec1}
============

The concept of integrated photonics was introduced in 1969 by Miller^[@CR1]^, which made a great contribution to the on-chip integrated optical interconnection. On-chip integrated optical interconnection is an emerging technique for large capacity data communications. The waveguide bends and crossings are one of the most critical components of optical interconnection. There are many waveguide bends and crossings designed by conventional approaches. A 5 μm radius 90-degree circular bend with 400 × 220 nm silicon strip waveguide was designed^[@CR2]^. Although the loss of the circular bend is 0.023 dB in FDTD simulation, the design method need a lot of complicated artificial adjustment experiments. A 90-degree partial Euler bends with an effective radius of 50 µm was fabricated on a silicon nitride photonic^[@CR3]^. The bend loss was within 0.2 dB, but it was only limited to the wavelength of 850 nm and had a large footprint. A 275 µm radius circular bend was designed using a silicon waveguide^[@CR4]^, which achieved very low loss. Unfortunately the length of the radius and waveguide is too long. Bogaerts et al.^[@CR5]^ proposed a 90-degree adiabatic bend for silicon waveguides. The total bend loss was 0.001 dB. However, the curve radius and circular wire section were up to 5 µm and 25 µm, respectively. Other waveguide bends with 2µm^[@CR6]^, 4µm^[@CR7]^, 10µm^[@CR8]^ radius were designed using a straight waveguide. Although the loss was below 0.3 dB loss, they had a large footprint and need a very long waveguide. Various studies have also been conducted on waveguide crossings. A low-loss waveguide crossing using the self-imaging properties of multimode interference (MMI) structures was designed^[@CR9]^, which achieved a loss of only 0.009 dB. But, the length is also fairly long (\> 140 µm in the case). A multimode-interference (MMI)-based crossing in high-index-contrast silicon wire waveguides was reported, which had an MMI crossing of \~ 0.4 dB insertion loss^[@CR10]^. However, its large 13 × 13µm^2^ footprint may be an issue for dense integration. Recently there have been reports attempting to further reduce the footprint with successful demonstrations with a 6 × 6µm^2^ ^[@CR11]--[@CR14]^. Note that − 1.7 dB insertion loss was achieved with only − 40 dB cross talk^[@CR11]^. Insertion loss was lower than 0.2 dB and crosstalk loss was below 40 dB in a broad bandwidth of 20 nm^[@CR12]^. Average insertion loss was 0.18 ± 0.03 dB and crosstalk was − 41 ± 2 dB^[@CR13]^. The waveguide crossing only operated at 1550 nm and 1310 nm, the results showed transmission insertion loss of − 0.028 ± 0.009 dB only for the 1550 nm device and − 0.017 ± 0.005 dB only for the 1310 nm device. Both crossings show crosstalk lower than − 37 dB^[@CR14]^. Another waveguide crossing showed an insertion loss of 1.82 dB as well as a crosstalk of \<  − 18 dB from 1,510 to 1600 nm^[@CR15]^. However, the footprint of the device was as compact as 21 × 21 μm^2^. A dual-mode waveguide crossing is proposed. The characterization results for the fabricated device show that low insertion loss below 1.5 dB and low crosstalk below − 18 dB can be achieved^[@CR16]^. But, it must need two waveguide crossing to achieve the function. Han et al.^[@CR17]^ used particle swarm optimization algorithm to design a SOI waveguide crossing with a footprint of only \~ 1 × 1 μm^2^. Although it has an insertion loss \< 0.175 dB and crosstalk around − 30 dB for the C-band, the designed method needs complex calculations and long-time simulation, and the number of iterations is set 500. A waveguide cross was designed that insertion loss and crosstalk were less than 0.6 dB and − 24 dB from 1530 to 1590 nm, respectively^[@CR18]^. However, the fabricated device still occupied a footprint of 4.8 × 4.8 μm^2^. Although traditional methods have been used to design devices with good performance, the footprint is always large that has a negative impact on dense integration, and the design of such devices rely heavily on the intuition and experience of the designer and cannot manually achieve a full-parameter space design. Therefore, we should use other methods to design waveguide bends and crossings. This method must be simple and feasible. It should greatly reduce the footprint of the device and achieve superior high performance.

In recent years, the inverse design method and DBS algorithm of ultra-compact silicon photonic devices have drawn more and more attentions^[@CR19]--[@CR25]^. The best-known institutions for inverse design research are the University of Utah and Stanford. The University of Utah achieved the inverse design of integrated nano-photonic polarization beam splitters by combining DBS algorithm and the current commercial software^[@CR19]^. The Stanford used alternating directions method of multipliers (ADMM) to design a wavelength demultiplexer of 1300 nm/1500 nm^[@CR20]^. For other optical device designs like flat optics and holograms, the inverse design also has a wide range of application^[@CR26]--[@CR29]^. Learning from their train of thought, Harbin Institute of Technology^[@CR30]--[@CR32]^, Huazhong University of Science and Technology^[@CR33]--[@CR35]^, National University of Defense Technology^[@CR36]--[@CR40]^ designed a variety of on-chip optical devices and achieved a good performance. One of the ways to realize inverse design is combining DBS algorithm and commercially available simulation software capable of implementing FDTD solutions, such as Lumerical FDTD Solutions^[@CR41]^. Inspired by them, we used inverse design and DBS algorithm to design waveguide bends and crossings.

In this paper, based on SOI platform, we successfully demonstrated a 90-degree waveguide bend, a waveguide crossing and a same direction waveguide crossing in the entire wavelength range of 1,400--1600 nm by using DBS algorithm. Their footprints are only 2.4 × 2.4 μm^2^, 2.4 × 2.4 μm^2^, 2.4 × 3.6μm^2^ respectively. The transmission efficiency of the 90-degree waveguide is \> 0.5 dB in the entire wavelength range of 1,440--1640 nm and up to 0.18 dB at some wavelengths. The waveguide crossing can provide insertion loss of less than 0.5 dB. The same direction waveguide crossing can provide the insertion loss of less than 0.56 dB. The measured crosstalk of the two crossings is lower than ‒19 dB and ‒21 dB respectively. Meanwhile, four kinds of optical interconnection modes are designed with the proposed bending and crossing, which realize the superior performance. Compared with traditional optical interconnection, the integration density and benefit various on-chip optical systems are greatly improved.

Results {#Sec2}
=======

Optical interconnection component {#Sec3}
---------------------------------

### A 90-degree waveguide bend {#Sec4}

Based on the "design method", we first design a 90-degree waveguide bend on a silicon-on-insulator (SOI) substrate with 220 nm silicon on top of 3 μm buried oxide. As shown in Fig. [1](#Fig1){ref-type="fig"}a--c, the computation domain is 2.4 × 2.4 μm^2^, where the optimization region is divided into 20 × 20 pixels and each pixel is a square of 120 nm × 120 nm with a central circular hole. The hole has a diameter of 90 nm and a depth of 220 nm, which can be fabricated with an electron-beam lithography (EBL) system (Vistec EBPG 5000 Plus).Figure 1Design and simulation results of the 90-degree waveguide bend. (**a**) Initial silicon slab before design optimization. (**b**) Final optimized structure. (**c**) Isometric view. (**d**) Simulated optical field distributions. (**e**) Simulated normalized transmission spectra. Here, we have plotted the electromagnetic energy density *U* = ϵ\|**E**\|^2^ + μ\|**H**\|^2^. An animated version of **d** is available in Supplementary Video [S1](#MOESM2){ref-type="media"}.

Each pixel can occupy two states: silicon or air. The white and red regions of the pixel denote the etched area where no silicon remains on top of the oxide and the unetched area where the 220 nm silicon remains, respectively. At the beginning, all the pixel states are chosen to be silicon, and the excitation the TE~0~ mode in the wavelength range of 1,400--1600 nm. Then the state of each pixel is altered one by one, and the figure-of-merit (FOM) is inspected. Here the FOM is defined as:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$ {\text{FOM}} = {\text{T}}_{{{\text{out}}\_{\text{port}}}} $$\end{document}$$where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\text{T}}_{{{\text{out}}\_{\text{port}}}}$$\end{document}$ is the transmission efficiency of the final structure at the wavelength of 1550 nm as the input light source is TE~0~ mode, which is calculated by 3D FDTD with 40 nm × 40 nm × 40 nm grid size and 8 layers PML. If the FOM is improved, the pixel state is retained. If not, the pixel state is reversed and the algorithm proceeds to the next pixel. One iteration ends after all the pixel states are inspected. Then the iterations continue until the FOM does not improve further.

We perform the device design and optimization using an eight-core desktop. For the 90-degree waveguide bend, it takes \~ 24 h on average to get the final results after 4 iterations. Given that fabrication variations can have a significant impact on real-world device performance, we characterize the same bends with a ± 10% variation in pillar diameter by simulation. The simulated optical field profile of the device when the continuous wave is launched from input port is shown in Fig. [1](#Fig1){ref-type="fig"}d. Figure [1](#Fig1){ref-type="fig"}e shows the calculated normalized transmission efficiency. The simulated normalized transmission efficiency of the bend is \< 0.5 dB in the entire range of 1,400--1600 nm and up to 0.18 dB at some wavelength. As we can show, a ± 10% variation in pillar diameter doesn\'t affect largely on the 90-degree waveguide bend. It is seen from the Table [1](#Tab1){ref-type="table"} that the waveguide bend we designed has very low loss losses and a relatively small footprint compared with other waveguide bends. The method we used is feasible and auto-actuated.Table 1Comparison of the waveguide bends.ReferencesBending area (μm^2^)Propagation losses (dB)WaveguideMethodNoteThis work2.4 × 2.40.18Etched siliconDBSAuto-actuated^[@CR2]^5 × 50.023\~ 10 μm long waveguideVariational optimizationHand-actuated^[@CR3]^50 × 500.2\~ 100 μm long waveguideNumerical analysisHand-actuated^[@CR4]^275 × 275\~ 0.1\~ 550 μm long waveguideBeam-propagationHand-actuated^[@CR5]^5 × 50.02725 μm long waveguideManual adjustmentHand-actuated^[@CR6]^2 × 20.08621 mm long waveguideManual adjustmentHand-actuated^[@CR7]^4 × 40.5\~ 8 μm long waveguideClothoid curvesHand-actuated^[@CR8]^10 × 100.02\~ 20 μm long waveguideMicron-scale coresHand-actuated

### Waveguide crossing {#Sec5}

A waveguide crossing The on-chip optical interconnection inevitably needs to address the waveguide cross-connect problem. Herein, we design a waveguide crossing with a footprint of only 2.4 × 2.4 μm^2^ using the same process as the bend, which has a quadrature symmetry. The structural schematic diagram of the device is described in Fig. [2](#Fig2){ref-type="fig"}a,b, the crossing has two input ports labeled as I1 and I2, and two output ports labeled as O1‒O2, respectively. The widths of the two input waveguides and the output waveguide are 480 nm.Figure 2Design and simulation results of the waveguide crossing. (**a**) Initial silicon slab before design optimization. (**b**) Final optimized structure. (**c**), (**d**) Simulated optical field distributions. (**e**), (**f**) Simulated normalized transmission spectra. Here, we have plotted the electromagnetic energy density *U* = ϵ\|**E**\|^2^ + μ\|**H**\|^2^. An animated version of (**c**) is available in Supplementary Video [S2](#MOESM3){ref-type="media"}.

This design region is also divided into 20 × 20 discrete pixels with minimum feature size of 90 nm. In the optimizing process, we use the TE~0~ mode in the wavelength range of 1,400--1600 nm as the only one light source placed at I1 port. The FOM of the crossing for inverse design is defined as: $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$ {\text{FOM}} = \left[ {{\text{T}}_{1} \left( {{\text{i}},{\text{j}} + 1} \right) > {\text{T}}_{1} \left( {{\text{i}},{\text{j}}} \right)} \right] \cap \left[ {{\text{T}}_{2} \left( {{\text{i}},{\text{j}} + 1} \right) \le {\text{T}}_{2} \left( {{\text{i}},{\text{j}}} \right)} \right] $$\end{document}$$where T~1~(i,j) and T~1~(i,j + 1) represent the average transmission of O1 port in the ith iterate of jth and (j + 1)th pixel, respectively. T~2~(i,j) and T~2~(i,j + 1) represent the average transmission of O2 port in the ith iterate of jth and (j + 1)th pixel, respectively. All the transmissions are calculated by 3D FDTD with 40 nm × 40 nm × 40 nm grid size and 8 layers PML. Actually, the first term in the right of Eq. ([2](#Equ2){ref-type=""}) is used to optimize average insertion loss (IL) and the second one is employed to minimize the crosstalk (CT). If FOM = 1, it means IL is improved and CT descend, the pixel state is retained. If not, the pixel state is reversed.

It takes about 36 h after 4 iterations to get the optimized crossing on an eight-core desktop. To investigate the manufacturing error of the crossing, we test the same crossing with a ± 10% variation in pillar diameter. Figure [2](#Fig2){ref-type="fig"}c,d shows the simulated optical. The simulated ILs and CTs spectra at O1and O2 ports are plotted in Fig. [2](#Fig2){ref-type="fig"}e,f when the continuous wave is launched from input ports I1 and I2, respectively. The simulated CTs for the crossing and the same crossings with ± 10% diameter are \< − 19 dB, \< −  18 dB, \< −  19 dB, respectively. The simulated ILs for all devices are \~ 0.5 dB from 1,400 to 1,600 nm wavelength range. A ± 10% variation in pillar diameter had nearly negligible impact on the IL. Due to a symmetric structure, the results obtained as the light source placed at I1 and I2 ports are consistent.

A same direction waveguide crossingHere, we design another same direction waveguide crossing. As shown in Fig. [3](#Fig3){ref-type="fig"}a,b, the inverse design region composed of 20 × 30 discrete pixels occupies a compact footprint of only 2.4 × 3.6 μm^2^. The crossing is optimized using the same process as the waveguide crossing. The FOM is also defined as Eq. ([2](#Equ2){ref-type=""}). It takes about 48 h after 5 iterations to get the final optimized structure on an eight-core desktop. Figure [3](#Fig3){ref-type="fig"}c,d shows the simulated optical and Fig. [3](#Fig3){ref-type="fig"}e,f shows the calculated ILs and CTs spectra when the continuous wave is launched from I1 and I2 ports, respectively. For the final crossing, the average IL is \~ 0.56 dB and the CT is \<  − 21 dB within the wavelength range from 1,400‒1600 nm. Even though the diameter of the crossing changes with a ± 10% variation, the IL and CT curves are still quite consistent. Table [2](#Tab2){ref-type="table"} summarizes results from some existing literature. It is seen that the waveguide crossings we designed have a small footprint. In contrast with other waveguide crossings, the measured insertion loss and crosstalk are relatively low.Figure 3Design and simulation results of the same direction waveguide crossing. (**a**) Initial silicon slab before design optimization. (**b**) Final optimized structure. (**c**), (**d**) Simulated optical field distributions. (**e**), (**f**) Simulated normalized transmission spectra. Here, we have plotted the electromagnetic energy density *U* = ϵ\|**E**\|^2^ + μ\|**H**\|^2^. An animated version of (**c**) is available in Supplementary Video [S3](#MOESM4){ref-type="media"}. Table 2Comparison of the waveguide crossings.ReferencesFootprint (μm^2^)IL (dB)CT (dB)MethodNoteThis work2.4 × 2.40.5− 19DBSAuto-actuated2.4 × 3.60.56− 21DBSAuto-actuated^[@CR9]^175.2 × 175.20.009\< − 37Self-imagingLarge footprint^[@CR10]^13 × 130.4− 30Manual adjustmentMultiple experiments^[@CR11]^6 × 61.7− 40Parabolic taperHand-actuated^[@CR12]^6 × 60.2− 40GA20 nm bandwidth^[@CR13]^6 × 60.18− 41PSO50 generations^[@CR14]^6 × 60.028, 0.017− 37Step by stepOnly operating at 1550 nm and 1310 nm^[@CR15]^21 × 211.82− 18Y-junctionLarge footprint^[@CR16]^2.4 × 2.41.5− 18Self-imagingRequiring two crossings^[@CR17]^1 × 10.175− 30PSO500 generations^[@CR18]^4.8 × 4.80.6− 24SymmetricStill large

Optical interconnection integrated structures {#Sec6}
---------------------------------------------

In principle, the optical signal can be delivered to any chip locations via any routing path provided the low loss bending and crossing are available. For a proof-of-concept demonstration, we arbitrarily design four optical interconnection integrated structures with the proposed waveguide bend and crossings.

The first structure consists of four waveguide crossings, which is named "2_2\_2_2" structure. The schematic diagram of this structure is shown in Fig. [4](#Fig4){ref-type="fig"}a, its footprint is 2.4 × 5.3 μm^2^, the span between the four crossings is 500 nm, and the input waveguides and the output waveguide are 480 nm. We define input ports I1--I4 and output ports O1‒O4, respectively. Figure [4](#Fig4){ref-type="fig"}b--e shows the simulated optical field when the continuous wave is launched from I1-I4 ports, respectively. The match IL and CT spectra are plotted in Fig. [4](#Fig4){ref-type="fig"}f--i. The result curves in Fig. [4](#Fig4){ref-type="fig"}f,h, and g,i are consistent due to the symmetric structure. The average IL is \~ 1 dB and the CTs are \<  − 18 dB within the wavelength range from 1,400 to 1600 nm.Figure 4Design and simulation results of the "2_2\_2_2" structure. (**a**) Schematic diagram of the structure. (**b**)--(**e**) Simulated optical field distributions. (**f**)--(**i**) Simulated normalized transmission spectra.

The second structure with a footprint of 8.2 × 8.2 μm^2^ consists of nine waveguide crossings, which is named "3_3\_3_3" structure. The relevant specific parameters of this structure are shown in Fig. [5](#Fig5){ref-type="fig"}a. We define input ports I1--I6 and output ports O1‒O6, respectively. Figure [5](#Fig5){ref-type="fig"}b--d and h--j show the simulated optical field when the continuous wave is launched from I1-I6 ports, respectively. The match IL and CT spectra are plotted in Fig. [5](#Fig5){ref-type="fig"}e--g and k--m. The result curves in Fig. [5](#Fig5){ref-type="fig"}e and k, f and l, g and m are consistent due to the symmetric structure. The average IL for all output ports is 2.5 dB from 1,400 to 1,600 nm. The CTs are measured to be \<  − 16 dB within the same wavelength range.Figure 5Design and simulation results of the "3_3\_3_3" structure. (**a**) Schematic diagram of the structure. (**b**)--(**d**) (**i**)--(**k**) Simulated optical field distributions. (**e**)--(**g**) (**l**)--(**n**) Simulated normalized transmission spectra.

The third structure with a footprint of 7.7 × 6.14 μm^2^ consists of three same direction waveguide crossings, which is named "4_4" structure. The labels of input ports and output ports and other relevant specific parameters of this structure are shown in Fig. [6](#Fig6){ref-type="fig"}a. Figure [6](#Fig6){ref-type="fig"}b--e shows the simulated optical field when the continuous wave is launched from I1-I4 ports, respectively. The match IL and CT spectra are plotted in Fig. [6](#Fig6){ref-type="fig"}f--i. The result curves in Fig. [6](#Fig6){ref-type="fig"}f and i, g and h are consistent due to the symmetric structure. When the input waves from I1 or I4 input ports, the losses come from the one times cross connection, and the average IL is \~ 0.6 dB and the CTs are \<  − 20 dB within the wavelength range from 1,400 to 1600 nm. When the input waves from I2 or I3 input ports, the losses come from the two times cross connection, and the average IL is \~ 1.2 dB from 1,400 to 1,600 nm, the CTs are measured to be \<  − 16 dB within the same wavelength range.Figure 6Design and simulation results of the "4_4" structure. (**a**) Schematic diagram of the structure. (**b**)--(**e**) Simulated optical field distributions. (**f**)--(**i**) Simulated normalized transmission spectra.

The fourth structure with a footprint of 10.12 × 10.12 μm^2^ consists of three waveguide bends and three waveguide crossing, which is named "3_3" structure. The relevant specific parameters of this structure are shown in Fig. [7](#Fig7){ref-type="fig"}a. We define input ports I1--I3 and output ports O1‒O3, respectively. Figure [7](#Fig7){ref-type="fig"}b--d shows the simulated optical field when the continuous wave is launched from I1-I3 ports, respectively. The match IL and CT spectra are plotted in Fig. [7](#Fig7){ref-type="fig"}e--g. The average ILs for O1, O2 and O3 output ports are \~ 2.5 dB, \~ 2 dB and \~ 3 dB from 1,500 to 1,580 nm, respectively. The CTs are for O1, O2 and O3 output ports measured to be \<  − 27 dB, \<  − 25 dB and \<  − 21 dB within the same wavelength range, respectively.Figure 7Design and simulation results of the "3_3" structure. (**a**) Schematic diagram of the structure. (**b**)--(**d**) Simulated optical field distributions. (**e**)--(**g**) Simulated normalized transmission spectra.

The above is a few examples of optical interconnect structures using the waveguide bend and crossings. It is obviously that the number of devices contributes to the losses of the structures. The more devices we use, the greater the losses. The minimum intervals between all devices is 500 nm making sure there is not any out-of-plane light coupling.

Discussion {#Sec7}
==========

In this paper, we use inverse design concept and DBS algorithm to design ultra-compact high efficiency and low crosstalk waveguide bend and crossings for optical interconnection based on SOI platform. The footprints of them are only 2.4 × 2.4 μm^2^, 2.4 × 2.4 μm^2^ and 2.4 × 3.6 μm^2^ respectively. The waveguide bend realizes the 90 degree bending of light, up to 0.18 dB transmission efficiency. The ILs and CTs of the waveguide crossings are less than 0.5 dB, 0.56 dB, − 19 dB and − 21 dB, respectively. Furthermore, we arbitrarily design four optical interconnection integrated structures with the proposed waveguide bend and crossings. All of them exhibit high performances. Actually, more optical interconnect structures can be achieved to deliver optical signal using them. The ultra-compact structures will significantly improve the integration density and benefit various on-chip optical systems. Meantime, our design thought and algorithm flow can also be applied to the research and design of other on-chip photonic devices, such as the optical power beam splitter, the mode division multiplexer, the wavelength division multiplexer, etc.

Design method {#Sec8}
-------------

The overall design idea of inverse design: Firstly, setting the target performance of the devices. Then, according to the set performance requirements, the devices are designed and optimized by using various optimization algorithms. Here, we use DBS algorithm to design the devices. The general process is as follows:Dividing pixels: The optimized region is divided into a certain number of pixels. The geometry of each pixel can be circle, square, or other complex shapes. The minimum feature size can be determined according to the fabrication capability. In this paper, each pixel is a square of 120 nm × 120 nm with a central circular hole. The hole has a diameter of 90 nm.Establishing the initial pixel states: Each pixel can occupy two states: silicon or air. For the initial pattern of the devices, all the pixel states are chosen to be silicon in this paper.Defining figure-of-merit (FOM): the FOM is used to see if the pixel state is witched. Different FOMs are defined according to different devices. In this paper, the FOM is defined as Eqs. ([1](#Equ1){ref-type=""}) or ([2](#Equ2){ref-type=""}).Witching the pixel state: the state of each pixel is altered one by one, and the FOM is inspected. The pixel state is retained if the FOM is improved. If not, the pixel state is reversed and the algorithm proceeds to the next pixel. One iteration ends after all the pixel states are inspected.Iterative refinement: The iterations continue until the FOM does not improve further.

Supplementary information
=========================

 {#Sec9}

Supplementary Information. Supplementary Video S1. Supplementary Video S2. Supplementary Video S3.

**Publisher\'s note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

These authors contributed equally: Zikang Li and Guofeng Li.

Supplementary information
=========================

is available for this paper at 10.1038/s41598-020-68936-w.

This work was supported by the National Natural Science Foundation of China (60907003, 61805278, 61661004), the Guangxi Science Foundation (2017GXNSFAA198227), the Foundation of NUDT (JC13-02-13, ZK17-03-01), the Hunan Provincial Natural Science Foundation of China (13JJ3001), and Program for New Century Excellent Talents in University (NCET) (NCET-12-0142).

Z.L. and G.L. contributed equally to this work. Z.L. and G.L. finished the algorithm and English writing of the paper together; Z.L. designed and performed the analysis, and prepared the figures/tables; J.Y., J.H., C.Y., Y.Q. and W.X. contributed to manuscript preparation; J.Y. and J.H. conceptualized and coordinated the study; J.Y. and Z.Z. supervised the project. H.H. participated in the revision of the paper.

All data generated or analyzed during this study are included in this article.

The authors declare no competing interests.
