#
Reduction of classification of a two-dimensional weak topological insulator

- real-space DMFT study -

###### Abstract

One of the remarkable interaction effects on topological insulators is the reduction of topological classification in free-fermion systems. We address this issue in a bilayer honeycomb lattice model by taking into account temperature effects on the reduction. Our analysis, based on the real-space dynamical mean field theory, elucidates the following results. (i) Even when the reduction occurs, the winding number defined by the Green’s function can take a nontrivial value at zero temperature. (ii) The winding number taking the nontrivial value becomes consistent with the absence of gapless edge modes due to Mott behaviors emerging only at the edges. (iii) Temperature effects can restore the gapless edge modes, provided that the energy scale of interactions is smaller than the bulk gap. In addition, we observe the topological edge Mott behavior only in some finite temperature region.

###### pacs:

71.27.+a, 73.20.-r, 71.10.Fd## I Introduction

After the discovery of topological insulators, topological aspects of free-fermion systems have been extensively studied.Hasan and Kane (2010); Qi and Zhang (2011) In these systems, the ground state wave function has topologically nontrivial properties which can be characterized with topological invariants calculated from the Bloch wave function. Nontrivial topology (, or a nontrivial value of topological invariant) in the bulk predicts gapless edge modes, which is known as the bulk-edge correspondence. The presence of these gapless edge modes is a source of characteristic magneto-electric responses. Furthermore, gapless edge modes in the topological superconductor are described as Majorana fermions whose experimental realization has been addressed for one-dimensional quantum wires, recently.Mourik et al. (2012); Rokhinson et al. (2012); Das et al. (2012)

In parallel with study of free-fermion systems, topological insulators for strongly correlated compounds have been proposed.Shitade et al. (2009); Chadov et al. (2010); Lin et al. (2010); Takimoto (2011); Yan et al. (2012); Weng et al. (2014) The first principles calculation has pointed out that can be a strong topological insulator.Takimoto (2011) Furthermore, correlation effects on topological insulators are expected to induce novel phenomena. These facts urge us to study correlation effects on topological insulators.

For analysis of topological insulators with electron interactions, a topological invariant defined in the bulk plays an important role. In non-interacting systems, this quantity can be calculated from the Bloch wave function. The topological invariant for free-fermion systems can be extended to interacting systems.Ishikawa and Matsuyama (1987); Haldane (2004); Volovik and Volovik (2009); Gurarie (2011); Essin and Gurarie (2011); Wang and Zhang (2012a, b) Even though the Bloch wave function is no longer well-defined in interacting systems, the single-particle Green’s function is still available. Thus, the topological invariant can be defined with the Green’s function, which means that the group structure of topological insulators (e.g., or ) is the same as that of the non-interacting case.

Interestingly, however, it is reported that the interactions can induce the reduction of topological classification of free-fermion systems.Fidkowski and Kitaev (2011); Turner et al. (2011); Yao and Ryu (2013); Qi (2013); Lu and Vishwanath (2012); Levin and Stern (2012); Hsieh et al. (2014); Wang et al. (2014); You and Xu (2014); Morimoto et al. (2015) Namely, at the non-interacting level, topological insulators/superconductors labeled by integers show gapless edge modes due to nontrivial properties, while gapless edge modes in some of these topological insulators/superconductors can be gapped out without symmetry breaking under electron interactions. This means that some of topological insulators/superconductors at the non-interacting level become topologically trivial in the presence of interactions, resulting in the reduction of topological classification. The reduction of topological classification has been mainly analyzed in terms of symmetry protection of edge modes at zero temperature.

Now, one may ask the following three questions. (i) How the topological invariant defined by the Green’s function behaves if the reduction of topological classification occurs? Does it still take a nontrivial value? (ii) If the topological invariant of the Green’s function is well-defined and takes a nontrivial value, how this nontrivial value in the bulk becomes consistent in the absence of the gapless edge modes? (iii) What are temperature effects on the reduction of topological classification? The analysis at finite temperatures is desired because all of experiments are carried out there.

