Two-wire Junction of Inequivalent Tomonaga-Luttinger Liquids

We develop two novel numerical schemes to study the conductance of the two-wire junction of inequivalent Tomonaga-Luttinger Liquids. In the first scheme we use the static current-current correlation function across the junction to extract the linear conductance through a relation that is derived via the bosonization method. In the second scheme we apply a bias and evaluate the time-dependent current across the junction to obtain the current-voltage characteristic. The conductance is then extracted from the small bias result within the linear response regime. Both schemes are based on the infinite size matrix product state to minimize the finite-size effects. Due to the lack of the translational invariance, we focus on a finite-size window containing the junction. For time-independent calculations, we use infinite boundary conditions to evaluate the correlations within the window. For time-dependent calculations, we use the window technique to evaluate the local currents within the window. The numerical results obtained by both schemes show excellent agreement with the analytical predictions.


I. INTRODUCTION
Transport properties of the strongly correlated quasione-dimensional (1D) quantum systems have been the subject of intensive investigation in recent years due to the potential applications in nanoelectronics. In these systems, electron-electron interaction has drastic effects and the Fermi liquid theory breaks down. Instead, the system is described by the Tomonaga-Luttinger liquid (TLL) theory, which is parameterized by a Luttinger parameter g [1][2][3]. Experimentally, TLL's characteristic behavior has been observed experimentally in a variety of quasi-1D systems [4][5][6][7][8]. In this work we focus on an important class of the quasi-1D transport problem: junctions of multiple TLL wires. The simplest setup is a twowire junction of equivalent TLL. It consists of two TLL wires with the same Luttinger parameter connected by a weak link. Theoretically, it is well known that in this case, the system renormalizes either to the single fully connected wire fixed point or to the two disconnected wires fixed point, depending on the sign of the interaction [9][10][11][12]. For three-wire Y junctions of equivalent TLL, more conductance fixed points begin to emerge [13][14][15]. From the perspective of the experiment, it is also important to study the influence of the contact as well as the Fermi liquid leads [16] and the multiwire junction with inequivalent TLLs [17][18][19][20].
In the above-mentioned studies, the bosonization method has been used extensively to draw important conclusions [9-11, 14, 15, 17, 18, 21-23]. However, in order to go beyond the perturbative regime to reach a * yjkao@phys.ntu.edu.tw † pcchen@phys.nthu.edu.tw comprehensive understanding of quasi-1D systems' transport properties, other methods are called for. Analytically, an exact solution method based on Bethe ansatz has been developed and applied successfully to several systems [12,24,25]. However, it is restricted to integrable models. On the other hand, many numerical methods have been developed. Within fermion representation, methods based on renormalization group equations [26][27][28] and functional renormalization group (fRG) technique [29,30] have been developed to evaluate the one-particle Green's function from which the linear conductance can be extracted. Numerical renormalization group (NRG) method [31], which is originally developed for the equilibrium properties of quantum impurity systems, has also been generalized to study the transport properties of nanodevices with noninteracting leads [32][33][34]. However, it is difficult to generalize the NRG based method to study a junction with interacting leads.
The method developed here is based on the density matrix renormalization group (DMRG) technique, which is a powerful and versatile numerical tool to study quasi-1D systems [35]. Over the years, several DMRG based methods have been developed to study the transport properties of quasi-1D systems. One of the earliest approaches uses a ring geometry and extracts the conductance from the current induced by the flux [36][37][38]. Another approach uses the linear response theory to relate the conductance and the dynamical correlation functions, which can be evaluated via the correction vector DMRG method [39] or the dynamical DMRG method [39][40][41]. In Refs. [18,42,43] a general method to extract the conductance tensor of the multiwire junction is proposed and is used to study the multiwire junction of equivalent and inequivalent TLLs. By using boundary conformal field theory, the conductance tensor is related to the static cor-arXiv:2105.03104v2 [cond-mat.str-el] 24 Dec 2021 relation function of a semi-infinite system. A conformal transformation is then used to connect the correlators of a finite system to a semi-infinite one, making it possible to extract the conductance tensor from the static currentcurrent correlator of a finite system. We note in passing that in this approach it is necessary to add a mirrorimage junction during the transformation. However, the Hamiltonian of the mirror-image junction can be rigorously derived only for the case of noninteracting wires. On the other hand, since the invention of time-dependent DMRG (tDMRG), various approaches have been developed to simulate the time-dependent current of the multiwire junction, from which the linear response and the full current-voltage characteristics can be obtained [44][45][46][47][48][49][50][51][52][53]. Recently, a relation between the static charge correlations and the linear has also been put forward [54]. Typically, a finite-size system with open boundary condition is simulated. However, this may lead to strong finite-size effects. In order to reduce the finite-size effects, modified boundary conditions have also been explored [45,47,50,55].
In this work, we develop two novel numerical schemes to study the linear and the nonlinear conductance of the multiwire junctions. Our main strategy is to always work with an infinite system to minimize the finite-size effects but only perform measurements within a finite size window to make the simulation feasible. This also removes the need to perform the conformal transformation to obtain an effective finite size system and the addition of the mirror-image junction. Specifically, we revisit the problem of the two-wire junction of inequivalent TLLs to benchmark our methods. It is known that in this case the conductance is determined by a single effective Luttinger parameter [17][18][19][20]. In the first scheme, we use the method recently proposed by some of us in Ref. [56] to calculate the static current-current correlation function of the two-wire junction and use the method proposed in Refs. [18,42,43] to extract the linear conductance. However, we find that the key relation between the conductance and the static current-current correlation function derived in Refs. [18,42,43] needs minor modification when the Luttinger parameters on two wires are different. While we only measure the correlation function within a finite size window, we show that a moderate window size already allows us to observe the asymptotic behavior. Our results agree excellently with the theoretical prediction and verify that the conductance is indeed governed by the effective Luttinger parameter.
To probe the nonlinear conductance, we incorporate the window technique developed in Ref. [57] to evaluate directly the local currents within a finite size window, after a source-drain bias is applied to the system. We show that after the transient time, we can obtain a very flat quasistationary current up to a time scale that is limited by the window size and the carrier velocity. This allows us to define an average current with very small error, from which the current-voltage characteristics can be obtained. Our results show a wide range of linear response regime, and the linear conductance can be extracted from a small bias calculation. Furthermore, we verify that the linear conductance obtained from the nonequilibrium setup is highly consistent with the results via static correlations.
The rest of the paper is organized as follows. In Sec. II, we set up the notation and define the Hamiltonian of the model. In Sec. III, we use the bosonization method to derive the modified key relationship between the conductance and the static current-current correlation function. In Sec. IV, we discuss the main ingredients of our numerical method and outline the steps. In Sec.V, we present our numerical results from the time-independent calculation, while in Sec.VI we present our numerical results from the time-dependent calculations. Finally we summarize in Sec.VII and discuss future directions.

