Dynamically encircling exceptional points in a three-mode waveguide system

Dynamically encircling exceptional points (EPs) in non-Hermitian systems has attracted considerable attention recently, but all previous studies focused on two-state systems, and the dynamics in more complex multi-state systems is yet to be investigated. Here we consider a three-mode non-Hermitian waveguide system possessing two EPs, and study the dynamical encircling of each single EP and both EPs, the latter of which is equivalent to the dynamical encircling of a third-order EP that has a cube-root behavior of eigenvalue perturbations. We find that the dynamics depends on the location of the starting point of the loop, instead of the order of the EP encircled. Compared with two-state systems, the dynamical processes in multi-state systems exhibit more non-adiabatic transitions owing to the more complex topological structures of energy surfaces. Our findings enrich the understanding of the physics of multi-state non-Hermitian systems and may lead to the design of new wave manipulation schemes. Exceptional points are singularities in the parameter space of a system where, under certain conditions, gain and loss can coalesce and give rise to exotic behaviour. The authors report on the dynamical encircling of exceptional points in a waveguide system having more than two states, increasing our understanding of the chiral phenomena occurring in non-Hermitian systems.

I n non-Hermitian systems, both the eigenvalues and eigenvectors can coalesce at a degeneracy point known as the exceptional point (EP) [1][2][3][4] , which exhibits interesting properties and has given rise to a variety of applications [5][6][7][8][9][10][11][12] . Recently, considerable attention has been paid to the topological structure of the energy surface (also known as Riemann surface) around the EP, which is typically investigated by designing a loop in the parameter space to encircle the EP 13 . In a two-state system, eigenmode exchange occurs when the EP is encircled stroboscopically (i.e., assuming the process is adiabatic) [14][15][16] . Stroboscopic encircling of a high-order EP can reveal the more complex topological structure of multi-state systems 17,18 . Unexpected physics and phenomena have also been revealed when EPs are encircled stroboscopically along non-homotopic loops in multistate systems 19 .
In contrast to the stroboscopic encircling of EPs which assumes an adiabatic process, dynamically encircling EPs in non-Hermitian systems exhibits more complex dynamics since adiabaticity may break down due to the presence of non-Hermiticity [20][21][22][23][24][25][26][27][28][29][30] . One of the first two experiments on dynamical encircling of an EP was performed in an optomechanical device with the Hamiltonian evolving in real time 7 . Another experiment used a microwave waveguide with specially designed boundaries that are used to mimic the evolution of Hamiltonian in space 24 . A chiral transmission behavior is found in the sense that different encircling directions result in different output states. The starting point of the loop lies near the region where two eigenmodes (i.e., symmetric and anti-symmetric modes) share the same loss and they can be used for asymmetric mode switching. When the starting point moves far away from such region, the eigenmodes are symmetry-broken modes, on the contrary, and the dynamics is found to be non-chiral 27 . Recent study also demonstrated that the EP can be dynamically encircled in antiparity-time (PT) symmetric systems and the dynamics is chiral with the starting point in the PT-broken phase 29 . As a result, symmetry-broken modes can also be used for asymmetric mode switching as long as the system is anti-PT-symmetric.
Although the dynamical encircling of EPs has drawn much attention, all the previous works consider systems consisting of only two states and the EP is a second-order EP. A natural question to ask is what the dynamics would be when an EP is dynamically encircled in multi-state systems. And what if the EP is a high-order EP? Multi-state systems possess more complex topological structures of energy surface so that the dynamical processes should be more complicated. Does the starting point still have an important role to determine the dynamics? These questions remain open.
In this paper, we address the above questions by studying a three-waveguide non-Hermitian system which possesses two EPs. We design a configuration to dynamically encircle these EPs. The field distributions and the amplitudes of instantaneous eigenmodes along the waveguiding direction are calculated to help understand the dynamics in the encircling process. In the following, we first study the energy surface of the proposed system in a two-parameter phase space. Then we investigate the dynamics when each EP is encircled individually to clarify the condition for the dynamics to be chiral or non-chiral. In particular, we investigate the dynamics when the loop encloses both EPs, and we will see that dynamically encircling two EPs in our system is equivalent to the dynamical encircling of a third-order EP that has a cube-root behavior of eigenvalue perturbations.