In this paper, we address these three questions by analyzing a correlated weak topological insulator in two dimensions. The topological structure is characterized by the winding number at the non-interacting level, which is also well-defined even in correlated systems. Our analysis based on the real-space dynamical mean field theory with continuous-time quantum Monte Carlo (R-DMFT+CTQMC) elucidates the following facts. First, the winding number can take a nontrivial value even when the reduction of topological classification occurs. Second, the system shows a Mott behavior (i.e., the divergence of the self-energy) only around the edge, signaling the reduction of topological classification. Third, temperature effects result in a gradual crossover in gapless edge modes. Namely, the behavior of edge modes changes with increasing temperature from zero temperature: starting from the absence of gapless excitations, gapless spin excitations emerge at slightly higher temperatures, which finally leads to single-particle gapless excitations in the higher temperature region.

The rest of this paper is organized as follows. The next section (Sec. II) is devoted to the topological invariant with the Green’s function as well as description of our model and method. In Sec. III, we show that interactions induce the reduction of topological classification from to and analyze the resulting properties in detail by using the R-DMFT. Short summary is given in Sec. IV.

## Ii Model and method

### ii.1 Lattice model

We study a bilayer honeycomb lattice model. The Hamiltonian reads

(1a) | |||||

(1b) | |||||

(1c) |

where the operator creates an electron at site and in layer and spin state. is intra-layer hopping between sites and , and the inter-layer hopping is assumed to be zero. If an electron hops in -direction, we set , otherwise ( is a real number). A sketch of the hopping is shown in Fig. 1; the brown (gray) bonds denote hopping with () respectively. The electrons interact via intra-layer on-site interaction and spin-exchange interaction . We consider the system under the open (periodic) boundary condition for - (-) direction respectively. The system is composed of rings along the -direction. denotes the -coordinate of these rings.

### ii.2 Numerical approach: R-DMFT

For systematic analysis of the bulk and the edge properties, we consider the system with a cylinder geometry. Namely, we impose the open (periodic) boundary condition for - (-) direction, respectively. This system is an inhomogeneous correlated system. The real-space dynamical mean field theory,Potthoff and Nolting (1999); Okamoto and Millis (2004); Helmes et al. (2008); Snoek et al. (2008) an extended version of the DMFT,Mul̈ler-Hartmann (1989); Metzner and Vollhardt (1989); Georges et al. (1996) enables us to treat it.

In this method, the problem is solved self-consistently. A brief outline of this approach is given in the following. In the first step, a one-dimensional subsystem of -th sites along the -axis is mapped onto an effective impurity model characterized by the cavity Green’s function . In the second step, applying an impurity solver to this effective model yields the self-energy . In the third step, the effective models are updated with obtained self-energy, and then we go back to the first step. The self-consistent equation is given by

(2a) | |||||

(2b) | |||||

with diagonal matrices | |||||

(2c) | |||||

(2d) | |||||

where is the Fourier transform of the hopping matrix and is the Green’s function for the lattice model. |

The effective impurity model is essentially the same as that of two-orbital systems. With path integral formalism, the effective model for -th sites are written as

(3a) | |||||

(3b) | |||||

with . To solve the effective model we employ the continuous-time quantum Monte CalroRubtsov et al. (2005); Werner et al. (2006); Werner and Millis (2006); Haule (2007); Gull et al. (2011) which is a powerful tool to analyze multi-orbital systems.

### ii.3 The winding number

For and , the model (1) having two-sublattices A and B realizes a weak topological insulator which can be regarded as an array of one-dimensional topological insulators with chiral symmetry along -direction in the momentum space. Accordingly, gapless edge modes are observed along zigzag edges, while no gapless edge mode is observed along armchair edges. The topological properties of this insulator are protected by the chiral symmetry and characterized by the winding number which is defined as Gurarie (2011); Essin and Gurarie (2011); Wang and Zhang (2012b); Manmana et al. (2012)

(4) |

where with is the full Green’s function for and under the torus geometry. In our bilayer system, is a -matrix. ’s are Pauli matrices acting on the sublattice space. The winding number is integral (i.e., ), which reflects the mathematical fact, with a positive integer .Essin and Gurarie (2011)

