\section{Introduction}\label{sec:intro}

The defining property of topological insulators (TIs)~\cite{footnote} is the band gap in the energy spectrum of the bulk material and gapless conducting boundary states. They are edge states in the case of two-dimensional (2D) systems or surfaces states in the case of three-dimensional ones~\cite{Haldane2017}. The paradigmatic cases which gave rise to the main concepts~\cite{Hasan2010,Qi2011} in this field are: ({\em i}) quantum Hall insulator (QHI) in 2D electron gas  which requires an external magnetic field to break the time-reversal invariance and whose edge states are {\em chiral} by allowing spin-unpolarized electron to propagate in only one direction; and ({\em ii})  quantum spin Hall insulator (QSHI)~\cite{Kane2005} which is time-reversal invariant, thereby requiring strong spin-orbit coupling effects instead of external magnetic field, and whose edge states appear in pairs with  different chirality and spin polarization. The last experimentally discovered member of 2D TI family is quantum anomalous Hall insulator (QAHI), which requires both nonzero magnetization to break the time reversal invariance and strong spin-orbit coupling effects, with its edge states allowing only one spin species to flow unidirectionally~\cite{Liu2016}.

\begin{figure}
	\includegraphics[width=8.5cm]{fig1.pdf}
	\caption{Schematic view of a two-terminal devices where an infinite ZGNR is attached to two macroscopic reservoirs with chemical potentials  $\mu_\mathrm{L}$ and $\mu_\mathrm{R}$ on the left and right, respectively, so that $\mu_\mathrm{L}-\mu_\mathrm{R} = eV_b$ is externally  applied dc bias voltage. In panel (a), the scattering region (blue shaded) in the middle of {\em finite length} $L = 30\sqrt{3}a$ and width $W=29a$ is Floquet TI generated by irradiating~\cite{Oka2009} segment of ZGNR by circularly polarized light of intensity $z$ and frequency $\omega$ [Eq.~\eqref{eq:pphase}]. In panel (b), the scattering region (shaded green)  is quantum Hall, quantum anomalous Hall or quantum spin Hall insulator  with their parameters tuned to produce the same topologically nontrivial band gap $\Delta_g$ [Fig.~\ref{fig:fig2}(a)] in the bulk of all such conventional time-independent TIs.}
	\label{fig:fig1}
\end{figure}

In theoretical analysis, edge states are found as eigenfunctions $\Psi_{k_x}(x,y)=e^{i k_x x}\psi(y)$ of the Hamiltonian of an infinite wire (periodic along the $x$-axis, so that eigenfunctions are labeled by the wavevector $k_x$) made of 2D TIs. The corresponding eigenenergies $\varepsilon(k_x)$ form subbands crossing the band gap~\cite{Halperin1982,MacDonald1984}. The width of the edge states is defined by the spatial region where the probability density is nonzero,  \mbox{$|\psi(y)|^2 \neq 0$}, while decaying  exponentially fast towards the bulk of the wire. Interestingly, their width~\cite{Prada2013,Chang2014} can also depend on the arrangement of atoms along the edge, such as in the case of graphene wires where edge states of QHI and QAHI or QSHI are narrower in the case of zigzag arrangement of carbon atoms along the edge than in the case of their armchair arrangement~\cite{Prada2013,Chang2014,Sheng2017}. In paradigmatic three-dimensional TIs like Bi$_2$Se$_3$, surface states actually have spatial extent of about $\sim 2$ nm~\cite{Chang2015}.

The zigzag edge, which is employed in devices in Fig.~\ref{fig:fig1}, can also introduce a kink in the subband of edge state~\cite{Sheng2017}, so that subband   intersects with the Fermi energy $E_F$ at $N_\mathrm{R}$ points with positive velocity and $N_\mathrm{L}$ points with negative velocity. However,  only the difference $N_\mathrm{R} - N_\mathrm{L}=|\mathcal{C}|$ is topologically protected according to the {\em bulk-boundary correspondence}~\cite{Hasan2010,Qi2011}. Here  $\mathcal{C}$ is an integer topological invariant (like the Chern number in the case of QHI and QAHI) associated with band structure in the bulk. This makes   electronic transport through edge states of infinite length {\em perfectly quantized in a robust way}~\cite{Buttiker1988}---the zero-temperature two-terminal conductance is $G(E_F)=G_\mathrm{Q} |\mathcal{C}|$  for $E_F$ swept through the bulk band gap and insensitive to both magnetic and nonmagnetic disorder in the case of QHI and QAHI~\cite{Sheng2017}, or only nonmagnetic disorder in the case of QSHI. Although infinite ballistic wires, including those with topologically trivial edge states~\cite{Allen2015,Zhu2017,Marmolejo-Tejada2018}, also exhibit integer $G(E_F)/ G_\mathrm{Q}$, this is easily disrupted by disorder introduced around the edges or even within the bulk~\cite{Marmolejo-Tejada2018}. Here $G_\mathrm{Q}=2e^2/h$ or $G_Q=e^2/h$ is the conductance quantum for spin-degenerate or spin-polarized edges states, respectively.