II. MODEL
We consider a two-wire junction which consists of two semi-infinite long TLL wires connected by a weak link as sketched in Fig. 1. To model such a junction, we start from two semi-infinite long wires with the Hamiltonian where c µ † i (c µ i ) with µ ∈ α, β is the creation (annihilation) operator at the site i of the wire µ andñ µ i ≡ c µ † i c µ i − 1 2 . Z ≥ denotes the set of non-negative integers. We note that the hopping strength is set to unity as the energy scale. Furthermore, H µ wire can also be expressed in the form of a translational invariant matrix product operator (MPO) [58] where When |U µ | < 2, the wire is in the TLL phase. At half filling, its Luttinger parameter g µ is determined by U µ through the relation while the carrier velocity reads We note that U µ = 0 corresponds to a noninteracting wire with g = 1, while positive and negative interaction correspond to g > 1 and g < 1, respectively. We form a two-wires junction by connecting a semiinfinite long wire with label α extending to the left and a semi-infinite long wire with label β extending to the right at site i = ±1/2 by a link of strength t . The Hamiltonian of the link reads: while the Hamiltonian of the junction reads where Z ≥ denotes the set of non-negative integers. In the form of the MPO the H junc is expressed as where W − 1 2 is associated with the H link as follows: In the following we denote the ground state of H junc by |Ψ junc . We define the current operator J(r) as where we define the current operator across the link as J link . Furthermore, the static current-current correlation function is defined as It is also convenient to define a r-dependent conductance It will be shown in Sec. III that G αβ (r) approaches the linear conductance G αβ as r goes to infinity, i.e., We note that Eq. (12) is different from the result in Refs. [18,43] except when v α = v β . In Sec. V we shall report numerical results that support Eqs. (12) and (13).
To probe the nonlinear response and determine the current-voltage characteristics, we turn on a bias at time t = 0 and apply a voltage of ±V /2 to the left and right wires, respectively. We then evaluate the time-dependent local current where |Ψ junc (V, t) = e i(Hjunc+H bias (V ))t |Ψ junc (15) and is the bias Hamiltonian. Formally we define the steady state current across the link as Numerically it is difficult to reach the infinite time limit.
In practice we expect that after a transient time the local current across the link will become quasistationary, from which we can define a time-averaged local current within a time window [t 1 , t 2 ]. Here we show explicitly the V dependence for clarity. If one can simulate accurately the local current up to a large enough time to contain a window such that that t 1 is larger than the transient time scale and t 2 −t 1 is large enough to obtain a good average, then the current-voltage characteristics can be accurately extracted. Furthermore, the linear conductance G αβ can be obtained by from a time-dependent calculation. We note in passing that this provides a consistency check by comparing with results from the time-independent calculations.