## Iii Results

First, we address the instability of gapless edge modes against interactions, which leads to the reduction of the topological classification from to . After that, we discuss the numerical results obtained by the R-DMFT+CTQMC.

### iii.1 classification

We show that gapless edge modes are destroyed in four copies of a weak topological insulator in two dimensions with chiral symmetry. This result indicates that the topological classification in free fermions changes to . We observe the instability of gapless edge modes against the interactions with the following two steps: (i) Adiabatically deforming the effective model of Eq. (1b) into a simple model; (ii) Introducing the interactions to destroy the gapless edge modes without symmetry breaking.

Let us start with the effective Hamiltonian of the system under the cylinder geometry,

(5) |

where describes hopping of electrons with the momentum along the -direction. denotes the four-dimensional identity matrix acting on the spin and the pseudospin space of and . In this basis the operator for chiral symmetry is given by ; . Zero modes are localized at the kink for arbitrary , where changes its sign. Thus, we can adiabatically deform the Hamiltonian, , into the following Hamiltonian

(6) | |||||

Electrons in this deformed model, Eq. (6), do not hop in the -direction. So, each decoupled chain along the -direction behaves as an ordinary one-dimensional topological insulator with chiral symmetry.

In this one-dimensional topological insulator, we can observe the instability of gapless edge modes in the following way:Fidkowski and Kitaev (2010) (a) Introducing reduces the system to the effective spin model showing a free spin of at the edge of each chain; (b) Introducing antiferromagnetic coupling gaps out these free spins.

Therefore we can gap out all of the gapless edge modes without symmetry breaking, which indicates the classification.

### iii.2 Numerical results

Here, we address the reduction of the topological classification by using the R-DMFT. We focus on the paramagnetic phase because no continuous symmetry is broken in two-dimensional systems at finite temperatures. Here, we set the parameter and . In the following we mainly focus on the case of . This is because the qualitative behaviors do not differ from other cases of interaction strength except for the temperature effects which we will discuss in Sec. III.2.2.

The main results are summarized in Fig. 2, which are obtained under the cylinder geometry. First, we elucidate how the winding number behaves when the system shows the reduction of the topological classification for free fermions. As seen in Fig. 2(a), switching on the interactions and does not change the winding number in the bulk; the winding number remains four for sufficiently weak interactions. We note that the gapless edge modes are destroyed in the region of as seen in Fig. 2(b) where the local density of states at zero energy and is plotted.

Second, we elucidate how the winding number taking the nontrivial value becomes consistent with the absence of gapless edge modes. We demonstrate that the system shows Mott behaviors only around the edges, making the winding number and the absence of gapless modes consistent. The Mott behaviors only around the edge are observed via a divergence of the self-energy and an abrupt change of double occupancy.

Finally, we elucidate that the region of shows difference from the trivial band insulator in the finite temperature region, although these two phases are topologically identical at zero temperature.

In the following, we explain the details. In Sec. III.2.1, we analyze the reduction of the topological classification in the bulk and the edges and elucidates the behaviors of the winding number. In Sec. III.2.2, we address the temperature effects on the reduction of the classification.

#### iii.2.1 The reduction and the winding number

In this section, we show that the winding number takes four, while the gapless edge modes are destroyed due to correlations. We then elucidate how the nontrivial value of the winding number becomes consistent with the absence of gapless edge modes. For this aim, we focus on the low temperature region , where the bulk electron gap is large enough than the temperature scale so that we can characterize the topological structure.

First, we discuss behaviors of local correlations as a function of interaction strength . With turning on the interactions, the double occupancy for the bulk () gradually changes [see Fig. 3(a)]. This behavior indicates that the bulk behaves as a renormalized band insulator. Correspondingly, the winding number remains four for [see Fig. 3(c)]. Further increasing the interactions, the system changes into the Mott insulating phase where the electrons are localized and form a dimer phase [see Fig. 3(a) and (b)]. Therefore, the system is no longer a renormalized band insulator, and the topological nontrivial structure is destroyed. We note that our results indicate a first order Mott transition because the hysteresis behavior is observed in the region of for .foo () Increasing temperature narrows the region of hysteresis behavior which is sandwiched by dashed blue and red lines in Fig. 2(a). This indicates that the first order transition changes to a crossover in high temperature region. In the above we have seen that if the interactions, , are weak, then, the system behaves as a renormalized band insulator showing the winding number .