Results
System under investigation. We consider a non-Hermitian system consisting of three waveguides with the cross section shown in Fig. 1a. Waveguide-2 and waveguide-3 have the same width W, and that of waveguide-1 is αW with α being a detuning parameter. The three waveguides have the same height H. Symbols g and D represent the gap distance between the corresponding waveguides. Two absorbers (see the yellow region) are attached to waveguide-2 and waveguide-3 to introduce loss, with their height denoted by h and width by w 2 and w 3 , respectively. Throughout this paper, we fix the parameters W = 8 mm, H = 4 mm, D = 3.2 mm, g = 1 mm, w 3 = 5 mm, and h = 1 mm. The working frequency is 10 GHz. The waveguides are made of yttrium iron garnet (YIG) with permittivity of 15.26, and that of the microwave absorber is assumed to be 4 + 20i. For simplicity, the substrate and the background are assumed to have the same refractive index as vacuum.
We first used COMSOL to calculate the effective mode index n eff , i.e., the eigenvalue of the eigenmodes, as a function of α and w 2 . Figure 1b, c plots the real part and imaginary part of the eigenvalues, respectively. The system supports three eigenmodes and their complex eigenvalues form three Riemann sheets distinguished by the red, blue, and gray color. There are two EPs labeled as EP-1 and EP-2 in the parameter space. They are second-order EPs, since EP-1 is the coalescence of the red and blue sheet, while EP-2 is the coalescence of the red and gray sheet (see Fig. 1c). For each EP, we use the term RPC (short for "real part coalescence") to represent the coalescence of Re[n eff ] between the adjacent Riemann sheets. The two white dashed lines in Fig. 1b represent RPC-1 and RPC-2, respectively. For the sake of the following analysis, we also use the term IPC (short for "imaginary part coalescence") to represent the coalescence of Im [n eff ] for each EP (see the two white dashed lines in Fig. 1c). Figure 1d shows the parameter space with the two EPs marked by stars. It is known that in a two-state system, the dynamical encircling of an EP exhibits a chiral transmission behavior when the starting point of the loop lies on the IPC where the two eigenmodes have the same loss 24 . As our system is a three-state system with two IPCs and a more complex energy surface, it is a good candidate to study whether the chiral behavior still exists in multi-state systems.
We consider two loops in the parameter space using the following formula: where L is the system length, and z is ranging from −L/2 to L/2. The other parameters (α 0 , Δα and Δw) are used to tune the loop, and we choose α 0 = 1.09, Δα = 0.1, and Δw = 4 mm for loop-1.
For loop-2, we choose α 0 = 0.937, Δα = 0.1, and Δw = 5.5 mm. With these choices, the starting point of the loop (i.e., z = −L/2 or L/2) lie exactly on the IPCs. A schematic of the system with the above parameters is shown in Fig. 2a (top view) and Fig. 2b (side view), where waveguide-3 is set to be transparent from the top view for better visibility. The variations of α and w 2 as a function of z are plotted in Fig. 2c, d, respectively. The transmission of electromagnetic waves through the system is equivalent to eigenstate evolving along the designed loop in the parameter space. Input waves from the left-hand side lead to counterclockwise loops while injections from right-hand side correspond to clockwise loops (see Fig. 2a-d). In potential experiment, the waveguides with irregular shapes can be fabricated by polishing commercially available YIG stripes using a hand polishing machine 27 .
Single EP encirclement starting from IPC. In this section, we will consider loop-1 and loop-2, both of which correspond to the dynamical encircling of a single EP. We first use loop-1 to study the dynamical encircling of EP-1 with the starting point on IPC-1. The length of the system is set to be 800 mm. Since we can consider both counter-clockwise and clockwise loops, and for each loop the input can be either one of the three eigenmodes, there are a total of six cases to be investigated. Throughout this paper, we sort the three eigenmodes at any point of the parameter space based on their real part of the eigenvalues. An eigenmode is called mode-1 when it has the largest Re[n eff ], while mode-3 has the smallest Re[n eff ]. We first study counter-clockwise loops. Figure 3a-c plots the power flow distributions along the waveguiding direction with mode-1, mode-2, and mode-3 as the input eigenmode, of which the energy is predominantly localized in  a Top view of the three-waveguide system with two parameters (i.e., the detuning parameter α and the absorber width w 2 ) varying continuously along the z-axis. Waveguide-3 is set to be transparent for better readability. Symbols W, g, and L denote, respectively, the width, the gap distance, and the length of the waveguides. b Side view of the system, where h and w 3 represent the height and width of the absorber attached to waveguide-3, and the three waveguides have the same height H. c The variation of α along z-axis for loop-1 and loop-2. d The variation of w 2 along z-axis for loop-1 and loop-2. Input from the left-and right-hand side of the system corresponds to counter-clockwise and clockwise loops, respectively d e 1 2 . a-c Numerically simulated power flow distributions along the z-axis for counter-clockwise loop-1 with a mode-1, b mode-2, and c mode-3 being the input eigenmode. The field distribution in waveguide-3 is plotted between that in waveguide-1 and waveguide-2, since the field strength outside the waveguides is very weak. d-f Calculated amplitudes of the three instantaneous eigenmodes along the z-axis for counter-clockwise loop-1 with d mode-1, e mode-2, and f mode-3 being the input eigenmode, where c r , c b , and c g represent, respectively, the coefficient of the eigenmode on the red, blue, and gray sheet. In each figure, the two insets show the power flow distributions at the input (left) and output (right) interface. g-i Trajectories of the state evolution on the real part of Riemann sheets (i.e., the real part of the effective mode index n eff as a function of the detuning parameter α and the absorber width w 2 ) for counter-clockwise loop-1 with g mode-1, h mode-2, and i mode-3 being the input eigenmode. j-l Same as (g-i) but on the imaginary part of Riemann sheets. The output state is always mode-2 for counter-clockwise loop-1. The term "NAT" denotes the non-adiabatic transition waveguide-1, waveguide-2, and waveguide-3, respectively. In our simulations, the desired eigenmodes at the input side can be excited efficiently by pre-calculating their eigenfields and then using them as boundary conditions (see Methods section). In the experiment, injecting signals from either waveguide will excite all the three eigenmodes. The desired eigenmode will still dominate the input state since its eigenfields are more similar to the injected signal and hence will couple predominantly with the injected signal. We plot the field distribution of waveguide-3 between that of waveguide-1 and waveguide-2, since the field strength outside the waveguides is very weak. The fields are captured in a x-z plane which is the central plane of each waveguide. As the power flow decays along the z-axis, the fields are normalized to have the same maximum intensity at each cross section in order to improve the readability. We find that regardless of which eigenmode is injected into the system, the output mode is always mode-2 with the power flow mainly localized in waveguide-2.
To investigate the details in the wave transmission process, we consider the instantaneous eigenfields obtained from COMSOL simulations at each cross section of the system (i.e., each z position with a particular configuration in the parameter space of α and w 2 ). The instantaneous electric/magnetic fields at each point can be written as a sum of the instantaneous eigenfields corresponding to the eigenmodes of the system configuration at that point, e.g., are the instantaneous electric eigenfields of eigenmode on the red, blue, and gray sheet, respectively, and c r , c b , and c g are the corresponding amplitude coefficients. These eigenfields are in fact right eigenvectors of the non-Hermitian system so that we use a superscript "R". We can use them to construct the corresponding left eigenvectors [e.g., E L r z ð Þ, E L b z ð Þ, and E L g z ð Þ] and then the amplitude coefficients can be calculated by projecting the instantaneous fields onto the left eigenfields (see Methods section for details). The results with mode-1, mode-2, and mode-3 being the input mode are shown in Fig. 3d-f, respectively, where c r , c b , and c g can represent, respectively, the proportion of the eigenmode on the red, blue, and gray sheet. We note that at the starting point, the coefficients for the non-excited states (e.g., c b and c g in Fig. 3d) are not zero. They are mainly due to reflections since the system has a finite length, but it will not affect the physics studied in this work. Using the extracted amplitude coefficients, we can draw the trajectory of the mode evolution on the Riemann sheets. The trajectory is marked on the red/blue/gray sheet when c r /c b /c g dominates the instantaneous state at each point in parameter space. Figure 3g-i plots the trajectory on the real part of Riemann sheets, and Fig. 3j-l shows that on the imaginary part of Riemann sheets. In the following, we will analyse the dynamics based on the results in Fig. 3. Figure 3d shows that when mode-1 (mainly localized in waveguide-1) is injected, the process is adiabatic. The end state is mode-2, which corresponds to a state flip after encircling EP-1. Non-adiabatic transitions (NATs) may occur when the EP is dynamically encircled in non-Hermitian systems because of the presence of gain and loss 22 . But in the considered process, the state always evolves on the red sheet which has the lowest loss (see the trajectory in Fig. 3g, j). Therefore, the state evolution process is stable. The profiles of the power flow at the input and output interface are shown in the inset of Fig. 3d. The variation in the eigenmode is evident, which is due to the dynamical encircling of EP-1.
When we inject mode-2 (mainly localized in waveguide-2) into the system, the state at first propagates on the blue sheet (see Fig. 3e, h, k), which is in fact not the Riemann sheet with the lowest loss. The mode evolution is not stable and after some time, a NAT occurs and the state jumps to the red sheet which exhibits the lowest loss (see Fig. 3h, k). After the NAT, the state evolution follows the adiabatic trajectory and the final state is mode-2, indicating that the state returns to itself after the loop because of the NAT. From the power flow distributions shown in Fig. 3b, this NAT is characterized by a transfer of the power flow from waveguide-2 to waveguide-1 (marked by the yellow dashed line). The crossing between the red curve and the blue dashed curve in Fig. 3e is also an evidence of the NAT.
The dynamical process of the case with mode-3 (mainly localized in waveguide-3) as the input is even more complex. We find that there are a total of two NATs, since the state at first evolves on the gray sheet which has the highest loss (see Fig. 3l). The first NAT corresponds to the state jump from the gray sheet to the blue sheet, while the second NAT corresponds to that from the blue sheet to the red sheet (see Fig. 3i, l). Once the state reaches the red sheet, the following evolution process is adiabatic and the final state is still mode-2. These two NATs can also be observed from the field distributions (see the two yellow dashed lines in Fig. 3c) and the amplitude coefficients (see the two crossings in Fig. 3f).
The above results indicate that the output of counter-clockwise loops is always mode-2. We can use the same method to study clockwise loops in which the signal is injected from the righthand side of the system. Figure 4a-c shows the calculated amplitudes of the instantaneous eigenmodes with mode-1, mode-2, and mode-3 being the input, respectively. The corresponding trajectories on the real part of Riemann sheets are plotted in Fig. 4d-f. We find that the output of clockwise loops is mode-1 for the three cases. Although the final state is the same, the dynamical behaviors in the encircling process are quite different. In brief, the encircling process is adiabatic when mode-2 is injected and the energy is gradually transferred from waveguide-2 to waveguide-1 (see Fig. 4b, e). When the input mode is mode-1 (see Fig. 4a, d) or mode-3 (see Fig. 4c, f), however, there are NATs since at the beginning the state does not evolve on the red sheet which has the lowest loss. The corresponding trajectories on the imaginary part of Riemann sheets are given in Supplementary  Fig. 1.
The chiral transmission behavior is now evident for loop-1 with the starting point on IPC-1, i.e., the final state of counterclockwise loops is always mode-2 while that of clockwise loops is mode-1. This is an extension of the chiral dynamics from twostate systems 24 to three-state systems, which implies that the chiral transmission behavior should also apply to systems with more than three states, as long as the starting point in the parameter space lies on IPCs where two eigenmodes share the lowest loss. The phenomenon is due to the fact that in a dynamical process in non-Hermitian systems, the state prefers to stay on the Riemann sheet with the lowest loss. When the state evolves on the sheet that does not have the lowest loss, NATs would occur after some delay time 23 and the state jumps to the sheet of lower loss. Provided the system is long enough (i.e., L is sufficiently large so that each NAT has enough time to occur), the state will eventually land on the sheet with the lowest loss, i.e., the red sheet in our system, when it approaches the end point of the loop. We can use this principle to study loop-1 with the starting point on IPC-1. Figure 1b indicates that the final state should be mode-2 and mode-1, respectively, for counterclockwise and clockwise loops, since the red sheet is not continuous on IPC-1. This can give an intuitive explanation of the chiral dynamics.
We can apply the same principle to study loop-2 with the starting point lying on IPC-2. Figure 1b indicates that the final state should be mode-3 for counter-clockwise loops whereas mode-2 for clockwise loops. Such prediction can be verified by the calculated amplitudes of the instantaneous eigenmodes for counter-clockwise loops (Fig. 5a-c) and clockwise loops (Fig. 5d-f) with L = 1200 mm. We note at the input and output interface, the power flow of mode-1, mode-2, and mode-3 is mainly localized in waveguide-3, waveguide-2, and waveguide-1, respectively. This is different from the eigenmode distributions for loop-1 (see the inset of Fig. 3d-f). The difference is due to the different widths of waveguide-1 (i.e., α 0 = 1.09 for the starting point of loop-1 and α 0 = 0.937 for that of loop-2) which can affect the energy distribution of the three eigenmodes. Figure 5a-f shows that for both loops, no matter which eigenmode is injected into the system, the final state is always on the lowest red sheet, leading to a chiral transmission behavior due to the topological structure of the energy surfaces around EP-2 (see Fig. 1b).
The chiral transmission behaviors of loop-1 and loop-2 are summarized in Table 1 for comparison. We emphasize that the chiral behavior is not only a result of the location of the starting point on IPCs, but also due to the fact that the eigenmodes on     2). a-c Calculated amplitudes of the three instantaneous eigenmodes along the z-axis for counter-clockwise loop-2 with a mode-1, b mode-2, and c mode-3 being the input eigenmode, where c r , c b , and c g represent, respectively, the coefficient of the eigenmode on the red, blue, and gray sheet. d-f Same as (a-c) but for clockwise loop-2. The dynamics is chiral in the sense that the output state is mode-3 for counter-clockwise loop-2 but mode-2 for clockwise loop-2. The term "NAT" denotes the non-adiabatic transition the IPC exhibit the lowest loss. In our system, IPC-1 and IPC-2 both meet this requirement (see Fig. 1c). But if a system possesses an IPC on which the eigenmodes do not have the lowest loss, then the dynamics may not be chiral since NATs to the lowestloss sheet would occur. When a multi-state system possesses multiple IPCs that exhibit the lowest loss, different IPCs will show different chiral dynamics, i.e., loops in the same encircling direction but enclosing different EPs can give different final states.
Single EP encirclement starting from RPC. We have demonstrated the chiral dynamics when the loops start from IPCs in the proposed system. For loops with a starting point far away from IPCs, the dynamics would not be chiral. To show this point, we consider a loop-3 as illustrated in Fig. 6a. The trajectory of loop-3 is exactly the same as that of loop-1, except that the starting point lies near RPC-1 (α = 1.09 and w 2 = 4 mm; see the circle in Fig. 6a). In this case, the three eigenmodes at the starting/end point have different losses. There is only one eigenmode exhibiting the lowest loss, which is in contrast to the case with a starting point on IPCs where two eigenmodes share the lowest loss. Therefore, counter-clockwise and clockwise loops should make no difference to the final state.
Numerical simulations were performed to verify the above prediction in a system with L = 800 mm. Figure 6b-d plots the trajectories of state evolution on the imaginary part of Riemann sheets for counter-clockwise loops, which are obtained by calculating the amplitudes of instantaneous eigenmodes in Fig. 6e-g, with mode-1 (mainly localized in waveguide-1), mode-2 (mainly localized in waveguide-2), and mode-3 (mainly localized in waveguide-3) being the input eigenmode, respectively. At the starting point, mode-1 exhibits the lowest loss since waveguide-1 is lossless, while the other two waveguides have absorbers attached. Figure 6b shows that the input mode lies on the lowest-loss red sheet. When z > 0 (i.e., after IPC-1), the state climbs up to the blue sheet which is no longer the lowest loss  Non-chiral dynamics for loop-3 that encloses exceptional point 1 (EP-1). a Parameter space of the system where loop-3 has the same trajectory as loop-1 except that the starting point lies near real part coalescence 1 (RPC-1). b-d Trajectories of the state evolution on the imaginary part of Riemann sheets (i.e., the imaginary part of the effective mode index n eff as a function of the detuning parameter α and the absorber width w 2 ) for counter-clockwise loop-3 with b mode-1, c mode-2, and d mode-3 being the input eigenmode. e-g Calculated amplitudes of the three instantaneous eigenmodes along the z-axis for counter-clockwise loop-3 with e mode-1, f mode-2, and g mode-3 being the input eigenmode, where c r , c b , and c g represent, respectively, the coefficient of the eigenmode on the red, blue, and gray sheet. h-j Same as (e-g) but for clockwise loop-3. The dynamics is non-chiral in the sense that the output state is mode-1 for both counter-clockwise and clockwise loop-3. The term "NAT" denotes the non-adiabatic transition sheet. As a result, a NAT to the lowest-loss red sheet occurs at z =~120 mm and the final state is mode-1 (also see Fig. 6e for IPC-1 and the NAT). The dynamical processes with mode-2 ( Fig. 6c) or mode-3 (Fig. 6d) as injections are different at the beginning of the trajectories, i.e., there is a NAT since at the starting point, the eigenmodes on the blue and gray sheet do not have the lowest loss. The following processes are the same, i.e., the state becomes unstable again after IPC-1 and jumps to the red sheet via a second NAT (also see Fig. 6f, g). The final state is always mode-1 since it exhibits the lowest loss. The dynamical processes of clockwise loops are similar and the final state is still mode-1 (see Fig. 6h-j, also see Supplementary Fig. 2). The consequence is a non-chiral transmission behavior (also see Table 1), i.e., different encircling directions result in the same state. The non-chiral dynamics was first observed in two-state systems 27 . Our results in this work indicate that the non-chiral dynamics would also exist in multi-state systems as along as the starting point lies away from IPCs, i.e., there is only one eigenmode exhibiting the lowest loss.
Double EPs encirclement starting from IPC. Generally speaking, a three-state system can support two types of third-order EPs based on the topology of Riemann sheets 31 . The first type of third-order EP has three perturbed eigenvalues around it, with three branches of cube-root expansion. In the vicinity of the second type of third-order EP, two eigenvalues are given by branches of square-root expansion, whereas the third one has a conventional Taylor expansion. In our three-state system, there are two second-order EPs in the parameter space. It was shown that encircling two second-order EPs in such a system will shift eigenmodes in a cyclic manner like encircling the first type of third-order EP mentioned above 18 . To show this point, we consider in Fig. 1d (also see Supplementary Fig. 3a) a loop-4 (α 0 = 1.09, Δα = 0.2, and Δw = 6 mm) that encloses both EP-1 and EP-2 with the starting point on IPC-1. The adiabatic trajectory of state evolution following loop-4 is shown in Supplementary  Fig. 3b, where the sequence of eigenmode variations is indeed the same as that of encircling the first type of third-order EP mentioned above 17,18,31 . In fact, the position of EP-1 and EP-2 can be tuned by changing the system parameters. The degeneracy will become the first type of third-order EP when EP-1 and EP-2 merge together. Therefore, we can use loop-4 to study the dynamical encircling of an equivalent third-order EP.
Let us now consider the dynamical encircling. The calculated amplitudes of instantaneous eigenmodes for counter-clockwise loops are shown in Fig. 7a-c, with mode-1, mode-2, and mode-3 as the injection mode, respectively. The corresponding dynamical trajectories on the real part of Riemann sheets are plotted in Fig. 7d-f. By the same way, the results for clockwise loops are shown in Fig. 7g-l. We find the mode conversion behavior is exactly the same as that of loop-1 (which encloses EP-1 only), i.e., the final state for counter-clockwise loops is mode-2 whereas that for clockwise loops is mode-1 (also see Table 1 and Supplementary Fig. 4). This means the dynamical encircling of a third-order EP with the starting point on the IPC also exhibits a chiral transmission behavior. It is actually the location of the starting point that determines the chiral dynamics, instead of the number or even the order of the encircled EP. A recent study showed that the dynamics can be chiral when the loop excludes any EP 26 , as long as the loop is in the direct vicinity of the EP and the starting point lies near the IPC with two eigenmodes sharing the lowest loss.
Although the consequences are the same, the dynamical processes of encircling a third-order EP are more complex. In brief, there are more NATs (see Fig. 7) than the case when a second-order EP is dynamically encircled (see Figs. 3,4), since in each process the state can cross over the RPC and IPC of both EPs. As a result, the state has more chances to evolve on the unstable sheets and the number of NATs increases. This also indicates that the number of NATs may increase with the increase in the order of the EP.