III. BOSONIZATION
In this section we first review known results in the literature, then we derive Eqs. (12) and (13) via the bosonization method. It is shown in Ref. [17] that the conductance G αβ of a junction of inequivalent TLLs depends only on an effective Luttinger parameter g e , where 1 g e = 1 2 Combined with the result for a junction of equivalent TLLs [9,10] one finds: When g e > 1 the conductance takes a universal value regardless the link strength of the junction. In contrast, when g e < 1 one has It is worth commenting on the case of g e = 1. In general, the effective theory is given by the boundary sine-Gordon theory with the marginal boundary interaction. For the free fermion (g µ = 1), the coefficient of the boundary interaction can be determined exactly and it is related to the exact conductance. This results in [9,10] In the presence of the interaction (g µ = 1), the effective theory is still the same if g e is still 1. While there is no reason to expect a strong renormalization of the coefficient of the marginal boundary interaction, the conductance should stay approximately the same as in the free fermion limit. The conductance is, however, eventually determined by the nonuniversal coefficient of the operator in the field theory. We hence suspect that the coefficient is renormalized weakly by the interaction, resulting in a weak change of the conductance. It will been shown later that our numerical results do support such a scenario. Unfortunately, there is no simple way to determine analytically the weak change of the coefficient and the conductance using the field theory.
In the literature, the relation between the conductance and the static current-current correlation function has been derived within the framework of boundary conformal field theory for the case of (i) multiple wires with the same Luttinger parameter and carrier velocity [42], (ii) multiple wires with the same Luttinger parameter but different carrier velocities [43], and (iii) multiple wires with different Luttinger parameters and different carrier velocities [18]. We note that the form of the velocity dependence in Ref. [43] and Ref. [18] is the same, and it falls back to the results in Ref. [42] if all the velocities are the same. However, as will be shown below, we find that this velocity dependence needs to be modified. On the other hand, in Ref. [18] it is shown that by rescaling the bosonic fields, the junction of two wires with different Luttinger parameters can be mapped to a junction with a single effective Luttinger parameter g e , in agreement with Ref. [17]. We confirm that this stands valid even when the two wires have different charge velocities, because the charge velocities can be absorbed by a proper rescaling.
We first follow the derivation in Refs. [18,43] to obtain the result for the case of v α = v β = 1. We then pay special attention to the case of inequivalent TLLs with different Fermi velocities and identify the proper way to rescale the equation. We start from the Kubo formula for the conductance Eq. (38) where the electric field is applied uniformly on the finite segment l 1 < −y, x < l 2 of the infinite system, and l = l 2 − l 1 is the length of the segment. Here T indicates imaginary-time ordering. When the two wires are inequivalent with two different Luttinger parameters where z ≡ τ + ix,z ≡ τ − ix, and the coefficients A αβ B are universal and determined by the conformally invariant boundary condition on the real axis. By using the identity where H(u) is the Heaviside step function, and Eq. (25) (with z 1 = τ − iy and z 2 = ix), Eq. (24) is reduced to This relates the conductance with the universal coefficient A αβ B for each conformally invariant boundary condition.
For a nonchiral (time-reversal invariant) junction, and thus the conductance is related to the asymptotic behavior of the current-current correlation function as for sufficiently large r. In fact, in the present problem of the junction of two wires, for generic values of the Luttinger parameters, we only need to consider the Dirichlet boundary condition and the Neumann boundary condition They are stable when g e > 1 and g e < 1, resulting in the conductance as in Eqs. (21) and (22), respectively. Because the asymptotic current-current correlation function is dominated by the subleading corrections in the latter case, in this paper we are mostly interested in the Dirichlet boundary condition (31) which represent maximally conducting junction realized for g e > 1. So far we have just followed Refs. [18,43]. Now let us resurrect the carrier velocity of each wire. When the velocity is v, the holomorphic/antiholomorphic variables z,z should be proportional to vτ ± ix. In other words, we can either For a single, uniform wire, these two formulations are equivalent under a rescaling by the factor of v. However, when multiple inequivalent wires are coupled, they are not equivalent. In our problem, the wires are coupled at a single junction, and both wires share the same time. Namely, the electron hops from one of the wires at time t should appear on the other side of the junction at the same time t. If one scales the time differently for the two wires, the electron transfer at the junction becomes nonlocal in the rescaled time variable. We should better avoid such a complication and instead rescale the spatial coordinate of each wire to define holomorphic/antiholomorphic variables: This means that the same holomorphic/antiholomorphic coordinate corresponds to different distances from the origin (junction). However, this does not pose a problem, as the two wires are connected only at the junction. With this rescaling, the current J remains unchanged, while the charge density ρ in each wire is multiplied by the velocity v µ . This is because the current measures the charge passing through a specific point per unit time, while the charge density measures the charge per unit length. Only the latter is renormalized by the rescaling of the length.
Thus, including the velocity, the current-current correlation function Eq. (25) appearing in the Kubo formula Eq. (24) reads (with τ 1 = τ and τ 2 = 0) Upon integration, one has Thus the conductance remains the same universal value independent of the velocities v α , v β . On the other hand, the real-space correlation function depends on the velocity factor. The equal-time correlation function of currents is given by setting τ = 0 in Eq. (35) as Setting x = r, y = −r, we find Therefore, we find This is to be compared with Eqs. (72) and (76) of Ref. [43] or Eq. (4.4) of Ref. [18] under time-reversal invariance and in the thermodynamic limit: That is, our result is different from the result in Ref. [43] by the factor of Two results are identical if and only if v α = v β . In the next section, we will demonstrate that the present results (38) and (40) are indeed consistent with numerical simulations.
We shall also emphasize that derivation above only consider the universal and dominant part of the correlation function. In general there are also nonuniversal parts which decay faster. That is why in the numerical simulation, the conductance is extracted from the asymptotic behavior of G αβ (r) as defined in Eq. (13).

