Dynamical conductivity of disordered quantum chains

Abstract We study the transport properties of a one-dimensional quantum system with disorder. We numerically compute the frequency dependence of the conductivity of a fermionic chain with nearest-neighbor interaction and a random chemical potential by using the Chebyshev matrix product state (CheMPS) method. As a benchmark, we investigate the noninteracting case first. Comparison with exact diagonalization and analytical solutions demonstrates that the results of CheMPS are reliable over a wide range of frequencies. We then calculate the dynamical conductivity spectra of the interacting system for various values of the interaction and disorder strengths. In the high-frequency regime, the conductivity decays as a power law, with an interaction-dependent exponent. This behavior is qualitatively consistent with the bosonized field theory predictions, although the numerical evaluation of the exponent shows deviations from the analytically expected values. We also compute the characteristic pinning frequency at which a peak in the conductivity appears. We confirm that it is directly related to the inverse of the localization length, even in the interacting case. We demonstrate that the localization length follows a power law of the disorder strength with an exponent dependent on the interaction, and find good quantitative agreement with the field theory predictions. In the low-frequency regime, we find a behavior consistent with the one of the noninteracting system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ^{2}(\ln \omega )^{2}$$\end{document}ω2(lnω)2 independently of the interaction. We discuss the consequences of our finding for experiments in cold atomic gases. Graphical abstract


Introduction
Disorder has profound effects on quantum systems, as directly evidenced in the celebrated Anderson localization [1]. As shown by Anderson, disorder can change the plane waves of free particles to exponentially localized states with spectacular consequences for transport. Anderson localization is relevant for a host of systems ranging from condensed matter to classical waves. In particular cold atomic gases, due to their remarkable controllability, have been instrumental in evidencing the localization of the wavefunction through seminal experiments of the groups of A. Aspect with laser speckles [2,3] and M. Inguscio with quasi-periodic potentials [4,5]. After sixty years since its proposal, Anderson localization still continues to present new challenges and mathematical developments [6].
The combination of disorder and interactions poses an additional layer of challenge, especially in the context of condensed matter physics. This complicated problem in thermalized systems was tackled by perturbative [7], or renormalization group (RG) techniques in one dimension [8] and two dimensions [9][10][11]. While the disorder basically decelerates particles, which leads to a reinforcement of interactions, it can also weaken them a e-mail: Thierry.Giamarchi@unige.ch (corresponding author) due to the exponentially small overlap between two localized states. This competition is highly nontrivial and constitutes an intensively studied topic, which we call the problems of localization of interacting particles. Another direction related to this problem is to study the thermalization and ergodicity in isolated quantum systems with disorder and interactions, which is known under the name of many-body localization [12][13][14][15].
Our target in this paper is one-dimensional systems, where strong quantum fluctuations lead to special states such as the Tomonaga-Luttinger liquid (TLL) characterized by correlations decaying in power law [16], and the effects of interactions are particularly strong. Furthermore, the disorder effects are also at their maximum, and even an infinitesimal disorder localizes all states for noninteracting particles. Thus, one expects a severe competition between disorder and interactions. Renormalization group (RG) analysis shows the existence of a localized-delocalized transition both for fermions and for bosons [8]. The localized phase for bosons, the Bose glass, persists in higher dimensions [17], and cold atomic systems again provide a controlled experimental access, confirming the Bose glass phase in biperiodic systems [18,19].
To characterize disordered systems, transport is an important property. RG can access DC conductivity down to the temperature related to the inverse local- (1) without disorder. The horizontal dashed line is K = 3/2, below which the system is localized by introducing the disorder ization scale [8]. Below this scale, more phenomenological calculations predict the Mott variable range hopping behavior in the thermalized case of localization of interacting particles [20,21] and zero DC conductivity in the isolated case of many-body localization [22]. AC transport also reflects the competition between disorder and interactions. Dynamical conductivity is exactly known for Anderson localization [23][24][25]. At high frequency, the behavior for the interacting particles can be extracted from the RG [8], and the low-frequency behavior has been investigated by approximate methods such as a variational approach [26]. Despite these efforts, no general methods are applicable to the full frequency range for dynamical conductivity. This situation is regrettable since cold atoms would be perfect systems to investigate such AC behavior of the conductivity with methods such as phase shaking of the optical lattice [27][28][29]. Indeed in biperiodic lattices, signatures of the localization such as a localization peak in the amplitude shaking of the optical lattice have been predicted [30] and observed [19].
In the present paper, we study the dynamical conductivity, i.e., the AC transport property as a function of frequency, in simple spinless fermion chains with nearest-neighbor interactions in a random chemical potential. We perform numerical calculations using a variant of a density matrix renormalization group (DMRG) method to compute the dynamical quantity of disordered quantum systems with good precision. We compare the obtained result for dynamical conductivity with the field theory and discuss its AC behavior over the full-frequency regime. Such fermionic systems can be mapped either to Ising anisotropic spin chains in a random magnetic field or to hard core bosonic chains in a random chemical potential. Thus, the system on which we focus is quite generic to demonstrate the applicability of our method and to study the physics of the dynamical conductivity of one-dimensional disordered quantum systems. This paper is organized as follows: In Sect. 2, we introduce the model of spinless fermions with a random chemical potential and nearest-neighbor interactions, which is a target of this paper. This model is con-nected to the XXZ spin chain with a random magnetic field. We also give the expression of dynamical conductivity using the Kubo formula. Section 3 explains the numerical technique called Chebyshev matrix product state (CheMPS), which we mainly utilize to investigate the interacting system. In Sect. 4, we describe the calculated results of the dynamical conductivity obtained by the numerics and its behavior in the high-and lowfrequency regimes. The numerical results are compared with analytical prediction from field theory. We summarize our results and discuss future problems in Sect. 5.