Discussion
In summary, we have revealed the dynamics when second-order EP and equivalent third-order EP are dynamically encircled in a three-state system (see Table 1 for a summary). From the viewpoint of the encircling process, there are more NATs when compared to that in two-state systems, since the Riemann sheets possess more RPCs and IPCs, via which the state has more chances to evolve on the unstable sheets, which is a necessary condition for NATs to occur. From the viewpoint of outcome, the dynamics is found to be chiral if the starting point lies on IPCs where two eigenmodes share the lowest loss (e.g., loop-1, loop-2, and loop-4). For starting point lying elsewhere (e.g., loop-3 with a starting point near RPCs), the dynamics is non-chiral. Since there are more than one IPC, starting points on different IPCs can result in different chiral dynamics. We have used loop-1 and loop-2 to demonstrate this point. Specifically, for loop-1, the outcome of counter-clockwise loops is mode-2 and that of clockwise loops is mode-1. For loop-2, meanwhile, the final state of counter-clockwise loops is mode-3, while that of clockwise loops is mode-2. This is a phenomenon unique to multi-state systems. We also showed that the chiral and non-chiral dynamics are in fact a result of the location of the starting point, instead of the number of EPs encircled and even the order of the encircled EP. The same outcome for loop-1 and loop-4 is an evidence to this point, although the dynamics when a high-order EP is dynamically encircled (i.e., loop-4) is more complex. The conclusion should also apply to systems with more than three states. In experimental realizations, one may only find approximate chiral/non-chiral dynamics, since the system has a finite length which will induce scatterings into the system. The new physics revealed in our three-state system may lead to the design of new wave manipulation schemes in multi-state systems.