IV. NUMERICAL METHOD
We note that the linear conductance and the stationary current, which are defined by Eq. (13) and Eq. (17), respectively, are well defined for an infinite system and at infinite time limit. However, due to the lack of translational invariance, in general it is difficult to simulate a two-wire junction in the thermodynamic limit. In order to extract the conductance from a finite size calculation, various approaches have been proposed. One approach is to map an infinite multiwire junction to a finite strip by a conformal mapping [42,43]. The mapping enables one to relate the linear conductance to the static current-current correlation function in a finite system. However, in this approach it is necessary to apply an ad hoc mirror boundary condition, which can be rigorously proved only for the case of non-interacting wires. Time-dependent DMRG has also been used to simulate the time-dependent local current. Typically a finite-size system with open boundary condition is studied. However, it has been show that finite-size effects can be severe [45]. On one hand, the current will completely reverse its direction after a finite period of time. On the other hand, the quasistationary current may have oscillation due to the presence of the boundary. This makes it difficult to obtain a well defined averaged current. It is known that these finite-size effects can be reduced by using damped boundary condition, which has been applied to the quantum dot systems connected to metallic leads [49] and fractional quantum Hall systems [50].
In this work we develop a framework to evaluate the static correlation functions and time-dependent currents for an infinite two-wire system. Before we outline the major steps, we briefly describe the key points of our approach. The main strategy of our approach is to work with an infinite system to minimize the finite-size effects but only perform measurements within a finite size window to make the simulation feasible. To proceed, we always assume that the wave function is in the form of an infinite matrix product state (iMPS). For translationally invariant systems the corresponding iMPS can be represented by only a few matrices. In contrast, in principle infinite many different matrices are needed for the iMPS representation of the ground state of H junc , since the translation invariance is broken. For impurity or interface systems, some of us propose an iMPS ansatz that only requires a finite number of matrices in Ref. [56]. The main assumption is that far away from the impurity or interface, the wave function should be almost the same as the bulk system without the impurity or interface. We hence consider a finite size window. Outside the window, the wave function is assumed to be the same as the corresponding bulk one. Inside the window, we need to optimize the wave function with an effective Hamiltonian, which is obtained by attaching infinite boundary conditions (IBCs) at the left and right boundary of the window. It is shown in Ref. [56] that the static correlation functions within the window can be evaluated accurately. Furthermore, we ensure that the window size is large enough such that the current-current correlation has reached its asymptotic behavior.
In this work, we further extend the method to the timedependent problem by applying the technique proposed in Ref. [57]. By using this technique, we can evaluate accurately the time-dependent local currents within a finitesize window. We show that the amplitude of the residual oscillation in the quasistationary region is extremely small, and an excellent time-averaged current can be obtained. In addition, while a small amount of current still reflects at the window boundary, the current never totally reverse its direction. Finally, it will be shown later that a longer quasistationary region can be reached by simply enlarging the window size. Now we are in a position to outline the major steps of our framework. We start from the iMPS ansatz with window size 2L for the two-wire junction: where A α i (t) = A α (t) and B β i (t) = B β (t) are siteindependent D × D matrices, M α i (t) are site-dependent D×D matrices, | n = | . . . n − 3 2 n − 1 2 n 1 2 n 3 2 · · · is the product basis, and D is the maximum virtual bond dimension. We assume a one site unit cell but it is straightforward to allow a multisite unit cell. When t ≤ 0, we assume |Ψ 2L junc (t ≤ 0) is the variational ground state of H junc . When t ≥ 0, we time evolve |Ψ 2L junc (0) with H junc + H bias (V ).
To find A α (t ≤ 0), B β (t ≤ 0), and M i (t ≤ 0) we first use infinite-size DMRG algorithms [59] to obtain the ground state |ψ µ wire of wire µ with Hamiltonian H µ wire as an iMPS in the mixed canonical form [35]: where λ is a D × D diagonal matrix while A µ and B µ are D × D matrices. Furthermore,Ã µ satisfy the left canonical form constraint µÃ µ †Ãµ = I, andB µ satisfy the right canonical form constraint µB µBµ † = I [35].
We then assume A α (t ≤ 0) =Ã α and B β (t ≤ 0) =B β .  is defined in Eq. (9), and W L ( W R ) represents the left (right) infinite boundary condition (IBC). We note that the left IBC W L is constructed from A α and W α , while the right IBC W R is constructed from B β and W β . The detail of the procedure to construct IBC can be found in Ref. [56]. We use the finite-size DMRG algorithm to obtain the ground state of H 2L junc as a finite MPS:

We next construct an effective 2L-sites Hamiltonian as
We then assume M i (t < 0) = M i . At this stage we have obtained |Ψ 2L junc (t ≤ 0) , from which we can evaluate the ground-state current-current correlation function J(−r)J(r) junc . We note that due to the left and right canonical conditions of the A and B matrices, the correlation functions within the window can be evaluated using only M matrices. We next calculate the position dependent conductance G αβ (r) and use the asymptotic behavior of G αβ (r) to estimate the conductance according to Eq. (13). This concludes the time-independent part of the calculation.
To find A α (t ≥ 0), B β (t ≥ 0), and M i (t ≥ 0), we first break the time-evolution operator into products of time-evolution operator with a small time step dt 1. At each time step, we use second-order Suzuki-Trotter decomposition to approximate the evolution operator as products of local gates. Due to the lack of translational invariance we perform the time evolution in three substeps: • Perform infinite-size TEBD update for A α (t) and B β (t) with bulk Hamiltonian H α wire and H β wire , respectively. We note that they are site independent (up to a unit cell) and standard infinite size TEBD update is used.
• Perform finite-size TEBD update for M i (t) with the effective Hamiltonian H 2L junc , except M ∓n L± 1 2 (t) at the boundary of the window.
At this stage we have obtained A α (t+dt), B β (t+dt), and M i (t + dt). We then evaluate J(r, t + dt) and iterate the procedure.
Some comments are in order. While the wave function is always represented by an iMPS, only finite numbers of matrices need to be updated at each time step. In principle the window can co-move with the wave front. Here we fix the window size for simplicity. There are two main factors that limit the time scale that one can simulate accurately the time-dependent current. First, TEBD update will eventually break down. Typically, by increasing the maximum virtual bond dimension, a larger time scale can be reached. Second, due to the finite virtual bond dimension, the reflection at the boundary cannot be complete removed. One can reduce the amount of reflection by increasing the maximum virtual bond dimension of the IBC or one can increase the window size to delay the time of reflection. In this work we always check how our results depend on maximum virtual bond dimension and window size to ensure the convergence of the results.