Model and physical quantities
Let us consider a spinless fermion system with a nearest-neighbor interaction where N is the number of sites,â † l (â l ) is the fermion creation (annihilation) operator at the site l, andn l ≡â † lâ l is the number operator. The random chemical potential h l on each site distributes uniformly in the finite interval h l ∈ [−W, W ]. Thus, W represents the strength of disorder. In view of the numerical solution of this model, we assume open boundary conditions, while the analytic solutions are usually performed with periodic boundary conditions. Without disorder, this system is particle-hole symmetric and half-filled. In this case, the system is known to be described at low energy by a TLL hamiltonian (see "Appendix A"). Since this model is solvable by Bethe ansatz, the TLL parameters can be exactly computed. For example, the parameter K controlling the decay of correlation functions [16] is given by: The Jordan-Wigner transformation: maps the fermionic model Eq.
(1) to a spin 1/2 chain with anisotropy along the z axis, which is called XXZ model, in a random magnetic field.
whereŜ x,y,z l is the spin-1/2 operator and S ± l = S x l ± iS y l . For simplicity, we employ the unit = c = 1. Additionally, the spins can be mapped onto hard core bosons by the mappingŜ + l =b † l andŜ z l =b † lb l − 1 2 . For Δ = 0, the Hamiltonian Eq. (1) represents free fermions with a random chemical potential. It is known that in one dimension, such a system is always exponentially localized for W > 0 [1]. The situation becomes complicated in the interacting case Δ = 0, but an analysis using field theory and bosonization is possible [8,16]. Such an analysis, using a RG procedure, shows the existence of a quantum phase transition between a localized and delocalized phases. The system is localized when K is smaller than 3/2 even for an infinitesimal W . The dependency of K on Δ [Eq. (2)] is shown in Fig. 1, and it can be seen that K < 3/2 corresponds to Δ > −1/2.
Besides the phase diagram itself the main physical quantity we will be computing in this paper is the frequency dependence of the conductivity. All calculations will be done at zero temperature. We use the Kubo formula relating the conductivity to current-current correlations. The current operator is: and the current-current retarded correlation function is written as: where ϑ(t) is the step function andĵ tot = N −1 l=1ĵ l . We denote byÔ(t) the usual Heisenberg time evolution The dynamical conductivity can be obtained by the Fourier transform of the retarded correlation function Eq. (6), where σ(ω) and σ (ω) represent the real and imaginary part of the dynamical conductivity, respectively. = 0 + is an infinitesimal convergence factor. Note that in the above expression we have not written the so-called diamagnetic term, which is purely imaginary, since we will concentrate on the real part of the conductivity σ(ω).
In the spectral representation, the real part of the conductivity is rewritten from Eq. (7) as where |Ψ 0 (|Ψ ν ) is the ground (ν-th excited) state and E 0 (E ν ) is its energy eigenvalue. We will use the expression Eq. (8) for the numerical evaluation of the conductivity.