Here, we show that even for weak interactions, the gapless edge modes are destroyed in spite of the nontrivial values of the winding number . Destruction of gapless electron excitations at edges can be seen in Figs. 4 (a) and (b). In these figures, the local density of states (LDOSs), with , for the bulk and edges are plotted. For the non-interacting case [Fig. 4(a)], the gap is observed in the bulk, while a sharp peak of LDOS is observed for at the edge, signaling the presence of gapless edge modes. Switching on the interaction does not qualitatively change the bulk gap. At the edge, however, the sharp peak at disappears [Fig. 4(b)]. A similar behavior is also observed in the collective excitations. Therefore, we can conclude that the introducing interactions with destroys the gapless edge modes observed in the non-interacting case. This indicates the reduction of the topological classification for free fermions.

In the above we have seen that the winding number defined in Eq. (4) takes four although the system does not show any gapless modes. At first glance, the nontrivial winding number in the bulk seems to contradict the absence of gapless edge modes. Here, we elucidate how these two results become consistent with each other. Mathematically, the system does not change its topological properties under continuous deformation, provided that the Green’s function at () does not satisfy and . Thus, when the bulk has the nontrivial topology, the following two scenarios are expected at the edge, which separates two distinct phases with different winding numbers. (a) One possibility is that the Green’s function diverges at edges and becomes singular because the topological invariant is not well-defined at edges. The physical meaning of this divergence is the emergence of zero energy excitations. This is indeed observed in the non-interacting case [see Fig. 4(a)], which is consistent with the bulk-edge correspondence in the non-interacting case. (b) The other possibility is that the Green function becomes zero (i.e., the divergence of the self-energy). The zeros of the Green’s function are induced only by divergence of the self-energy (i.e., the zeros are absent for the non-interacting case). In the region of and in Fig. 2(a), the scenario (b) holds. In the following, we demonstrate that the system shows the Mott behavior only around the edge, which keeps the winding number and thus makes the absence of gapless edge modes consistent. A signal of the Mott behavior emerging only around the edge is observed in the interaction dependence of the local correlations [Figs. 3(a) and (b)]. Once the interactions are turned on, the double occupancy and the spin correlation suddenly drop only at , while these local correlations gradually change in the bulk. In order to observe the Mott behavior more clearly, we plot the imaginary part of self-energy at , , and . This figure indicates that the self-energy diverges at but not at and . Therefore, we conclude that the system shows the Mott behavior only around the edge, whereby the winding number and the absence of gapless edge modes become consistent.

#### iii.2.2 Temperature effects

In this section, we address the third question; What are temperature effects on the reduction of topological classification? In the previous section, we have discussed the reduction. The point is that if the reduction occurs the zeros of Green’s function emerge at the edge due to consistency with the nontrivial values of winding number. Therefore, there is a chance to recover gapless modes; if the singularity of the self-energy is suppressed, the gapless mode would show up again. Finite temperature effects indeed recover the gapless modes.

For , we can observe that increasing temperature suppresses the singularity of the self-energy and restores the gapless electron excitations; Fig. 4(c) indicates that the gapless edge modes are restored with increasing temperature, while the edge modes are destroyed for . Accordingly, the double occupancy at the edge gradually increases with increasing temperature.

Furthermore, for , the temperature effects induce richer properties. Increasing temperature changes statistical properties of gapless excitations if the energy scale of interactions is smaller than the bulk gap. Namely, with increasing temperature, the bosonic excitations are restored first, after that, the fermionic excitations are restored. Suppose that the edge modes are gapped out by interactions in the phase with the winding number . Then, with increasing temperature a topological edge Mott state emerges, which is accompanied by gapless edge modes only in the spin excitations. This is because the phase with is decoupled to two-copies of the weak topological insulator with , each of which has gapless spin excitations as edge modes, if temperature effects are dominant than the spin exchange interactions but smaller than the on-site interaction .