V. TIME-INDEPENDENT RESULTS
In this section we present our results of the static current-current correlation function and the linear conductance. We first focus on the case of a two-wire junction with an effective Luttinger parameter g e > 1. In this case it is expected that the system renormalizes to a single fully connected wire with a linear conductance G th = g e e 2 /h, regardless of the Luttinger parameter of each individual wire and the link strength t . Specifically we consider three combinations of g α and g β that give rise to g e = 1.2: (i) g α = 0.8 and g β = 2.4, (ii) g α = 1.0 and g β = 1.5, and (iii) g α = 1.10, and g β = 1.32.
In the left column of Fig. 2 we plot the static currentcurrent correlation function J α (y)J β (x) as a function of −y/v α + x/v β on a log-log scale, where the variable x and y are restricted to the corresponding lines in the 10 0 10 1 10 2 right column. This choice of the scaling variable is motivated by the theoretical analysis in Sec. III. Here we set link strength t = 0.7, maximum virtual bond dimension D = 800, and window size 2L = 400. We observe that the correlation function quickly converges to Eq. (38) (dotted black line) with the condition (31), regardless of the parameters used, in excellent agreement with the theoretical prediction. Furthermore, we also perform simulation with t = 0.9 and we find that the results are indistinguishable from the results shown above. We note in passing that this also verifies that the window size and maximum virtual bond dimension are large enough to capture the asymptotic behavior. where v µ is the carrier velocity.
Let us make a more quantitative comparison between the theory and the numerical data, in terms of the conductance estimated using Eq. (40). In Fig. 3(a) we show the distance dependent conductance G(r) (in units of e 2 /h) as a function of r. We observe that for all combinations of g α,β and t , G(r) quickly approaches the theoretical prediction (dashed black line). To investigate how G(r) approaches its asymptotic value, in Fig. 3(b) we show G th − G(r) vs r on a log-log scale. We find that asymptotically it decays nearly as a power law. We note that the bump at large distance is due to the finite window effect, which can be eliminated by enlarging the window size.
We next study the case of a two-wire junction with an effective Luttinger parameter g e < 1. In this case it is expected that the system renormalizes to two disconnected wires fixed point with G th = 0. Specifically, we use two combinations of g α and g β : (i) g α = 0.6, g β = 1.2 and (ii) g α = 0.72, g β = 0.9 that give rise to g e = 0.8. To probe the disconnected wires behavior we use two values of link strength t = 0.3 and t = 0.5. Smaller t is used here to ensure that the correlation function can reach its asymptotic behavior within the window. In Fig. 3(c) we show t −2 G(r) as a function of r on a log-log scale. We observe a nonuniversal behavior and the value of G(r) depends on g α,β and t . However, asymptotically G(r) always decays as a power law. While our data is limited by the window size, we expect that lim r→∞ G(r) = 0, consistent with the broken wire interpretation.
Finally we investigate the case of g e = 1. Here we use three combinations of g α,β that give rise to g e = 1: (a) g α = 0.75, g β = 1.5, (b) g α = 0.9, g β = 1.125, and (c) g α = g β = 1. In Fig. 4 we plot G αβ (r)/G th as a function of r. The ratio approaches one regardless of the parameters used, but we also find a tiny deviation when the wires are interacting. (d) .002 vs t for the case of g α = g β = 1.2. The vertical dashed and dash-dot lines corresponds to t1 and t2, respectively, where [t1, t2] is the time interval over which the time average of the current is taken.