Methods
We explain the numerical method that we use to treat the interacting systems with disorder. In Sect. 4, we discuss the dynamical conductivity of such systems by comparing the results from this numerical method with those from a field theory for the low energy limit of this model. The details of the field theory are described in "Appendix A".
To tackle the problem of low-dimensional interacting quantum systems, numerical methods utilizing matrix product states, such as DMRG [31,32], are very effective. The dynamical quantities such as conductivity and Green's function can be computed by performing a realtime evolution with, for example, time-evolving block decimation, after obtaining the ground state by DMRG, and such techniques have been widely used [33][34][35]. The spectral functions are calculable through the Fourier transformation of the temporal correlation functions. However, in the real-time evolution of matrix product states, the entanglement of the systems grows exponentially and the achievable time interval is limited. The acquired frequency resolution in this way was not sufficient for our purpose. Therefore, we used a numerical method which calculate the spectral functions directly in the frequency space.
In particular, to perform our calculation of the conductivity, we focus on the method CheMPS [36]. This is a combination of DMRG and the kernel polynomial method [36,37], a method to evaluate the spectral function where |Ψ 0 is the ground state and E 0 is its energy. We assume that the spectra have nonzero weight in ω ∈ [0, Ω], which can be mapped to the interval ω ∈ [−1 + s , 1 − s ] by redefining the energy scale as where s is a small safety factor. We take s = 0.0125 in this study.
The Hamiltonian is then mapped to: and the spectral function Eq. (9) becomes Using the Chebyshev polynomials we can expand the spectral function as The Chebyshev moments are represented as μ n = Ψ 0 |Ô 1 |t n , where |t n = T n (H )Ô 2 |Ψ 0 are Chebyshev vectors. The recurrence equations are useful to evaluate the coefficients μ n numerically. In the numerical calculations, the expansion of Eq. (14) is performed up to some finite-order M , and we multiply the weight μ n by the Jackson damping factor [37] M + 1 (16) to smoothen the spectrum. Therefore, the spectral function is numerically obtained as: In this study, we calculate the ground state |Ψ 0 using DMRG and then obtain the matrix product state representation of |t n by the recurrence equations Eq. (15). The system size is N = 250, the energy width is Ω = 6, the bond dimension of matrix product representation is M B = 64, and the order of expansion is M = 200.
The dynamical conductivity is calculated from the current-current correlation function as given by Eq. (7). However, in the low energy region ω/J 1, the 1/ω factor in the right-hand side of Eq. (7) enhances the numerical error. Hence, to avoid this problem we employ the polarization-current correlation function instead of the current-current one. The polarization operator is defined asP and it is related to the current operator through the time derivative The polarization-current correlation function becomes and its time derivative is the current-current correlation function Performing the integral by part in Eq. (7), we obtain Thus, the dynamical conductivity is represented as which can be evaluated by CheMPS withÔ 1 =P and O 2 =ĵ tot in Eq. (9).

Results
We now examine the results obtained by the numerical method described in Sect. 3. Before dealing with the interacting case, which is the focus of our study, let us first discuss the dynamical conductivity in the noninteracting case, which corresponds to Anderson localization, where both simpler numerical solutions and analytical approaches are available. Through the comparison of the results from CheMPS with from simpler methods, we can provide a benchmark for the reliability and applicability of CheMPS.