Methods
Numerical simulations. The effective mode index (e.g., in Fig. 1b, c) and eigenfields distributions (e.g., in Fig. 3d-f) of the eigenmodes were calculated using the eigenmode analysis study in electromagnetic waves module of COMSOL. For the simulation of the wave transmissions (e.g., in Fig. 3a-c), we first solved the eigenmode field distributions at the input and output interfaces using boundary mode analysis study in electromagnetic waves module of COMOSL. The obtained eigenfields were then used as boundary conditions and the wave transmission in the system was calculated using the frequency domain study of COMSOL.

Data availability
The data that support the findings of this work are available from the corresponding author upon reasonable request.     2). a-c Calculated amplitudes of the three instantaneous eigenmodes along the z-axis for counter-clockwise loop-4 with a mode-1, b mode-2, and c mode-3 being the input eigenmode, where c r , c b , and c g represent, respectively, the coefficient of the eigenmode on the red, blue and gray sheet. d-f Trajectories of the state evolution on the real part of Riemann sheets (i.e., the real part of the effective mode index n eff as a function of the detuning parameter α and the absorber width w 2 ) for counterclockwise loop-4 with d mode-1, e mode-2, and f mode-3 being the input eigenmode. g-l Same as (a-f) except for clockwise loop-4. The output state is mode-2 for counter-clockwise loop-4, whereas mode-1 for clockwise loop-4. The term "NAT" denotes the non-adiabatic transition

Code availability
All the simulations in this work were performed using COMSOL (www.comsol.com). The source code is available from the corresponding author upon reasonable request.