VI. TIME-DEPENDENT RESULTS
In this section we present our time-dependent results. We first benchmark our method with a junction of two equivalent wires with g µ = 1.2. The hopping strength t between wires is set to 0.8. We use second order Suzuki-Trotter expansion with dt = 0.002 to perform the time evolution. The window size 2L is set to 100 and the maximum virtual bond dimension D is 200. In Fig. 5(a), we plot J(r, t) V (in units of e 2 /h) on the r − t plane after a small bias V = 0.002 is applied to the system at t = 0. It is expected that the current will first appear in a location where the voltage is reversed. Indeed we observe that the current emerges at the junction ( r = 0 ) once we turn on the bias. After that, the fronts propagate in both directions and a light cone is formed. In Fig. 5(a) we also plot r = ±v µ t as dashed lines, where v µ is the carrier velocity. It is evident that the slope of the light cone agrees well with the carrier velocity. While the window technique can minimize the reflection at the window boundary, partial reflection still occurs due to the finite virtual bond dimension. We find that the front hits the window around t ≈ L/v µ ≈ 30. However, the partial reflection does not appear until around t ≈ 37.
We next turn our attention to the time-dependent current across the link. In Fig. 6(a) we plot J link (t) V /V (in units of e 2 /h) as a function of time. We observe a fast increase from zero once we turn on the bias, followed by the transiently decaying oscillations. This oscillation is bias dependent and due to the backward scattering and Andreev-type reflection at the junction [60][61][62][63]. After that, the current becomes quasistationary until t ≈ 70. At this time the fronts of the partially reflected current reach r = 0 and unphysical oscillations start to emerge. However, the current across the link never reverses its direction in our simulation. This is in contrast to the simulations with open boundary condition, in which the current will change direction periodically. It is worth mentioning that the amplitude of residual oscillations in the quasistationary region is extremely small. This is because we work with an infinite system and the finite-size induced oscillation is eliminated. To obtain the averaged current we identify a time interval (t 1 , t 2 ) as follows: We first set t 2 to be a time that is slightly before the reflected front reaches the center. We then move t 1 away from t 2 until the estimated error is minimized. In Fig. 6(a) we draw t 1 and t 2 as vertical dashed and dashed-dotted lines, respectively. From the time-averaged current we findJ link /V e 2 h = 1.200 ± 0.004, which agrees excellently with the expected results of G th = e 2 h g e (= 1.2 e 2 h ). In general, there are two simulation parameters which determine the accuracy of the time-averaged current J link . The window size L determines the maximal time scale before the unphysical oscillations appears. The maximum virtual bond dimension D determines the amount of reflection as well as the quality of time evolution. To investigate how our numerical results depend on L and D, we first fix D and run the simulations with various L as shown in Fig. 6(b). We observe that results from different L almost collapse into a single curve before the unphysical oscillation occurs. Next we fix L = 50 and run the simulations with various D. As shown in Fig. 6(c), we find that both the amplitude of the residual oscillations in the quasistationary region and the amplitude of the unphysical oscillation decrease as D increases. This is consistent with the picture that quality of timeevolution increases as D increases. In Table I, we list our numerical results ofJ link /V e 2 h and the estimated error, which are obtained with various D. We find that results obtained with D = 40 already agree well with the theoretical prediction and the estimated error is already small. Furthermore, the error continues to decrease as D increases. It is worth mentioning that this is one of the advantages of the method proposed in this work. For a similar simulation in Ref. [47], the tDMRG results are not convergent until D ≥ 300.
We note that while the finite-size induced oscillations are eliminated in the quasistationary regime, J link (t) V /V e 2 h still has residual oscillations. We find that, however, the amplitude of the residual oscillation is G αβ (t = 0.9) g α = 0.8, g β = 2. 4 1.1955 ± 0.0008 1.1975 ± 0.0013 g α = 1.0, g β = 1. 5 1.1845 ± 0.0020 1.2028 ± 0.0034 g α = 1.1, g β = 1. 32 1.1783 ± 0.0017 1.1998 ± 0.0076 almost independent of the bias. Consequently, for larger bias the residual oscillations for J link (t) V /V e 2 h are invisible. In contrast for a very small bias the residual oscillation becomes visible, leading to a larger error in J link /V e 2 h . In order to ensure that the system is in the linear response region with minimized error, in the following we will use V = 0.002 when evaluating the linear conductance. Now we are in a position to study the case of two inequivalent wires with an effective g e > 1. Three combinations of g α,β with g e = 1.2 are considered: (a) g α = 0.8, g β = 2.4, (b) g α = 1.0, g β = 1.5, and (c) g α = 1.1, g β = 1.32, where we use the convention of g α ≤ g β . Two hopping strengths t = 0.7 and t = 0.9 are used. In Fig. 5(b) we show the time-dependent local current on the r − t plane. Due to the different propagation velocities, the light cone is asymmetric. Since a smaller Luttinger parameter implies larger carrier velocity, the current reflected by the left boundary will come back first. In Fig. 7, we show J link (t) V /V e 2 h as a function of time. We observe that in all cases it converges to the expected result of G th = e 2 h g e (= 1.2 e 2 h ), which is denoted by the dotted horizontal line. By averaging the current between the dashed and the dash dot vertical lines, we obtain the conductance as summarized in Table II. The results are highly consistent with the theoretical prediction as well as the results from the current-current correlation function as presented in Sec. V.
To probe the nonlinear response, we set D = 200, L = 50 and run the simulations with various V . We use g α = g β = 1.2 and g α = 0.8, g β = 2.4, both with g e = 1.2. In Fig. 6 h as a function of time. A long period of quasistationary region is observed for all cases except when V is close to 1. This is due to the finite-window effect and can be eliminated by enlarging the window. In Fig. 8 we show the time-averaged current J link in units of e 2 /h as a function of V on a log-log scale. Two hopping strengths t = 0.7 and t = 0.9 are considered. We also plot a straight line that corresponds to the linear responseJ = G th V = e 2 h g e V . We observe a universal linear response up to V ≈ 0.2 while at large bias the deviation becomes substantial. As has been pointed out in Ref. [64], the deviation at large bias is due to the finite band width of the model.
We now investigate the case of g e < 1. We consider two combinations (a) g α = 0.6, g β = 1.2 and (b) g α = 0.72, g β = 0.9, which give rise to g e = 0.8. It is expected that in this case the the junction will flow to the broken wire fixed point with zero conductance. In order to observe the broken wire behavior within the time window, we use smaller link strength, t = 0.3 and 0.5. From the perturbation theory we expect that the current might scale as (t ) −2 . In Fig. 9 we plot (t ) −2 J link (t) V /V e 2 h on a loglog plot. At short time, we observe a universal rise of the rescaled current. After that, it decays as a power law but the decay exponent is parameter dependent. These results indicate that at large time one has G αβ = 0, which is in agreement with the theoretical predictions and the results of the time-independent method.
Finally we investigate the case of g e = 1. In Fig. 10 10 −1 10 0 10 1 t we show our results of the conductance obtained by the time-dependent calculation as a function of link strength. For junction of equivalent noninteracting wires (g µ = 1), the results agree excellently with the theoretical prediction which is shown as the dotted line. For junction of inequivalent wires (g α = g β ), the results are very close to the case of noninteracting wires. However, we do observe a small but systematic deviation. The results agree with our theoretical conjecture that only a weak change of the conductance is expected in this case.