Noninteracting case
Let us first revisit the noninteracting case Δ = 0. Without interactions, Eq. (1) is a tight binding system of free spinless fermions at half-filling.
If one linearizes the dispersion relation around the Fermi wave number k = ±π/2 in the continuum limit. This model reduces to a Dirac model with a random chemical potential for which an exact analytical solution was obtained by Berezinskii [23] and the conductivity could be computed analytically [23][24][25]38].
The expression for the dynamical conductivity is given as [25]: where Q n and R n are the solution of the recurrence equations with the boundary condition In practice, starting from the initial condition Q n = 0 and R n = 0 for large enough n, we can obtain numerically Q n−1 , . . . , Q 0 and R n−1 , . . . , R 0 from Eqs. (19) and (20) and finally normalize the sequence {Q 0 , . . . , Q n } and {R 0 , . . . , R n } so as to satisfy the condition Eq. (21). Here, τ i and σ 0 are fitting parameters. Figure 2 shows the dynamical conductivity σ(ω) calculated by the above procedure. In the high-frequency region, the dynamical conductivity decays as a power law σ(ω) ∝ ω −2 . In the low-frequency region, the predicted analytical behavior is σ(ω) ∝ ω 2 (ln ω) 2 [20,23,38]. As can be seen from Fig. 2, such a behavior fits much better the data than σ(ω) ∝ ω 2 .
To compare with the CheMPS calculations, we also perform exact diagonalization directly on the lattice model Eq. (1) with Δ = 0. We evaluate the dynamical conductivity by In this equation, we have introduced a nonzero broadening δ 0 to the delta function in Eq.  The dynamical conductivity of the lattice-free fermion model is well fitted by the curve obtained from Eq. (18), confirming that the lattice model correctly captures the σ(ω) ∝ ω 2 (ln ω) 2 behavior in the low-frequency region and σ(ω) ∝ ω −2 in the high-frequency region. Given the finite bandwidth J of the lattice model, above ω/J > 1, σ(ω) decays exponentially and deviates from the curve obtained from Eq. (18).
In addition to the asymptotic behaviors at small and large frequency, we can extract the pinning frequency ω p and the peak value σ(ω p ) from the dynamical conductivity curves shown in Fig. 3a. The quantities ω p and σ(ω p ) are plotted in Fig. 3b as a function of W . We observe that ω p and σ(ω p ) are scaled with the disorder strength W as This is the expected behavior since the pinning frequency is directly related to the inverse localization length, which for the noninteracting case scales as the mean free path in one dimension and thus as ξ ∝ W −2 , as can be seen from Eqs. (A16) and (22). The results for the noninteracting case also serve as a benchmark of the CheMPS method in the next section.