Thus, it has been considered that the key experimental signature of topology in 2D condensed matter is conductance quantization in transport through edge states, which persists even in the presence of disorder as long as it does not break underlying symmetries of the topological phase or generates energy scales larger than the bulk band gap. However, for QHI, QAHI and QSHI of finite length, the zero-temperature longitudinal  conductance $G=I/V_b$, also denoted as `two-terminal' since current $I$ and small bias voltage $V_b$ are measured between the same normal metal (NM) leads,   oscillates in Fig.~\ref{fig:fig2} just below the quantized plateau at $2e^2/h$ while remaining very close to it. We use zigzag graphene nanoribbon (ZGNR) within which 2D TI of finite length [Fig.~\ref{fig:fig1}(b)] is established  using sufficiently large external magnetic field~\cite{Cresti2016}, or additional terms of the Haldane~\cite{Haldane1988} or the Kane-Mele~\cite{Kane2005} models, to generate QHI, QAHI and QSHI, respectively. Their parameters are tuned so that all three examples of conventional time-independent 2D TIs in Fig.~\ref{fig:fig2}(a) have identical topologically nontrivial bulk band gap $\Delta_g$. Even though $G(E_F)$, for $E_F$ swept through bulk band gap $\Delta_g$, is not perfectly quantized in Fig.~\ref{fig:fig2}(a), its oscillations zoomed in Figs.~\ref{fig:fig2}(b)--(g) are insensitive to nonmagnetic edge disorder (ED) introduced in the form of edge vacancies [illustrated in Fig.~\ref{fig:fig4}].

It is worth mentioning that imperfectly quantized two-terminal $G(E_F)$ was observed in early experiments on QSHI~\cite{Roth2009}, provoking a lively search for exotic many-body inelastic effects~\cite{Maciejko2009,Tanaka2011,Budich2012,Vayrynen2013,Kainaris2014,Vayrynen2018,Novelli2019} which can circumvent band-topology constraints and introduce  backscattering of electrons as they propagate through edge states. On the other hand, Fig.~\ref{fig:fig2} demonstrates that imperfectly quantized $G(E_F)$ can be due to a much simpler mechanism---backscattering at the NM-lead/TI-region interface.

Lacking perfectly quantized two-terminal conductance $G(E_F)$ as the experimental signature of 2D TI phase of finite length, one can resort to direct imaging of spatial profiles of local current density that should confirm electronic flux confined to a narrow region defined by the edge states. Continuous experimental advances have made this possible, such as by using superconducting interferometry in Josephson junction setup~\cite{Allen2015,Zhu2017} or scanning superconducting quantum interference device (SQUID)~\cite{Nowack2013}. In the latter case, one images  magnetic field produced by the current from  which one can reconstruct the local current density with  $\sim \mu$m spatial resolution~\cite{Nowack2013}. Even higher resolution, with reconstructed images having spatial resolution of $\sim 10$ nm, has been achieved by using scanning tip based on electronic spin of a diamond nitrogen-vacancy (NV)  centers~\cite{Tetienne2017,Chang2017,Ku2020}. Particularly {\em intriguing questions} that such images can answer is how electron flux transitions from topologically trivial NM lead present in every experimental device into the region of 2D TI of finite length where the flux is confined within narrow edge currents, as well as how processes at the NM-lead/2D-TI interface affect the total current and the corresponding conductance.

\begin{figure*}
	\includegraphics[scale=1.1]{fig2.pdf}
	\caption{(a) The zero-temperature two-terminal conductance vs. the Fermi energy $E_F$ of device in Fig.~\ref{fig:fig1}(b) where central scattering region hosts {\em finite-length} conventional 2D  time-independent TIs, such as QHI, QAHI and QSHI. The TIs are defined on pristine or edge disordered (denoted by ED) ZGNR due to vacancies illustrated in Figs.~\ref{fig:fig4}(d), ~\ref{fig:fig4}(f) and ~\ref{fig:fig4}(h).  The zoom in of conductance values within the rectangle in panel (a) is shown in: (b)--(d) for pristine ZGNR; and (e)--(g) for edge-disordered ZGNR.  The two NM leads, from which electrons are injected into the topologically protected edge states of finite length with the corresponding local current density profiles shown in Fig.~\ref{fig:fig4}(c)--(h), are also made of ZGNR of the same width as the scattering region [Fig.~\ref{fig:fig1}(b)]. The gap in the bulk of all three 2D TIs is tuned to $\Delta_g = 0.54 \gamma$ and marked in panel (a).}
	\label{fig:fig2}
\end{figure*}

Imaging of local current density could also offer new avenue for resolving a {\em crucial issue} for recently conjectured new class of 2D TIs---the so-called Floquet TI~\cite{Oka2009,Lindner2011,Oka2019,Rudner2020}---which is the connection between the Floquet quasi-energy
spectrum and experimentally measurable dc transport properties. The Floquet TI phase arises in 2D electron systems driven out of equilibrium by strong light-matter interaction.  For example,  graphene~\cite{Oka2009,Lindner2011,Oka2019,Rudner2020}, as well as other 2D materials with honeycomb lattice structure like transition-metal dichalcogenides~\cite{Huaman2019}, subject to a spatially uniform and circularly polarized  light are predicted to transmute into Floquet TI with quasi-energy spectrum~\cite{Shirley1965,Sambe1973}. Its multiple gaps share~\cite{Lindner2011} the same topological properties as the band  gap of QAHI described by the Haldane model~\cite{Haldane1988}. This means that the laser induced band gaps, such as $\Delta_0$ in Fig.~\ref{fig:fig3}(a) emerging at the charge neutral point (CNP) of graphene and $\Delta_1$ away from CNP, are crossed by subbands of chiral edge states~\cite{Perez-Piskunow2014,Perez-Piskunow2014a}. The eigenfunctions of these subbands  decay exponentially towards the bulk with a decay length that depends only on the ratio of the laser frequency and its intensity.

The $\Delta_1$ gaps are called dynamical gaps~\cite{Syzranov2008} and they occur at energy $\hbar \omega/2$ above/below the CNP. They can be reached using experimentally accessible parameters. For example, the very recent experiment~\cite{McIver2020} has been interpreted in terms of creation of a transient Floquet TI by driving graphene flake by \mbox{$500$~fs} laser pulse at a frequency of \mbox{$\omega=46$~THz}, so that the photon energy is \mbox{$\hbar \omega \approx 191$ meV} and its wavelength is \mbox{$\lambda \approx 6.5$ $\mu$m}. However, the experiment of Ref.~\cite{McIver2020} did not observe either quantized longitudinal or transverse (Hall) conductance. Instead, they found that at the peak laser pulse fluence the transverse conductance within $\Delta_0$ gap saturated at plateau around \mbox{$G_{xy}=(1.8 \pm 0.4)e^2/h$}, while no such plateau of $G_{xy}$ was observed within $\Delta_1$ gap.

The calculations of two-terminal [as in Figs.~\ref{fig:fig2} and  ~\ref{fig:fig3}(b)] or multi-terminal conductance typically employ  the Landauer-B\"{u}ttiker setup~\cite{Buttiker1988,Imry1999} depicted in Fig.~\ref{fig:fig1} where finite-size
scattering region---time-dependent due to light irradiation in Fig.~\ref{fig:fig1}(a) or conventional time-independent in Fig.~\ref{fig:fig1}(b)---is attached to semi-infinite NM leads terminating at infinity into the macroscopic particle reservoirs. This is highly appropriate for Floquet TI since time-dependent potential applied in experiments~\cite{McIver2020} is confined to a finite region, either because of  a finite laser spot or the screening inside metallic contacts. On the conceptual side, such setup ensures well-defined asymptotic states and their occupation far away from the irradiated region, thereby evading technical difficulties when using time-dependent leads or reservoirs~\cite{Gaury2014}. It also ensures {\em continuous energy spectrum} of the whole system which plays a key role in both  the Landauer-B\"{u}ttiker and Kubo~\cite{Sato2019} formulation of quantum transport because it  effectively introduces dissipation at infinity and thereby steady-state transport~\cite{Perfetto2010}, while not requiring~\cite{Esin2018} to explicitly model many-body inelastic scattering processes responsible for dissipation~\cite{Imry1999}.

\begin{figure}
	\includegraphics[width=\linewidth]{fig3.pdf}
	\caption{(a) Quasi-energy spectrum $\xi_\mathrm{QE} (k_x)$ for an infinite ZGNR that is irradiated by circularly polarized monochromatic laser light of frequency \mbox{$\hbar \omega = 3 \gamma$} and intensity \mbox{$z=0.5$} over its whole length. The spectrum is obtained by diagonalizing the corresponding Floquet Hamiltonian [Eq.~\eqref{eq:floquet_ham}] truncated to $-N_\mathrm{ph}<n<N_\mathrm{ph}$ Floquet replicas where $N_\mathrm{ph}=7$ is chosen. The yellow shaded region marks the topological gap $\Delta_0$ around $\xi_\mathrm{QE} = 0$ corresponding to CNP, while the red shaded region marks the dynamical  topological gap $\Delta_1$  around  $\xi_\mathrm{QE} = \pm\hbar\omega/2$. (b) The zero-temperature two-terminal conductance vs. $E_F$ (computed using $N_\mathrm{ph} = 7$) of two-terminal device in Fig.~\ref{fig:fig1}(a) whose scattering region is Floquet TI of finite length due to irradiation by circularly polarized light. The pristine irradiated ZGNR is marked by FTI and irradiated edge-disordered ZGNR is marked by FTI-ED. The conductance of an infinite nonirradiated (NIR) pristine ZGNR is also show as a reference. (c) Total DOS for the same device marked by FTI in panel (b). (d) Convergence of lead currents $I_\mathrm{L}$ and $I_\mathrm{R}$ vs. $N_\mathrm{ph}$ at \mbox{$E_F = \hbar\omega/2$ }. }
	\label{fig:fig3}
\end{figure}

However, for the same two-terminal Landauer-B\"{u}ttiker setup with irradiated scattering region a plethora of conflicting theoretical conclusions have been reached~\cite{Rudner2020}. For example, Refs.~\cite{Kitagawa2011,Gu2011} predict quantization of longitudinal dc conductance within a few percent of $2e^2/h$, while Ref.~\cite{Kundu2014} finds its anomalous suppression. To recover the quantized value, Ref.~\cite{Kundu2013} proposed an {\em ad hoc}  summation procedure over different energies in the lead. Without utilizing such ``Floquet sum rule''~\cite{Farrell2016,Yap2017,Kundu2020}, both Refs.~\cite{FoaTorres2014,Farrell2016} confirm nonquantized $G<2e^2/h$ within $\Delta_0$ gap and $G<4e^2/h$ within $\Delta_1$ gap which, however, are largely insensitive to disorder like vacancies or on-site impurities. The precise quantization could be disrupted by dc component of pumping current~\cite{Brouwer1998,Moskalets2002}, which appears~\cite{FoaTorres2014}  even at zero bias voltage  due to time-dependent potential in the Hamiltonian whenever the left-right symmetry of the device is broken statically or dynamically~\cite{FoaTorres2005,Bajpai2019}.

The absence of quantization is explained~\cite{Esin2018,Rudner2020,FoaTorres2014,Farrell2016} by the mismatch between nonirradiated electronic states in the NM leads and edge  states within the gaps of the Floquet TI. The mismatch between states in topologically trivial NM leads and TI scattering region exists also in conventional time-independent TI devices, but without significant disruption of quantized conductance in Fig.~\ref{fig:fig2}. However, specific to  Floquet TIs is possibility of Floquet replicas to couple to bulk bands~\cite{FoaTorres2014,Fedorova2019}. That is, although edge states within the gap $\Delta_0$ are  primarily built from states near the CNP of nonirradiated graphene, they also contain harmonic components near $\pm n\hbar \omega$ which open possibility for electronic photon-assisted tunneling into or out of states in the NM leads whose energies are far away from the CNP. Thus, engineering the density of state of the leads, in order to connect Floquet TI and macroscopic reservoirs through a narrow band of filter states, can recover longitudinal dc conductance within a few percent of $2e^2/h$~\cite{Esin2018}.

The issue of experimentally detectable quantized conductance can be examined without resorting to time-independent Floquet formalism, that is, by performing direct time-dependent quantum transport simulations~\cite{Gaury2014}. Due to high computational demand, such calculations are rarely pursued, but some
attempts yield longitudinal conductance reaching close to quantized value after sufficiently long time~\cite{Fruchart2016}. This then poses a question on the accuracy of truncation procedure that is inevitably done to reduce infinite matrices in the  Floquet formalism where artifacts~\cite{Mahfouzi2012} can be introduced. One  such artifact is dc current which is not conserved (i.e., different in the left and right lead)~\cite{Kitagawa2011,Wang2003}, or insufficient number of Floquet replicas is retained for achieving converged results.

In this study, we employ {\em charge-conserving} solution~\cite{Mahfouzi2012} for the Floquet-nonequilibrium Green functions (Floquet-NEGF)~\cite{Kitagawa2011,Wang2003,Mahfouzi2012} which ensures that dc current in the left (L) and the right (R) lead  are identical at each level of truncation of matrices in the Floquet formalism, i.e., the number of of ``photon'' excitations $N_\mathrm{ph}$ retained. As an overtire, Fig.~\ref{fig:fig3}(d) demonstrates $|I_\mathrm{L}| \equiv |I_\mathrm{R}|$ at each $N_\mathrm{ph}$, as well as that dc component of current converges at  $N_\mathrm{ph}=7$. Nevertheless, the conductance in Fig.~\ref{fig:fig3}(b) remains nonquantized in both $\Delta_0$ and $\Delta_1$ gaps. We then proceed to compare spatial profiles of local current density in conventional and Floquet TIs in Fig.~\ref{fig:fig4} which offer detailed microscopic insight on how electrons propagate from one to another carbon atom as they transition from topologically trivial NM leads into the TI region, or within the TI region with possible edge or bulk vacancies introduced as disorder.

The paper is organized as follows. Section~\ref{sec:mm_a} describes the Hamiltonian of Floquet TI defined on ZGNR, as well as charge-conserving Floquet-NEGF from Ref.~\cite{Mahfouzi2012}, which is extended here to nonzero bias voltage and to obtain local current density. The same ZGNR is used in Sec.~\ref{sec:mm_b} to define  Hamiltonians for the conventional time-independent QHI, QAHI and QSHI, where we also provide steady-state NEGF expressions for  local current density in these systems. Section~\ref{sec:results_a} presents results for two-terminal conductance of these four TIs, and Sec.~\ref{sec:results_b} compares spatial profiles of local current density as it flows from the NM leads into those four TIs.  In Sec.~\ref{sec:exp}  we discuss experimental schemes to quantify  bulk vs. edge contributions to total current within Floquet TI using either a nanopore~\cite{Heerema2016,Chang2014,Chang2012} drilled in the interior of irradiated ZGNR, whose effect on the conductance is also explicitly calculated, or magnetic field imaging via diamond NV centers~\cite{Ku2020}. We conclude in Sec.~\ref{sec:conclusions}.

 \section{Models and Methods}\label{sec:models}

 \subsection{Hamiltonian and quantum transport formalism for Floquet TI}\label{sec:mm_a}
The semi-infinite leads and the scattering region in Fig.~\ref{fig:fig1} combined constitute, prior to introducing light or external magnetic field or spin-orbit coupling into the scattering region, an infinite homogeneous ZGNR described by the nearest-neighbor tight-binding Hamiltonian
\begin{equation}\label{eq:gr_ham_nir}
\hat{H}_\mathrm{ZGNR} = - \sum_{\braket{ij}} \gamma_{ij} \hat{c}^\dagger_i\hat{c}_j,
\end{equation}
Here $\braket{ij}$ indicates the sum over the nearest-neighbor sites; $\hat{c}_i^\dagger$($\hat{c}_j$) creates (annihilates) an electron on site $i$ of the honeycomb lattice hosting a single $p_z$-orbital $\langle \mathbf {r}| i \rangle=\pi(\mathbf{r}-\mathbf{R}_i)$; and \mbox{$\gamma_{ij} = \gamma = 2.7$ eV} is the  nearest-neighbor hopping from site $i$ to $j$. The width of the ZGNR is chosen as $W= 29a$, where $a$ is the distance between two  nearest-neighbor carbon atoms in graphene. It is well-known that, in general, TIs thinner than twice the width of their boundary states will experience hybridization of those states and opening of a topologically trivial mini gap~\cite{Hasan2010,Qi2011,Chang2015} at the crossing point. For Floquet TI studied in Fig.~\ref{fig:fig3}(a) this would happen if $W \le 14a$, so that our choice of $W$ evades such size artifacts. This is also ensured in the cases for QHI, QAHI and QSHI in Fig.~\ref{fig:fig4} where ZGNR is always wider than the width of edge currents. The ZGNR terminates at infinity into the macroscopic reservoirs of electrons  whose chemical potentials are $\mu_\mathrm{L} = E_{F} + eV_b/2$ and $\mu_\mathrm{R} = E_{F} - eV_b/2$ for $E_F$ as the Fermi energy and $V_b$ as the applied dc bias voltage. Note that zero-temperature two-terminal conductance $G(E_F)$ of an infinite homogeneous ZGNR described by Hamiltonian in Eq.~\eqref{eq:gr_ham_nir} is plotted for comparison in Fig.~\ref{fig:fig3}(b) and labeled as nonirradiated (NIR).

In the case of Floquet TI, circularly polarized monochromatic laser light irradiates the scattering region (shaded blue) of finite length $L = 30\sqrt{3}a$ in Fig.~\ref{fig:fig1}(a). The electromagnetic field of light is introduced into the Hamiltonian through the vector potential \mbox{$\mathbf{A}(t) =  A_0(\bold{e}_x \cos\omega t  + \bold{e}_y \sin\omega t )$}, where $\bold{e}_x$ ($\bold{e}_y$) is the unit vector along the $+x$-axis ($+y$-axis). The corresponding electric field generated by $\mathbf{A}(t)$ is \mbox{$\bold{E}(t) = -\partial \bold{A}(t)/\partial t$}. We neglect the relativistic magnetic field of light, so that electronic spin degree of freedom maintains its degeneracy and it is excluded from our analysis. The vector potential
modifies the Hamiltonian in Eq.~\eqref{eq:gr_ham_nir} via the standard Peierls substitution~\cite{Li2020}
\begin{equation}\label{eq:pphase}
 \hat{c}^\dagger_i \hat{c}_j \longmapsto \exp\bigg[ i2z(\bold{e}_x \cos\omega t  + \bold{e}_y \sin\omega t)\cdot \bold{r}_{ij}\bigg]\hat{c}^\dagger_i \hat{c}_j,
\end{equation}
which is rigorously proven~\cite{Panati2003} to be sufficient to capture the leading order effects due to the presence of the vector potential $\mathbf{A}(t)$. Here $z = eaA_0/2\hbar$ is a dimensionless measure of intensity of the circularly polarized light; $\omega$ is the frequency and $\bold{r}_{ij}$ is the position vector connecting site $i$ with site $j$. The new Hamiltonian $\hat{H}(t)$ with time-dependent hopping  between sites $i$ and $j$, $\gamma_{ij}(t) = \gamma \exp\bigg[ i2z(\bold{e}_x \cos\omega t  + \bold{e}_y \sin\omega t)\cdot \bold{r}_{ij}\bigg]$, is time-periodic, $\hat{H}(t+T)=\hat{H}(t)$, with period $T = 2\pi/\omega$.

Any solution of the Schr\"{o}dinger equation, $i\hbar \partial \Psi(t)/\partial t = \hat{H}(t) \Psi(t)$, with time-periodic Hamiltonian $\hat{H}(t)=\hat{H}(t+T)$  can be expressed as a linear combination, $\Psi(t)=\sum_\alpha c_\alpha \phi_\alpha^\mathrm{F}(t)$, of the so-called Floquet functions~\cite{Shirley1965,Sambe1973}
\begin{equation}\label{eq:eigenstates}
\phi_\alpha^\mathrm{F}(t) = e^{-i \xi_\mathrm{QE}^\alpha t/\hbar} u_\alpha(t).
\end{equation}
Here  $\xi_\mathrm{QE}^\alpha$ is the Floquet {\em quasi-energy} and \mbox{$u_\alpha(t+T)=u_\alpha(t)$} are periodic functions which can, therefore, be expanded into a Fourier series
\begin{equation}\label{eq:ufunction}
u_\alpha(\mathbf{r},t)=\sum_{n=-\infty}^{\infty} e^{i n \omega t} u_n^\alpha(\mathbf{r}).
\end{equation}
The time-periodic Hamiltonian $\hat{H}(t)=\hat{H}(t+T)$ can also be expanded into a Fourier series
\begin{equation}\label{eq:ham_fourier}
\hat{H}(t) = \sum_{n=-\infty}^{\infty} \hat{H}_n e^{i n \omega t},
\end{equation}
where $\hat{H}_n$ is given in  terms of the Bessel functions $J_{m}(z)$ of the first kind
\begin{subequations}\label{eq:jmcs}
\begin{eqnarray}\label{eq:jms}
\exp(iz\sin x) = \sum_{m=-\infty}^{\infty} J_{m}(z)e^{i m x},
\end{eqnarray}
\begin{eqnarray}\label{eq:jmc}
\exp(iz\cos x) = \sum_{m=-\infty}^{\infty} i^m J_{m}(z)e^{i m x}.
\end{eqnarray}
\end{subequations}
Using the matrix representation of the Fourier coefficients $\bold{H}_n$ in Eq.~\eqref{eq:ham_fourier} in the basis of orbitals $|i\rangle$, we construct the  Floquet Hamiltonian~\cite{Shirley1965,Sambe1973}
\begin{equation}\label{eq:floquet_ham}
\check{\bold{H}}_\mathrm{F} = \begin{pmatrix}
\ddots & \vdots & \vdots & \vdots & \iddots\\
\cdots & \bold{H}_0 & \bold{H}_1 & \bold{H}_2 & \cdots \\
\cdots & \bold{H}_{-1} & \bold{H}_0 & \bold{H}_1 & \cdots \\
\cdots & \bold{H}_{-2} & \bold{H}_{-1}& \bold{H}_0 & \cdots \\
\iddots & \vdots & \vdots & \vdots & \ddots
\end{pmatrix}.
\end{equation}
which is time-independent but infinite.

The time-dependent NEGF formalism~\cite{Gaury2014} operates with two fundamental quantities~\cite{Stefanucci2013}---the retarded $\mathbf{G}^r(t,t')$ and the lesser $\mathbf{G}^<(t,t')$ Green functions (GF)---which describe the density of available quantum states and how electrons occupy those states in nonequilibrium, respectively. They depend on two times, but solutions can be sought in other representations, such as the double-time-Fourier-transformed~\cite{Mahfouzi2012,Wang2003} GFs, ${\bf G}^{r,<}(E,E')$. In the case of periodic time-dependent Hamiltonian, they must take the form~\cite{Martinez2003}
\begin{equation}
{\bf G}^{r,<}(E,E')={\bf G}^{r,<}(E,E+n \hbar \omega)={\bf G}^{r,<}_n(E),
\end{equation}
in accord with the Floquet theorem~\cite{Shirley1965,Sambe1973}. The coupling of energies $E$ and $E+ n\hbar\omega$ ($n$ is integer) indicate  ``multiphoton'' exchange processes. In the absence of many-body (electron-electron or electron-boson)  interactions, currents can be expressed using solely the Floquet-retarded-GF $\check{\bold{G}}^r(E)$
\begin{equation}\label{eq:floquet_GF}
[E + \check{\bold{\Omega}} - \check{\bold{H}}_\mathrm{F} - \check{\bold{\Sigma}}^r(E) ]\check{\bold{G}}^r(E) = \check{\bold{1}},
\end{equation}
which is composed of ${\bf G}^r_n(E)$ submatrices along the diagonal. Here we use notation
\begin{equation}\label{eq:omega_mtx}
\check{\bold{\Omega}} = \begin{pmatrix}
\ddots & \vdots & \vdots & \vdots & \iddots\\
\cdots & -\hbar\omega\bold{1} & 0 & 0 & \cdots \\
\cdots & 0 & 0 & 0 & \cdots \\
\cdots & 0 & 0 & \hbar\omega\bold{1}& \cdots \\
\iddots & \vdots & \vdots & \vdots & \ddots
\end{pmatrix},
\end{equation}
and $\check{\bold{\Sigma}}^r(E)$  is the retarded Floquet self-energy matrix
\begin{equation}\label{eq:floquet_self}
\check{\bold{\Sigma}}^r(E) = \begin{pmatrix}
\ddots & \vdots & \vdots & \vdots & \iddots\\
\cdots & \bold{\Sigma}^r(E-\hbar\omega) & 0 & 0 & \cdots \\
\cdots & 0 & \bold{\Sigma}^r(E)  & 0 & \cdots \\
\cdots & 0 & 0 & \bold{\Sigma}^r(E+\hbar\omega)  & \cdots \\
\iddots & \vdots & \vdots & \vdots & \ddots
\end{pmatrix},
\end{equation}
composed of the usual self-energies of the  leads~\cite{Velev2004}, $\bold{\Sigma}^r(E) = \sum_{p=\mathrm{L,R}}\bold{\Sigma}_p^r(E)$, on the diagonal. All matrices labeled as $\check{\mathbf{O}}$ are representations of operators acting in the Floquet-Sambe~\cite{Sambe1973} space, $\mathcal{H}_\mathrm{F} =  \mathcal{H}_T \otimes \mathcal{H}_e$, where $\mathcal{H}_e$ is the Hilbert space of electronic states spanned by localized orbitals $|i \rangle$ and $\mathcal{H}_T$ is the Hilbert space of periodic functions with period $T=2\pi/\omega$ spanned by orthonormal Fourier vectors $\langle t|n \rangle = \exp(i n \omega t)$.

The charge current $I_{p}(t)$ in the lead $p=\mathrm{L,R}$ is time-dependent  due to Eq.~\eqref{eq:pphase}, and it also has periodicity   $T=2\pi/\omega$ like the Hamiltonian itself. The dc component of current, either due to pumping by time-dependent potential~\cite{Brouwer1998,Moskalets2002,FoaTorres2005,Bajpai2019}
or due to applied bias voltage $V_b$ or both, is given by
\begin{equation}
I_{p}  = \frac{1}{T}\int_{t}^{t+T}  I_p(t') dt'.
\end{equation}
Such dc component, or time-averaged current over one period $T$, that is injected into the lead $p$ is obtained from the  following NEGF expression~\cite{Mahfouzi2012}
\begin{equation}\label{eq:floq_chj}
I_{p} = \frac{e}{2N_\mathrm{ph}}\int\limits_{-\infty}^{+\infty}dE\, \mathrm{Tr}[\check{\bold{\Gamma}}_p\check{\bold{f}}_p\check{\bold{G}}^r\check{\bold{\Gamma}}\check{\bold{G}}^a - \sum_{p=\mathrm{L,R}}\check{\bold{\Gamma}}_p\check{\bold{G}}^r\check{\bold{\Gamma}}_\alpha\check{\bold{f}}_\alpha\check{\bold{G}}^a ].
\end{equation}
In our convention, $I_p>0$ indicates that charge current is flowing into the lead. Here $\check{\bold{f}}_p$ is the Floquet Fermi matrix
\begin{equation}\label{eq:fermi_mtx}
\check{\bold{f}}_p(E) = \begin{pmatrix}
\ddots & \vdots & \vdots & \vdots & \iddots\\
\cdots & f_p(E-\hbar\omega)\bold{1} & 0 & 0 & \cdots \\
\cdots & 0 & f_p(E) & 0 & \cdots \\
\cdots & 0 & 0 & f_p(E+\hbar\omega)\bold{1}& \cdots \\
\iddots & \vdots & \vdots & \vdots & \ddots
\end{pmatrix},
\end{equation}
where $f_p(E)$ is the Fermi function of the macroscopic particle reservoir attached to  lead $p$; \mbox{$\check{\bold{\Gamma}}_p(E) = i [\check{\bold{\Sigma}}_p^r(E) - (\check{\bold{\Sigma}}_p^r(E))^\dagger]$} is the Floquet level broadening matrix; \mbox{$\check{\bold{\Gamma}}(E) = \sum_{p=\mathrm{L,R}} \check{\bold{\Gamma}}_p(E)$}; the Floquet-advanced-GF is defined as $\check{\bold{G}}^a(E) = [\check{\bold{G}}^r(E)]^\dagger$; and $\mathbf{1}$ is the unit matrix in $\mathcal{H}_e$ space. We note that Eq.~\eqref{eq:floq_chj} is generalization of the expression for charge current in Ref.~\cite{Mahfouzi2012} to include the applied bias voltage $V_b$. The linear-response two-terminal conductance  is then given by
\begin{equation} \label{eq:conductance}
G = \frac{I_\mathrm{R}}{V_b},
\end{equation}
for small applied bias voltage $eV_b \ll E_F$.

While the space $\mathcal{H}_e$ is finite-dimensional, with dimension equal to the number of sites $N_e$ within the scattering region, the space $\mathcal{H}_T$ is infinite-dimensional and has to be truncated using $|n| \le N_\mathrm{ph}$. For truncation we employ the following convergence criterion
\begin{equation}\label{eq:convergence}
\bigg| \frac{I_p(N_\mathrm{ph}) - I_p(N_\mathrm{ph}-1)}{I_p(N_\mathrm{ph}-1)}\bigg|\times 100  < \delta,
\end{equation}
where $\delta$ is the convergence tolerance. Since the operators acting in $\mathcal{H}_e$ are represented by matrices of dimension $N_e \times N_e$, the operators $\check{\mathbf{O}}$ acting on the truncated Floquet-Sambe space $\mathcal{H}_\mathrm{F}$ are represented by matrices of dimension $(2N_\mathrm{ph}+1)N_e \times (2N_\mathrm{ph}+1)N_e$. Note that the trace in Eq.~\eqref{eq:floq_chj}, $\mathrm{Tr} \equiv  \mathrm{Tr}_e \mathrm{Tr}_T$, is summing over contributions from different subspaces of $\mathcal{H}_T$ so that the denominator includes $2N_\mathrm{ph}$ to avoid double counting. The part of the trace operating in $\mathcal{H}_T$ space ensures that at each chosen truncation $N_\mathrm{ph}$ of Floquet replicas charge current is conserved, $I_\mathrm{L} = - I_\mathrm{R}$, unlike other types of  solutions~\cite{Wang2003,Kitagawa2011} of the Floquet-NEGF equations where current conservation is ensured only in the limit $N_\mathrm{ph} \rightarrow \infty$.

The bond current operator~\cite{Nikolic2006} between sites $i$ and $j$ is time-dependent due to Eq.~\eqref{eq:pphase} and it is given by~\cite{Gaury2014}
\begin{equation}\label{eq:bc_oper}
\begin{split}
\bold{J}_{ij}(t) &= \frac{e}{i\hbar} [\gamma_{ij}(t)\hat{c}_i^\dagger\hat{c}_j - \gamma_{ji}(t)\hat{c}_j^\dagger\hat{c}_i] \\
&= \sum_{n=-\infty}^{\infty} \bold{J}_{n}^{ij}e^{in\omega t}.
\end{split}
\end{equation}
We define the Floquet bond current matrix as
\begin{equation}\label{eq:floquet_curr}
\check{\bold{J}}_{ij} = \begin{pmatrix}
\ddots & \vdots & \vdots & \vdots & \iddots\\
\cdots & \bold{J}_0^{ij} & \bold{J}_{-1}^{ij} & \bold{J}_{-2}^{ij} & \cdots \\
\cdots & \bold{J}_{1}^{ij} & \bold{J}_0^{ij} & \bold{J}_{-1}^{ij} & \cdots \\
\cdots & \bold{J}_{2}^{ij} & \bold{J}_{1}^{ij}& \bold{J}_0^{ij} & \cdots \\
\iddots & \vdots & \vdots & \vdots & \ddots
\end{pmatrix},
\end{equation}
which yields nonequilibrium part~\cite{Nikolic2006} of dc bond (or local) charge current flowing between site $i$ and $j$ as
\begin{equation}\label{eq:floquet_bc}
J_{ij}^\mathrm{neq} =  \frac{1}{2\pi i}\sum_{n=-N_\mathrm{ph}}^{N_\mathrm{ph}} \int_{E_F+n\hbar\omega-eV_b/2}^{E_F+n\hbar\omega + eV_b/2} dE\, \mathrm{Tr}[\check{\bold{G}}^<(E)\check{\bold{J}}_{ij}],
\end{equation}
where $\check{\bold{G}}^<(E) = \sum_{p=\mathrm{L,R}} i\check{\bold{G}}^r(E)\check{\bold{\Gamma}}_p(E)\check{\bold{f}}_p(E)\check{\bold{G}}^a(E)$.

\subsection{Hamiltonian and quantum transport formalism for QHI, QAHI and QSHI}\label{sec:mm_b}

The two-terminal setup in Fig.~\ref{fig:fig1}(b) hosts one of the three conventional time-independent TIs as the scattering region (shaded green) of finite length $L = 30\sqrt{3}a$. The QHI is realized by applying an external time-independent magnetic field perpendicular to ZGNR. The magnetic field is described  by a static vector potential \mbox{$\bold{A} = (By,0,0)$} in the Landau gauge, which is then included into the Hamiltonian in Eq.~\eqref{eq:gr_ham_nir} via the Peierls substitution~\cite{Li2020,Panati2003}
\begin{equation}\label{eq:qh}
 \hat{c}^\dagger_i \hat{c}_j \longmapsto \exp\bigg[ i \frac{\beta}{a_0^2}(x_i-x_j)(y_i+y_j)\bigg]\hat{c}^\dagger_i \hat{c}_j.
\end{equation}
Here $(x_i,y_i)$ indicates the position vector of a carbon atom at site $i$, and $\beta = eBa_0^2/\sqrt{3}\hbar \approx 0.07$ is a dimensionless measure of the magnetic field strength.

The QAHI~\cite{Liu2016} is described by the Haldane model~\cite{Haldane2017,Haldane1988} on the honeycomb lattice
\begin{eqnarray}\label{eq:qah}
\hat{H}_\mathrm{QAHI} & = & \sum_{\braket{ij}}- \gamma_{ij}  \hat{c}^\dagger_i\hat{c}_j + \sum_{\braket{\braket{ij}}} \widetilde{\gamma}_{ij}\hat{c}^\dagger_i\hat{c}_j  \nonumber \\
&& + \sum_{i \in A} m \hat{c}_i^\dagger\hat{c}_i + \sum_{i \in B} (-m) \hat{c}_i^\dagger\hat{c}_i.
\end{eqnarray}
Here $\braket{\braket{ij}}$ indicates the sum over the next-nearest-neighbor sites, and $\widetilde{\gamma}_{ij} = - \widetilde{\gamma}_{ji} = i \widetilde{\gamma}$ where we use $\tilde{\gamma}=0.14 \gamma$. The last two mass terms on the right hand side have different sign on the triangular sublattices A and B of the honeycomb lattice, where $m=0.2 \gamma$ specifies the ``mass'' term. Note that circularly polarized light employed in Eq.~\eqref{eq:pphase} is mandatory for Floquet TI to mimic QAHI phase of the Haldane model. In contrast, linearly polarized light, which is made of equal superposition of clockwise and anticlockwise circular polarizations, does not break time-reversal symmetry and cannot lead to Haldane ``mass'' term.

Finally, the QSHI is described by the Kane-Mele model~\cite{Kane2005}
\begin{equation}\label{eq:kane_mele_ham}
\hat{H}_\mathrm{QSHI} = \sum_{\braket{ij}}-\gamma_{ij} \bold{c}_i^\dagger \bold{c}_j + \sum_{\braket{\braket{ij}}} it_\mathrm{SO} \bold{c}_i^\dagger \bm{\sigma} \cdot (\bold{d}_{kj}\times \bold{d}_{ik} )\bold{c}_j,
\end{equation}
whose edge states crossing the topological nontrivial band gap are both chiral and spin-polarized~\cite{Hasan2010,Qi2011}. Here $\bold{c}_i^\dagger = (\hat{c}_{i\uparrow}^\dagger,\hat{c}_{j\downarrow}^\dagger)$ is a row vector of creation operators $\hat{c}^\dagger_{i\sigma}$ that create an electron on site $i$ with spin $\sigma=\uparrow, \downarrow$; $\bold{c}_i$ is the corresponding column vector of annihilation operators; $\bold{d}_{ik}$ is the unit vector pointing from site $k$ to $i$; $\bm{\sigma} = (\hat{\sigma}_x, \hat{\sigma}_y, \hat{\sigma}_z)$ is the vector of the Pauli matrices; and $t_\mathrm{SO}$ is the strength of the intrinsic spin-orbit coupling~\cite{Kane2005,Chang2014}.

The zero-temperature two-terminal conductance $G(E_F)=G_Q \mathcal{T}(E)$ of the setup in Fig.~\ref{fig:fig1}(b) is calculated using the Landauer transmission function~\cite{Imry1999,Stefanucci2013}
\begin{equation}
\mathcal{T}(E) = \mathrm{Tr}[\bm{\Gamma}_\mathrm{R}(E)\bold{G}^r(E)\bm{\Gamma}_\mathrm{L}(E)\bold{G}^a(E)],
\end{equation}
where the conductance quantum is $G_Q=2e^2/h$ for QHI and QAHI and  $G_Q=e^2/h$ for QSHI. Here the retarded GF of the scattering region is given by \mbox{$\bold{G}^r(E) = [E - \bold{H} - \bm{\Sigma}^r(E)]^{-1}$}; the advanced GF is $\bold{G}^a(E) = [\bold{G}^r(E)]^\dagger$; and $\bold{\Gamma}_p(E) = i [\bm{\Sigma}^r_p(E) - \bm{\Sigma}^a_p(E)]$ are the level-broadening matrices. To compute the nonequilibrium bond current between sites $i$ and $j$ we use~\cite{Mahfouzi2013}
\begin{equation}\label{eq:loc_j_ss}
J_{ij}^\mathrm{neq} = \frac{eV_b}{2\pi} \Tr[\bold{G}^r(E_F)\bm{\Gamma}_\mathrm{L}(E_F)\bold{G}^a(E_F) \bold{J}_{ij}],
\end{equation}
where $\bold{J}_{ij}$ is the bond current operator in Eq.~\eqref{eq:bc_oper} but with time-independent hopping $\gamma_{ij}$.

\begin{figure*}
	\includegraphics[width=\linewidth]{fig4.pdf}
	\caption{Spatial profiles of local current density in two-terminal devices of Fig.~\ref{fig:fig1} where the scattering region (dotted rectangle) of finite length is: (a) irradiated pristine ZGNR hosting Floquet TI; (b) irradiated edge-disordered ZGNR; (c) pristine QHI; (d) edge-disordered QHI; (e) pristine QAHI; (f) edge-disordered QAHI; (g) pristine QSHI; and (h) edge-disordered QSHI.  In panels (a) and (b)  we use $\hbar\omega = 3\gamma$, $z=0.5$, $N_\mathrm{ph} = 7$, and $E_F = \hbar \omega /2$ corresponding to the middle of $\Delta_1$ gap in Fig.~\ref{fig:fig3}(a). In panels (c)--(h), \mbox{$E_F = 0.2 \gamma$}, and in (g) and (h) \mbox{$t_\mathrm{SO} = 0.1 \gamma$}. The black solid arrows are guide to the eye to indicate the spatial region with large flux  of local current density.}
	\label{fig:fig4}
\end{figure*}

\section{Results and Discussion}\label{sec:results}
\subsection{Conductance within the topological gap: FTI vs. conventional TIs }\label{sec:results_a}

By diagonalizing $\check{\mathbf{H}}_\mathrm{F}$ in Eq.~\eqref{eq:floquet_ham} for an infinite ZGNR that is periodic along the $x$-axis and irradiated by circularly polarized light over its {\em whole} length, we obtain the quasi-energy spectrum $\xi_\mathrm{QE}(k_x)$ shown in Fig.~\ref{fig:fig3}(a). The chiral edge states crossing the  light-induced gap $\Delta_0$ at $\xi_\mathrm{QE} = 0$ (shaded yellow) and $\Delta_1$ at \mbox{$\xi_\mathrm{QE} = \pm \hbar\omega/2$} (shaded red) suggest naively that upon applying small bias voltage the zero-temperature linear-response two-terminal conductance in Eq.~\eqref{eq:conductance} should be quantized: \mbox{$G(E_F)=2e^2/h$} for $E_F$ within $\Delta_0$ gap; and \mbox{$G(E_F)=4e^2/h$} for $E_F$ within $\Delta_1$ gap due to one or two spin-degenerate conduction channels provided by the edge states, respectively. This is in analogy with chiral edge states of conventional time-independent TIs and their quantized conductance in Fig.~\ref{fig:fig2}. In contrast, the average conductance in Fig.~\ref{fig:fig3}(b) is $G(E_F) \approx 0.73 \times 2e^2/h$ within $\Delta_0$ gap and $G(E_F) \approx 1.87  \times 2e^2/h$ within $\Delta_1$ gap. We emphasize that these results are not an artifact of truncation of the Floquet Hamiltonian $\check{\mathbf{H}}_\mathrm{F}$ in Eq.~\eqref{eq:floquet_ham} since current in the L  and R lead in Fig.~\ref{fig:fig3}(d) converge  at $N_\mathrm{ph}=7$ using $\delta = 1\%$ criterion in Eq.~\eqref{eq:convergence}. Also,  our Floquet-NEGF formalism~\cite{Mahfouzi2012} ensures $|I_\mathrm{L}| \equiv | I_\mathrm{R}|$ in  Fig.~\ref{fig:fig3}(d) at each chosen $N_\mathrm{ph}$.

We additionally plot the total density  of states (DOS) $D(E) = \sum_j D_j(E)$ in Fig.~\ref{fig:fig3}(c) which is nonzero within the gaps $\Delta_0$ and $\Delta_1$ due to contributions from the local DOS (LDOS) $D_j(E)$ originating [Fig.~\ref{fig:fig6}(a)] from both edges and bulk of ZGNR. The LDOS is extracted from the Floquet-retarded-GF in Eq.~\eqref{eq:floquet_GF} using
\begin{equation}\label{eq:ldos}
D_j (E) = \frac{i}{2\pi} \bra{j} \mathrm{Tr}_T [\check{\bold{G}}^r(E) - \check{\bold{G}}^a(E)] \ket{j},
\end{equation}
where $\mathrm{Tr}_T$ is the partial trace over states in $\mathcal{H}_T$. The issue of positivity of DOS and LDOS obtained from the Floquet-retarded-GF  has been discussed extensively~\cite{Rudner2020,Uhrig2019}.

Even though Floquet TI in irradiated ZGNR does not exhibit quantized conductance plateau in Fig.~\ref{fig:fig3}(b), its conductance is largely insensitive to ED. For example, $G(E_F)$ is reduced by $\sim 2$\% within $\Delta_0$ gap and by $\sim 15$\% within $\Delta_1$ gap upon introducing edge vacancies. Nevertheless, this is still less robust than conventional time-independent TIs whose conductance within the topologically nontrivial band gap is completely insensitive to ED, as shown in Figs.~\ref{fig:fig2}(e)--(g). The disorder is introduced by removing carbon atoms on the top and bottom edges of the scattering region, as illustrated in Fig.~\ref{fig:fig4}(b), while imposing the following conditions: (\emph{i}) ED introduced in NIR ZGNR leads to complete conductance suppression  $G(E_F) \rightarrow 0$ within the same energy interval defined by $\Delta_0$ gap; ({\em ii}) ED preserves the left-right symmetry of the device, so that charge pumping is absent when the ED ZGNR is irradiated with circularly polarized  light~\cite{FoaTorres2005,Bajpai2019,FoaTorres2014} in the absence of dc bias voltage $V_b=0$.

Note that in the case of vacancies at the edges of QSHI, our tight-binding Hamiltonian in Eq.~\eqref{eq:kane_mele_ham} does not capture possibility of formation of a localized magnetic moment at the vacancy site which requires first-principles Hamiltonians~\cite{Vannucci2020}. This opens a possibility of backscattering involving spin flip which will disrupt~\cite{Vannucci2020} (nearly) quantized conductance in Fig.~\ref{fig:fig2}(g).

\subsection{Spatial profiles of local current density: Floquet TI vs. QHI, QAHI and QSHI}\label{sec:results_b}

The spatial profiles of local current density, i.e., bond current  $J_{ij}^\mathrm{neq}$ defined in Eqs.~\eqref{eq:floquet_bc} and ~\eqref{eq:loc_j_ss} for Floquet TI and conventional time-independent TIs, respectively, allows us to visualize how electrons transition from topologically trivial NM leads into chiral edge states within the TI region. Figures~\ref{fig:fig4}(c), ~\ref{fig:fig4}(e) and ~\ref{fig:fig4}(g) shows that bulk states contribute to current density within the leads, but current density becomes confined to narrow flux near the  edges of QHI, QAHI and QSHI. The width of the flux corresponds to spatial extent of the edge state. As expected due to chirality of edge states, current flows only along the top edge in QHI and QAHI, while in QSHI it flows on both the top and bottom edges~\cite{Nowack2013}. This is because boundaries of QSHI host a pair of spin-polarized edge states~\cite{Kane2005}, so that on the top edge  electrons with spin $\sigma=\uparrow$ move from left to right while at the bottom edge electrons with spin $\sigma=\downarrow$ move from left to right. Upon introducing ED in Figs.~\ref{fig:fig4}(d), ~\ref{fig:fig4}(f) and ~\ref{fig:fig4}(h), topological protection and quantized transport through edge states manifest by local current density circumventing the disorder since any backscattering would require to cross over to  the other edge which is forbidden due to the absence of bulk states~\cite{Buttiker1988}.

In contrast, local current density is nonzero within the whole Floquet TI in Fig.~\ref{fig:fig4}(a), with larger flux near the edges [Fig.~\ref{fig:fig5}(a)]. Upon introducing ED, the edge flux circumvents the disorder but due to general nonlocality of quantum transport bulk flux is also reduced [Fig.~\ref{fig:fig5}(b)] which explains slight reduction of conductance in Fig.~\ref{fig:fig3}(b) within gaps $\Delta_0$ and $\Delta_1$.

\begin{figure}
	\includegraphics[width=\linewidth]{fig5.pdf}
	\caption{Spatial profile of local current density in Figs.~\ref{fig:fig4}(a) and ~\ref{fig:fig4}(b) over the transverse cross section within: (a) pristine Floquet TI; or (b) Floquet TI with edge disorder. The position of the transverse cross section is marked by dashed vertical line in Figs.~\ref{fig:fig4}(a) and ~\ref{fig:fig4}(b), respectively. The horizontal dashed line in both panels marks the extent of the edge state.}
	\label{fig:fig5}
\end{figure}

Interestingly, SQUID-based imaging of QSHI made from HgTe quantum wells has found that gate tuning of bulk conductivity can lead to transport regime where edge and bulk local current densities {\em coexist}~\cite{Nowack2013}. The trace of local current density is scanned by detecting its magnetic field produced according to the  Biot-Savart law, which is possible even through the top gate employed to tune the carrier density. In this regime, experimental images were analyzed to quantify contribution of edge and bulk local currents to the total current. We perform similar analysis in Fig.~\ref{fig:fig5} which shows that in pristine Floquet TI from Fig.~\ref{fig:fig4}(a), edge current contributes 44\% and bulk current contributes 56\% to the total current over the transverse cross section [marked by dashed line in Fig.~\ref{fig:fig4}(a)]. Conversely, in the presence of edge disorder in  Fig.~\ref{fig:fig4}(b), edge current contributes 52\% and bulk current contributes 48\% to the total current over  the same transverse cross section.

\subsection{Proposed experimental schemes for probing edge vs. bulk transport within Floquet TI: Graphene nanopore and magnetic field imaging}\label{sec:exp}

\begin{figure}
	\includegraphics[width=8.5cm]{fig6.pdf}
	\caption{(a) The LDOS [Eq.~\eqref{eq:ldos}] evaluated at $E = \hbar\omega/2$ in the center of $\Delta_1$ gap in Fig.~\ref{fig:fig3}(a) for irradiated pristine ZGNR. (b) The LDOS evaluated at $E=\hbar\omega/2$  for irradiated ZGNR with a nanopore drilled in the interior of the nanoribbon. (c) Time-averaged local bond current [Eq.~\eqref{eq:floquet_bc}] in the same irradiated ZGNR with a nanopore as in (b). (d) Zero-temperature two-terminal conductance $G(E_F)$ of irradiated ZGNR with a nanopore (orange line) vs. conductance of irradiated pristine ZGNR (blue line) within the gap $\Delta_1$ in Fig.~\ref{fig:fig3}(a). The former is reduced by $\sim 28\%$ with respect to the latter. In all panels we use $\hbar\omega = 3\gamma$, $z=0.5$ and $N_\mathrm{ph} = 7$.}
	\label{fig:fig6}
\end{figure}

The spatial profiles of local current density of conventional time-independent TIs in Figs.~\ref{fig:fig4}(c)--(h) indicate that any disorder introduced in the interior of ZGNR will have no effect on the two-terminal conductance in Fig.~\ref{fig:fig2}. This was explicitly demonstrated in Ref.~\cite{Chang2014} for the case of QSHI (based on graphene plus heavy adatoms). Therefore, we propose to employ a nanopore in the ZGNR interior as the simplest technique that can detect the presence of bulk current density in Figs.~\ref{fig:fig4}(a) and ~\ref{fig:fig4}(b) in the case of Floquet TI. We introduce nanopore in Figs.~\ref{fig:fig6}(b) and ~\ref{fig:fig6}(c) in  such a way that it preserves the left-right symmetry of the device in order to avoid any charge pumping by time-dependent potential of light~\cite{FoaTorres2005,Bajpai2019,FoaTorres2014}. In experiments, nanopores are routinely drilled, without disrupting the surrounding honeycomb lattice of graphene, for applications like DNA sequencing~\cite{Heerema2016}, and they could also be deployed to block phonon transport in thermoelectric applications~\cite{Chang2014,Chang2012}.  Figures~\ref{fig:fig6}(a) and ~\ref{fig:fig6}(b) confirms that nanopore does not impair high LDOS [Eq.~\eqref{eq:ldos}] near the edges of the Floquet TI, which correspond to chiral edge states from Fig.~\ref{fig:fig3}(a).  Figure~\ref{fig:fig6}(c) shows that local transport in the presence of nanopore utilizes {\em both} left-to-right moving chiral edge states and bulk states. Since electrons flowing through the bulk  are backscattered by the nanopore, presence of nanopore reduces conductance by about $\simeq 28$\% in Fig.~\ref{fig:fig6}(d) within  the gap $\Delta_1$.

A more detailed probing of edge vs. bulk transport in  $\sim\mu$m-sized devices, such as those employed in recent experiments~\cite{McIver2020} to convert graphene into Floquet TI, could be achieved using diamond NV centers. The device can be fabricated on a diamond containing high-density, near-surface NV  ensemble~\cite{Tetienne2017,Ku2020}, along with a graphite top gate separated by hexagonal boron nitride to tune the carrier density~\cite{Ku2020}. The spin state of NV centers, which serve as the sensor of magnetic field produced by local current density, can be optically initialized and readout via imaging the NV photoluminescence onto a camera. Such a setup has the advantages of being able to operate over a wide range of temperatures, from cryogenic to room temperature (e.g., the experiment in Ref.~\cite{McIver2020} was done at 80 K); it can be readily integrated with an optical cryostat necessary for experiments involving THz radiation; and it has less stringent vibrational requirement compared to scanning setups. We note that THz radiation is far detuned from any of the NV orbital/spin transitions and hence it will not affect NV centers at all.

In Ref.~\cite{McIver2020}, a constant DC bias generates a current \mbox{$I \simeq 125$ $\mu$A}, and THz pulses drive the system into Floquet TI state for about \mbox{$3$ ps} at \mbox{$\sim 210$ kHz} repetition rate. Hence, one has a time-averaged typical current density \mbox{$\bar{J}_{\mathrm{F}}\sim 80$ pA/$\mu$m} in a 1-$\mu$m-wide device. This corresponds to a typical stray magnetic field \mbox{$\mu_0 \bar{J}\sim 0.1$ nT}, where $\mu_0$ is the permeability of free space. While it is a small field, its measurement is attainable with existing NV sensing technologies. For example, a single NV can sense $\sim$ nT field with a Hahn-echo sequence over 100 seconds signal averaging at room temperature~\cite{Maze2008}. Detection of \mbox{$\sim 0.1$ nT} field is attainable in combination with entanglement-assisted repetitive readout~\cite{Neumann2010,Lovchinsky2016}, which is available for NV ensemble, as well as with enhanced coherence at cryogenic temperatures and with dynamical decoupling sequences~\cite{BarGill2013}. One can measure the differential current density \mbox{$\Delta J(x,y)\equiv J_{\mathrm{FTI}}(x,y)-J_{\mathrm{NIR}}(x,y)$}, where $J_{\mathrm{FTI}}(x,y)$ [$J_{\mathrm{NIR}}(x,y)$] is current density within the Floquet TI [nonirradiated normal phase], by pulsing on the THz radiation during one free precession time of the Hahn-echo and keeping the THz drive off during the other free precession. The current density $J_{\mathrm{NIR}}(x,y)$ can be measured separately in a Hahn-echo measurement without any THz pulses to enable one to extract $J_{\mathrm{FTI}}(x,y)$. Diffraction-limited optical imaging of magnetic field has resolution $\sim$ 400 nm~\cite{Ku2020}, which is enough to resolve edge currents separated by a width of 1 $\mu$m. With further improvement in spatial resolution, we anticipate that $\sim$ 10 nm resolution can be achieved by using the Fourier gradient imaging~\cite{Arai2015}.

\section{Conclusions}\label{sec:conclusions}

In conclusion, using steady-state NEGF formalism applied to two-terminal Landauer-B\"{u}ttiker setup [Fig.~\ref{fig:fig1}(b)] with scattering region consisting of conventional time-independent TIs---such as QHI, QAHI and QSHI defined on graphene nanoribbon in order to generate chiral edge states of {\em finite} length---we demonstrate that their conductance is never perfectly quantized [Fig.~\ref{fig:fig2}]. This is  due to backscattering at the NM-lead/2D-TI interface. Nevertheless, it remains very close to perfect plateau at $2e^2/h$ within the topologically nontrivial band gap, and it is completely insensitive to edge disorder. The spatial profiles of local current density visualize how electrons flow from bulk states within topologically trivial NM leads into the narrow flux defined by edge states within the TI region, while circumventing any edge disorder within the TI region.

In contrast, when the scattering region is converted into the Floquet TI by irradiating graphene nanoribbon by circularly polarized light, conductance within light-induced topologically nontrivial band gaps is not quantized, but it changes little with edge disorder. This results confirm previous findings in the literature~\cite{FoaTorres2014,Farrell2016} while ensuring proper convergence and charge current conservation in the solution of Floquet-NEGF equations~\cite{Mahfouzi2012}. Furthermore, we employ such charge-conserving  Floquet-NEGF formalism to compute spatial profiles of local current density. They are higher along the edges [Fig.~\ref{fig:fig5}(a)], following high LDOS near the edges  [Fig.~\ref{fig:fig6}(a)], but they remains nonzero also in the interior of the Floquet TI [Fig.~\ref{fig:fig4}(a)]. Such spatial profiles make it also possible to refine  previous qualitative estimates of edge vs. bulk contribution to current through Floquet TI~\cite{Esin2018} with precise measure from Figs.~\ref{fig:fig4}(a), ~\ref{fig:fig4}(b) and ~\ref{fig:fig5} which show that edge currents and bulk currents contribute nearly equally to the total current. Thus, observing quantized transport in Floquet TI would require to minimize coupling to bulk states.

We propose a very simple experimental technique to detect presence or absence of bulk states in quantum transport through Floquet TI---conductance measurements under laser irradiation should be performed using uniform graphene flake, as in the very recent experiments~\cite{McIver2020}, as well as repeated after a nanopore~\cite{Heerema2016} is drilled in the interior of the flake. If local current density is nonzero in the bulk, it will be scattered by the nanopore which leads to $\simeq 28$\% reduction [Fig.~\ref{fig:fig6}(d)] of the two-terminal conductance when compared to graphene nanoribbon without the nanopore. Finally, we delineate more sophisticated experimental schemes for direct imaging~\cite{Ku2020} of magnetic field produced by edge and bulk local current densities based on diamond NV centers whose orbital/spin transitions are far detuned from THz radiation employed~\cite{McIver2020} in recent experiments to convert graphene into Floquet TI.

\begin{acknowledgments}
We thank L. E. F. Foa Torres, J. Gong and L. Vannucci  for insightful discussions. This research was supported by the US National Science Foundation (NSF) under Grant No. ECCS 1922689.
\end{acknowledgments}


\section{Introduction}

The discovery of superconductivity in doped nickel oxides Nd$_{0.8}$Sr$_{0.2}$NiO$_2$~\cite{Li_2019,Sawatzky19} has attracted intensive interests
both in experiment~\cite{Hepting_2020,Q_Li_2020,Zhou_2020,Fu_arXiv,Lee_2020,D_Li_2020,Zeng_arXiv,Goodge_arXiv,BX_Wang_2020,Q_Gu_arXiv,Osada_2020} and theory~\cite{Hepting_2020,Botana_2020,Sakakibara_2020,Hirsch_2019,
Nomura_2019,Hirayama_2020,Gao_arXiv,Singh_arXiv,
Jiang_2020,Ryee_2020,HuZhang_2020,GM_Zhang_2020,Z_Liu_2020,
Wu_2020,Been_arXiv,Lang_arXiv,Leonov_2020,Leonov_arXiv_2,
Werner_2019,Petocchi_arXiv,Y_Gu_2020,Liang_Si_2020,Lechermann_2020,Lechermann_arXiv,Karp_2020,Kitatani_arXiv,Y_Wang_arXiv
Zhang_2020,LH_Hu_2019,Chang_arXiv,
Z_Wang_arXiv,
P_Jiang_2019,Choi_2020,
Geisler_2020,He_2020,Bernardini_2020c,
Talantsev_2020,T_Zhou_2020,Bernardini_2020,
Bernardini_2020b,Olevano_2020,Choi_arXiv,Adhikary_arXiv,Nica_arXiv},
because the nickelate might be an analog of the well known high-$T_{\rm c}$ superconductor, cuprates.
Recently, the doping dependence has been explored both theoretically~\cite{Kitatani_arXiv} and experimentally~\cite{D_Li_2020,Zeng_arXiv},
and the presence of the superconducting dome has been confirmed~\cite{D_Li_2020,Zeng_arXiv}.
The maximum superconducting transition temperature $T_{\rm c}$ is about 15 K, not very high compared to that of the high-$T_{\rm c}$ cuprates.
However, because the Bardeen-Cooper-Schrieffer (BCS) phonon mechanism cannot explain the observed $T_{\rm c}$~\cite{Nomura_2019}, the superconducting mechanism is most likely unconventional, in which the electron correlations play an important role~\cite{Sakakibara_2020,Wu_2020,Kitatani_arXiv,Adhikary_arXiv}.
A recent observation of $d$-wave like superconducting gap also supports this scenario~\cite{Q_Gu_arXiv}.
Here, a natural question arises: is there any possibility to realize $T_{\rm c}$ as high as the cuprates in nickelates?

In the cuprates, the superconductivity emerges by doping carriers into the antiferromagnetic Mott insulator having a large magnetic exchange coupling $J$ ($\sim$130 meV)~\cite{Lee_Nagaosa_Wen_2006}.
One of the reasons for the large $J$ is because the cuprates belong to the charge-transfer type in the Zaanen-Sawatzky-Allen diagram~\cite{Zaanen_1985}, and the charge-transfer energy $\Delta_{dp}$ (the energy difference between the copper $3d$ and oxygen 2$p$ orbitals) is small among transition metal oxides.
Although the mechanism of the high-$T_{\rm c}$ superconductivity is highly controversial, the large $J$ is a plausible factor in enhancing the $d$-wave superconductivity in the cuprates~\cite{Ogata_2008}.
This large value of $J$ is certainly a characteristic feature of the cuprates, which makes the cuprates very different from other transition metal oxides.

On the other hand, in the case of the nickelate NdNiO$_2$, $\Delta_{dp}$ is larger than that of the cuprates~\cite{Lee_2004}.
Thus, naively, we expect smaller $J$ for nickelates.
Indeed, a recent experimental estimate using the Raman spectroscopy gives $J = 25$ meV~\cite{Fu_arXiv}.
However, it should be noted that the origin of small $J$ in NdNiO$_2$ may be ascribed to another notable difference from the cuprates, namely, NdNiO$_2$ is not a Mott-insulator due to the self-doping effect.
In NdNiO$_2$, orbitals in the Nd layer form extra Fermi pockets on top of the large Fermi surface formed by the Ni 3$d_{x^2-y^2}$\xspace orbital, and the Ni 3$d_{x^2-y^2}$\xspace orbital is hole-doped,
i.e., the filling of the Ni $3d$ orbitals deviates from $d^9$~\cite{Lee_2004,Botana_2020,Sakakibara_2020,Gao_arXiv,Liang_Si_2020,Olevano_2020}.
The self-doping naturally explains the absence of Mott-insulating behavior in NdNiO$_2$.
Although it has been shown that the Ni 3$d_{x^2-y^2}$\xspace orbital forms a two-dimensional strongly-correlated system~\cite{Nomura_2019}, $J$ at the $d^9$ configuration with half-filled $d_{x^2-y^2}$\xspace orbital is masked by the self-doping.
The experimental estimate should be understood as the $J$ value including the effect of the self-doping, not the $J$ value at the ideal $d^9$ configuration.
One of the reasons for the controversy in theory about the size of $J$~\cite{Jiang_2020,Ryee_2020,HuZhang_2020,GM_Zhang_2020,Z_Liu_2020,Wu_2020,Been_arXiv,Lang_arXiv,Leonov_2020,Leonov_arXiv_2} is ascribed to the ambiguity in calculating $J$ (whether we calculate $J$ at $d^9$ filling or $J$ including the self-doping effect).
In any case, it is a non-trivial problem whether we can justify the mapping onto a simple spin model to understand the property of NdNiO$_2$.
This fact makes NdNiO$_2$ an imperfect analog of the cuprates.

Recently, there was a proposal to design cuprate-analog nickelates without the complication of the self-doping~\cite{Hirayama_2020} \footnote{See also Refs.~\cite{Bernardini_2020b,Nica_arXiv} for other attempts to find nickelate superconductors.}.
Since NdNiO$_2$ is a layered material, one can systematically propose nickelate family materials by changing the composition of the ``block-layer" \cite{Tokura_1990} between NiO$_2$ layers.
Proposed dynamically stable nickelates have smaller Fermi pockets of the block-layer orbitals than NdNiO$_2$.
In some materials, the self-doping is completely suppressed, and the ideal $d^9$ system with half-filled 3$d_{x^2-y^2}$\xspace orbital is realized.
An {\it ab initio} estimate of Hubbard $U$ using the constrained random-phase approximation (cRPA)~\cite{Aryasetiawan_2004} shows that the correlation strength $U/t$ ($t$: nearest-neighbor hopping) is comparable to that of cuprates~\cite{Hirayama_2020}.
Therefore, once such nickelates are synthesized, the mother compounds will be a Mott insulator similarly to the cuprates, and the effective model becomes the Heisenberg model, which gets rid of the ambiguity in calculating $J$.

In this paper, we study the strength of $J$ in the two ideal $d^9$ nickelates, which are free from the self-doping (see Sec~\ref{sec_materials} for the details of the materials).
We estimate the $J$ value by the following three methods~\cite{Lichtenstein_2013}.
First, we start from a single-orbital Hubbard model derived in Ref.~\onlinecite{Hirayama_2020} and then evaluate $J$ by the expansion in terms of $t/U$. Second, we perform an energy mapping between the classical Heisenberg model and the total energy of different magnetic configurations calculated by the LDA+$U$ (LDA: local density approximation) method.
Third, we employ a scheme based on the so-called local force theorem.
Hereafter, we simply call these three methods ``strong-coupling expansion", ``energy mapping method", and ``local force approach", respectively.
We show that the three independent estimates show reasonable agreement and conclude that the $d^9$ nickelates have sizeable $J$ (about 100 meV),
which is not far smaller than that of the cuprates.
Therefore, the proposed $d^9$ nickelates provide an interesting playground to explore the cuprate-analog high-$T_{\rm c}$ superconductivity.

The paper is organized as follows.
In Sec.~\ref{sec_materials}, we introduce two ideal $d^9$ nickelates, RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$: a cation with the valence of $2.5+$) and discuss the advantage over NdNiO$_2$. In Sec.~\ref{sec_methods}, we explain the three methods employed in the present study, and we show the results in Sec.~\ref{sec_results}. Section~\ref{sec_summary} is devoted to the summary.

\begin{figure*}[tb]
\vspace{0.0cm}
\begin{center}
\includegraphics[width=0.98\textwidth]{fig_structure_band.pdf}
\caption{
Crystal structure of (a) RbCa$_2$NiO$_3$\xspace and (c) $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$: a cation with the valence of $2.5+$) and the paramagnetic DFT band structure [(b) RbCa$_2$NiO$_3$\xspace and (d) $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace)].
The blue dotted curves are the Wannier band dispersion of the Ni 3$d_{x^2-y^2}$\xspace single-orbital Hamiltonian.
In (b) and (d), the consistent ${\bf k}$ path is employed: $(0,0,0)$ $\rightarrow$ $(\pi/a,0,0)$ $\rightarrow$ $(\pi/a,\pi/a,0)$ $\rightarrow$  $(0,0,0)$ $\rightarrow$ $(0,0,\pi/c)$ $\rightarrow$ $(\pi/a,0,\pi/c)$ $\rightarrow$ $(\pi/a,\pi/a,\pi/c)$ $\rightarrow$  $(0,0,\pi/c)$ (The symbols are different because the primitive cells of RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace are tetragonal and bace-centered tetragonal, respectively).
}
\label{Fig_crystal_structure}
\end{center}
\end{figure*}

\section{Materials: $d^9$ nickelates}
\label{sec_materials}

In Ref.~\cite{Hirayama_2020}, various layered nickelates have been systematically proposed.
They are classified into ``1213", ``1214", ``H$_2$", and ``G" families, depending on the composition and the type of the block-layer~\cite{Tokura_1990}.
Among the four families, the compounds without the self-doping exist in the 1213 and G families.
We here take RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$: a cation with the valence of $2.5+$) for a representative of the ideal $d^9$ nickelates belonging to 1213 and G families, respectively (see Figs.~\ref{Fig_crystal_structure}(a) and (c) for the crystal structure).
In the following, we employ Ba$_{0.5}$La$_{0.5}$\xspace as $A$.
The phonon calculations have shown that both RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace) are dynamically stable~\cite{Hirayama_2020}.
We take the crystal structure optimized in Ref.~\cite{Hirayama_2020}, and perform density-functional theory (DFT) calculations~\footnote{Here, we ignore the interface effect~\cite{Geisler_2020,He_2020,Bernardini_2020c} and consider the bulk property.
Note that the thickness of the film reaches around 10 nm and there are several tens of NiO$_2$ layers in the sample~\cite{Lee_2020}.}.
Figs.~\ref{Fig_crystal_structure}(b) and \ref{Fig_crystal_structure}(d) show the paramagnetic DFT band structure for RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), respectively.
As is shown in Ref.~\cite{Hirayama_2020}, only the Ni 3$d_{x^2-y^2}$\xspace orbital crosses the Fermi level.
As far as the topology of the band structure is concerned, these systems are more similar to the cuprates than NdNiO$_2$.

The advantages of studying these nickelates rather than NdNiO$_2$ are as follows.
First, it is still controversial whether the role of Nd-layer (block-layer) orbitals is essential or not.
If the hybridization between Ni $3d$ and Nd-layer orbitals is substantial, the Nd-layer orbitals are not only a charge reservoir, but they might give Kondo-like physics~\cite{Sawatzky19,Hepting_2020,GM_Zhang_2020,Z_Wang_arXiv,Y_Gu_2020}.
In the cases of the $d^9$ nickelates, RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), the block-layer orbitals do not show up at the Fermi level, and this controversy can be avoided.
We can also exclude the possible role of the $4f$ orbitals with localized moments proposed in Refs.~\cite{P_Jiang_2019,Choi_2020}.

Another controversial issue for NdNiO$_2$ is to which orbitals the doped holes go ($d^9 \underline{L}$ vs. $d^8$, where $\underline{L}$ denotes a hole in a ligand oxygen).
In the case of the cuprates (charge-transfer insulator), the holes are doped into the oxygen $2p$ orbitals.
On the other hand, the nickelates have larger $\Delta_{dp}$ and are classified as Mott-Hubbard type~\cite{Hepting_2020,Jiang_2020,Nomura_2019,HuZhang_2020,Fu_arXiv,Goodge_arXiv}.
Because there is nonzero hybridization between Ni 3$d_{x^2-y^2}$\xspace and O $2p$ orbitals, some of the holes should be doped into oxygen $2p$ orbitals~\cite{Hirsch_2019,Karp_2020,Lang_arXiv}. However, the amount should be smaller than that of the cuprates.

When the system is Mott-Hubbard type and the holes mainly reside in the Ni $3d$ orbitals, another issue arises: which model is more appropriate, the single-orbital or multi-orbital model?
In other words, whether the doped $d^8$ configuration favors high-spin state or low-spin state.
If the crystal field splitting between Ni 3$d_{x^2-y^2}$\xspace and the other $3d$ orbitals is much larger than the Hund's coupling, holes stay within the Ni 3$d_{x^2-y^2}$\xspace orbital, and the single-orbital model is justified.
On this issue, several studies insist that Ni $3d$ multi-orbital nature cannot be
ignored~\cite{Jiang_2020,Zhang_2020,Werner_2019,Petocchi_arXiv,LH_Hu_2019,Lechermann_2020,Lechermann_arXiv,Y_Wang_arXiv,Chang_arXiv,Choi_arXiv}.
To resolve this issue, we certainly need more experimental evidences.
In the cases of RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), compared to NdNiO$_2$, the Ni 3$d_{x^2-y^2}$\xspace orbital is more isolated in energy space from the other $3d$ orbitals [see Figs.~\ref{Fig_crystal_structure}(b) and \ref{Fig_crystal_structure}(d)]:
In the case of NdNiO$_2$, due to the dispersion along the $k_z$ direction, the position of the Ni $d_{3z^2-r^2}$\xspace band becomes close to the Fermi level on the $k_z = \pi/c$ plane; however, such $k_z$ dependence is much weaker in RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace).
Considering also the above-mentioned absence of the complication from the self-doping, in this study, we adopt the single-orbital Hubbard model as a minimal model for RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace).
In the absence of the carrier doping, we can further map onto the spin model with the exchange coupling $J$.

\section{Methods}
\label{sec_methods}

Here, we introduce three different methods to estimate $J$
(see eg., Ref.~\cite{Lichtenstein_2013} for the ideas behind the three methods).
We employ the following convention for the spin Hamiltonian: ${\mathcal H} = \sum_{\langle i,j \rangle } J_{ij} {\bf S}_i \cdot {\bf S}_j$, where $\langle i,j \rangle$ is the bond consisting of sites $i$ and $j$, and ${\bf S}_i$ is the spin-1/2 operator at site $i$.
$J$ stands for the nearest-neighbor $J_{ij}$ interaction in the NiO$_2$ layer.

\subsection{Strong-coupling expansion}
\label{sec_method_superexchange}

When the single-orbital Hubbard model is a good description,
the magnetic interactions in the Mott insulating region can be obtained by strong-coupling perturbation expansion.
The strong-coupling expansion becomes valid in the region $U \gtrsim W$ with the bandwidth $W$ (in the square lattice $W =8t$)~\cite{Otsuki_2019}. RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace) with $U/t = $9.522 and 10.637, respectively~\cite{Hirayama_2020}, satisfy the condition $U > W$.

In the strong-coupling expansion, the superexchange interaction $J_{\rm s}$ (with $t^4$-order correction term) and cyclic ring-exchange interaction $J_{\rm c}$ are given by $J_{\rm s} = 4 t^2 / U - 24 t^4 / U^3$ and $J_{\rm c} = 80 t^4 / U^3$, respectively~\cite{Takahashi_1977,MacDonald_1988,Delannoy_2009}.
If we effectively take into account the effect of the ring-exchange interaction in the nearest-neighbor interaction $J$, the $J$ value becomes
\begin{eqnarray}
J = J_{\rm s} - 2J_{\rm c} S^2 = \frac{4 t^2} {U} -  \frac{64 t^4 }{U^3}
\end{eqnarray}
with $S=1/2$.

\subsection{Energy mapping method}
\label{sec_method_LDA+U}

Within the LDA+$U$~\cite{Anisimov_1991,Anisimov_1993,Liechtenstein_1995,Cococcioni_2012}, we perform the magnetic calculations. Here, $U$ is introduced into the Ni 3$d$ orbital subspace.
We employ $2\times2\times1$ supercell consisting of four conventional cells.
We simulate two different magnetic solutions: one is N\'eel type [$(\pi/a,\pi/a,0)$ antiferromagnetic order] and the other is stripe type [$(\pi/a,0,0)$ antiferromagnetic order].
We calculate the energy difference $\Delta E$ between the two antiferromagnetic solutions.
When we assume the two-dimensional classical spin-1/2 Heisenberg model up to next-nearest-neighbor magnetic interaction $J'$, $\Delta E$ per formula unit is given by $\Delta E = J/2 - J' \simeq J/2$.
We estimate $J$ with this equation.

\subsection{Local force approach}
\label{sec_method_Liechtenstein}

Based on the N\'eel-type solutions of the LDA+$U$ calculations, we estimate $J$ and $J'$ using the local force theorem~\cite{Lichtenstein_2013}.
The local force approach estimates the magnetic interactions from the small energy change induced by the infinitesimal spin-rotation from the magnetic solutions (N\'eel-type in the present case).
We employ the so-called Lichtenstein formula, which is recently developed in the low-energy Hamiltonian with the Wannier orbitals~\cite{Korotin_2015,Nomoto_2020_1,Nomoto_2020_2}, given by,
\begin{align}
    (-1)^PJ_{ij}= 4T\sum_{\omega_n}{\rm Tr}[G_{ij}(\omega_n) M_j G_{ji}(\omega_n) M_i],\label{eq:licht}
\end{align}
where $\omega_n=  (2n+1) \pi T$ denotes the Matsubara frequency. Here, we set $P=0$ (1) when the spins at $i$ and $j$ sites are aligned parallel (anti-parallel) to each other.
The Green's function $G_{ij}$ is defined by $G_{ij}^{-1}(i\omega_n)=(i\omega_n+\mu)\delta_{ij}-{\mathcal H}^0_{ij}$,
where ${\mathcal H}^0_{ij}$ is the hopping matrix of the Wannier tight-binding model, and $\mu$ is the chemical potential.
Note that ${\mathcal H}^0_{ij}$ is a $N_{{\rm orb}_i}\times N_{{\rm orb}_j}$ matrix, where $N_{{\rm orb}_i}$ is the number of Wannier orbitals at $i$-site including the spin index.
In the case of collinear magnets, one may write ${\mathcal H}^0_{ii}$ as ${\mathcal H}^0_{ii}=  \varepsilon_{i}\otimes\sigma_0 + m_i\otimes\sigma_z$.
Then, $M_i$ is defined by $M_i=m_i \otimes \sigma_x$ and proportional to the exchange splitting at $i$-site $m_i$.
Here, we have neglected the spin-dependent hopping term of $M_i$ (see Ref.~\onlinecite{Nomoto_2020_1} for details~\footnote{Note that the $J_{ij}$ value in this paper is defined to be eight times as large as that in Ref. ~\onlinecite{Nomoto_2020_1}}).

\subsection{Comparison among the three methods}
\label{sec_method_comparison}
The strong-coupling expansion gives local (in real-space) $J$, the energy mapping method sees the energy difference between the global and local minima of the magnetic solutions, and the local force method sees the low-energy excitations around the global minimum. These $J$ are complementary to each other, and hence we employ all the three methods.
When the Coulomb repulsion is much larger than the bandwidth and the mapping to the Heisenberg model becomes valid, these three methods see the same $J$.
As we will show in Sec.~\ref{sec_results}, the three results agree reasonably well as expected from a Mott insulating behavior of the proposed $d^9$ nickelates.

\subsection{Calculation conditions}
\label{sec_method_condition}

The DFT band structure calculations are performed using {\textsc{Quantum ESPRESSO}}\xspace~\cite{QE-2017}.
We employ Perdew-Burke-Ernzerhof (PBE)~\cite{Perdew_1996} norm-conserving pseudopotentials downloaded from PseudoDojo~\cite{Setten_2018} [the pseudopotentials are based on ONCVPSP (Optimized Norm-Conserving Vanderbilt PSeudopotential)~\cite{Hamann_2013}].

The energy comparison between the N\'eel- and stripe-type antiferromagnetic solutions is performed using $9\times 9 \times 7$ and $9\times 9 \times 3$ {\bf k}-mesh for RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), respectively.
We treat Ba$_{0.5}$La$_{0.5}$\xspace by the virtual crystal approximation.
The energy cutoff is set to be 100 Ry for the Kohn-Sham wave functions, and 400 Ry for the electron charge density.

For the estimate of $J$ based on the local force approach, we first construct the maximally localized Wannier functions~\cite{Marzari_1997,Souza_2001} for the  N\'eel-type antiferromagnetic band structure using RESPACK~\cite{Nakamura_arXiv,RESPACK_URL}.
For RbCa$_2$NiO$_3$\xspace, we use $5\times 5 \times 5$ {\bf k}-mesh for the construction of Wannier orbitals.
We put Ni $d$, O $p$, Ca $d$, and interstitial-$s$ (located at the interstitial positions surrounded by Ni$^{+}$, Ca$^{2+}$, and Rb$^{+}$ cations) projections.
The interstitial orbitals are stabilized because they feel attractions from the surrounded cations~\cite{Nomura_2019}.
Then, we obtain 104 orbital (per spin) tight-binding Hamiltonian.
For $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), we employ $5\times 5 \times 3$ {\bf k}-mesh for constructing Wannier orbitals.
We derive 232 orbital (per spin) tight-binding Hamiltonian using the projections of Ni $d$, O $p$, Br $p$, $A$ $d$, and interstitial-$s$ (located at the interstitial positions surrounded by Ni$^{+}$, $A^{2.5+}$, and Br$^{-}$ ions) orbitals.

In the calculation of Eq.~\eqref{eq:licht},
we employ $16\times16\times16$ $\bf k$-mesh and set the inverse temperature $\beta= 200$~eV$^{-1}$ for both cases.
We have confirmed that the difference of $J_{ij}$ values at $\beta= 200$ and 400~eV$^{-1}$ is less than 1 \%.
We use the intermediate representation basis for the Matsubara frequency summation~\cite{Shinaoka_2017,Chikano_2019,Li_Shinaoka_2020}, and set the cutoff parameter $\Lambda=10^5$, which is sufficiently larger than $W\beta$ where $W$ is the band width.

\begin{figure*}[tb]
\vspace{0.0cm}
\begin{center}
\includegraphics[width=0.99\textwidth]{fig_LDA+U.pdf}
\caption{
N\'eel-type antiferromagnetic band structure (red curves) and orbital-resolved density of states (per formula unit, per spin) for (a) RbCa$_2$NiO$_3$\xspace and (b) $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), calculated with LDA+$U$ method ($\overline{U} = 3$ eV).
The blue dotted curves are band dispersion calculated from the Wannier tight-binding Hamiltonian.
The symbols for the high-symmetry ${\bf k}$ points with the underlines are defined based on $2\times2\times1$ supercell consisting of four conventional cells.
The origin of the energy axis is set to be the middle of the gap.
The orbital-resolved density of states is calculated from the Wannier tight-binding Hamiltonian.
``I-$s$" stands for the interstitial-$s$ orbitals
(see Sec.~\ref{sec_method_condition} for the details of the projections used in the Wannier construction).
(c) The energy difference $\Delta E$ per formula unit between N\'eel- and stripe-type antiferromagnetic solutions.
The N\'eel-type solutions always show lower energy.
}
\label{Fig_LDA+U}
\end{center}
\end{figure*}

\section{\mbox{\boldmath$J$} in $d^9$ nickelates}
\label{sec_results}

\begin{figure*}[tb]
\vspace{0.0cm}
\begin{center}
\includegraphics[width=0.99\textwidth]{fig_J.pdf}
\caption{
Estimated exchange coupling $J$ for (a) RbCa$_2$NiO$_3$\xspace and (b) $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace).
$\overline{U}$ is the Hubbard interaction in the LDA+$U$ calculation (the Coulomb repulsion between the Ni 3$d$ orbitals), which we distinguish from the Hubbard $U$ in the single-orbital Hubbard model used in the strong-coupling expansion (the Coulomb repulsion between the Wannier orbitals made from the Ni 3$d_{x^2-y^2}$ orbital with O $2p$ tails). See text for details.
}
\label{Fig_J}
\end{center}
\end{figure*}

In the previous study~\cite{Hirayama_2020}, the effective single-orbital Hamiltonians for RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace) are constructed using maximally-localized Wannier functions~\cite{Marzari_1997,Souza_2001} and cRPA~\cite{Aryasetiawan_2004}.
The derived nearest-neighbor hopping and Hubbard parameters are
$t = -0.352$ eV, $U = 3.347$ eV for RbCa$_2$NiO$_3$\xspace,
and
$t = -0.337$ eV, $U = 3.586$ eV for $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace).
Then, the strong-coupling expansion described in Sec.~\ref{sec_method_superexchange} gives
$J = 122$ meV and $J=109$ meV for RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), respectively
(see Appendix~\ref{Appendix_dp} for the estimate from three-orbital $d$-$p$ model).

Figures~\ref{Fig_LDA+U}(a) and \ref{Fig_LDA+U}(b) show the band structure calculated by the LDA+$U$ method for the
N\'eel-type antiferromagnetic state.
While the Hubbard $U$ in the single-orbital Hubbard model is the Coulomb repulsion between the Wannier orbitals made from the Ni 3$d_{x^2-y^2}$ orbital with O $2p$ tails,
the $U$ interaction in the LDA+$U$ calculation is the Coulomb repulsion between the Ni 3$d$ orbitals.
To make the difference clearer, we call $U$ in the LDA+$U$ calculation $\overline{U}$.
In Figs.~\ref{Fig_LDA+U}(a) and \ref{Fig_LDA+U}(b), we have used $\overline{U}$ = 3 eV.

In contrast to the case of the LDA+$U$ calculation for NdNiO$_2$, where the system stays metal even in the presence of antiferromagnetic order~\cite{Botana_2020,HuZhang_2020,Z_Liu_2020}, both systems become insulating.
The top of the valence band has mainly Ni $3d$ character, in agreement with the classification into the Mott-Hubbard type insulator.
We see that both systems are insulating even at smaller $\overline{U}$ (= 1 eV).
For all the $\overline{U}$ region we studied (1-5 eV), there exists well defined spin-1/2 Ni spin moment.
The results suggest that, if these $d^9$ nickelates are synthesized, they become antiferromagnetic Mott insulator as in the cuprates.

Figure~\ref{Fig_LDA+U}(c) shows the energy difference $\Delta E$ per formula unit between the N\'eel- and stripe-type antiferromagnetic solutions.
$\Delta E$ decreases as $\overline{U}$ increases, which is a natural behavior given that $\Delta E$ is governed by $J$ and the origin of $J$ is the superexchange interaction.

In Figs.~\ref{Fig_LDA+U}(a) and \ref{Fig_LDA+U}(b), the band dispersions obtained by the Wannier tight-binding Hamiltonian, which are used in the local force approach, are also shown.
The Wannier bands well reproduce the LDA+$U$ magnetic band dispersions.

From $\Delta E$ in Fig.~\ref{Fig_LDA+U}(c), we perform the order estimate of $J$ by the energy mapping method with assuming $J'/J = 0.05$ (Sec.~\ref{sec_method_LDA+U}) \footnote{We do not pay special attention to the precise value of the ratio $J'/J$ because we are only interested in the order estimate of $J$.}.
Then $J$ is given by $J = \Delta E / 0.45 $.
We also estimate $J$ using the local force approach (Sec.~\ref{sec_method_Liechtenstein}).

These results on top of the $J$ value estimated by the strong-coupling expansion (see above) are summarized in Figs.~\ref{Fig_J}(a) and \ref{Fig_J}(b) for RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), respectively.
The $J$ value in the energy mapping method changes from about 140 meV ($\overline{U}=1$ eV) to 60 meV ($\overline{U}=5$ eV).
The local force approach gives $J \simeq 70$-80 meV.
These estimates give the same order of $J$ as the strong-coupling expansion results [$J = 122$ meV and $J=109$ meV for RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace), respectively].

Although the energy mapping method and local force approach are based on the same LDA+$U$ calculations,
we see that there is a discrepancy between the two results at small $\overline{U}$ values (although the difference is no more than two times).
It should be noted that the former method sees the global change of the energy between the completely different magnetic patterns,
whereas the latter approach only sees the local landscape around the N\'eel-type solutions,
as described in Sec.~\ref{sec_method_comparison}.
For larger $\overline{U}$, the agreement between these two results becomes better as is expected: The system can be mapped to the classical spin model with a constant $J$ regardless of the assumed magnetic structure in the local force approach.

Overall, all the three estimates of $J$ lie within 60-140 meV, and we conclude that the $d^9$ nickelates have a sizable $J$ of the order of 100 meV.
The agreement in the order estimate of $J$ among three independent methods shows that
RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace) are indeed Mott insulators with the effective model being the Heisenberg model, and the magnetic exchange coupling $J$ is governed by the superexchange interaction [if the materials were, for example, weakly correlated, the three methods would not agree well (see Sec.~\ref{sec_method_comparison})].

Finally, we compare the $J$ value with that of the cuprates.
In the cuprates, the magnitude of $J$ is intensively studied by Raman spectroscopy in the early stage~\cite{Lyons_1988_1,Lyons_1988_2,Sugai_1988}.
The $J$ value for La$_2$CuO$_4$ is estimated to be about 130 meV~\cite{Singh_1989}.
Systematic investigations have shown that the material dependence of $J$ in the cuprates family is weak~\cite{Sulewski_1990,Tokura_1990_raman}.
The numerical study on the $d$-$p$ model has also derived $J$ as large as about 130 meV~\cite{Hybertsen_1990}.
Compared to the $J$ value of 130 meV for the cuprates, our estimate based on the $d$-$p$ model giving 90-100 meV (see Appendix~\ref{Appendix_dp}) is small, which is consistent with the fact that $\Delta_{dp}$ is larger in the nickelates.
However, we note that the $J$ value of about 100 meV is still significantly large, and the $d^9$ nickelates would serve as interesting cuprate-analog materials.

\section{Summary}
\label{sec_summary}

One of the remarkable features of the high $T_{\rm c}$ cuprates is the large exchange coupling $J$, whose size is as large as 130 meV~\cite{Lee_Nagaosa_Wen_2006}. In the present study, we have evaluated the size of $J$ for $d^9$ nickelates from first principles.
While the cuprates having small $\Delta_{dp}$ belong to the charge-transfer type in the Zaanen-Sawatzky-Allen diagram~\cite{Zaanen_1985}, nickelates with larger $\Delta_{dp}$ belong to the Mott-Hubbard type.
To answer how large $J$ can be expected in the Mott-Hubbard insulating $d^9$ nickelates, we studied RbCa$_2$NiO$_3$\xspace and $A_{2}$NiO$_{2}$Br$_2$\xspace ($A$ = Ba$_{0.5}$La$_{0.5}$\xspace),
which were recently proposed theoretically and shown to be free from the self-doping in Ref.~\onlinecite{Hirayama_2020}.
By means of the strong-coupling expansion, energy mapping method, and local force approach, we found that $J$ for these nickelates is as large as 100 meV, which is not far smaller than that of the cuprates.
This result suggests that the $d^9$ nickelates and cuprates share a notable common feature in the Mott insulating phase, although the former and latter belong to the Mott-Hubbard and charge-transfer regime, respectively.

Finally, we note that the proposed $d^9$ nickelates might give rare examples of realizing the square-lattice Hubbard model with sizeable $J$ in real materials.
Recent numerical studies show that the phase diagram of the doped Hubbard model is under severe competition between the stripe state with charge/spin modulation and $d$-wave superconductivity~\cite{Zheng_2017,Darmawan_2018,Ohgoe_2020,HC_Jiang_2019}.
Therefore, once synthesized, the $d^9$ nickelates will serve as a valuable test-bed system to understand the superconductivity in the Hubbard-like model.
They are also an important reference to understand the superconducting mechanism in the cuprates, because they would tell us whether the charge-transfer nature in the cuprates is essential in the high-$T_{\rm c}$ superconductivity or not.

\begin{acknowledgments}
We acknowledge the financial support by JSPS KAKENHI
Grant No. 16H06345 (YN, MH, and RA), 17K14336 (YN), 18H01158 (YN), 19K14654 (TN), 19H05825 (RA), 20K14390 (MH), and 20K14423 (YN).
This work was supported by MEXT as ``Program for Promoting Researches on the Supercomputer Fugaku" (Basic Science for Emergence and Functionality in Quantum Matter).
A part of the calculations was performed at Supercomputer Center, Institute for Solid State Physics, University of Tokyo.
\end{acknowledgments}