To show this, we compute LDOSs and local spin excitations at the edge for . In Figs. 5 (a) and (b), the spectral weight of the LDOS , and spin excitations for at the edge () are plotted. With increasing temperature, the spectral weights at increases, which indicates the restoration of gapless edge modes. Here, the dashed (solid) lines denote crossover temperature where the single-particle excitations (the spin excitations) become gapless, respectively. Namely, in the region sandwiched by dashed and solid green lines, the system shows gapless modes only in spin (i.e., bosonic) excitations. Further increasing interactions we can observe the gapless fermionic excitations. Here, we show the details of the data. The LDOS (the local spin excitation) is plotted in Figs. 5 (c) and (d) [(e) and (f)]. The data for and , [see Figs. 5 (c) and (e)] indicate that the spin excitations show a coherence peak for while electron excitations are still gapped at this temperature. Absence of gapless modes in the low temperature region is more clearly observed for and , where the restoration of gapless modes is also observed [see Figs. 5 (d) and (f)]. Performing similar calculations we end up with Figs. 5 (a) and (b) where we can find the emergence of the topological edge Mott state only at finite temperatures.

## Iv summary

In this paper, we have analyzed the reduction of topological classification in the two-dimensional correlated systems by using the R-DMFT+CTQMC. In particular, we have elucidated the following behaviors. (i) The topological invariant (4) takes a nontrivial value even when the reduction occurs. (ii) The absence of zero modes becomes consistent with the finite winding number due to the emergence of the zeros of the Green’s function only at the edge. The system shows the Mott behaviors only around the edge, resulting in the singularity of the self-energy only at the edge. (iii) Increasing temperature removes the zeros of the Green’s function emerging only around the edges and restores gapless electron excitations around the edge. Besides, we have found that the temperature effects can change statistical properties of the restored gapless edge modes.

The result (i) indicates that calculating the winding number in Eq. (4) is insufficient to describe the reduction in the bulk, although it is well-defined even in the presence of electron correlation. For the direct observation of the reduction in the bulk, it is recently proposed to calculate the angle of partition functions on unoriented space-time.Shapourian et al. (2016); Shiozaki et al. (2016) The numerical simulation of it is left for our future work.

## V acknowledgements

This work is partly supported by a Grand-in-Aid for Scientific Research on Innovative Areas (JSPS KAKENHI Grant No. JP15H05855) and also JSPS KAKENHI (No. 16K05501). The numerical calculations were per- formed on supercomputer at the ISSP in the University of Tokyo, and the SR16000 at YITP in Kyoto University.

## References

- Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
- Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
- Mourik et al. (2012) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
- Rokhinson et al. (2012) L. P. Rokhinson, X. Liu, and J. K. Furdyna, Nature Physics 8, 795 (2012).
- Das et al. (2012) A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman, Nature Physics 8, 887 (2012).
- Shitade et al. (2009) A. Shitade, H. Katsura, J. Kuneš, X.-L. Qi, S.-C. Zhang, and N. Nagaosa, Phys. Rev. Lett. 102, 256403 (2009).
- Chadov et al. (2010) S. Chadov, X. Qi, J. Kübler, G. H. Fecher, C. Felser, and S. C. Zhang, Nat. Mater. 9, 541 (2010).
- Lin et al. (2010) H. Lin, L. A. Wray, Y. Xia, S. Xu, S. Jia, R. J. Cava, A. Bansil, and M. Z. Hasan, Nat. Mater. 9, 546 (2010).
- Takimoto (2011) T. Takimoto, Journal of the Physical Society of Japan 80, 123710 (2011), http://dx.doi.org/10.1143/JPSJ.80.123710 .
- Yan et al. (2012) B. Yan, L. Müchler, X.-L. Qi, S.-C. Zhang, and C. Felser, Phys. Rev. B 85, 165125 (2012).
- Weng et al. (2014) H. Weng, J. Zhao, Z. Wang, Z. Fang, and X. Dai, Phys. Rev. Lett. 112, 016403 (2014).
- Ishikawa and Matsuyama (1987) K. Ishikawa and T. Matsuyama, Nuclear Physics B 280, 523 (1987).
- Haldane (2004) F. D. M. Haldane, Phys. Rev. Lett. 93, 206602 (2004).
- Volovik and Volovik (2009) G. E. Volovik and G. Volovik, (2009).
- Gurarie (2011) V. Gurarie, Phys. Rev. B 83, 085426 (2011).
- Essin and Gurarie (2011) A. M. Essin and V. Gurarie, Phys. Rev. B 84, 125132 (2011).
- Wang and Zhang (2012a) Z. Wang and S.-C. Zhang, Phys. Rev. X 2, 031008 (2012a).
- Wang and Zhang (2012b) Z. Wang and S.-C. Zhang, Phys. Rev. B 86, 165116 (2012b).
- Fidkowski and Kitaev (2011) L. Fidkowski and A. Kitaev, Phys. Rev. B 83, 075103 (2011).
- Turner et al. (2011) A. M. Turner, F. Pollmann, and E. Berg, Phys. Rev. B 83, 075102 (2011).
- Yao and Ryu (2013) H. Yao and S. Ryu, Phys. Rev. B 88, 064507 (2013).
- Qi (2013) X.-L. Qi, New J. Phys. 15, 065002 (2013).
- Lu and Vishwanath (2012) Y.-M. Lu and A. Vishwanath, Phys. Rev. B 86, 125119 (2012).
- Levin and Stern (2012) M. Levin and A. Stern, Phys. Rev. B 86, 115131 (2012).
- Hsieh et al. (2014) C.-T. Hsieh, T. Morimoto, and S. Ryu, Phys. Rev. B 90, 245111 (2014).
- Wang et al. (2014) C. Wang, A. C. Potter, and T. Senthil, Science 343, 629 (2014), http://science.sciencemag.org/content/343/6171/629.full.pdf .
- You and Xu (2014) Y.-Z. You and C. Xu, Phys. Rev. B 90, 245120 (2014).
- Morimoto et al. (2015) T. Morimoto, A. Furusaki, and C. Mudry, Phys. Rev. B 92, 125104 (2015).
- Potthoff and Nolting (1999) M. Potthoff and W. Nolting, Phys. Rev. B 59, 2549 (1999).
- Okamoto and Millis (2004) S. Okamoto and A. J. Millis, Phys. Rev. B 70, 241104 (2004).
- Helmes et al. (2008) R. W. Helmes, T. A. Costi, and A. Rosch, Phys. Rev. Lett. 100, 056403 (2008).
- Snoek et al. (2008) M. Snoek, I. Titvinidze, C. Tőke, K. Byczuk, and W. Hofstetter, New Journal of Physics 10, 093008 (2008).
- Mul̈ler-Hartmann (1989) E. Mul̈ler-Hartmann, Zeitschrift fur̈ Physik B Condensed Matter 74, 507 (1989).
- Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
- Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- Rubtsov et al. (2005) A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
- Werner et al. (2006) P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
- Werner and Millis (2006) P. Werner and A. J. Millis, Phys. Rev. B 74, 155107 (2006).
- Haule (2007) K. Haule, Phys. Rev. B 75, 155113 (2007).
- Gull et al. (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
- Manmana et al. (2012) S. R. Manmana, A. M. Essin, R. M. Noack, and V. Gurarie, Phys. Rev. B 86, 205119 (2012).
- Fidkowski and Kitaev (2010) L. Fidkowski and A. Kitaev, Phys. Rev. B 81, 134509 (2010).
- (43) In two-dimensional systems, spatial fluctuations is important. In our approach, the spatial fluctuation is neglected. Hence, to obtain accurate results concerning the order of the transition, one should perform simulations by taking into account spatial fluctuation .
- Shapourian et al. (2016) H. Shapourian, K. Shiozaki, and S. Ryu, arXiv preprint arXiv:1607.03896 (2016).
- Shiozaki et al. (2016) K. Shiozaki, H. Shapourian, and S. Ryu, arXiv preprint arXiv:1609.05970 (2016).