Interacting case
Let us now turn to the interacting disordered systems. In this case, since the dimension of the Hilbert space grows exponentially by increasing the system size, ED is not a viable method any more, and we employ the CheMPS method described in Sect. 3 as a numerical approach as well as the bosonization and variational replica approach as analytical methods to calculate the dynamical conductivity.
To benchmark the method, let us first compare the results of CheMPS for system size N = 250 with Δ = 0 and disorders W/J = 0.1, 0.2, 0.4, and 0.8 with the ones obtained in the previous section with ED. The comparison is shown in Fig. 4. As can be seen from the To check the finite size effect, we also compare the same CheMPS data with the ED result for N = 3200 (the same data as Fig. 3) in Fig. 4b. While the agreement is good for large W , there is a deviation for small W . However, in the high-energy region, the deviation is just a multiplication of a constant factor, and the power of the decay does not change. This confirms that our numerical approach properly captures the behavior of the dynamical conductivity in both low-and highfrequency regions at zero temperature. Hence CheMPS is a very useful technique for dealing with interacting systems, for which no other quantitative method is available. In Fig. 5a, we show the numerically calculated dynamical conductivity for the cases of Δ = 0, 0.2, and 0.4. We can see that the decay power of the highfrequency region changes as we increase the interaction Δ. To discuss the high-frequency region and to compare the results with the bosonized field theory, a relatively small disorder strength is desirable, and we adopt W/J = 0.1 here. We fit the data in the high-frequency region by a power law ∝ ω −μ and plot μ as a function of Δ in Fig. 5c. We confirm that the variance of the data for dynamical conductivity shown in Fig. 5a is negligibly small by comparing the results for three bins of the various disorder configurations over which the average is taken. The variations of the data in Fig. 5c mainly arise from the power law fitting, and we estimate the error bars from the fittings in several frequency intervals on the high-frequency regime.
Let us compare the numerical results with the prediction from the field theory on the continuum model [8] (see "Appendix A"). The behavior in the highfrequency regime remains a power law, while the exponent is indeed modified by introducing the interaction.
Note that the precise analytical form is more complicated than a simple power law since the TLL parameter K is renormalized and depends on the scale. However, far from the transition point K = 3/2, the simple power law decay that neglects this renormalization becomes a good approximation [8,16]. The modification of the exponent by the interaction reflects the renormalization of the scattering on the disorder and the power law behavior of the susceptibility of the charge-charge correlations in TLL. The analytical prediction for the exponent is μ = 4 − 2K [8,16], which naturally reproduces the exponent μ = 2 for the noninteracting case K = 1. The parameter K takes the value of K < 1 for repulsive interactions and K > 1 for attractive ones. As seen in Fig. 5c, the numerically evaluated μ has a similar global trend as the expectations from the field theory (dashed line), which is obtained by substituting the Bethe ansatz evaluation of K for TLL in the XXZ chain Eq. (2) into μ = 4 − 2K. In particular, in the region of −0.2 Δ 0.5, the numerically calculated μ agrees quantitatively well with the analytical prediction. Figure 5b shows the plotting of (ω/J) 2 σ(ω) for the same data as Fig. 5a. The plotting of the data for Δ = 0 is almost horizontal, which indicates the decay follows σ(ω) ∝ ω −2 . We can see that the decay power μ increases as Δ becomes larger.
However, there exist surprising quantitative differences between the numerical results and the field theoretical predictions in the attractive regime (Δ −0.2) and the strong repulsive regime (Δ 0.5). The origin of these discrepancies is not clear at the present stage. On the repulsive side, a plateau-like structure appears in the regime Δ 0.5 for the numerical results, which is incompatible with the exponent predicted by the continuum model. Several reasons are conceivable for this behavior. One possibility is the effect of irrelevant operators neglected in the field theoretical treatment. In particular, the scaling dimension of the irrelevant operator cos(4φ) representing the umklapp processes on the lattice lowers as Δ is increased, and it finally becomes marginal in the Δ → 1 limit. Another possibility is an error of the numerical extraction of μ. As seen from Fig. 5a, the localization frequency scale, i.e., the pinning frequency characterized by the peak of the dynamical conductivity, shifts to higher frequency as Δ is increased. Hence, the region where the curve is fitted by the power law becomes narrower, and the estimation becomes less precise. On the attractive side, if we approach the localization-delocalization transition The solid lines represent the power-law fitting in the high energy region. The frequency dependence of the conductivity is well represented by an interaction-dependent exponent at high frequencies. b The plotting of (ω/J) 2 σ(ω) for the same data as (a). c Exponent, as a function of the interaction Δ, of the decay of the conductivity with frequency in the high-frequency region. This exponent results from fits σ(ω) ∝ ω −µ similar to the ones of (a). The disorder strength is W/J = 0.1. The dashed line is the theoretical value of the exponent from the field theory μ = 4 − 2K and the Bethe-ansatz value of the TLL parameter K (see text). For repulsive interactions Δ > 0, there is a reasonable agreement up to Δ ∼ 0.5 beyond which there is a plateau like behavior not expected by the field theory. On the attractive side Δ < 0 strong deviations are observed even for relatively small attraction point K = 3/2 (Δ = −0.5), the renormalization flow is closer to the separatrix and the direction of the flow is not parallel to the disorder axis. Thus, the effect of renormalization of the parameter K should become more important. To elucidate the reasons for this discrepancy between the numerics and the field theory is a very challenging problem and we leave it for a future study.
Let us now turn to the disorder dependence of the pinning frequency and the conductivity at the peak for the interacting system. In Fig. 7a, we show the dynamical conductivity calculated for the interaction Δ = 0.5 and various disorder strengths W/J = 0.1, 0.2, 0.4, 0.8 (Fig. 6).
Similarly than for the noninteracting case [see Fig. 3a], the pinning frequency ω p increases and the peak height σ(ω p ) decreases as the disorder strength W is increased. We plot ω p and σ(ω p ) as a function of W in Fig. 3b. The pinning frequency and the peak height are well fitted as: Note that the data points of ω p and σ(ω p ) for W/J = 0.1 deviate from the fitting curves. We attribute it to the large finite size effect in the case of small W , as mentioned in the benchmark result (Fig. 4b).
This scaling is in good agreement with the analytical predictions that ω p ∝ ξ −1 and σ(ω p ) ∝ ξ [see Eq. (A16) in "Appendix A"]. Using the dependence of K on interactions Eq. (2), we obtain K = 3/4 for Δ = 0.5. This leads to ξ ∝ W −4/3 by using the formula Eq. (A14) in excellent agreement with the numerical results. As can be seen both from the numerics and the above formula, repulsive interactions lead to a shorter localization length than for the noninteracting case.
Finally let us discuss the behavior in the lowfrequency regime. In order to get a sizable range of the frequency interval below the pinning peak, and to prevent finite size effects from playing a major role, we take a relatively large value of disorder W/J = 0.8. The dynamical conductivity calculated for Δ = 0, 0.3, 0.5, and 0.8 is shown in Fig. 7a. While there is a clear dependence of the exponent of the power law decay on the interaction Δ in the high-frequency regime, the behavior remains similar on the low-frequency side. One can see it more clearly in Fig. 7b, where the curves have been rescaled by the value of the pinning frequency and conductivity at the peak.
The fitting curves by σ(ω) ∝ ω 2 (ln ω) 2 and σ(ω) ∝ ω 2 are also shown, and the former fitting looks better than the latter. However, the present interval of data fitting below the peak is just one decade of frequency, and the acquisition of the data over a wider range is desirable for a more precise analysis.
On a qualitative side, one would indeed expect to recover at low enough frequencies essentially the noninteracting behavior. Indeed, frequencies lower than the pinning frequency correspond to probing scales large compared to the localization length. At that scales, since the particles are exponentially localized the effect of interactions should drastically decrease, leading back to the noninteracting behavior. On a more quantitative level, it is difficult to make an unambiguous compari- son with existing analytical formulas, since the various calculations of the low-frequency behavior suffer from their own limitations. The variational calculation [26] is unable to capture the logarithmic correction to the ω 2 behavior. Calculations containing the logarithmic correction rely either on an extreme classical limit K → 0 [39], which is far from the values of K reached here or an instanton expansion [21,40]. Given the rescaling of the curves in Fig. 7 and the fact that the different curves broadly superimpose, this suggests that the prefactor of the ω 2 (ln ω) 2 term varies as: from the analytical predictions ω p ∝ ξ −1 and σ(ω p ) ∝ ξ [see Eq. (A16) in "Appendix A"].