VII. SUMMARY AND DISCUSSION
In summary, we have developed a novel numerical framework to study the transport properties of two-wire junctions with inequivalent TLL wires, based on a finite window embedded in an infinite wire. Within this frame-work, we implemented two schemes. In the first scheme, the linear conductance is extracted from the asymptotic value of the static current-current correlation function, through a relation which is derived via the bosonization method. It is worth mentioning that the relation derived in this work corrected a subtle error in Ref. [43]. Our result agrees with Ref. [43] if and only two TLL wires have the same charge velocity, but it is generically different for two inequivalent wires. In the second scheme, we evaluate the time-dependent local current across the junction after a bias is applied. The current-voltage characteristic is then extracted by averaging the local current across the junction in the quasistationary region. In particular, the linear conductance is estimated by applying a small bias.
The main advantage of our schemes is to always work with an infinite system to eliminate the finite-size effects but only perform measurements within a finite window to make the calculation tractable. We benchmark our schemes against known theoretical results. For the timeindependent calculations, we show that the asymptotic behavior of the current-current correlation function can be obtained with a moderate window size. Furthermore, the linear conductance extracted agrees excellently with the theoretical prediction. For the time-dependent calculations, we show that with moderate window size and maximum virtual bond dimension, a long quasistationary region can be reached. By averaging the current within the quasistationary region one can accurately determine the current-voltage characteristics. Furthermore, we obtain the linear conductance by applying a small bias and the results agree excellently with the theoretical prediction as well as the results via the time-independent method. It is worth pointing out that computationally it is less demanding to use the time-dependent method to reach the same accuracy. This is because for the timeindependent method, it is essential to have high precision results of the correlation function at large distance to obtain its asymptotic behavior. However, the values of large distance correlations are quite small, making it more demanding to calculate accurately. On the other hand, for the time-dependent method accurate results can be obtained as long as a long quasistationary region with small residual oscillations can be reached.
Some comments are now in order. First, it is straightforward to generalize both schemes to study the transport properties of complex multiwire junctions. For complex geometry, one can study the Y junction with three TLL wires. By changing the enclosing magnetic flux and the Luttinger parameters, many conductance fixed points can be reached. It is also possible to include spin degrees of freedom and study the conductance of nanostructures. Going beyond linear response, it is interesting to study the single impurity Anderson model, where Kondo physics is important. In summary, our schemes allow one to determine accurately the correlation functions to a very large distance, simulate the dynamics to a very large time scale, and the results are free of finite size effects. We believe that these schemes can become important tools to study the transport properties of strongly interacting nanoscopic systems.