Discussions and summary
In this paper, we have numerically computed the frequency dependence of the conductivity in onedimensional spinless fermions with a nearest-neighbor interaction and a random chemical potential. This problem is equivalent to XXZ spin chains in a random magnetic field along the z axis and to hard core bosons with nearest-neighbor interactions and a random chemical potential. Using the CheMPS method, a variant of DMRG, adapted to the case we study, we have numerically calculated the dynamical conductivity over a broad range of frequencies, interactions, and disorder strength. We have benchmarked our method by comparisons with the noninteracting case. Since analytical approaches and numerical exact diagonalizations are applicable for the noninteracting systems, these results have been compared with those from CheMPS and we have confirmed that our method does capture the frequency dependence of the conductivity in a broad range of frequencies.
We have then investigated the dynamical conductivity of the interacting systems with CheMPS. In the high-frequency regime, the conductivity follows a power law σ(ω) ∝ ω −μ . We have calculated μ for various interaction Δ, and found it agrees quantitatively with the expectation from the field theory μ = 4−2K and the K value from BA in −0.2 Δ 0.5. However, there exists a deviation in the attractive and strongly repulsive regions. On the attractive side Δ −0.2 (K 1.15), the estimated exponent μ seems larger than expected from the TLL determination. On the repulsive side, a plateau-like structure appears in the region Δ 0.5 (0.5 ≤ K 0.75). We leave the clarification of the reasons for this deviation for future studies. We have also evaluated the localization length ξ from the pinning frequency ω p ∝ ξ −1 as a function of the disorder strength and found a reasonable agreement with the expected behavior of the localization length as determined by RG: ξ ∝ (1/W 2 ) 1/(3−2K) and which is now dependent on both interactions and disorder. In the low-frequency regime, we have performed numerical calculations for a large disorder W/J = 0.8. The scaling of dynamical conductivity is compatible with an ω 2 (ln ω) 2 behavior similar to the one of the noninteracting case but with a prefactor varying as (ξ(Δ)/ξ(Δ = 0)) 3 .
The low-frequency behavior is difficult to access, and although the numerics is indeed compatible with the ω 2 (ln ω) 2 behavior to ascertain the existence and power of the log correction data over a wider range of frequency is needed, another challenge for future studies.
In addition to pushing the numerical investigations of the conductivity, it would of course be extremely interesting to test the results obtained in our study in cold atom experiments. The existence of a peak as a response to shaking, similar to the pinning peak in the conductivity discussed here, was indeed observed for bosons in a quasiperiodic potential and has been used as a key signature of the existence of the Bose glass in these systems. However, due to the large inhomogeneities arising from the quadratic trapping potential, testing the power law behavior was not practically feasible. The existence of fermionic systems with disorder in quantum microscopes makes it possible to observe the features computed here more easily. In preforming a comparison with experiments, we need to note the following points: (i) the experimental system size realizable for the moment is still relatively small, typically of the order of 20 to 50 atoms per chain; (ii) the temperature is still relatively high in fermionic systems. These conditions should not drastically affect the properties in the high-frequency regime, roughly for ω > T , but of course will essentially change the low-frequency behavior of the dynamical conductivity. The low-frequency regime is ultimately connected to the question of variable range hopping and many-body localization. Addressing these issues via experiments and by the extension of numerical techniques to finite temperature is clearly a considerable challenge.
Acknowledgements The authors appreciate fruitful discussions with Michele Filippone. S. T. is supported by JSPS KAKENHI Grant No. JP21K03412 and JST CREST Grant No. JPMJCR19T3, Japan. This work is partly supported by the Swiss National Science Foundation under Division II. This paper is dedicated to Alain Aspect, who in addition to his numerous contributions in atomic physics and quantum information is responsible, via key experiments, to having made, the subject of disorder a very hot topic in cold atomic systems. Thank you Alain for being a constant source of inspiration and stimulation, also in this field of "dirty" quantum systems.

Author contributions
Both authors contributed equally to the idea behind the work, the analytical calculations and the writing of the manuscript. ST devised and implemented the numerical CheMPS methods used for the numerical calculations.

Funding Information Open access funding provided by University of Geneva
Data Availability Statement This manuscript has associated data in a data repository. [Authors' comment The data are available for download at the address https://doi. org/10.26037/yareta:4hpihszjk5h7pegfpvdq36rtra [DOI will be provided in the published paper] and will be kept for 15 years.]

Conflict of interest
The authors declare that they have no known competing financial interests that could have influenced the work reported in this paper Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A Field theory solution
As a framework to discuss the numerical solution, let us give here a short summary of the field theory approach to this problem [8]. The most convenient way to deal with the interactions is to use the bosonization representation [16].
Through the Fourier transformation of the fermionic operators (the lattice constant is set to unity): the disorder term in the Hamiltonian becomeŝ For fermionic systems, Eq. (1) with a Fermi momentum k F = π/2 (half-filled band) corresponds to the nonmagnetized case in spin chains. For such systems, special care should be exerted to treat the disorder and the situation is slightly more involved than starting from an incommensurate filling. As is well known in one dimension, the system has low energy excitations at momenta q 0 and q 2k F [16]. It is thus useful to separate the disorder into two slowly varying fields centered around these momenta, which are called forward and backward scattering, respectively, (A3) Note that in the case of 2k F = π, the backward scattering leads to a real field h b (x). Since h f (f ) and h b (x) involve components of h(q) which are different in q, their cross averages are zero and they can be considered as two independent random fields with essentially delta function correlations upon average. The fermionic field ψ(x) = k e ikxâ k can be represented in terms of right-and left-movers [16] and the disorder term is recast intô by noting h * b (x) = h b (x). The current operator becomes: where v F is the Fermi velocity. The forward scattering part of the disorder can be eliminated by the gauge This transformation does not affect the current operator, but changes the backscattering term to One thus see that the backscattering part of disorder is replaced by a complex field ξ( with correlations ξ(x)ξ(y) = 0, where D b ∝ W 2 . This reflects the breaking of the particle-hole symmetry, which is caused by one realization of the random chemical potential. Note that a system with perfect particle hole symmetry, such as a random bond system, would lead to a different fixed point, namely the random singlet phase [41]. We employ the usual representation of the fermion operators in terms of collective fields related to density and current, respectively [16], and connect them in the spin language [Eq. (4)] to the two angles needed to dictate the direction of the spin vector S + l (−1) l e iθ(x) + e iθ(x) cos(2φ(x)), S z l −∇φ(x) π + (−1) l cos(2φ(x)), where x is the position of l-th spin, and the prefactors are omitted. Then, the Hamiltonian is written as [8,16] Here u is the velocity of excitations (u = v F for the noninteracting case) and K is the dimensionless parameter controlling the decay of correlations. We have neglected the irrelevant operator cos(4φ(x)) which appears for the special case of nonmagnetized spin chains, or equivalently half-filled fermion chains. This operator is irrelevant for K > 1/2. In this representation, the current operator is the time derivative of the field φ(x) [16], and thus, the conductivity Eq. (7) is simply related to the Green's function of the field φ(x) by where σ is in units of e/h and = 0 + is a convergence factor. The disorder term in Eq. (A11) can be eliminated for a Gaussian disorder by using the replica trick [8,16]. This leads to an action of the form where α, ν = 1, 2, . . . , n are the replica indices and the limit n → 0 must be taken. The disorder term is relevant for K < 3/2 even for an infinitesimal strength of the disorder D b . For K > 3/2 (Δ < −0.5), there is a separatrix, with a transition of the Berezinskii-Kosterlitz-Thouless universality class between the localized and a delocalized TLL phase [8].
The localization length can be captured by an RG evaluation of this problem. Far from the transition point K = 3/2, the localization length behaves as [8] ξ In addition to the RG method which can only give access to the physical quantities, either in the regime where the disorder is irrelevant or for scales smaller than the localization length, physical observables can be computed in the localized phase using a replica Gaussian variational approach [26]. This approximate approach leads to a low-frequency behavior of the real part of the conductivity as: where ξ is the localization length. It does not include the logarithmic correction, which is known to occur for noninteracting particles σ(ω) ∝ ω 2 (ln ω) 2 [23,25]. This is clearly an artifact of the variational approach, and a logarithmic correction is also a priori expected for interacting systems as suggested by semiclassical calculations K → 0 [39] and instanton expansions [21,40]. Equation (A15) implies that the pinning frequency ω p and the corresponding value of the dynamical conductivity σ(ω p ) scale as in agreement with the predictions from the RG calculations as well. We use this scaling in comparison with the numerical results.